ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/echarge.c
Revision: 115
Committed: Thu May 7 14:35:26 2009 UTC (11 years, 3 months ago) by gilbertke
File size: 6133 byte(s)
Log Message:
added charge-charge routines
Line File contents
1 #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 printf("\nCharge-Charge Interactions : Diele = %7.3f\n",units.dielec);
32 printf(" At1 At2 q1 q2 Rik Eqq\n");
33 }
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 printf("QQ: (%-3d)- (%-3d) %-8.3f %-8.3f %-8.3f = %8.4f\n", i, j,charge[i], charge[j],rik,e);
69 }
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 }