ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/eangle.c
Revision: 85
Committed: Mon Jan 12 15:05:44 2009 UTC (10 years, 9 months ago) by gilbertke
File size: 41041 byte(s)
Log Message:
more cleanup and code removal
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5
6 #include "energies.h"
7 #include "angles.h"
8 #include "derivs.h"
9 #include "hess.h"
10 #include "utility.h"
11 #include "field.h"
12
13 void eangle2a(int);
14 void eangle2b(int);
15 void eangle3(int, double (*)[3]);
16 //double acos(double);
17
18
19 EXTERN struct t_minim_values {
20 int iprint, ndc, nconst;
21 float dielc;
22 } minim_values;
23
24 void eangle()
25 {
26 int i, ia, ib, ic, id;
27 double e, angle, dot, cosine, factor;
28 double dt, dt2, dt3, dt4;
29 double xia, yia, zia;
30 double xib, yib, zib;
31 double xic, yic, zic;
32 double xid, yid, zid;
33 double xab, yab, zab, rab2;
34 double xcb, ycb, zcb, rcb2;
35 double xad, yad, zad;
36 double xbd, ybd, zbd;
37 double xcd,ycd,zcd;
38 double xip,yip,zip;
39 double xap,yap,zap,rap2;
40 double xcp,ycp,zcp,rcp2;
41 double xt,yt,zt,rt2,delta;
42
43
44
45 if (minim_values.iprint)
46 {
47 fprintf(pcmlogfile,"\nAngle Terms \n");
48 fprintf(pcmlogfile," At1 At2 At3 Angle Thet0 Tconst Ebend\n");
49 }
50
51 for (i=0; i < angles.nang; i++)
52 {
53 ia = angles.i13[i][0];
54 ib = angles.i13[i][1];
55 ic = angles.i13[i][2];
56 id = angles.i13[i][3];
57 if (atom[ia].use || atom[ib].use || atom[ic].use)
58 {
59 xia = atom[ia].x;
60 yia = atom[ia].y;
61 zia = atom[ia].z;
62 xib = atom[ib].x;
63 yib = atom[ib].y;
64 zib = atom[ib].z;
65 xic = atom[ic].x;
66 yic = atom[ic].y;
67 zic = atom[ic].z;
68 if ( (angles.angin[i] == FALSE) || field.type == MMFF94 )
69 {
70 xab = xia - xib;
71 yab = yia - yib;
72 zab = zia - zib;
73 rab2 = xab*xab + yab*yab + zab*zab;
74 xcb = xic - xib;
75 ycb = yic - yib;
76 zcb = zic - zib;
77 rcb2 = xcb*xcb + ycb*ycb + zcb*zcb;
78 if (rab2 != 0.0 && rcb2 != 0.00)
79 {
80 dot = xab*xcb + yab*ycb + zab*zcb;
81 cosine = dot / sqrt(rab2*rcb2);
82 if (cosine > 1.0)
83 cosine = 1.0;
84 if (cosine < -1.0)
85 cosine = -1.0;
86 angle = radian*acos(cosine);
87
88 if (angles.angtype[i] == HARMONIC)
89 {
90 dt = angle - angles.anat[i];
91 dt2 = dt*dt;
92 dt3 = dt2*dt;
93 dt4 = dt2*dt2;
94 e = units.angunit*angles.acon[i]*dt2
95 *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
96 } else
97 {
98 factor = angle/radian;
99 e = units.angunit*angles.fcon[i]*
100 (angles.c0[i] + angles.c1[i]*cosine + angles.c2[i]*cos(2*factor)
101 + angles.c3[i]*cos(3*factor) + angles.c4[i]*cos(4*factor)
102 + angles.c5[i]*cos(5*factor) + angles.c6[i]*cos(6*factor));
103 }
104 energies.ebend += e;
105 atom[ia].energy += e;
106 atom[ib].energy += e;
107 atom[ic].energy += e;
108 if (minim_values.iprint)
109 {
110 if (angles.angtype[i] == HARMONIC)
111 fprintf(pcmlogfile,"Angle 1-3: %2s(%-3d)- %2s(%-3d)- %2s(%-3d) %-8.3f %-8.3f %-8.4f = %-8.4f\n",
112 atom[ia].name,ia,atom[ib].name, ib,atom[ic].name, ic, angle, angles.anat[i],angles.acon[i], e);
113 else
114 fprintf(pcmlogfile,"Angle 1-3: %2s(%-3d)- %2s(%-3d)- %2s(%-3d) %-8.3f %-8.4f = %-8.4f\n",
115 atom[ia].name,ia,atom[ib].name, ib,atom[ic].name, ic, angle, angles.fcon[i], e);
116 }
117 }
118 } else
119 {
120 xid = atom[id].x;
121 yid = atom[id].y;
122 zid = atom[id].z;
123 xad = xia - xid;
124 yad = yia - yid;
125 zad = zia - zid;
126 xbd = xib - xid;
127 ybd = yib - yid;
128 zbd = zib - zid;
129 xcd = xic - xid;
130 ycd = yic - yid;
131 zcd = zic - zid;
132 xt = yad*zcd - zad*ycd;
133 yt = zad*xcd - xad*zcd;
134 zt = xad*ycd - yad*xcd;
135 rt2 = xt*xt + yt*yt + zt*zt;
136 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
137 xip = xib + xt*delta;
138 yip = yib + yt*delta;
139 zip = zib + zt*delta;
140 xap = xia - xip;
141 yap = yia - yip;
142 zap = zia - zip;
143 rap2 = xap*xap + yap*yap + zap*zap;
144 xcp = xic - xip;
145 ycp = yic - yip;
146 zcp = zic - zip;
147 rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
148 if (rap2 != 0.0 && rcp2 != 0.0)
149 {
150 dot = xap*xcp + yap*ycp + zap*zcp;
151 cosine = dot/ sqrt(rap2*rcp2);
152 if (cosine < -1.0)
153 cosine = -1.0;
154 if (cosine > 1.0)
155 cosine = 1.0;
156 angle = radian*acos(cosine);
157 dt = angle - angles.anat[i];
158 dt2 = dt*dt;
159 dt3 = dt2*dt;
160 dt4 = dt2*dt2;
161 e = units.angunit*angles.acon[i]*dt2
162 *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
163 energies.ebend += e;
164 atom[ia].energy += e;
165 atom[ib].energy += e;
166 atom[ic].energy += e;
167 if (minim_values.iprint)
168 fprintf(pcmlogfile,"Angle InP: %2s(%-3d)- %2s(%-3d)- %2s(%-3d) %-8.3f %-8.3f %-8.4f = %-8.4f\n",
169 atom[ia].name, ia, atom[ib].name, ib, atom[ic].name, ic, angle, angles.anat[i], angles.acon[i], e);
170
171 }
172 }
173 }
174 }
175
176 }
177 // ===============================================
178 void eangle1()
179 {
180 int i, ia, ib, ic, id;
181 double temp1, temp2;
182 double e, angle, dot, cosine, factor, sine;
183 double dt, dt2, dt3, dt4;
184 double deddt,terma,termc,term;
185 double xia, yia, zia;
186 double xib, yib, zib;
187 double xic, yic, zic;
188 double xid, yid, zid;
189 double xab, yab, zab, rab2;
190 double xcb, ycb, zcb, rcb2;
191 double xp,yp,zp,rp;
192 double xad, yad, zad;
193 double xbd, ybd, zbd;
194 double xcd,ycd,zcd;
195 double xip,yip,zip;
196 double xap,yap,zap,rap2;
197 double xcp,ycp,zcp,rcp2;
198 double xt,yt,zt,rt2,ptrt2;
199 double xm,ym,zm,rm,delta,delta2;
200 double dedxia,dedyia,dedzia;
201 double dedxib,dedyib,dedzib;
202 double dedxic,dedyic,dedzic;
203 double dedxid,dedyid,dedzid;
204 double dedxip,dedyip,dedzip;
205 double dpdxia,dpdyia,dpdzia;
206 double dpdxic,dpdyic,dpdzic;
207 double dtemp;
208
209 energies.ebend = 0.0F;
210 for (i=0; i <= natom; i++)
211 {
212 deriv.dea[i][0] = 0.0;
213 deriv.dea[i][1] = 0.0;
214 deriv.dea[i][2] = 0.0;
215 }
216 for (i=0; i < angles.nang; i++)
217 {
218 ia = angles.i13[i][0];
219 ib = angles.i13[i][1];
220 ic = angles.i13[i][2];
221 id = angles.i13[i][3];
222 if (atom[ia].use || atom[ib].use || atom[ic].use)
223 {
224 xia = atom[ia].x;
225 yia = atom[ia].y;
226 zia = atom[ia].z;
227 xib = atom[ib].x;
228 yib = atom[ib].y;
229 zib = atom[ib].z;
230 xic = atom[ic].x;
231 yic = atom[ic].y;
232 zic = atom[ic].z;
233 if ( (angles.angin[i] == FALSE) || field.type == MMFF94 )
234 {
235 xab = xia - xib;
236 yab = yia - yib;
237 zab = zia - zib;
238 rab2 = xab*xab + yab*yab + zab*zab;
239 xcb = xic - xib;
240 ycb = yic - yib;
241 zcb = zic - zib;
242 rcb2 = xcb*xcb + ycb*ycb + zcb*zcb;
243 xp = ycb*zab - zcb*yab;
244 yp = zcb*xab - xcb*zab;
245 zp = xcb*yab - ycb*xab;
246 rp = sqrt(xp*xp + yp*yp + zp*zp);
247 if (rp != 0.0 )
248 {
249 dot = xab*xcb + yab*ycb + zab*zcb;
250 cosine = dot / sqrt(rab2*rcb2);
251 if (cosine > 1.0)
252 cosine = 1.0;
253 if (cosine < -1.0)
254 cosine = -1.0;
255 angle = radian*acos(cosine);
256 if (angles.angtype[i] == HARMONIC)
257 {
258 dt = angle - angles.anat[i];
259 dt2 = dt*dt;
260 dt3 = dt2*dt;
261 dt4 = dt2*dt2;
262
263 e = units.angunit*angles.acon[i]*dt2
264 *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
265 deddt = units.angunit*angles.acon[i]* dt * radian
266 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
267 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
268 } else
269 {
270 factor = angle/radian;
271 sine = sin(factor);
272 e = units.angunit*angles.fcon[i]*
273 (angles.c0[i] + angles.c1[i]*cosine + angles.c2[i]*cos(2*factor)
274 + angles.c3[i]*cos(3*factor) + angles.c4[i]*cos(4*factor)
275 + angles.c5[i]*cos(5*factor) + angles.c6[i]*cos(6*factor));
276 deddt = -units.angunit*angles.fcon[i]*
277 (angles.c1[i]*sine + 2*angles.c2[i]*sin(2*factor)
278 + 3*angles.c3[i]*sin(3*factor) + 4*angles.c4[i]*sin(4*factor)
279 + 5*angles.c5[i]*sin(5*factor) + 6*angles.c6[i]*sin(6*factor));
280 }
281
282 /* terma = -deddt / (rab2*rp);
283 termc = deddt / (rcb2*rp);
284 dedxia = terma * (yab*zp-zab*yp);
285 dedyia = terma * (zab*xp-xab*zp);
286 dedzia = terma * (xab*yp-yab*xp);
287 dedxic = termc * (ycb*zp-zcb*yp);
288 dedyic = termc * (zcb*xp-xcb*zp);
289 dedzic = termc * (xcb*yp-ycb*xp);
290 dedxib = -dedxia - dedxic;
291 dedyib = -dedyia - dedyic;
292 dedzib = -dedzia - dedzic; */
293
294 dtemp = dot*dot/(rab2*rcb2);
295 if (fabs(dtemp) >= 1.0)
296 dtemp = 0.9999*dtemp/(fabs(dtemp));
297 terma = sqrt(1.0-dtemp);
298 termc = sqrt(rab2*rcb2);
299 dedxia = -deddt*( xcb/termc - (xab*cosine)/rab2)/terma;
300 dedyia = -deddt*( ycb/termc - (yab*cosine)/rab2)/terma;
301 dedzia = -deddt*( zcb/termc - (zab*cosine)/rab2)/terma;
302
303 dedxic = -deddt*( xab/termc - (xcb*cosine)/rcb2)/terma;
304 dedyic = -deddt*( yab/termc - (ycb*cosine)/rcb2)/terma;
305 dedzic = -deddt*( zab/termc - (zcb*cosine)/rcb2)/terma;
306
307 dedxib = -deddt*( (-xab-xcb)/termc + (dot*(xcb*rab2 + xab*rcb2))/(rab2*rcb2*termc))/terma;
308 dedyib = -deddt*( (-yab-ycb)/termc + (dot*(ycb*rab2 + yab*rcb2))/(rab2*rcb2*termc))/terma;
309 dedzib = -deddt*( (-zab-zcb)/termc + (dot*(zcb*rab2 + zab*rcb2))/(rab2*rcb2*termc))/terma;
310
311
312 deriv.dea[ia][0] += dedxia;
313 deriv.dea[ia][1] += dedyia;
314 deriv.dea[ia][2] += dedzia;
315
316 deriv.dea[ib][0] += dedxib;
317 deriv.dea[ib][1] += dedyib;
318 deriv.dea[ib][2] += dedzib;
319
320 deriv.dea[ic][0] += dedxic;
321 deriv.dea[ic][1] += dedyic;
322 deriv.dea[ic][2] += dedzic;
323
324 virial.virx += xab*dedxia + xcb*dedxic;
325 virial.viry += yab*dedyia + ycb*dedyic;
326 virial.virz += zab*dedzia + zcb*dedzic;
327
328 energies.ebend += e;
329 }
330 } else
331 {
332 xid = atom[id].x;
333 yid = atom[id].y;
334 zid = atom[id].z;
335 xad = xia - xid;
336 yad = yia - yid;
337 zad = zia - zid;
338 xbd = xib - xid;
339 ybd = yib - yid;
340 zbd = zib - zid;
341 xcd = xic - xid;
342 ycd = yic - yid;
343 zcd = zic - zid;
344 xt = yad*zcd - zad*ycd;
345 yt = zad*xcd - xad*zcd;
346 zt = xad*ycd - yad*xcd;
347 rt2 = xt*xt + yt*yt + zt*zt;
348 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
349 xip = xib + xt*delta;
350 yip = yib + yt*delta;
351 zip = zib + zt*delta;
352 xap = xia - xip;
353 yap = yia - yip;
354 zap = zia - zip;
355 rap2 = xap*xap + yap*yap + zap*zap;
356 xcp = xic - xip;
357 ycp = yic - yip;
358 zcp = zic - zip;
359 rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
360 xm = ycp*zap - zcp*yap;
361 ym = zcp*xap - xcp*zap;
362 zm = xcp*yap - ycp*xap;
363 rm = sqrt(xm*xm + ym*ym + zm*zm);
364 if (rm != 0.0)
365 {
366 dot = xap*xcp + yap*ycp + zap*zcp;
367 cosine = dot/ sqrt(rap2*rcp2);
368 if (cosine < -1.0)
369 cosine = 1.0;
370 if (cosine > 1.0)
371 cosine = 1.0;
372 temp1 = acos(cosine);
373 temp2 = radian;
374 angle = temp1*temp2;
375 dt = angle - angles.anat[i];
376 dt2 = dt*dt;
377 dt3 = dt2*dt;
378 dt4 = dt2*dt2;
379 e = units.angunit*angles.acon[i]*dt2
380 *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
381
382 deddt = units.angunit *angles.acon[i]* dt
383 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
384 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
385
386 deddt = deddt * radian;
387 terma = -deddt / (rap2*rm);
388 termc = deddt / (rcp2*rm);
389 dedxia = terma * (yap*zm-zap*ym);
390 dedyia = terma * (zap*xm-xap*zm);
391 dedzia = terma * (xap*ym-yap*xm);
392 dedxic = termc * (ycp*zm-zcp*ym);
393 dedyic = termc * (zcp*xm-xcp*zm);
394 dedzic = termc * (xcp*ym-ycp*xm);
395 dedxip = -dedxia - dedxic;
396 dedyip = -dedyia - dedyic;
397 dedzip = -dedzia - dedzic;
398
399 delta2 = 2.0 * delta;
400 ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
401 term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
402 dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
403 term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
404 dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
405 term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
406 dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
407 term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
408 dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
409 term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
410 dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
411 term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
412 dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
413
414 dedxia = dedxia + dpdxia;
415 dedyia = dedyia + dpdyia;
416 dedzia = dedzia + dpdzia;
417 dedxib = dedxip;
418 dedyib = dedyip;
419 dedzib = dedzip;
420 dedxic = dedxic + dpdxic;
421 dedyic = dedyic + dpdyic;
422 dedzic = dedzic + dpdzic;
423 dedxid = -dedxia - dedxib - dedxic;
424 dedyid = -dedyia - dedyib - dedyic;
425 dedzid = -dedzia - dedzib - dedzic;
426
427 energies.ebend += e;
428 deriv.dea[ia][0] += dedxia;
429 deriv.dea[ia][1] += dedyia;
430 deriv.dea[ia][2] += dedzia;
431
432 deriv.dea[ib][0] += dedxib;
433 deriv.dea[ib][1] += dedyib;
434 deriv.dea[ib][2] += dedzib;
435
436 deriv.dea[ic][0] += dedxic;
437 deriv.dea[ic][1] += dedyic;
438 deriv.dea[ic][2] += dedzic;
439
440 deriv.dea[id][0] += dedxid;
441 deriv.dea[id][1] += dedyid;
442 deriv.dea[id][2] += dedzid;
443
444 virial.virx += xad*dedxia + xbd*dedxib + xcd*dedxic;
445 virial.viry += yad*dedyia + ybd*dedyib + ycd*dedyic;
446 virial.virz += zad*dedzia + zbd*dedzib + zcd*dedzic;
447 }
448 }
449 }
450 }
451 }
452
453
454 void eangle2(int i)
455 {
456 int j,k,ia,ib,ic,id;
457 double old,eps;
458 double **d0; // d0[MAXATOM][3];
459
460 d0 = dmatrix(0,natom,0,3);
461
462 eangle2a(i);
463
464 eps = 1.0e-7;
465 for (k=0; k < angles.nang; k++)
466 {
467 if ( (angles.angin[k] == TRUE) && (field.type != MMFF94) )
468 {
469 ia = angles.i13[k][0];
470 ib = angles.i13[k][1];
471 ic = angles.i13[k][2];
472 id = angles.i13[k][3];
473 if (i == ia || i == ib || i == ic || i == id)
474 {
475 eangle2b(k);
476 for (j=0; j < 3; j++)
477 {
478 d0[ia][j] = deriv.dea[ia][j];
479 d0[ib][j] = deriv.dea[ib][j];
480 d0[ic][j] = deriv.dea[ic][j];
481 d0[id][j] = deriv.dea[id][j];
482 }
483 // numerical x derivative
484 old = atom[i].x;
485 atom[i].x += eps;
486 eangle2b(k);
487 atom[i].x = old;
488 for (j=0; j < 3; j++)
489 {
490 hess.hessx[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
491 hess.hessx[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
492 hess.hessx[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
493 hess.hessx[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
494 }
495 // numerical y derivative
496 old = atom[i].y;
497 atom[i].y += eps;
498 eangle2b(k);
499 atom[i].y = old;
500 for (j=0; j < 3; j++)
501 {
502 hess.hessy[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
503 hess.hessy[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
504 hess.hessy[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
505 hess.hessy[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
506 }
507 // numerical z derivative
508 old = atom[i].z;
509 atom[i].z += eps;
510 eangle2b(k);
511 atom[i].z = old;
512 for (j=0; j < 3; j++)
513 {
514 hess.hessz[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
515 hess.hessz[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
516 hess.hessz[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
517 hess.hessz[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
518 }
519 }
520 }
521 }
522 free_dmatrix(d0 ,0,natom,0,3);
523 }
524 // =================================================
525 void eangle2a(int iatom)
526 {
527 int i,ia,ib,ic;
528 double angle,dot,cosine,factor,sine;
529 double dt,dt2,dt3,dt4;
530 double xia,yia,zia;
531 double xib,yib,zib;
532 double xic,yic,zic;
533 double xab,yab,zab,rab;
534 double xcb,ycb,zcb,rcb;
535 double xp,yp,zp,rp,rp2;
536 double xrab,yrab,zrab,rab2;
537 double xrcb,yrcb,zrcb,rcb2;
538 double xabp,yabp,zabp;
539 double xcbp,ycbp,zcbp;
540 double deddt,d2eddt2,terma,termc;
541 double ddtdxia,ddtdyia,ddtdzia;
542 double ddtdxib,ddtdyib,ddtdzib;
543 double ddtdxic,ddtdyic,ddtdzic;
544 double dxiaxia,dxiayia,dxiazia;
545 double dxibxib,dxibyib,dxibzib;
546 double dxicxic,dxicyic,dxiczic;
547 double dyiayia,dyiazia,dziazia;
548 double dyibyib,dyibzib,dzibzib;
549 double dyicyic,dyiczic,dziczic;
550 double dxibxia,dxibyia,dxibzia;
551 double dyibxia,dyibyia,dyibzia;
552 double dzibxia,dzibyia,dzibzia;
553 double dxibxic,dxibyic,dxibzic;
554 double dyibxic,dyibyic,dyibzic;
555 double dzibxic,dzibyic,dzibzic;
556 double dxiaxic,dxiayic,dxiazic;
557 double dyiaxic,dyiayic,dyiazic;
558 double dziaxic,dziayic,dziazic;
559
560 for (i=0; i < angles.nang; i++)
561 {
562 ia = angles.i13[i][0];
563 ib = angles.i13[i][1];
564 ic = angles.i13[i][2];
565 if (iatom == ia || iatom == ib || iatom == ic)
566 {
567 xia = atom[ia].x;
568 yia = atom[ia].y;
569 zia = atom[ia].z;
570 xib = atom[ib].x;
571 yib = atom[ib].y;
572 zib = atom[ib].z;
573 xic = atom[ic].x;
574 yic = atom[ic].y;
575 zic = atom[ic].z;
576
577 if ( (angles.angin[i] == FALSE) || field.type == MMFF94 )
578 {
579 xab = xia - xib;
580 yab = yia - yib;
581 zab = zia - zib;
582 rab = sqrt(xab*xab + yab*yab + zab*zab);
583 xcb = xic - xib;
584 ycb = yic - yib;
585 zcb = zic - zib;
586 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
587 xp = ycb*zab - zcb*yab;
588 yp = zcb*xab - xcb*zab;
589 zp = xcb*yab - ycb*xab;
590 rp = sqrt(xp*xp + yp*yp + zp*zp);
591 if (rp != 0.0)
592 {
593 dot = xab*xcb + yab*ycb + zab*zcb;
594 cosine = dot / (rab*rcb);
595 if (cosine < -1.0)
596 cosine = -1.0;
597 if (cosine > 1.0)
598 cosine = 1.0;
599 angle = radian * acos(cosine);
600
601 if (angles.angtype[i] == HARMONIC)
602 {
603 dt = angle - angles.anat[i];
604 dt2 = dt * dt;
605 dt3 = dt2 * dt;
606 dt4 = dt3 * dt;
607 deddt = units.angunit * angles.acon[i] * dt
608 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
609 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
610 d2eddt2 = units.angunit * angles.acon[i]
611 * (2.0 + 6.0*units.cang*dt + 12.0*units.qang*dt2
612 + 20.0*units.pang*dt3 + 30.0*units.sang*dt4);
613 terma = -radian/ (rab*rab*rp);
614 termc = radian / (rcb*rcb*rp);
615 } else
616 {
617 factor = angle/radian;
618 sine = sin(factor);
619 deddt = -units.angunit*angles.fcon[i]*
620 (angles.c1[i]*sine + 2*angles.c2[i]*sin(2*factor)
621 + 3*angles.c3[i]*sin(3*factor) + 4*angles.c4[i]*sin(4*factor)
622 + 5*angles.c5[i]*sin(5*factor) + 6*angles.c6[i]*sin(6*factor));
623 d2eddt2 = -units.angunit*angles.fcon[i]*
624 (angles.c1[i]*cosine + 4*angles.c2[i]*cos(2*factor)
625 + 9*angles.c3[i]*cos(3*factor) + 16*angles.c4[i]*cos(4*factor)
626 + 25*angles.c5[i]*cos(5*factor) + 36*angles.c6[i]*cos(6*factor));
627 terma = -1.0 / (rab*rab*rp);
628 termc = 1.0 / (rcb*rcb*rp);
629 }
630 ddtdxia = terma * (yab*zp-zab*yp);
631 ddtdyia = terma * (zab*xp-xab*zp);
632 ddtdzia = terma * (xab*yp-yab*xp);
633 ddtdxic = termc * (ycb*zp-zcb*yp);
634 ddtdyic = termc * (zcb*xp-xcb*zp);
635 ddtdzic = termc * (xcb*yp-ycb*xp);
636 ddtdxib = -ddtdxia - ddtdxic;
637 ddtdyib = -ddtdyia - ddtdyic;
638 ddtdzib = -ddtdzia - ddtdzic;
639 rab2 = 2.0 / (rab*rab);
640 xrab = xab * rab2;
641 yrab = yab * rab2;
642 zrab = zab * rab2;
643 rcb2 = 2.0 / (rcb*rcb);
644 xrcb = xcb * rcb2;
645 yrcb = ycb * rcb2;
646 zrcb = zcb * rcb2;
647 rp2 = 1.0 / (rp*rp);
648 xabp = (yab*zp-zab*yp) * rp2;
649 yabp = (zab*xp-xab*zp) * rp2;
650 zabp = (xab*yp-yab*xp) * rp2;
651 xcbp = (ycb*zp-zcb*yp) * rp2;
652 ycbp = (zcb*xp-xcb*zp) * rp2;
653 zcbp = (xcb*yp-ycb*xp) * rp2;
654
655 dxiaxia = terma*(xab*xcb-dot) + ddtdxia*(xcbp-xrab);
656 dxiayia = terma*(zp+yab*xcb) + ddtdxia*(ycbp-yrab);
657 dxiazia = terma*(zab*xcb-yp) + ddtdxia*(zcbp-zrab);
658 dyiayia = terma*(yab*ycb-dot) + ddtdyia*(ycbp-yrab);
659 dyiazia = terma*(xp+zab*ycb) + ddtdyia*(zcbp-zrab);
660 dziazia = terma*(zab*zcb-dot) + ddtdzia*(zcbp-zrab);
661 dxicxic = termc*(dot-xab*xcb) - ddtdxic*(xabp+xrcb);
662 dxicyic = termc*(zp-ycb*xab) - ddtdxic*(yabp+yrcb);
663 dxiczic = -termc*(yp+zcb*xab) - ddtdxic*(zabp+zrcb);
664 dyicyic = termc*(dot-yab*ycb) - ddtdyic*(yabp+yrcb);
665 dyiczic = termc*(xp-zcb*yab) - ddtdyic*(zabp+zrcb);
666 dziczic = termc*(dot-zab*zcb) - ddtdzic*(zabp+zrcb);
667 dxiaxic = terma*(yab*yab+zab*zab) - ddtdxia*xabp;
668 dxiayic = -terma*xab*yab - ddtdxia*yabp;
669 dxiazic = -terma*xab*zab - ddtdxia*zabp;
670 dyiaxic = -terma*xab*yab - ddtdyia*xabp;
671 dyiayic = terma*(xab*xab+zab*zab) - ddtdyia*yabp;
672 dyiazic = -terma*yab*zab - ddtdyia*zabp;
673 dziaxic = -terma*xab*zab - ddtdzia*xabp;
674 dziayic = -terma*yab*zab - ddtdzia*yabp;
675 dziazic = terma*(xab*xab+yab*yab) - ddtdzia*zabp;
676
677 dxibxia = -dxiaxia - dxiaxic;
678 dxibyia = -dxiayia - dyiaxic;
679 dxibzia = -dxiazia - dziaxic;
680 dyibxia = -dxiayia - dxiayic;
681 dyibyia = -dyiayia - dyiayic;
682 dyibzia = -dyiazia - dziayic;
683 dzibxia = -dxiazia - dxiazic;
684 dzibyia = -dyiazia - dyiazic;
685 dzibzia = -dziazia - dziazic;
686 dxibxic = -dxicxic - dxiaxic;
687 dxibyic = -dxicyic - dxiayic;
688 dxibzic = -dxiczic - dxiazic;
689 dyibxic = -dxicyic - dyiaxic;
690 dyibyic = -dyicyic - dyiayic;
691 dyibzic = -dyiczic - dyiazic;
692 dzibxic = -dxiczic - dziaxic;
693 dzibyic = -dyiczic - dziayic;
694 dzibzic = -dziczic - dziazic;
695 dxibxib = -dxibxia - dxibxic;
696 dxibyib = -dxibyia - dxibyic;
697 dxibzib = -dxibzia - dxibzic;
698 dyibyib = -dyibyia - dyibyic;
699 dyibzib = -dyibzia - dyibzic;
700 dzibzib = -dzibzia - dzibzic;
701
702 if (ia == iatom)
703 {
704 hess.hessx[ia][0] += deddt*dxiaxia + d2eddt2*ddtdxia*ddtdxia;
705 hess.hessx[ia][1] += deddt*dxiayia + d2eddt2*ddtdxia*ddtdyia;
706 hess.hessx[ia][2] += deddt*dxiazia + d2eddt2*ddtdxia*ddtdzia;
707 hess.hessy[ia][0] += deddt*dxiayia + d2eddt2*ddtdyia*ddtdxia;
708 hess.hessy[ia][1] += deddt*dyiayia + d2eddt2*ddtdyia*ddtdyia;
709 hess.hessy[ia][2] += deddt*dyiazia + d2eddt2*ddtdyia*ddtdzia;
710 hess.hessz[ia][0] += deddt*dxiazia + d2eddt2*ddtdzia*ddtdxia;
711 hess.hessz[ia][1] += deddt*dyiazia + d2eddt2*ddtdzia*ddtdyia;
712 hess.hessz[ia][2] += deddt*dziazia + d2eddt2*ddtdzia*ddtdzia;
713 hess.hessx[ib][0] += deddt*dxibxia + d2eddt2*ddtdxia*ddtdxib;
714 hess.hessx[ib][1] += deddt*dyibxia + d2eddt2*ddtdxia*ddtdyib;
715 hess.hessx[ib][2] += deddt*dzibxia + d2eddt2*ddtdxia*ddtdzib;
716 hess.hessy[ib][0] += deddt*dxibyia + d2eddt2*ddtdyia*ddtdxib;
717 hess.hessy[ib][1] += deddt*dyibyia + d2eddt2*ddtdyia*ddtdyib;
718 hess.hessy[ib][2] += deddt*dzibyia + d2eddt2*ddtdyia*ddtdzib;
719 hess.hessz[ib][0] += deddt*dxibzia + d2eddt2*ddtdzia*ddtdxib;
720 hess.hessz[ib][1] += deddt*dyibzia + d2eddt2*ddtdzia*ddtdyib;
721 hess.hessz[ib][2] += deddt*dzibzia + d2eddt2*ddtdzia*ddtdzib;
722 hess.hessx[ic][0] += deddt*dxiaxic + d2eddt2*ddtdxia*ddtdxic;
723 hess.hessx[ic][1] += deddt*dxiayic + d2eddt2*ddtdxia*ddtdyic;
724 hess.hessx[ic][2] += deddt*dxiazic + d2eddt2*ddtdxia*ddtdzic;
725 hess.hessy[ic][0] += deddt*dyiaxic + d2eddt2*ddtdyia*ddtdxic;
726 hess.hessy[ic][1] += deddt*dyiayic + d2eddt2*ddtdyia*ddtdyic;
727 hess.hessy[ic][2] += deddt*dyiazic + d2eddt2*ddtdyia*ddtdzic;
728 hess.hessz[ic][0] += deddt*dziaxic + d2eddt2*ddtdzia*ddtdxic;
729 hess.hessz[ic][1] += deddt*dziayic + d2eddt2*ddtdzia*ddtdyic;
730 hess.hessz[ic][2] += deddt*dziazic + d2eddt2*ddtdzia*ddtdzic;
731 } else if (ib == iatom)
732 {
733 hess.hessx[ib][0] += deddt*dxibxib + d2eddt2*ddtdxib*ddtdxib;
734 hess.hessx[ib][1] += deddt*dxibyib + d2eddt2*ddtdxib*ddtdyib;
735 hess.hessx[ib][2] += deddt*dxibzib + d2eddt2*ddtdxib*ddtdzib;
736 hess.hessy[ib][0] += deddt*dxibyib + d2eddt2*ddtdyib*ddtdxib;
737 hess.hessy[ib][1] += deddt*dyibyib + d2eddt2*ddtdyib*ddtdyib;
738 hess.hessy[ib][2] += deddt*dyibzib + d2eddt2*ddtdyib*ddtdzib;
739 hess.hessz[ib][0] += deddt*dxibzib + d2eddt2*ddtdzib*ddtdxib;
740 hess.hessz[ib][1] += deddt*dyibzib + d2eddt2*ddtdzib*ddtdyib;
741 hess.hessz[ib][2] += deddt*dzibzib + d2eddt2*ddtdzib*ddtdzib;
742 hess.hessx[ia][0] += deddt*dxibxia + d2eddt2*ddtdxib*ddtdxia;
743 hess.hessx[ia][1] += deddt*dxibyia + d2eddt2*ddtdxib*ddtdyia;
744 hess.hessx[ia][2] += deddt*dxibzia + d2eddt2*ddtdxib*ddtdzia;
745 hess.hessy[ia][0] += deddt*dyibxia + d2eddt2*ddtdyib*ddtdxia;
746 hess.hessy[ia][1] += deddt*dyibyia + d2eddt2*ddtdyib*ddtdyia;
747 hess.hessy[ia][2] += deddt*dyibzia + d2eddt2*ddtdyib*ddtdzia;
748 hess.hessz[ia][0] += deddt*dzibxia + d2eddt2*ddtdzib*ddtdxia;
749 hess.hessz[ia][1] += deddt*dzibyia + d2eddt2*ddtdzib*ddtdyia;
750 hess.hessz[ia][2] += deddt*dzibzia + d2eddt2*ddtdzib*ddtdzia;
751 hess.hessx[ic][0] += deddt*dxibxic + d2eddt2*ddtdxib*ddtdxic;
752 hess.hessx[ic][1] += deddt*dxibyic + d2eddt2*ddtdxib*ddtdyic;
753 hess.hessx[ic][2] += deddt*dxibzic + d2eddt2*ddtdxib*ddtdzic;
754 hess.hessy[ic][0] += deddt*dyibxic + d2eddt2*ddtdyib*ddtdxic;
755 hess.hessy[ic][1] += deddt*dyibyic + d2eddt2*ddtdyib*ddtdyic;
756 hess.hessy[ic][2] += deddt*dyibzic + d2eddt2*ddtdyib*ddtdzic;
757 hess.hessz[ic][0] += deddt*dzibxic + d2eddt2*ddtdzib*ddtdxic;
758 hess.hessz[ic][1] += deddt*dzibyic + d2eddt2*ddtdzib*ddtdyic;
759 hess.hessz[ic][2] += deddt*dzibzic + d2eddt2*ddtdzib*ddtdzic;
760 }else if (ic == iatom)
761 {
762 hess.hessx[ic][0] += deddt*dxicxic + d2eddt2*ddtdxic*ddtdxic;
763 hess.hessx[ic][1] += deddt*dxicyic + d2eddt2*ddtdxic*ddtdyic;
764 hess.hessx[ic][2] += deddt*dxiczic + d2eddt2*ddtdxic*ddtdzic;
765 hess.hessy[ic][0] += deddt*dxicyic + d2eddt2*ddtdyic*ddtdxic;
766 hess.hessy[ic][1] += deddt*dyicyic + d2eddt2*ddtdyic*ddtdyic;
767 hess.hessy[ic][2] += deddt*dyiczic + d2eddt2*ddtdyic*ddtdzic;
768 hess.hessz[ic][0] += deddt*dxiczic + d2eddt2*ddtdzic*ddtdxic;
769 hess.hessz[ic][1] += deddt*dyiczic + d2eddt2*ddtdzic*ddtdyic;
770 hess.hessz[ic][2] += deddt*dziczic + d2eddt2*ddtdzic*ddtdzic;
771 hess.hessx[ib][0] += deddt*dxibxic + d2eddt2*ddtdxic*ddtdxib;
772 hess.hessx[ib][1] += deddt*dyibxic + d2eddt2*ddtdxic*ddtdyib;
773 hess.hessx[ib][2] += deddt*dzibxic + d2eddt2*ddtdxic*ddtdzib;
774 hess.hessy[ib][0] += deddt*dxibyic + d2eddt2*ddtdyic*ddtdxib;
775 hess.hessy[ib][1] += deddt*dyibyic + d2eddt2*ddtdyic*ddtdyib;
776 hess.hessy[ib][2] += deddt*dzibyic + d2eddt2*ddtdyic*ddtdzib;
777 hess.hessz[ib][0] += deddt*dxibzic + d2eddt2*ddtdzic*ddtdxib;
778 hess.hessz[ib][1] += deddt*dyibzic + d2eddt2*ddtdzic*ddtdyib;
779 hess.hessz[ib][2] += deddt*dzibzic + d2eddt2*ddtdzic*ddtdzib;
780 hess.hessx[ia][0] += deddt*dxiaxic + d2eddt2*ddtdxic*ddtdxia;
781 hess.hessx[ia][1] += deddt*dyiaxic + d2eddt2*ddtdxic*ddtdyia;
782 hess.hessx[ia][2] += deddt*dziaxic + d2eddt2*ddtdxic*ddtdzia;
783 hess.hessy[ia][0] += deddt*dxiayic + d2eddt2*ddtdyic*ddtdxia;
784 hess.hessy[ia][1] += deddt*dyiayic + d2eddt2*ddtdyic*ddtdyia;
785 hess.hessy[ia][2] += deddt*dziayic + d2eddt2*ddtdyic*ddtdzia;
786 hess.hessz[ia][0] += deddt*dxiazic + d2eddt2*ddtdzic*ddtdxia;
787 hess.hessz[ia][1] += deddt*dyiazic + d2eddt2*ddtdzic*ddtdyia;
788 hess.hessz[ia][2] += deddt*dziazic + d2eddt2*ddtdzic*ddtdzia;
789 }
790 }
791 }
792 }
793 }
794 }
795 // ================================================================
796 void eangle2b(int i)
797 {
798 int ia,ib,ic,id;
799 double angle,dot,cosine;
800 double dt,dt2,dt3,dt4;
801 double deddt,terma,termc,term;
802 double xia,yia,zia;
803 double xib,yib,zib;
804 double xic,yic,zic;
805 double xid,yid,zid;
806 double xad,yad,zad;
807 double xbd,ybd,zbd;
808 double xcd,ycd,zcd;
809 double xip,yip,zip;
810 double xap,yap,zap,rap2;
811 double xcp,ycp,zcp,rcp2;
812 double xt,yt,zt,rt2,ptrt2;
813 double xm,ym,zm,rm,delta,delta2;
814 double dedxia,dedyia,dedzia;
815 double dedxib,dedyib,dedzib;
816 double dedxic,dedyic,dedzic;
817 double dedxid,dedyid,dedzid;
818 double dedxip,dedyip,dedzip;
819 double dpdxia,dpdyia,dpdzia;
820 double dpdxic,dpdyic,dpdzic;
821
822 ia = angles.i13[i][0];
823 ib = angles.i13[i][1];
824 ic = angles.i13[i][2];
825 id = angles.i13[i][3];
826 xia = atom[ia].x;
827 yia = atom[ia].y;
828 zia = atom[ia].z;
829 xib = atom[ib].x;
830 yib = atom[ib].y;
831 zib = atom[ib].z;
832 xic = atom[ic].x;
833 yic = atom[ic].y;
834 zic = atom[ic].z;
835 xid = atom[id].x;
836 yid = atom[id].y;
837 zid = atom[id].z;
838
839 deriv.dea[ia][0] = 0.0;
840 deriv.dea[ia][1] = 0.0;
841 deriv.dea[ia][2] = 0.0;
842 deriv.dea[ib][0] = 0.0;
843 deriv.dea[ib][1] = 0.0;
844 deriv.dea[ib][2] = 0.0;
845 deriv.dea[ic][0] = 0.0;
846 deriv.dea[ic][1] = 0.0;
847 deriv.dea[ic][2] = 0.0;
848 deriv.dea[id][0] = 0.0;
849 deriv.dea[id][1] = 0.0;
850 deriv.dea[id][2] = 0.0;
851
852 xad = xia - xid;
853 yad = yia - yid;
854 zad = zia - zid;
855 xbd = xib - xid;
856 ybd = yib - yid;
857 zbd = zib - zid;
858 xcd = xic - xid;
859 ycd = yic - yid;
860 zcd = zic - zid;
861 xt = yad*zcd - zad*ycd;
862 yt = zad*xcd - xad*zcd;
863 zt = xad*ycd - yad*xcd;
864 rt2 = xt*xt + yt*yt + zt*zt;
865 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
866 xip = xib + xt*delta;
867 yip = yib + yt*delta;
868 zip = zib + zt*delta;
869 xap = xia - xip;
870 yap = yia - yip;
871 zap = zia - zip;
872 rap2 = xap*xap + yap*yap + zap*zap;
873 xcp = xic - xip;
874 ycp = yic - yip;
875 zcp = zic - zip;
876 rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
877 xm = ycp*zap - zcp*yap;
878 ym = zcp*xap - xcp*zap;
879 zm = xcp*yap - ycp*xap;
880 rm = sqrt(xm*xm + ym*ym + zm*zm);
881 if (rm != 0.0)
882 {
883 dot = xap*xcp + yap*ycp + zap*zcp;
884 cosine = dot / sqrt(rap2*rcp2);
885 if (cosine < -1.0)
886 cosine = -1.0;
887 if (cosine > 1.0)
888 cosine = 1.0;
889 angle = radian * acos(cosine);
890
891 dt = angle - angles.anat[i];
892 dt2 = dt * dt;
893 dt3 = dt2 * dt;
894 dt4 = dt2 * dt2;
895 deddt = units.angunit * angles.acon[i] * dt
896 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
897 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
898
899 deddt = deddt * radian;
900 terma = -deddt / (rap2*rm);
901 termc = deddt / (rcp2*rm);
902 dedxia = terma * (yap*zm-zap*ym);
903 dedyia = terma * (zap*xm-xap*zm);
904 dedzia = terma * (xap*ym-yap*xm);
905 dedxic = termc * (ycp*zm-zcp*ym);
906 dedyic = termc * (zcp*xm-xcp*zm);
907 dedzic = termc * (xcp*ym-ycp*xm);
908 dedxip = -dedxia - dedxic;
909 dedyip = -dedyia - dedyic;
910 dedzip = -dedzia - dedzic;
911
912 delta2 = 2.0 * delta;
913 ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
914 term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
915 dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
916 term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
917 dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
918 term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
919 dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
920 term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
921 dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
922 term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
923 dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
924 term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
925 dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
926 dedxia = dedxia + dpdxia;
927 dedyia = dedyia + dpdyia;
928 dedzia = dedzia + dpdzia;
929 dedxib = dedxip;
930 dedyib = dedyip;
931 dedzib = dedzip;
932 dedxic = dedxic + dpdxic;
933 dedyic = dedyic + dpdyic;
934 dedzic = dedzic + dpdzic;
935 dedxid = -dedxia - dedxib - dedxic;
936 dedyid = -dedyia - dedyib - dedyic;
937 dedzid = -dedzia - dedzib - dedzic;
938
939 deriv.dea[ia][0] = dedxia;
940 deriv.dea[ia][1] = dedyia;
941 deriv.dea[ia][2] = dedzia;
942 deriv.dea[ib][0] = dedxib;
943 deriv.dea[ib][1] = dedyib;
944 deriv.dea[ib][2] = dedzib;
945 deriv.dea[ic][0] = dedxic;
946 deriv.dea[ic][1] = dedyic;
947 deriv.dea[ic][2] = dedzic;
948 deriv.dea[id][0] = dedxid;
949 deriv.dea[id][1] = dedyid;
950 deriv.dea[id][2] = dedzid;
951 }
952 }
953