ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/elj.c
Revision: 119
Committed: Tue Jun 16 21:36:00 2009 UTC (11 years, 10 months ago) by wdelano
File size: 8196 byte(s)
Log Message:
printf( to fprintf(pcmlogfile,
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "utility.h"
5
6 void elj(int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,double *evdw,double *e14);
7 void elj1(int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,double *evdw,double *e14,
8 double **devdw,double **de14);
9 void elj2(int iatom,int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,
10 float **hessx,float **hessy,float **hessz);
11
12 EXTERN struct t_minim_values {
13 int iprint, ndc, nconst;
14 float dielc;
15 } minim_values;
16
17 void elj(int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,double *evdw,double *e14)
18 {
19 int i,ia, ib,j, ita, itb;
20 double xi,yi,zi,xr,yr,zr;
21 double xk,yk,zk;
22 double e,p2,p6,p12,rv,eps;
23 double rik2;
24 double cutoff;
25
26 cutoff = vdwcut*vdwcut;
27
28 *evdw = 0.0;
29 *e14 = 0.0;
30
31 if (minim_values.iprint)
32 {
33 fprintf(pcmlogfile,"\nVDW Terms - Lennard-Jones Potential\n");
34 fprintf(pcmlogfile," At1 At2 Rik Radius Eps Evdw\n");
35 }
36 for (i=1; i < natom; i++)
37 {
38 for(j=i+1; j <= natom; j++)
39 {
40 if (skip[i][j] != i)
41 {
42 if (use[i] || use[j])
43 {
44 ia = i;
45 ib = j;
46 ita = type[ia];
47 itb = type[ib];
48
49 xi = x[ia];
50 yi = y[ia];
51 zi = z[ia];
52 xk = x[ib];
53 yk = y[ib];
54 zk = z[ib];
55
56 xr = xi - xk;
57 yr = yi - yk;
58 zr = zi - zk;
59 rik2 = xr*xr + yr*yr + zr*zr;
60 if ( rik2 < cutoff)
61 {
62 rv = vrad[ita][itb];
63 eps = veps[ita][itb];
64
65 if (skip[i][j] == -i)
66 eps /= units.v14scale;
67 p2 = (rv*rv)/rik2;
68 p6 = p2*p2*p2;
69 p12 = p6*p6;
70 e = eps*(p12 - 2.0*p6);
71 if (skip[i][j] == -i)
72 *e14 += e;
73 else
74 *evdw += e;
75 if (minim_values.iprint)
76 fprintf(pcmlogfile,"VDW: %-4d - %-4d %-8.3f %-8.3f %-8.3f = %8.4f\n",ia, ib, sqrt(rik2),rv,eps,e);
77 }
78 }
79 }
80 }
81 }
82 }
83 // =====================================
84 void elj1(int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,double *evdw,double *e14,
85 double **devdw,double **de14)
86 {
87 int i,ia, ib, j, ita,itb;
88 double xi,yi,zi,xr,yr,zr;
89 double xk,yk,zk;
90 double e,p2,p6,p12,rv,eps;
91 double rik2, rik;
92 double dedx,dedy,dedz,de;
93 double cutoff;
94
95 cutoff = vdwcut*vdwcut;
96 *evdw = 0.0;
97 *e14 = 0.0;
98
99 for (i=1; i <= natom; i++)
100 {
101 devdw[i][0] = 0.0;
102 devdw[i][1] = 0.0;
103 devdw[i][2] = 0.0;
104 de14[i][0] = 0.0;
105 de14[i][1] = 0.0;
106 de14[i][2] = 0.0;
107 }
108
109
110 for (i=1; i < natom; i++)
111 {
112 ia = i;
113 ita = type[ia];
114 xi = x[ia];
115 yi = y[ia];
116 zi = z[ia];
117 for(j=i+1; j <= natom; j++)
118 {
119 if (use[i] || use[j])
120 {
121 if (skip[i][j] != i)
122 {
123 ib = j;
124 itb = type[j];
125 xk = x[ib];
126 yk = y[ib];
127 zk = z[ib];
128
129 xr = xi - xk;
130 yr = yi - yk;
131 zr = zi - zk;
132 rik2 = xr*xr + yr*yr + zr*zr;
133 if (rik2 < cutoff)
134 {
135 rv = vrad[ita][itb];
136 eps = veps[ita][itb];
137 if (skip[i][j] == -i)
138 eps /= units.v14scale;
139 rik = sqrt(rik2);
140 p2 = (rv*rv)/rik2;
141 p6 = p2*p2*p2;
142 p12 = p6*p6;
143 e = eps*(p12 - 2.0*p6);
144 de = eps * (p12 - p6) * (-12.0/rik);
145 de = de / rik;
146 dedx = de * xr;
147 dedy = de * yr;
148 dedz = de * zr;
149
150 if (skip[i][j] == -i)
151 {
152 *e14 += e;
153 de14[ia][0] += dedx;
154 de14[ia][1] += dedy;
155 de14[ia][2] += dedz;
156 de14[ib][0] -= dedx;
157 de14[ib][1] -= dedy;
158 de14[ib][2] -= dedz;
159 } else
160 {
161 *evdw += e;
162 devdw[ia][0] += dedx;
163 devdw[ia][1] += dedy;
164 devdw[ia][2] += dedz;
165 devdw[ib][0] -= dedx;
166 devdw[ib][1] -= dedy;
167 devdw[ib][2] -= dedz;
168 }
169 }
170 }
171 }
172 }
173 }
174 }
175 // ===============================
176 void elj2(int iatom,int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,
177 float **hessx,float **hessy,float **hessz)
178 {
179 int ia, ib,j, k, ita, itb;
180 double xi,yi,zi,xr,yr,zr;
181 double xk,yk,zk;
182 double e,p2,p6,p12,rv,eps;
183 double rik2, rik;
184 double de,d2e;
185 double d2edx,d2edy,d2edz,term[3][3];
186 double cutoff;
187
188 cutoff = vdwcut*vdwcut;
189
190 skip[iatom][iatom] = iatom;
191 ia = iatom;
192 ita = type[ia];
193
194 xi = x[ia];
195 yi = y[ia];
196 zi = z[ia];
197 for (j=1; j <= natom; j++)
198 {
199 ib = j;
200 if (skip[iatom][j] != iatom)
201 {
202 itb = type[ib];
203
204 xk = x[ib];
205 yk = y[ib];
206 zk = z[ib];
207
208 xr = xi - xk;
209 yr = yi - yk;
210 zr = zi - zk;
211 rik2 = xr*xr + yr*yr + zr*zr;
212 if (rik2 < cutoff)
213 {
214 rv = vrad[ita][itb];
215 eps = veps[ita][itb];
216 if (skip[iatom][j] == -iatom)
217 eps /= units.v14scale;
218 rik = sqrt(rik2);
219 p2 = (rv*rv)/rik2;
220 p6 = p2*p2*p2;
221 p12 = p6*p6;
222 e = eps*(p12 - 2.0*p6);
223 de = eps * (p12 - p6) * (-12.0/rik);
224 d2e = eps * (13.0*p12 - 7.0*p6) * (12.0/rik2);
225 de = de / rik;
226 d2e = (d2e-de) / rik2;
227 d2edx = d2e * xr;
228 d2edy = d2e * yr;
229 d2edz = d2e * zr;
230 term[0][0] = d2edx*xr + de;
231 term[1][0] = d2edx*yr;
232 term[2][0] = d2edx*zr;
233 term[0][1] = term[1][0];
234 term[1][1] = d2edy*yr + de;
235 term[2][1] = d2edy*zr;
236 term[0][2] = term[2][0];
237 term[1][2] = term[2][1];
238 term[2][2] = d2edz*zr + de;
239
240 if (ia == iatom)
241 {
242 for (k=0; k < 3; k++)
243 {
244 hessx[ia][k] += term[k][0];
245 hessy[ia][k] += term[k][1];
246 hessz[ia][k] += term[k][2];
247 hessx[ib][k] -= term[k][0];
248 hessy[ib][k] -= term[k][1];
249 hessz[ib][k] -= term[k][2];
250 }
251 } else
252 {
253 for (k=0; k < 3; k++)
254 {
255 hessx[ia][k] -= term[k][0];
256 hessy[ia][k] -= term[k][1];
257 hessz[ia][k] -= term[k][2];
258 hessx[ib][k] += term[k][0];
259 hessy[ib][k] += term[k][1];
260 hessz[ib][k] += term[k][2];
261 }
262 }
263 }
264 }
265 }
266 }