ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/ebond.c
Revision: 27
Committed: Tue Jul 8 19:10:25 2008 UTC (11 years, 5 months ago) by tjod
File size: 12371 byte(s)
Log Message:
move mengine src to reflect new freemol directory structure

Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5
6 #include "energies.h"
7 #include "bonds_ff.h"
8 #include "derivs.h"
9 #include "hess.h"
10 #include "field.h"
11 #include "fix.h"
12
13 EXTERN struct t_units {
14 double bndunit, cbnd, qbnd;
15 double angunit, cang, qang, pang, sang, aaunit;
16 double stbnunit, ureyunit, torsunit, storunit, v14scale;
17 double aterm, bterm, cterm, dielec, chgscale;
18 } units;
19
20 EXTERN struct t_minim_values {
21 int iprint, ndc, nconst;
22 float dielc;
23 } minim_values;
24
25 int find_bond(int,int);
26 int find_fixed_bond(int,int);
27 void bnd_corr(int,int, int,int,double *);
28 double dihdrl(int,int,int,int);
29 double deltaks(double,double);
30
31 void ebond()
32 {
33 /* compute stretching energy */
34
35 int i, it, kt, itype, ktype, classi,classk;
36 double xr, yr, zr, rik, rik2, bcorr;
37 double dt, dt2, e;
38
39 energies.estr = 0.0F;
40
41 if (minim_values.iprint)
42 {
43 fprintf(pcmoutfile,"\nBond Terms \n");
44 fprintf(pcmoutfile," At1 At2 R BLen Bconst Eb\n");
45 }
46
47 for (i=0; i < bonds_ff.nbnd; i++)
48 {
49 it = bonds_ff.i12[i][0];
50 kt = bonds_ff.i12[i][1];
51 itype = atom[it].type;
52 ktype = atom[kt].type;
53 classi = atom[it].tclass;
54 classk = atom[kt].tclass;
55 bcorr = 0.0;
56 if ( atom[it].use || atom[kt].use )
57 {
58 xr = atom[it].x - atom[kt].x;
59 yr = atom[it].y - atom[kt].y;
60 zr = atom[it].z - atom[kt].z;
61 rik2 = xr*xr + yr*yr + zr*zr;
62 rik = sqrt(rik2);
63 if (field.type == MM3 && ( (itype ==1 && ktype == 6) || (itype == 6 && ktype ==1) ||
64 (classi == 1 && classk == 6) || (classi == 6 && classk == 1) ) )
65 bnd_corr(it,kt, itype, ktype, &bcorr);
66
67 dt = rik - (bonds_ff.bl[i] - bcorr);
68 dt2 = dt*dt;
69 e = units.bndunit *bonds_ff.bk[i]*dt2*
70 (1.0 + units.cbnd*dt + units.qbnd*dt2);
71 energies.estr += e;
72 atom[it].energy += e;
73 atom[kt].energy += e;
74 if (minim_values.iprint)
75 fprintf(pcmoutfile,"Bond: %2s(%-3d) - %2s(%-3d) %-8.3f %-8.3f %-8.3f = %-8.4f\n",atom[it].name,it,atom[kt].name
76 ,kt, rik, bonds_ff.bl[i]-bcorr, bonds_ff.bk[i],e);
77 }
78 }
79 // fixed distances
80 for (i=0; i < fixdis.nfxstr; i++)
81 {
82 it = fixdis.ifxstr[i][0];
83 kt = fixdis.ifxstr[i][1];
84 itype = atom[it].type;
85 ktype = atom[kt].type;
86 bcorr = 0.0;
87 if ( atom[it].use || atom[kt].use )
88 {
89 xr = atom[it].x - atom[kt].x;
90 yr = atom[it].y - atom[kt].y;
91 zr = atom[it].z - atom[kt].z;
92 rik2 = xr*xr + yr*yr + zr*zr;
93 rik = sqrt(rik2);
94 dt = rik - fixdis.fxstrd[i];
95
96 dt2 = dt*dt;
97 e = units.bndunit *fixdis.fxstrc[i]*dt2;
98
99 energies.estr += e;
100 atom[it].energy += e;
101 atom[kt].energy += e;
102 if (minim_values.iprint)
103 fprintf(pcmoutfile,"Fixed Dist: %2s(%-3d) - %2s(%-3d) %-8.3f %-8.3f %-8.3f = %-8.4f\n",atom[it].name,it,atom[kt].name
104 ,kt, rik, fixdis.fxstrd[i], fixdis.fxstrc[i],e);
105 }
106 }
107 }
108
109 void ebond1()
110 {
111 /* compute stretching energy and first derivatives */
112 int i, it, kt, itype, ktype, classi,classk;
113 double xr, yr, zr, rik, rik2, bcorr;
114 double dt, dt2, e, deddt;
115 double de,dedx,dedy,dedz;
116
117 energies.estr = 0.0F;
118 for (i=0; i <= natom; i++)
119 {
120 deriv.deb[i][0] = 0.0;
121 deriv.deb[i][1] = 0.0;
122 deriv.deb[i][2] = 0.0;
123 }
124
125 for (i=0; i < bonds_ff.nbnd; i++)
126 {
127 it = bonds_ff.i12[i][0];
128 kt = bonds_ff.i12[i][1];
129 itype = atom[it].type;
130 ktype = atom[kt].type;
131 classi = atom[it].tclass;
132 classk = atom[kt].tclass;
133 bcorr = 0.0;
134 if ( atom[it].use || atom[kt].use )
135 {
136 xr = atom[it].x - atom[kt].x;
137 yr = atom[it].y - atom[kt].y;
138 zr = atom[it].z - atom[kt].z;
139 rik2 = xr*xr + yr*yr + zr*zr;
140 rik = sqrt(rik2);
141
142 if (field.type == MM3 && ( (itype ==1 && ktype == 6) || (itype == 6 && ktype ==1) ||
143 (classi == 1 && classk == 6) || (classi == 6 && classk == 1) ) )
144 bnd_corr(it,kt, itype, ktype, &bcorr);
145 dt = rik - (bonds_ff.bl[i] - bcorr);
146
147 dt2 = dt*dt;
148 e = units.bndunit *bonds_ff.bk[i]*dt2*
149 (1.0 + units.cbnd*dt + units.qbnd*dt2);
150 deddt = 2.0 * units.bndunit * bonds_ff.bk[i] * dt
151 * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
152 if (rik == 0.0)
153 de = 0.0;
154 else
155 de = deddt/rik;
156 dedx = de*xr;
157 dedy = de*yr;
158 dedz = de*zr;
159 energies.estr += e;
160 deriv.deb[it][0] += dedx;
161 deriv.deb[it][1] += dedy;
162 deriv.deb[it][2] += dedz;
163 deriv.deb[kt][0] -= dedx;
164 deriv.deb[kt][1] -= dedy;
165 deriv.deb[kt][2] -= dedz;
166
167 virial.virx += xr*dedx;
168 virial.viry += yr*dedy;
169 virial.virz += zr*dedz;
170 }
171 }
172 // fixed distances
173 for (i=0; i < fixdis.nfxstr; i++)
174 {
175 it = fixdis.ifxstr[i][0];
176 kt = fixdis.ifxstr[i][1];
177 if ( atom[it].use || atom[kt].use )
178 {
179 xr = atom[it].x - atom[kt].x;
180 yr = atom[it].y - atom[kt].y;
181 zr = atom[it].z - atom[kt].z;
182 rik2 = xr*xr + yr*yr + zr*zr;
183 rik = sqrt(rik2);
184 dt = rik - fixdis.fxstrd[i];
185
186 dt2 = dt*dt;
187 e = units.bndunit *fixdis.fxstrc[i]*dt2;
188 deddt = 2.0 * units.bndunit * fixdis.fxstrc[i] * dt;
189 if (rik == 0.0)
190 de = 0.0;
191 else
192 de = deddt/rik;
193 dedx = de*xr;
194 dedy = de*yr;
195 dedz = de*zr;
196 energies.estr += e;
197 deriv.deb[it][0] += dedx;
198 deriv.deb[it][1] += dedy;
199 deriv.deb[it][2] += dedz;
200 deriv.deb[kt][0] -= dedx;
201 deriv.deb[kt][1] -= dedy;
202 deriv.deb[kt][2] -= dedz;
203
204 virial.virx += xr*dedx;
205 virial.viry += yr*dedy;
206 virial.virz += zr*dedz;
207 }
208 }
209 }
210
211
212 void ebond2(int ia)
213 {
214 int j,k,m,ibond;
215 double dt,dt2;
216 double xr,yr,zr,rik,rik2,deddt,d2eddt2;
217 double de,term,termx,termy,termz,d2e[3][3];
218
219 for (m=0; m < MAXIAT; m++)
220 {
221 if (atom[ia].iat[m]!= 0 && atom[ia].bo[m] != 9)
222 {
223 ibond = find_bond(ia, atom[ia].iat[m]);
224 if (bonds_ff.i12[ibond][0] == ia)
225 k = bonds_ff.i12[ibond][1];
226 else
227 k = bonds_ff.i12[ibond][0];
228
229 xr = atom[ia].x - atom[k].x;
230 yr = atom[ia].y - atom[k].y;
231 zr = atom[ia].z - atom[k].z;
232 rik2 = xr*xr + yr*yr + zr*zr;
233 rik = sqrt(rik2);
234 dt = rik - bonds_ff.bl[ibond];
235 dt2 = dt * dt;
236
237 deddt = 2.0 * units.bndunit * bonds_ff.bk[ibond] * dt
238 * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
239 d2eddt2 = 2.0 * units.bndunit * bonds_ff.bk[ibond]
240 * (1.0+3.0*units.cbnd*dt+6.0*units.qbnd*dt2);
241
242 if (rik2 == 0.0)
243 {
244 de = 0.0;
245 term = 0.0;
246 }else
247 {
248 de = deddt / rik;
249 term = (d2eddt2-de) / rik2;
250 }
251
252 termx = term * xr;
253 termy = term * yr;
254 termz = term * zr;
255 d2e[0][0] = termx*xr + de;
256 d2e[1][0] = termx*yr;
257 d2e[2][0] = termx*zr;
258 d2e[0][1] = d2e[1][0];
259 d2e[1][1] = termy*yr + de;
260 d2e[2][1] = termy*zr;
261 d2e[0][2] = d2e[2][0];
262 d2e[1][2] = d2e[2][1];
263 d2e[2][2] = termz*zr + de;
264
265 for (j=0; j < 3; j++)
266 {
267 hess.hessx[ia][j] += d2e[j][0];
268 hess.hessy[ia][j] += d2e[j][1];
269 hess.hessz[ia][j] += d2e[j][2];
270 hess.hessx[k][j] -= d2e[j][0];
271 hess.hessy[k][j] -= d2e[j][1];
272 hess.hessz[k][j] -= d2e[j][2];
273 }
274 }
275 }
276 // fixed distances
277 for (m=0; m < MAXIAT; m++)
278 {
279 if (atom[ia].iat[m]!= 0 && atom[ia].bo[m] != 9)
280 {
281 ibond = find_fixed_bond(ia, atom[ia].iat[m]);
282 if (ibond > -1)
283 {
284 if (fixdis.ifxstr[ibond][0] == ia)
285 k = fixdis.ifxstr[ibond][1];
286 else
287 k = fixdis.ifxstr[ibond][0];
288
289 xr = atom[ia].x - atom[k].x;
290 yr = atom[ia].y - atom[k].y;
291 zr = atom[ia].z - atom[k].z;
292 rik2 = xr*xr + yr*yr + zr*zr;
293 rik = sqrt(rik2);
294 dt = rik - fixdis.fxstrd[ibond];
295 dt2 = dt * dt;
296
297 deddt = 2.0 * units.bndunit * fixdis.fxstrc[ibond] * dt;
298 d2eddt2 = 2.0 * units.bndunit * fixdis.fxstrc[ibond];
299
300 if (rik2 == 0.0)
301 {
302 de = 0.0;
303 term = 0.0;
304 }else
305 {
306 de = deddt / rik;
307 term = (d2eddt2-de) / rik2;
308 }
309
310 termx = term * xr;
311 termy = term * yr;
312 termz = term * zr;
313 d2e[0][0] = termx*xr + de;
314 d2e[1][0] = termx*yr;
315 d2e[2][0] = termx*zr;
316 d2e[0][1] = d2e[1][0];
317 d2e[1][1] = termy*yr + de;
318 d2e[2][1] = termy*zr;
319 d2e[0][2] = d2e[2][0];
320 d2e[1][2] = d2e[2][1];
321 d2e[2][2] = termz*zr + de;
322
323 for (j=0; j < 3; j++)
324 {
325 hess.hessx[ia][j] += d2e[j][0];
326 hess.hessy[ia][j] += d2e[j][1];
327 hess.hessz[ia][j] += d2e[j][2];
328 hess.hessx[k][j] -= d2e[j][0];
329 hess.hessy[k][j] -= d2e[j][1];
330 hess.hessz[k][j] -= d2e[j][2];
331 }
332 }
333 }
334 }
335
336 }
337
338 /* ------------------------------------------ */
339 int find_fixed_bond(int ia, int ib)
340 {
341 int i;
342 int iz;
343
344 iz = -1;
345 if (fixdis.nfxstr <= 0)
346 return (iz);
347 for (i=0; i < fixdis.nfxstr; i++)
348 {
349 if (ia == fixdis.ifxstr[i][0] && ib == fixdis.ifxstr[i][1])
350 return(i);
351 if (ib == fixdis.ifxstr[i][0] && ia == fixdis.ifxstr[i][1])
352 return(i);
353 }
354 return(iz);
355 }
356 /* ------------------------------------------- */
357 void bnd_corr(int ia, int ib, int itype, int ktype, double *bcorr)
358 {
359 // two terms - first O-ia-ib-R
360 // then X-O-ia-ib
361 // it = carbon, kt = oxygen
362
363 int i, j,k, iatt, jatt, latt,iatype, jatype, iox;
364 int it, kt;
365 double angle;
366
367 *bcorr = 0.0;
368
369 if (itype == 1)
370 {
371 it = ia;
372 kt = ib;
373 } else
374 {
375 it = ib;
376 kt = ia;
377 }
378
379 for (i=0; i < MAXIAT; i++)
380 {
381 if (atom[kt].iat[i] != 0 && atom[kt].iat[i] != it)
382 {
383 iatt = atom[kt].iat[i];
384 iatype = atom[iatt].type;
385 for(j=0; j < MAXIAT; j++)
386 {
387 iox = 0;
388 if (atom[it].iat[j] != 0 && atom[it].iat[j] != kt)
389 {
390 jatt = atom[it].iat[j];
391 jatype = atom[jatt].type;
392
393 if (itype == 6 && jatype == 6) iox = 1;
394 if (ktype == 6 && iatype == 6) iox = 2;
395 if (iox != 0)
396 {
397 angle = dihdrl(iatt,kt,it,jatt)/radian;
398 *bcorr += 0.5* 0.01*(1.0-cos(2.0*angle));
399 if (iox == 1)
400 {
401 for (k = 0 ; k < MAXIAT; k++)
402 {
403 if (atom[jatt].iat[k] != 0 && atom[jatt].iat[k] != it)
404 {
405 latt = atom[jatt].iat[k];
406 angle = dihdrl(kt,it,jatt,latt)/radian;
407 *bcorr += 0.5*0.01*(1.0 - cos(2.0*angle))*(-1.30)+0.005;
408 break;
409 }
410 }
411 }
412 }
413 }
414 }
415 }
416 }
417 }