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 (11 years, 9 months ago) by gilbertke
File size: 41041 byte(s)
Log Message:
more cleanup and code removal
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     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