ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/esolv.c
(Generate patch)
# Line 12 | Line 12
12      double *shct,*asolv,*rsolv,*vsolv,*gpol,*rborn;
13      } solvent;
14  
15 < void esolv(int natom,double chrgcut,int **skip,int *mmx_type,int *use,double *charge,double *x,double *y,double *z,double *esolv);
16 < void esolv1(int natom,double chrgcut,int **skip,int *mmx_type,int *use,double *charge,double *x,double *y,double *z,double *esolv,double **desolv,double *drb);
17 < void esolv2(int ia,int natom,double chrgcut,int *mmx_type,double *charge,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz);
18 < void born1(int natom,int *mmx_type,int **skip,double *x,double *y,double *z,double **desolv,double *drb);
19 < void born(int natom,int *mmx_type,int **skip, double *x,double *y,double *z);
15 > void esolv(int natom,double chrgcut,int **skip,int *use,double *charge,double *x,double *y,double *z,double *esolv);
16 > void esolv1(int natom,double chrgcut,int **skip,int *use,double *charge,double *x,double *y,double *z,double *esolv,double **desolv,double *drb);
17 > void esolv2(int ia,int natom,double chrgcut,double *charge,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz);
18 > void born1(int natom,int **skip,double *x,double *y,double *z,double **desolv,double *drb);
19 > void born(int natom,int **skip, double *x,double *y,double *z);
20   // ===============================
21 < void esolv(int natom,double chrgcut,int **skip,int *mmx_type,int *use, double *charge,double *x,double *y,double *z,double *esolv)
21 > void esolv(int natom,double chrgcut,int **skip,int *use, double *charge,double *x,double *y,double *z,double *esolv)
22   {
23      int i,j;
24      double e,ai,ri,rb;
# Line 32 | Line 32
32      cc = 4.80298*4.80298*14.39418;
33      cutoff = chrgcut*chrgcut;
34  
35 <    born(natom,mmx_type,skip,x,y,z);
35 >    born(natom,skip,x,y,z);
36      term = 4.0*PI;
37      for (i=1; i <= natom; i++)
38      {
39      if (mmx_type[i] != 20)
40      {
39          ai = solvent.asolv[i];
40          ri = solvent.rsolv[i];
41          rb = solvent.rborn[i];
# Line 47 | Line 45
45              e = ai*term*(ri+probe)*(ri+probe)*term1*term1*term1;
46              *esolv += e;
47          }
48 <      }
48 >
49      }
50   // polarization term
51      f = -cc*(1.0/solvent.EPSin - 1.0/solvent.EPSsolv);
52  
53      for (i=1; i <= natom; i++)
54      {
55 <      if (charge[i] != 0.0 && mmx_type[i] != 20)
55 >      if (charge[i] != 0.0)
56        {
57           for(j=i; j <= natom; j++)
58           {
59             if (use[i] || use[j])
60             {
61 <              if (charge[j] != 0.0 && mmx_type[j] != 20)
61 >              if (charge[j] != 0.0)
62                {
63                   xi = x[i];
64                   yi = y[i];
# Line 85 | Line 83
83      }    
84   }
85   // =======================
86 < void esolv1(int natom,double chrgcut,int **skip,int *mmx_type,int *use,double *charge,double *x,double *y,double *z,double *esolv,double **desolv,double *drb)
86 > void esolv1(int natom,double chrgcut,int **skip,int *use,double *charge,double *x,double *y,double *z,double *esolv,double **desolv,double *drb)
87   {
88      int i,j;
89      double e,de,derb,ai,ri,rb;
# Line 107 | Line 105
105         desolv[i][2] = 0.0;
106     }
107  
108 <    born(natom,mmx_type,skip,x,y,z);
108 >    born(natom,skip,x,y,z);
109      term = 4.0*PI;
110      for (i=1; i <= natom; i++)
111      {
114      if (mmx_type[i] != 20)
115      {
112          ai = solvent.asolv[i];
113          ri = solvent.rsolv[i];
114          rb = solvent.rborn[i];
# Line 123 | Line 119
119              *esolv += e;
120              drb[i] -= 6.0*e/rb;
121          }
126      }
122      }
123   // polarization term
124      f = -cc*(1.0/solvent.EPSin - 1.0/solvent.EPSsolv);
125  
126      for (i=1; i <= natom; i++)
127      {
128 <      if (fabs(charge[i]) > 0.001 && mmx_type[i] != 20)
128 >      if (fabs(charge[i]) > 0.001)
129        {
130           xi = x[i];
131           yi = y[i];
# Line 140 | Line 135
135           {
136             if (use[i] || use[j])
137             {
138 <              if (fabs(charge[j]) > 0.001 && mmx_type[j] != 20)
138 >              if (fabs(charge[j]) > 0.001)
139                {
140                   xr = xi - x[j];
141                   yr = yi - y[j];
# Line 189 | Line 184
184        }
185      }
186   //
187 <    born1(natom,mmx_type,skip,x,y,z,desolv,drb);
187 >    born1(natom,skip,x,y,z,desolv,drb);
188   }
189   /* ===================================================== */
190 < void born1(int natom,int *mmx_type,int **skip,double *x,double *y,double *z,double **desolv,double *drb)
190 > void born1(int natom,int **skip,double *x,double *y,double *z,double **desolv,double *drb)
191   {
192      int i,k;
193      double xi,yi,zi,cc,rvdw;
# Line 217 | Line 212
212          pip5 = PI*solvent.p5;
213          for (i=1; i <= natom; i++)
214          {
220          if (mmx_type[i] != 20)
221          {
215              skip[i][i] = i;
216              xi = x[i];
217              yi = y[i];
# Line 226 | Line 219
219              gpi = 2.0*solvent.rborn[i]*solvent.rborn[i]/cc;
220              for (k = 1; k <= natom; k++)
221              {
222 <                if (skip[i][k] != i  && mmx_type[k] != 20)
222 >                if (skip[i][k] != i)
223                  {
224                      xr = x[k] - xi;
225                      yr = y[k] - yi;
# Line 263 | Line 256
256                      
257                  }
258              }
259 <          }
259 >          
260          }
261      } else if (solvent.type == HCT)
262      {
263          for (i=1; i <= natom; i++)
264          {
272          if (mmx_type[i] != 20)
273          {
265              xi = x[i];
266              yi = y[i];
267              zi = z[i];
# Line 278 | Line 269
269              rb2 = solvent.rborn[i]*solvent.rborn[i];
270              for (k=1; k <= natom; k++)
271              {
272 <                if (i != k  && mmx_type[k] != 20)
272 >                if (i != k)
273                  {
274                      xr = x[k] - xi;
275                      yr = y[k] - yi;
# Line 315 | Line 306
306                      }
307                  }
308              }
309 <          }
309 >
310          }
311      }
312   }
313   /*    ============================================== */
314 < void esolv2(int ia,int natom,double chrgcut,int *mmx_type,double *charge,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz)
314 > void esolv2(int ia,int natom,double chrgcut,double *charge,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz)
315   {
316      int jj,k;
317      double probe,cc,cutoff;
# Line 337 | Line 328
328      dwater = 78.3;
329      cc = 4.80298*4.80298*14.39418;
330  
331 <    if (fabs(charge[ia]) > 0.001 || mmx_type[ia] == 20)
331 >    if (fabs(charge[ia]) > 0.001)
332         return;
333          
334      cutoff = chrgcut*chrgcut;
# Line 348 | Line 339
339  
340      for (k = 1; k <= natom; k++)
341      {
342 <        if (ia != k && charge[k] != 0.0 && mmx_type[k] != 20)
342 >        if (ia != k && charge[k] != 0.0)
343          {
344              xr = xi - x[k];
345              yr = yi - y[k];

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines