ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/echarge.c
Revision: 119
Committed: Tue Jun 16 21:36:00 2009 UTC (11 years, 2 months ago) by wdelano
File size: 6169 byte(s)
Log Message:
printf( to fprintf(pcmlogfile,
Line User Rev File contents
1 gilbertke 115 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "utility.h"
5    
6     EXTERN struct t_minim_values {
7     int iprint, ndc, nconst;
8     float dielc;
9     } minim_values;
10    
11     void echarge(int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu);
12     void echarge1(int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu,double **deqq);
13     void echarge2(int i,int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,float **hessx,float **hessy,float **hessz);
14    
15     void echarge(int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu)
16     {
17     int i, j;
18     double xi, yi, zi, xk,yk,zk;
19     double xr, yr, zr;
20     double e, rik, rik2,cc;
21     double rdielc;
22     double cutoff;
23    
24     cutoff = chrgcut*chrgcut;
25     *eu = 0.0;
26     cc = 4.80298*4.80298*14.39418;
27     rdielc = 1.0/units.dielec;
28    
29     if (minim_values.iprint)
30     {
31 wdelano 119 fprintf(pcmlogfile,"\nCharge-Charge Interactions : Diele = %7.3f\n",units.dielec);
32     fprintf(pcmlogfile," At1 At2 q1 q2 Rik Eqq\n");
33 gilbertke 115 }
34     for (i=1; i < natom; i++)
35     {
36     if (charge[i] != 0.0 )
37     {
38     xi = x[i];
39     yi = y[i];
40     zi = z[i];
41    
42     for(j=i+1; j <= natom; j++)
43     {
44     if (use[i] || use[j])
45     {
46     if ( skip[i][j] != i)
47     {
48     if (charge[j] != 0.0)
49     {
50    
51     xk = x[j];
52     yk = y[j];
53     zk = z[j];
54    
55     xr = xi - xk;
56     yr = yi - yk;
57     zr = zi - zk;
58     rik2 = xr*xr + yr*yr + zr*zr;
59     if (rik2 < cutoff)
60     {
61     rik = sqrt(rik2);
62     e = cc*rdielc*charge[i]*charge[j]/rik;
63     if (skip[i][j] == -i)
64     e /= units.chgscale;
65     *eu += e;
66    
67     if (minim_values.iprint)
68 wdelano 119 fprintf(pcmlogfile,"QQ: (%-3d)- (%-3d) %-8.3f %-8.3f %-8.3f = %8.4f\n", i, j,charge[i], charge[j],rik,e);
69 gilbertke 115 }
70     }
71     }
72     }
73     }
74     }
75     }
76     }
77     // ====================================================
78     void echarge1(int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu,double **deqq)
79     {
80     int i,j;
81     double xi, yi, zi;
82     double xr, yr, zr;
83     double xk,yk,zk;
84     double e, rik, cc, rik2;
85     double de, dedx, dedy, dedz;
86     double rdielc;
87     double cutoff;
88    
89     cutoff = chrgcut*chrgcut;
90     cc = 4.80298*4.80298*14.39418;
91     rdielc = 1.0/units.dielec;
92    
93     for (i=1; i < natom; i++)
94     {
95     if (charge[i] != 0.0 )
96     {
97     xi = x[i];
98     yi = y[i];
99     zi = z[i];
100    
101     for(j=i+1; j <= natom; j++)
102     {
103     if (use[i] || use[j])
104     {
105     if ( skip[i][j] != i)
106     {
107     if (charge[j] != 0.0)
108     {
109    
110     xk = x[j];
111     yk = y[j];
112     zk = z[j];
113     xr = xi - xk;
114     yr = yi - yk;
115     zr = zi - zk;
116     rik2 = xr*xr +yr*yr + zr*zr;
117     if (rik2 < cutoff)
118     {
119     rik = sqrt( rik2);
120     e = cc*rdielc*charge[i]*charge[j]/rik;
121     if (skip[i][j] == -i)
122     e /= units.chgscale;
123    
124     de = -cc*rdielc*charge[i]*charge[j]/rik2;
125    
126     de = de/rik;
127     dedx = de*xr;
128     dedy = de*yr;
129     dedz = de*zr;
130    
131     deqq[i][0] += dedx;
132     deqq[i][1] += dedy;
133     deqq[i][2] += dedz;
134    
135     deqq[j][0] -= dedx;
136     deqq[j][1] -= dedy;
137     deqq[j][2] -= dedz;
138     *eu += e;
139     }
140     }
141     }
142     }
143     }
144     }
145     }
146     }
147     // ==========================================
148     void echarge2(int i,int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,float **hessx,float **hessy,float **hessz)
149     {
150     int j,k;
151     double fik,de,d2e,d2edx,d2edy,d2edz;
152     double xi,yi,zi,xr,yr,zr,term[3][3];
153     double xk,yk,zk;
154     double r,r2;
155     double cc, rdielc;
156     double cutoff;
157    
158     if (charge[i] == 0.0)
159     return;
160    
161     cutoff = chrgcut*chrgcut;
162    
163     cc = 4.80298*4.80298*14.39418;
164     rdielc = 1.0/units.dielec;
165    
166     skip[i][i] = i;
167    
168     xi = x[i];
169     yi = y[i];
170     zi = z[i];
171    
172    
173     for (k=1; k <= natom; k++)
174     {
175     if (i != k && charge[k] != 0.0 && (skip[i][k] != i) )
176     {
177     xk = x[k];
178     yk = y[k];
179     zk = z[k];
180    
181     xr = xi - xk;
182     yr = yi - yk;
183     zr = zi - zk;
184     r2 = xr*xr + yr*yr + zr*zr;
185     if (r2 < cutoff)
186     {
187     r = sqrt(r2);
188     fik = cc*rdielc*charge[i]*charge[k];
189     if (skip[i][k] == -i)
190     fik /= units.chgscale;
191     de = -fik/r2;
192     d2e = -2.0*de/r;
193     de = de / r;
194     d2e = (d2e-de) / r2;
195     d2edx = d2e * xr;
196     d2edy = d2e * yr;
197     d2edz = d2e * zr;
198     term[0][0]= d2edx*xr + de;
199     term[1][0] = d2edx*yr;
200     term[2][0] = d2edx*zr;
201     term[0][1] = term[1][0];
202     term[1][1] = d2edy*yr + de;
203     term[2][1] = d2edy*zr;
204     term[0][2] = term[2][0];
205     term[1][2] = term[2][1];
206     term[2][2] = d2edz*zr + de;
207    
208     for (j=0; j < 3; j++)
209     {
210     hessx[i][j] += term[j][0];
211     hessy[i][j] += term[j][1];
212     hessz[i][j] += term[j][2];
213     hessx[k][j] -= term[j][0];
214     hessy[k][j] -= term[j][1];
215     hessz[k][j] -= term[j][2];
216     }
217     }
218     }
219     }
220     }