ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/eangle.c
Revision: 99
Committed: Mon Jan 19 04:31:29 2009 UTC (10 years, 10 months ago) by wdelano
File size: 40643 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
Line User Rev File contents
1 wdelano 58 #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     energies.ebend += e;
325     }
326     } else
327     {
328     xid = atom[id].x;
329     yid = atom[id].y;
330     zid = atom[id].z;
331     xad = xia - xid;
332     yad = yia - yid;
333     zad = zia - zid;
334     xbd = xib - xid;
335     ybd = yib - yid;
336     zbd = zib - zid;
337     xcd = xic - xid;
338     ycd = yic - yid;
339     zcd = zic - zid;
340     xt = yad*zcd - zad*ycd;
341     yt = zad*xcd - xad*zcd;
342     zt = xad*ycd - yad*xcd;
343     rt2 = xt*xt + yt*yt + zt*zt;
344     delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
345     xip = xib + xt*delta;
346     yip = yib + yt*delta;
347     zip = zib + zt*delta;
348     xap = xia - xip;
349     yap = yia - yip;
350     zap = zia - zip;
351     rap2 = xap*xap + yap*yap + zap*zap;
352     xcp = xic - xip;
353     ycp = yic - yip;
354     zcp = zic - zip;
355     rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
356     xm = ycp*zap - zcp*yap;
357     ym = zcp*xap - xcp*zap;
358     zm = xcp*yap - ycp*xap;
359     rm = sqrt(xm*xm + ym*ym + zm*zm);
360     if (rm != 0.0)
361     {
362     dot = xap*xcp + yap*ycp + zap*zcp;
363     cosine = dot/ sqrt(rap2*rcp2);
364     if (cosine < -1.0)
365     cosine = 1.0;
366     if (cosine > 1.0)
367     cosine = 1.0;
368     temp1 = acos(cosine);
369     temp2 = radian;
370     angle = temp1*temp2;
371     dt = angle - angles.anat[i];
372     dt2 = dt*dt;
373     dt3 = dt2*dt;
374     dt4 = dt2*dt2;
375     e = units.angunit*angles.acon[i]*dt2
376     *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
377    
378     deddt = units.angunit *angles.acon[i]* dt
379     * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
380     + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
381    
382     deddt = deddt * radian;
383     terma = -deddt / (rap2*rm);
384     termc = deddt / (rcp2*rm);
385     dedxia = terma * (yap*zm-zap*ym);
386     dedyia = terma * (zap*xm-xap*zm);
387     dedzia = terma * (xap*ym-yap*xm);
388     dedxic = termc * (ycp*zm-zcp*ym);
389     dedyic = termc * (zcp*xm-xcp*zm);
390     dedzic = termc * (xcp*ym-ycp*xm);
391     dedxip = -dedxia - dedxic;
392     dedyip = -dedyia - dedyic;
393     dedzip = -dedzia - dedzic;
394    
395     delta2 = 2.0 * delta;
396     ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
397     term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
398     dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
399     term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
400     dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
401     term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
402     dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
403     term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
404     dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
405     term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
406     dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
407     term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
408     dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
409    
410     dedxia = dedxia + dpdxia;
411     dedyia = dedyia + dpdyia;
412     dedzia = dedzia + dpdzia;
413     dedxib = dedxip;
414     dedyib = dedyip;
415     dedzib = dedzip;
416     dedxic = dedxic + dpdxic;
417     dedyic = dedyic + dpdyic;
418     dedzic = dedzic + dpdzic;
419     dedxid = -dedxia - dedxib - dedxic;
420     dedyid = -dedyia - dedyib - dedyic;
421     dedzid = -dedzia - dedzib - dedzic;
422    
423     energies.ebend += e;
424     deriv.dea[ia][0] += dedxia;
425     deriv.dea[ia][1] += dedyia;
426     deriv.dea[ia][2] += dedzia;
427    
428     deriv.dea[ib][0] += dedxib;
429     deriv.dea[ib][1] += dedyib;
430     deriv.dea[ib][2] += dedzib;
431    
432     deriv.dea[ic][0] += dedxic;
433     deriv.dea[ic][1] += dedyic;
434     deriv.dea[ic][2] += dedzic;
435    
436     deriv.dea[id][0] += dedxid;
437     deriv.dea[id][1] += dedyid;
438     deriv.dea[id][2] += dedzid;
439    
440     }
441     }
442     }
443     }
444     }
445    
446    
447     void eangle2(int i)
448     {
449     int j,k,ia,ib,ic,id;
450     double old,eps;
451     double **d0; // d0[MAXATOM][3];
452    
453     d0 = dmatrix(0,natom,0,3);
454    
455     eangle2a(i);
456    
457     eps = 1.0e-7;
458     for (k=0; k < angles.nang; k++)
459     {
460     if ( (angles.angin[k] == TRUE) && (field.type != MMFF94) )
461     {
462     ia = angles.i13[k][0];
463     ib = angles.i13[k][1];
464     ic = angles.i13[k][2];
465     id = angles.i13[k][3];
466     if (i == ia || i == ib || i == ic || i == id)
467     {
468     eangle2b(k);
469     for (j=0; j < 3; j++)
470     {
471     d0[ia][j] = deriv.dea[ia][j];
472     d0[ib][j] = deriv.dea[ib][j];
473     d0[ic][j] = deriv.dea[ic][j];
474     d0[id][j] = deriv.dea[id][j];
475     }
476     // numerical x derivative
477     old = atom[i].x;
478     atom[i].x += eps;
479     eangle2b(k);
480     atom[i].x = old;
481     for (j=0; j < 3; j++)
482     {
483     hess.hessx[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
484     hess.hessx[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
485     hess.hessx[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
486     hess.hessx[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
487     }
488     // numerical y derivative
489     old = atom[i].y;
490     atom[i].y += eps;
491     eangle2b(k);
492     atom[i].y = old;
493     for (j=0; j < 3; j++)
494     {
495     hess.hessy[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
496     hess.hessy[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
497     hess.hessy[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
498     hess.hessy[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
499     }
500     // numerical z derivative
501     old = atom[i].z;
502     atom[i].z += eps;
503     eangle2b(k);
504     atom[i].z = old;
505     for (j=0; j < 3; j++)
506     {
507     hess.hessz[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
508     hess.hessz[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
509     hess.hessz[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
510     hess.hessz[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
511     }
512     }
513     }
514     }
515     free_dmatrix(d0 ,0,natom,0,3);
516     }
517     // =================================================
518     void eangle2a(int iatom)
519     {
520     int i,ia,ib,ic;
521     double angle,dot,cosine,factor,sine;
522     double dt,dt2,dt3,dt4;
523     double xia,yia,zia;
524     double xib,yib,zib;
525     double xic,yic,zic;
526     double xab,yab,zab,rab;
527     double xcb,ycb,zcb,rcb;
528     double xp,yp,zp,rp,rp2;
529     double xrab,yrab,zrab,rab2;
530     double xrcb,yrcb,zrcb,rcb2;
531     double xabp,yabp,zabp;
532     double xcbp,ycbp,zcbp;
533     double deddt,d2eddt2,terma,termc;
534     double ddtdxia,ddtdyia,ddtdzia;
535     double ddtdxib,ddtdyib,ddtdzib;
536     double ddtdxic,ddtdyic,ddtdzic;
537     double dxiaxia,dxiayia,dxiazia;
538     double dxibxib,dxibyib,dxibzib;
539     double dxicxic,dxicyic,dxiczic;
540     double dyiayia,dyiazia,dziazia;
541     double dyibyib,dyibzib,dzibzib;
542     double dyicyic,dyiczic,dziczic;
543     double dxibxia,dxibyia,dxibzia;
544     double dyibxia,dyibyia,dyibzia;
545     double dzibxia,dzibyia,dzibzia;
546     double dxibxic,dxibyic,dxibzic;
547     double dyibxic,dyibyic,dyibzic;
548     double dzibxic,dzibyic,dzibzic;
549     double dxiaxic,dxiayic,dxiazic;
550     double dyiaxic,dyiayic,dyiazic;
551     double dziaxic,dziayic,dziazic;
552    
553     for (i=0; i < angles.nang; i++)
554     {
555     ia = angles.i13[i][0];
556     ib = angles.i13[i][1];
557     ic = angles.i13[i][2];
558     if (iatom == ia || iatom == ib || iatom == ic)
559     {
560     xia = atom[ia].x;
561     yia = atom[ia].y;
562     zia = atom[ia].z;
563     xib = atom[ib].x;
564     yib = atom[ib].y;
565     zib = atom[ib].z;
566     xic = atom[ic].x;
567     yic = atom[ic].y;
568     zic = atom[ic].z;
569    
570     if ( (angles.angin[i] == FALSE) || field.type == MMFF94 )
571     {
572     xab = xia - xib;
573     yab = yia - yib;
574     zab = zia - zib;
575     rab = sqrt(xab*xab + yab*yab + zab*zab);
576     xcb = xic - xib;
577     ycb = yic - yib;
578     zcb = zic - zib;
579     rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
580     xp = ycb*zab - zcb*yab;
581     yp = zcb*xab - xcb*zab;
582     zp = xcb*yab - ycb*xab;
583     rp = sqrt(xp*xp + yp*yp + zp*zp);
584     if (rp != 0.0)
585     {
586     dot = xab*xcb + yab*ycb + zab*zcb;
587     cosine = dot / (rab*rcb);
588     if (cosine < -1.0)
589     cosine = -1.0;
590     if (cosine > 1.0)
591     cosine = 1.0;
592     angle = radian * acos(cosine);
593    
594     if (angles.angtype[i] == HARMONIC)
595     {
596     dt = angle - angles.anat[i];
597     dt2 = dt * dt;
598     dt3 = dt2 * dt;
599     dt4 = dt3 * dt;
600     deddt = units.angunit * angles.acon[i] * dt
601     * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
602     + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
603     d2eddt2 = units.angunit * angles.acon[i]
604     * (2.0 + 6.0*units.cang*dt + 12.0*units.qang*dt2
605     + 20.0*units.pang*dt3 + 30.0*units.sang*dt4);
606     terma = -radian/ (rab*rab*rp);
607     termc = radian / (rcb*rcb*rp);
608     } else
609     {
610     factor = angle/radian;
611     sine = sin(factor);
612     deddt = -units.angunit*angles.fcon[i]*
613     (angles.c1[i]*sine + 2*angles.c2[i]*sin(2*factor)
614     + 3*angles.c3[i]*sin(3*factor) + 4*angles.c4[i]*sin(4*factor)
615     + 5*angles.c5[i]*sin(5*factor) + 6*angles.c6[i]*sin(6*factor));
616     d2eddt2 = -units.angunit*angles.fcon[i]*
617     (angles.c1[i]*cosine + 4*angles.c2[i]*cos(2*factor)
618     + 9*angles.c3[i]*cos(3*factor) + 16*angles.c4[i]*cos(4*factor)
619     + 25*angles.c5[i]*cos(5*factor) + 36*angles.c6[i]*cos(6*factor));
620     terma = -1.0 / (rab*rab*rp);
621     termc = 1.0 / (rcb*rcb*rp);
622     }
623     ddtdxia = terma * (yab*zp-zab*yp);
624     ddtdyia = terma * (zab*xp-xab*zp);
625     ddtdzia = terma * (xab*yp-yab*xp);
626     ddtdxic = termc * (ycb*zp-zcb*yp);
627     ddtdyic = termc * (zcb*xp-xcb*zp);
628     ddtdzic = termc * (xcb*yp-ycb*xp);
629     ddtdxib = -ddtdxia - ddtdxic;
630     ddtdyib = -ddtdyia - ddtdyic;
631     ddtdzib = -ddtdzia - ddtdzic;
632     rab2 = 2.0 / (rab*rab);
633     xrab = xab * rab2;
634     yrab = yab * rab2;
635     zrab = zab * rab2;
636     rcb2 = 2.0 / (rcb*rcb);
637     xrcb = xcb * rcb2;
638     yrcb = ycb * rcb2;
639     zrcb = zcb * rcb2;
640     rp2 = 1.0 / (rp*rp);
641     xabp = (yab*zp-zab*yp) * rp2;
642     yabp = (zab*xp-xab*zp) * rp2;
643     zabp = (xab*yp-yab*xp) * rp2;
644     xcbp = (ycb*zp-zcb*yp) * rp2;
645     ycbp = (zcb*xp-xcb*zp) * rp2;
646     zcbp = (xcb*yp-ycb*xp) * rp2;
647    
648     dxiaxia = terma*(xab*xcb-dot) + ddtdxia*(xcbp-xrab);
649     dxiayia = terma*(zp+yab*xcb) + ddtdxia*(ycbp-yrab);
650     dxiazia = terma*(zab*xcb-yp) + ddtdxia*(zcbp-zrab);
651     dyiayia = terma*(yab*ycb-dot) + ddtdyia*(ycbp-yrab);
652     dyiazia = terma*(xp+zab*ycb) + ddtdyia*(zcbp-zrab);
653     dziazia = terma*(zab*zcb-dot) + ddtdzia*(zcbp-zrab);
654     dxicxic = termc*(dot-xab*xcb) - ddtdxic*(xabp+xrcb);
655     dxicyic = termc*(zp-ycb*xab) - ddtdxic*(yabp+yrcb);
656     dxiczic = -termc*(yp+zcb*xab) - ddtdxic*(zabp+zrcb);
657     dyicyic = termc*(dot-yab*ycb) - ddtdyic*(yabp+yrcb);
658     dyiczic = termc*(xp-zcb*yab) - ddtdyic*(zabp+zrcb);
659     dziczic = termc*(dot-zab*zcb) - ddtdzic*(zabp+zrcb);
660     dxiaxic = terma*(yab*yab+zab*zab) - ddtdxia*xabp;
661     dxiayic = -terma*xab*yab - ddtdxia*yabp;
662     dxiazic = -terma*xab*zab - ddtdxia*zabp;
663     dyiaxic = -terma*xab*yab - ddtdyia*xabp;
664     dyiayic = terma*(xab*xab+zab*zab) - ddtdyia*yabp;
665     dyiazic = -terma*yab*zab - ddtdyia*zabp;
666     dziaxic = -terma*xab*zab - ddtdzia*xabp;
667     dziayic = -terma*yab*zab - ddtdzia*yabp;
668     dziazic = terma*(xab*xab+yab*yab) - ddtdzia*zabp;
669    
670     dxibxia = -dxiaxia - dxiaxic;
671     dxibyia = -dxiayia - dyiaxic;
672     dxibzia = -dxiazia - dziaxic;
673     dyibxia = -dxiayia - dxiayic;
674     dyibyia = -dyiayia - dyiayic;
675     dyibzia = -dyiazia - dziayic;
676     dzibxia = -dxiazia - dxiazic;
677     dzibyia = -dyiazia - dyiazic;
678     dzibzia = -dziazia - dziazic;
679     dxibxic = -dxicxic - dxiaxic;
680     dxibyic = -dxicyic - dxiayic;
681     dxibzic = -dxiczic - dxiazic;
682     dyibxic = -dxicyic - dyiaxic;
683     dyibyic = -dyicyic - dyiayic;
684     dyibzic = -dyiczic - dyiazic;
685     dzibxic = -dxiczic - dziaxic;
686     dzibyic = -dyiczic - dziayic;
687     dzibzic = -dziczic - dziazic;
688     dxibxib = -dxibxia - dxibxic;
689     dxibyib = -dxibyia - dxibyic;
690     dxibzib = -dxibzia - dxibzic;
691     dyibyib = -dyibyia - dyibyic;
692     dyibzib = -dyibzia - dyibzic;
693     dzibzib = -dzibzia - dzibzic;
694    
695     if (ia == iatom)
696     {
697     hess.hessx[ia][0] += deddt*dxiaxia + d2eddt2*ddtdxia*ddtdxia;
698     hess.hessx[ia][1] += deddt*dxiayia + d2eddt2*ddtdxia*ddtdyia;
699     hess.hessx[ia][2] += deddt*dxiazia + d2eddt2*ddtdxia*ddtdzia;
700     hess.hessy[ia][0] += deddt*dxiayia + d2eddt2*ddtdyia*ddtdxia;
701     hess.hessy[ia][1] += deddt*dyiayia + d2eddt2*ddtdyia*ddtdyia;
702     hess.hessy[ia][2] += deddt*dyiazia + d2eddt2*ddtdyia*ddtdzia;
703     hess.hessz[ia][0] += deddt*dxiazia + d2eddt2*ddtdzia*ddtdxia;
704     hess.hessz[ia][1] += deddt*dyiazia + d2eddt2*ddtdzia*ddtdyia;
705     hess.hessz[ia][2] += deddt*dziazia + d2eddt2*ddtdzia*ddtdzia;
706     hess.hessx[ib][0] += deddt*dxibxia + d2eddt2*ddtdxia*ddtdxib;
707     hess.hessx[ib][1] += deddt*dyibxia + d2eddt2*ddtdxia*ddtdyib;
708     hess.hessx[ib][2] += deddt*dzibxia + d2eddt2*ddtdxia*ddtdzib;
709     hess.hessy[ib][0] += deddt*dxibyia + d2eddt2*ddtdyia*ddtdxib;
710     hess.hessy[ib][1] += deddt*dyibyia + d2eddt2*ddtdyia*ddtdyib;
711     hess.hessy[ib][2] += deddt*dzibyia + d2eddt2*ddtdyia*ddtdzib;
712     hess.hessz[ib][0] += deddt*dxibzia + d2eddt2*ddtdzia*ddtdxib;
713     hess.hessz[ib][1] += deddt*dyibzia + d2eddt2*ddtdzia*ddtdyib;
714     hess.hessz[ib][2] += deddt*dzibzia + d2eddt2*ddtdzia*ddtdzib;
715     hess.hessx[ic][0] += deddt*dxiaxic + d2eddt2*ddtdxia*ddtdxic;
716     hess.hessx[ic][1] += deddt*dxiayic + d2eddt2*ddtdxia*ddtdyic;
717     hess.hessx[ic][2] += deddt*dxiazic + d2eddt2*ddtdxia*ddtdzic;
718     hess.hessy[ic][0] += deddt*dyiaxic + d2eddt2*ddtdyia*ddtdxic;
719     hess.hessy[ic][1] += deddt*dyiayic + d2eddt2*ddtdyia*ddtdyic;
720     hess.hessy[ic][2] += deddt*dyiazic + d2eddt2*ddtdyia*ddtdzic;
721     hess.hessz[ic][0] += deddt*dziaxic + d2eddt2*ddtdzia*ddtdxic;
722     hess.hessz[ic][1] += deddt*dziayic + d2eddt2*ddtdzia*ddtdyic;
723     hess.hessz[ic][2] += deddt*dziazic + d2eddt2*ddtdzia*ddtdzic;
724     } else if (ib == iatom)
725     {
726     hess.hessx[ib][0] += deddt*dxibxib + d2eddt2*ddtdxib*ddtdxib;
727     hess.hessx[ib][1] += deddt*dxibyib + d2eddt2*ddtdxib*ddtdyib;
728     hess.hessx[ib][2] += deddt*dxibzib + d2eddt2*ddtdxib*ddtdzib;
729     hess.hessy[ib][0] += deddt*dxibyib + d2eddt2*ddtdyib*ddtdxib;
730     hess.hessy[ib][1] += deddt*dyibyib + d2eddt2*ddtdyib*ddtdyib;
731     hess.hessy[ib][2] += deddt*dyibzib + d2eddt2*ddtdyib*ddtdzib;
732     hess.hessz[ib][0] += deddt*dxibzib + d2eddt2*ddtdzib*ddtdxib;
733     hess.hessz[ib][1] += deddt*dyibzib + d2eddt2*ddtdzib*ddtdyib;
734     hess.hessz[ib][2] += deddt*dzibzib + d2eddt2*ddtdzib*ddtdzib;
735     hess.hessx[ia][0] += deddt*dxibxia + d2eddt2*ddtdxib*ddtdxia;
736     hess.hessx[ia][1] += deddt*dxibyia + d2eddt2*ddtdxib*ddtdyia;
737     hess.hessx[ia][2] += deddt*dxibzia + d2eddt2*ddtdxib*ddtdzia;
738     hess.hessy[ia][0] += deddt*dyibxia + d2eddt2*ddtdyib*ddtdxia;
739     hess.hessy[ia][1] += deddt*dyibyia + d2eddt2*ddtdyib*ddtdyia;
740     hess.hessy[ia][2] += deddt*dyibzia + d2eddt2*ddtdyib*ddtdzia;
741     hess.hessz[ia][0] += deddt*dzibxia + d2eddt2*ddtdzib*ddtdxia;
742     hess.hessz[ia][1] += deddt*dzibyia + d2eddt2*ddtdzib*ddtdyia;
743     hess.hessz[ia][2] += deddt*dzibzia + d2eddt2*ddtdzib*ddtdzia;
744     hess.hessx[ic][0] += deddt*dxibxic + d2eddt2*ddtdxib*ddtdxic;
745     hess.hessx[ic][1] += deddt*dxibyic + d2eddt2*ddtdxib*ddtdyic;
746     hess.hessx[ic][2] += deddt*dxibzic + d2eddt2*ddtdxib*ddtdzic;
747     hess.hessy[ic][0] += deddt*dyibxic + d2eddt2*ddtdyib*ddtdxic;
748     hess.hessy[ic][1] += deddt*dyibyic + d2eddt2*ddtdyib*ddtdyic;
749     hess.hessy[ic][2] += deddt*dyibzic + d2eddt2*ddtdyib*ddtdzic;
750     hess.hessz[ic][0] += deddt*dzibxic + d2eddt2*ddtdzib*ddtdxic;
751     hess.hessz[ic][1] += deddt*dzibyic + d2eddt2*ddtdzib*ddtdyic;
752     hess.hessz[ic][2] += deddt*dzibzic + d2eddt2*ddtdzib*ddtdzic;
753     }else if (ic == iatom)
754     {
755     hess.hessx[ic][0] += deddt*dxicxic + d2eddt2*ddtdxic*ddtdxic;
756     hess.hessx[ic][1] += deddt*dxicyic + d2eddt2*ddtdxic*ddtdyic;
757     hess.hessx[ic][2] += deddt*dxiczic + d2eddt2*ddtdxic*ddtdzic;
758     hess.hessy[ic][0] += deddt*dxicyic + d2eddt2*ddtdyic*ddtdxic;
759     hess.hessy[ic][1] += deddt*dyicyic + d2eddt2*ddtdyic*ddtdyic;
760     hess.hessy[ic][2] += deddt*dyiczic + d2eddt2*ddtdyic*ddtdzic;
761     hess.hessz[ic][0] += deddt*dxiczic + d2eddt2*ddtdzic*ddtdxic;
762     hess.hessz[ic][1] += deddt*dyiczic + d2eddt2*ddtdzic*ddtdyic;
763     hess.hessz[ic][2] += deddt*dziczic + d2eddt2*ddtdzic*ddtdzic;
764     hess.hessx[ib][0] += deddt*dxibxic + d2eddt2*ddtdxic*ddtdxib;
765     hess.hessx[ib][1] += deddt*dyibxic + d2eddt2*ddtdxic*ddtdyib;
766     hess.hessx[ib][2] += deddt*dzibxic + d2eddt2*ddtdxic*ddtdzib;
767     hess.hessy[ib][0] += deddt*dxibyic + d2eddt2*ddtdyic*ddtdxib;
768     hess.hessy[ib][1] += deddt*dyibyic + d2eddt2*ddtdyic*ddtdyib;
769     hess.hessy[ib][2] += deddt*dzibyic + d2eddt2*ddtdyic*ddtdzib;
770     hess.hessz[ib][0] += deddt*dxibzic + d2eddt2*ddtdzic*ddtdxib;
771     hess.hessz[ib][1] += deddt*dyibzic + d2eddt2*ddtdzic*ddtdyib;
772     hess.hessz[ib][2] += deddt*dzibzic + d2eddt2*ddtdzic*ddtdzib;
773     hess.hessx[ia][0] += deddt*dxiaxic + d2eddt2*ddtdxic*ddtdxia;
774     hess.hessx[ia][1] += deddt*dyiaxic + d2eddt2*ddtdxic*ddtdyia;
775     hess.hessx[ia][2] += deddt*dziaxic + d2eddt2*ddtdxic*ddtdzia;
776     hess.hessy[ia][0] += deddt*dxiayic + d2eddt2*ddtdyic*ddtdxia;
777     hess.hessy[ia][1] += deddt*dyiayic + d2eddt2*ddtdyic*ddtdyia;
778     hess.hessy[ia][2] += deddt*dziayic + d2eddt2*ddtdyic*ddtdzia;
779     hess.hessz[ia][0] += deddt*dxiazic + d2eddt2*ddtdzic*ddtdxia;
780     hess.hessz[ia][1] += deddt*dyiazic + d2eddt2*ddtdzic*ddtdyia;
781     hess.hessz[ia][2] += deddt*dziazic + d2eddt2*ddtdzic*ddtdzia;
782     }
783     }
784     }
785     }
786     }
787     }
788     // ================================================================
789     void eangle2b(int i)
790     {
791     int ia,ib,ic,id;
792     double angle,dot,cosine;
793     double dt,dt2,dt3,dt4;
794     double deddt,terma,termc,term;
795     double xia,yia,zia;
796     double xib,yib,zib;
797     double xic,yic,zic;
798     double xid,yid,zid;
799     double xad,yad,zad;
800     double xbd,ybd,zbd;
801     double xcd,ycd,zcd;
802     double xip,yip,zip;
803     double xap,yap,zap,rap2;
804     double xcp,ycp,zcp,rcp2;
805     double xt,yt,zt,rt2,ptrt2;
806     double xm,ym,zm,rm,delta,delta2;
807     double dedxia,dedyia,dedzia;
808     double dedxib,dedyib,dedzib;
809     double dedxic,dedyic,dedzic;
810     double dedxid,dedyid,dedzid;
811     double dedxip,dedyip,dedzip;
812     double dpdxia,dpdyia,dpdzia;
813     double dpdxic,dpdyic,dpdzic;
814    
815     ia = angles.i13[i][0];
816     ib = angles.i13[i][1];
817     ic = angles.i13[i][2];
818     id = angles.i13[i][3];
819     xia = atom[ia].x;
820     yia = atom[ia].y;
821     zia = atom[ia].z;
822     xib = atom[ib].x;
823     yib = atom[ib].y;
824     zib = atom[ib].z;
825     xic = atom[ic].x;
826     yic = atom[ic].y;
827     zic = atom[ic].z;
828     xid = atom[id].x;
829     yid = atom[id].y;
830     zid = atom[id].z;
831    
832     deriv.dea[ia][0] = 0.0;
833     deriv.dea[ia][1] = 0.0;
834     deriv.dea[ia][2] = 0.0;
835     deriv.dea[ib][0] = 0.0;
836     deriv.dea[ib][1] = 0.0;
837     deriv.dea[ib][2] = 0.0;
838     deriv.dea[ic][0] = 0.0;
839     deriv.dea[ic][1] = 0.0;
840     deriv.dea[ic][2] = 0.0;
841     deriv.dea[id][0] = 0.0;
842     deriv.dea[id][1] = 0.0;
843     deriv.dea[id][2] = 0.0;
844    
845     xad = xia - xid;
846     yad = yia - yid;
847     zad = zia - zid;
848     xbd = xib - xid;
849     ybd = yib - yid;
850     zbd = zib - zid;
851     xcd = xic - xid;
852     ycd = yic - yid;
853     zcd = zic - zid;
854     xt = yad*zcd - zad*ycd;
855     yt = zad*xcd - xad*zcd;
856     zt = xad*ycd - yad*xcd;
857     rt2 = xt*xt + yt*yt + zt*zt;
858     delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
859     xip = xib + xt*delta;
860     yip = yib + yt*delta;
861     zip = zib + zt*delta;
862     xap = xia - xip;
863     yap = yia - yip;
864     zap = zia - zip;
865     rap2 = xap*xap + yap*yap + zap*zap;
866     xcp = xic - xip;
867     ycp = yic - yip;
868     zcp = zic - zip;
869     rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
870     xm = ycp*zap - zcp*yap;
871     ym = zcp*xap - xcp*zap;
872     zm = xcp*yap - ycp*xap;
873     rm = sqrt(xm*xm + ym*ym + zm*zm);
874     if (rm != 0.0)
875     {
876     dot = xap*xcp + yap*ycp + zap*zcp;
877     cosine = dot / sqrt(rap2*rcp2);
878     if (cosine < -1.0)
879     cosine = -1.0;
880     if (cosine > 1.0)
881     cosine = 1.0;
882     angle = radian * acos(cosine);
883    
884     dt = angle - angles.anat[i];
885     dt2 = dt * dt;
886     dt3 = dt2 * dt;
887     dt4 = dt2 * dt2;
888     deddt = units.angunit * angles.acon[i] * dt
889     * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
890     + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
891    
892     deddt = deddt * radian;
893     terma = -deddt / (rap2*rm);
894     termc = deddt / (rcp2*rm);
895     dedxia = terma * (yap*zm-zap*ym);
896     dedyia = terma * (zap*xm-xap*zm);
897     dedzia = terma * (xap*ym-yap*xm);
898     dedxic = termc * (ycp*zm-zcp*ym);
899     dedyic = termc * (zcp*xm-xcp*zm);
900     dedzic = termc * (xcp*ym-ycp*xm);
901     dedxip = -dedxia - dedxic;
902     dedyip = -dedyia - dedyic;
903     dedzip = -dedzia - dedzic;
904    
905     delta2 = 2.0 * delta;
906     ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
907     term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
908     dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
909     term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
910     dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
911     term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
912     dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
913     term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
914     dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
915     term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
916     dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
917     term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
918     dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
919     dedxia = dedxia + dpdxia;
920     dedyia = dedyia + dpdyia;
921     dedzia = dedzia + dpdzia;
922     dedxib = dedxip;
923     dedyib = dedyip;
924     dedzib = dedzip;
925     dedxic = dedxic + dpdxic;
926     dedyic = dedyic + dpdyic;
927     dedzic = dedzic + dpdzic;
928     dedxid = -dedxia - dedxib - dedxic;
929     dedyid = -dedyia - dedyib - dedyic;
930     dedzid = -dedzia - dedzib - dedzic;
931    
932     deriv.dea[ia][0] = dedxia;
933     deriv.dea[ia][1] = dedyia;
934     deriv.dea[ia][2] = dedzia;
935     deriv.dea[ib][0] = dedxib;
936     deriv.dea[ib][1] = dedyib;
937     deriv.dea[ib][2] = dedzib;
938     deriv.dea[ic][0] = dedxic;
939     deriv.dea[ic][1] = dedyic;
940     deriv.dea[ic][2] = dedzic;
941     deriv.dea[id][0] = dedxid;
942     deriv.dea[id][1] = dedyid;
943     deriv.dea[id][2] = dedzid;
944     }
945     }
946