ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/eangle.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 11 months ago) by tjod
File size: 57397 byte(s)
Log Message:
test

Line User Rev File contents
1 tjod 3 #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 "gmmx.h"
12     #include "field.h"
13     #include "fix.h"
14    
15     void eangle2a(int);
16     void eangle2a_fixed(int);
17     void eangle2b(int);
18     void eangle3(int, double (*)[3]);
19     //double acos(double);
20    
21     EXTERN struct t_units {
22     double bndunit, cbnd, qbnd;
23     double angunit, cang, qang, pang, sang, aaunit;
24     double stbnunit, ureyunit, torsunit, storunit, v14scale;
25     double aterm, bterm, cterm, dielec, chgscale;
26     } units;
27    
28     EXTERN struct t_minim_values {
29     int iprint, ndc, nconst;
30     float dielc;
31     } minim_values;
32    
33     void eangle()
34     {
35     int i, ia, ib, ic, id;
36     double e, angle, dot, cosine, factor;
37     double dt, dt2, dt3, dt4;
38     double xia, yia, zia;
39     double xib, yib, zib;
40     double xic, yic, zic;
41     double xid, yid, zid;
42     double xab, yab, zab, rab2;
43     double xcb, ycb, zcb, rcb2;
44     double xad, yad, zad;
45     double xbd, ybd, zbd;
46     double xcd,ycd,zcd;
47     double xip,yip,zip;
48     double xap,yap,zap,rap2;
49     double xcp,ycp,zcp,rcp2;
50     double xt,yt,zt,rt2,delta;
51    
52    
53    
54     if (minim_values.iprint)
55     {
56     fprintf(pcmoutfile,"\nAngle Terms \n");
57     fprintf(pcmoutfile," At1 At2 At3 Angle Thet0 Tconst Ebend\n");
58     }
59     gmmx.hybrida = FALSE;
60    
61     for (i=0; i < angles.nang; i++)
62     {
63     ia = angles.i13[i][0];
64     ib = angles.i13[i][1];
65     ic = angles.i13[i][2];
66     id = angles.i13[i][3];
67     if (atom[ia].use || atom[ib].use || atom[ic].use)
68     {
69     xia = atom[ia].x;
70     yia = atom[ia].y;
71     zia = atom[ia].z;
72     xib = atom[ib].x;
73     yib = atom[ib].y;
74     zib = atom[ib].z;
75     xic = atom[ic].x;
76     yic = atom[ic].y;
77     zic = atom[ic].z;
78     if ( (angles.angin[i] == FALSE) || field.type == MMFF94 )
79     {
80     xab = xia - xib;
81     yab = yia - yib;
82     zab = zia - zib;
83     rab2 = xab*xab + yab*yab + zab*zab;
84     xcb = xic - xib;
85     ycb = yic - yib;
86     zcb = zic - zib;
87     rcb2 = xcb*xcb + ycb*ycb + zcb*zcb;
88     if (rab2 != 0.0 && rcb2 != 0.00)
89     {
90     dot = xab*xcb + yab*ycb + zab*zcb;
91     cosine = dot / sqrt(rab2*rcb2);
92     if (cosine > 1.0)
93     cosine = 1.0;
94     if (cosine < -1.0)
95     cosine = -1.0;
96     angle = radian*acos(cosine);
97    
98     if (angle < 0.75*angles.anat[i])
99     if ( atom[ib].mmx_type != 80)
100     gmmx.hybrida = TRUE;
101    
102     if (angles.angtype[i] == HARMONIC)
103     {
104     dt = angle - angles.anat[i];
105     dt2 = dt*dt;
106     dt3 = dt2*dt;
107     dt4 = dt2*dt2;
108     e = units.angunit*angles.acon[i]*dt2
109     *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
110     } else
111     {
112     factor = angle/radian;
113     e = units.angunit*angles.fcon[i]*
114     (angles.c0[i] + angles.c1[i]*cosine + angles.c2[i]*cos(2*factor)
115     + angles.c3[i]*cos(3*factor) + angles.c4[i]*cos(4*factor)
116     + angles.c5[i]*cos(5*factor) + angles.c6[i]*cos(6*factor));
117     }
118     energies.ebend += e;
119     atom[ia].energy += e;
120     atom[ib].energy += e;
121     atom[ic].energy += e;
122     if (minim_values.iprint)
123     {
124     if (angles.angtype[i] == HARMONIC)
125     fprintf(pcmoutfile,"Angle 1-3: %2s(%-3d)- %2s(%-3d)- %2s(%-3d) %-8.3f %-8.3f %-8.4f = %-8.4f\n",
126     atom[ia].name,ia,atom[ib].name, ib,atom[ic].name, ic, angle, angles.anat[i],angles.acon[i], e);
127     else
128     fprintf(pcmoutfile,"Angle 1-3: %2s(%-3d)- %2s(%-3d)- %2s(%-3d) %-8.3f %-8.4f = %-8.4f\n",
129     atom[ia].name,ia,atom[ib].name, ib,atom[ic].name, ic, angle, angles.fcon[i], e);
130     }
131     }
132     } else
133     {
134     xid = atom[id].x;
135     yid = atom[id].y;
136     zid = atom[id].z;
137     xad = xia - xid;
138     yad = yia - yid;
139     zad = zia - zid;
140     xbd = xib - xid;
141     ybd = yib - yid;
142     zbd = zib - zid;
143     xcd = xic - xid;
144     ycd = yic - yid;
145     zcd = zic - zid;
146     xt = yad*zcd - zad*ycd;
147     yt = zad*xcd - xad*zcd;
148     zt = xad*ycd - yad*xcd;
149     rt2 = xt*xt + yt*yt + zt*zt;
150     delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
151     xip = xib + xt*delta;
152     yip = yib + yt*delta;
153     zip = zib + zt*delta;
154     xap = xia - xip;
155     yap = yia - yip;
156     zap = zia - zip;
157     rap2 = xap*xap + yap*yap + zap*zap;
158     xcp = xic - xip;
159     ycp = yic - yip;
160     zcp = zic - zip;
161     rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
162     if (rap2 != 0.0 && rcp2 != 0.0)
163     {
164     dot = xap*xcp + yap*ycp + zap*zcp;
165     cosine = dot/ sqrt(rap2*rcp2);
166     if (cosine < -1.0)
167     cosine = -1.0;
168     if (cosine > 1.0)
169     cosine = 1.0;
170     angle = radian*acos(cosine);
171     dt = angle - angles.anat[i];
172     dt2 = dt*dt;
173     dt3 = dt2*dt;
174     dt4 = dt2*dt2;
175     e = units.angunit*angles.acon[i]*dt2
176     *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
177     energies.ebend += e;
178     atom[ia].energy += e;
179     atom[ib].energy += e;
180     atom[ic].energy += e;
181     if (minim_values.iprint)
182     fprintf(pcmoutfile,"Angle InP: %2s(%-3d)- %2s(%-3d)- %2s(%-3d) %-8.3f %-8.3f %-8.4f = %-8.4f\n",
183     atom[ia].name, ia, atom[ib].name, ib, atom[ic].name, ic, angle, angles.anat[i], angles.acon[i], e);
184    
185     }
186     }
187     }
188     }
189    
190     }
191    
192     void eangle1()
193     {
194     int i, ia, ib, ic, id;
195     double temp1, temp2;
196     double e, angle, dot, cosine, factor, sine;
197     double dt, dt2, dt3, dt4;
198     double deddt,terma,termc,term;
199     double xia, yia, zia;
200     double xib, yib, zib;
201     double xic, yic, zic;
202     double xid, yid, zid;
203     double xab, yab, zab, rab2;
204     double xcb, ycb, zcb, rcb2;
205     double xp,yp,zp,rp;
206     double xad, yad, zad;
207     double xbd, ybd, zbd;
208     double xcd,ycd,zcd;
209     double xip,yip,zip;
210     double xap,yap,zap,rap2;
211     double xcp,ycp,zcp,rcp2;
212     double xt,yt,zt,rt2,ptrt2;
213     double xm,ym,zm,rm,delta,delta2;
214     double dedxia,dedyia,dedzia;
215     double dedxib,dedyib,dedzib;
216     double dedxic,dedyic,dedzic;
217     double dedxid,dedyid,dedzid;
218     double dedxip,dedyip,dedzip;
219     double dpdxia,dpdyia,dpdzia;
220     double dpdxic,dpdyic,dpdzic;
221     double dtemp;
222    
223     energies.ebend = 0.0F;
224     for (i=0; i <= natom; i++)
225     {
226     deriv.dea[i][0] = 0.0;
227     deriv.dea[i][1] = 0.0;
228     deriv.dea[i][2] = 0.0;
229     }
230     gmmx.hybrida = FALSE;
231     for (i=0; i < angles.nang; i++)
232     {
233     ia = angles.i13[i][0];
234     ib = angles.i13[i][1];
235     ic = angles.i13[i][2];
236     id = angles.i13[i][3];
237     if (atom[ia].use || atom[ib].use || atom[ic].use)
238     {
239     xia = atom[ia].x;
240     yia = atom[ia].y;
241     zia = atom[ia].z;
242     xib = atom[ib].x;
243     yib = atom[ib].y;
244     zib = atom[ib].z;
245     xic = atom[ic].x;
246     yic = atom[ic].y;
247     zic = atom[ic].z;
248     if ( (angles.angin[i] == FALSE) || field.type == MMFF94 )
249     {
250     xab = xia - xib;
251     yab = yia - yib;
252     zab = zia - zib;
253     rab2 = xab*xab + yab*yab + zab*zab;
254     xcb = xic - xib;
255     ycb = yic - yib;
256     zcb = zic - zib;
257     rcb2 = xcb*xcb + ycb*ycb + zcb*zcb;
258     xp = ycb*zab - zcb*yab;
259     yp = zcb*xab - xcb*zab;
260     zp = xcb*yab - ycb*xab;
261     rp = sqrt(xp*xp + yp*yp + zp*zp);
262     if (rp != 0.0 )
263     {
264     dot = xab*xcb + yab*ycb + zab*zcb;
265     cosine = dot / sqrt(rab2*rcb2);
266     if (cosine > 1.0)
267     cosine = 1.0;
268     if (cosine < -1.0)
269     cosine = -1.0;
270     angle = radian*acos(cosine);
271     if (angle < 0.75*angles.anat[i])
272     if ( atom[ib].mmx_type != 80)
273     gmmx.hybrida = TRUE;
274     if (angles.angtype[i] == HARMONIC)
275     {
276     dt = angle - angles.anat[i];
277     dt2 = dt*dt;
278     dt3 = dt2*dt;
279     dt4 = dt2*dt2;
280    
281     e = units.angunit*angles.acon[i]*dt2
282     *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
283     deddt = units.angunit*angles.acon[i]* dt * radian
284     * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
285     + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
286     } else
287     {
288     factor = angle/radian;
289     sine = sin(factor);
290     e = units.angunit*angles.fcon[i]*
291     (angles.c0[i] + angles.c1[i]*cosine + angles.c2[i]*cos(2*factor)
292     + angles.c3[i]*cos(3*factor) + angles.c4[i]*cos(4*factor)
293     + angles.c5[i]*cos(5*factor) + angles.c6[i]*cos(6*factor));
294     deddt = -units.angunit*angles.fcon[i]*
295     (angles.c1[i]*sine + 2*angles.c2[i]*sin(2*factor)
296     + 3*angles.c3[i]*sin(3*factor) + 4*angles.c4[i]*sin(4*factor)
297     + 5*angles.c5[i]*sin(5*factor) + 6*angles.c6[i]*sin(6*factor));
298     }
299    
300     /* terma = -deddt / (rab2*rp);
301     termc = deddt / (rcb2*rp);
302     dedxia = terma * (yab*zp-zab*yp);
303     dedyia = terma * (zab*xp-xab*zp);
304     dedzia = terma * (xab*yp-yab*xp);
305     dedxic = termc * (ycb*zp-zcb*yp);
306     dedyic = termc * (zcb*xp-xcb*zp);
307     dedzic = termc * (xcb*yp-ycb*xp);
308     dedxib = -dedxia - dedxic;
309     dedyib = -dedyia - dedyic;
310     dedzib = -dedzia - dedzic; */
311    
312     dtemp = dot*dot/(rab2*rcb2);
313     if (fabs(dtemp) >= 1.0)
314     dtemp = 0.9999*dtemp/(fabs(dtemp));
315     terma = sqrt(1.0-dtemp);
316     termc = sqrt(rab2*rcb2);
317     dedxia = -deddt*( xcb/termc - (xab*cosine)/rab2)/terma;
318     dedyia = -deddt*( ycb/termc - (yab*cosine)/rab2)/terma;
319     dedzia = -deddt*( zcb/termc - (zab*cosine)/rab2)/terma;
320    
321     dedxic = -deddt*( xab/termc - (xcb*cosine)/rcb2)/terma;
322     dedyic = -deddt*( yab/termc - (ycb*cosine)/rcb2)/terma;
323     dedzic = -deddt*( zab/termc - (zcb*cosine)/rcb2)/terma;
324    
325     dedxib = -deddt*( (-xab-xcb)/termc + (dot*(xcb*rab2 + xab*rcb2))/(rab2*rcb2*termc))/terma;
326     dedyib = -deddt*( (-yab-ycb)/termc + (dot*(ycb*rab2 + yab*rcb2))/(rab2*rcb2*termc))/terma;
327     dedzib = -deddt*( (-zab-zcb)/termc + (dot*(zcb*rab2 + zab*rcb2))/(rab2*rcb2*termc))/terma;
328    
329    
330     deriv.dea[ia][0] += dedxia;
331     deriv.dea[ia][1] += dedyia;
332     deriv.dea[ia][2] += dedzia;
333    
334     deriv.dea[ib][0] += dedxib;
335     deriv.dea[ib][1] += dedyib;
336     deriv.dea[ib][2] += dedzib;
337    
338     deriv.dea[ic][0] += dedxic;
339     deriv.dea[ic][1] += dedyic;
340     deriv.dea[ic][2] += dedzic;
341    
342     virial.virx += xab*dedxia + xcb*dedxic;
343     virial.viry += yab*dedyia + ycb*dedyic;
344     virial.virz += zab*dedzia + zcb*dedzic;
345    
346     energies.ebend += e;
347     }
348     } else
349     {
350     xid = atom[id].x;
351     yid = atom[id].y;
352     zid = atom[id].z;
353     xad = xia - xid;
354     yad = yia - yid;
355     zad = zia - zid;
356     xbd = xib - xid;
357     ybd = yib - yid;
358     zbd = zib - zid;
359     xcd = xic - xid;
360     ycd = yic - yid;
361     zcd = zic - zid;
362     xt = yad*zcd - zad*ycd;
363     yt = zad*xcd - xad*zcd;
364     zt = xad*ycd - yad*xcd;
365     rt2 = xt*xt + yt*yt + zt*zt;
366     delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
367     xip = xib + xt*delta;
368     yip = yib + yt*delta;
369     zip = zib + zt*delta;
370     xap = xia - xip;
371     yap = yia - yip;
372     zap = zia - zip;
373     rap2 = xap*xap + yap*yap + zap*zap;
374     xcp = xic - xip;
375     ycp = yic - yip;
376     zcp = zic - zip;
377     rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
378     xm = ycp*zap - zcp*yap;
379     ym = zcp*xap - xcp*zap;
380     zm = xcp*yap - ycp*xap;
381     rm = sqrt(xm*xm + ym*ym + zm*zm);
382     if (rm != 0.0)
383     {
384     dot = xap*xcp + yap*ycp + zap*zcp;
385     cosine = dot/ sqrt(rap2*rcp2);
386     if (cosine < -1.0)
387     cosine = 1.0;
388     if (cosine > 1.0)
389     cosine = 1.0;
390     temp1 = acos(cosine);
391     temp2 = radian;
392     angle = temp1*temp2;
393     dt = angle - angles.anat[i];
394     dt2 = dt*dt;
395     dt3 = dt2*dt;
396     dt4 = dt2*dt2;
397     e = units.angunit*angles.acon[i]*dt2
398     *(1.0 + units.cang*dt + units.qang*dt2 + units.pang*dt3 + units.sang*dt4);
399    
400     deddt = units.angunit *angles.acon[i]* dt
401     * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
402     + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
403    
404     deddt = deddt * radian;
405     terma = -deddt / (rap2*rm);
406     termc = deddt / (rcp2*rm);
407     dedxia = terma * (yap*zm-zap*ym);
408     dedyia = terma * (zap*xm-xap*zm);
409     dedzia = terma * (xap*ym-yap*xm);
410     dedxic = termc * (ycp*zm-zcp*ym);
411     dedyic = termc * (zcp*xm-xcp*zm);
412     dedzic = termc * (xcp*ym-ycp*xm);
413     dedxip = -dedxia - dedxic;
414     dedyip = -dedyia - dedyic;
415     dedzip = -dedzia - dedzic;
416    
417     delta2 = 2.0 * delta;
418     ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
419     term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
420     dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
421     term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
422     dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
423     term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
424     dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
425     term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
426     dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
427     term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
428     dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
429     term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
430     dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
431    
432     dedxia = dedxia + dpdxia;
433     dedyia = dedyia + dpdyia;
434     dedzia = dedzia + dpdzia;
435     dedxib = dedxip;
436     dedyib = dedyip;
437     dedzib = dedzip;
438     dedxic = dedxic + dpdxic;
439     dedyic = dedyic + dpdyic;
440     dedzic = dedzic + dpdzic;
441     dedxid = -dedxia - dedxib - dedxic;
442     dedyid = -dedyia - dedyib - dedyic;
443     dedzid = -dedzia - dedzib - dedzic;
444    
445     energies.ebend += e;
446     deriv.dea[ia][0] += dedxia;
447     deriv.dea[ia][1] += dedyia;
448     deriv.dea[ia][2] += dedzia;
449    
450     deriv.dea[ib][0] += dedxib;
451     deriv.dea[ib][1] += dedyib;
452     deriv.dea[ib][2] += dedzib;
453    
454     deriv.dea[ic][0] += dedxic;
455     deriv.dea[ic][1] += dedyic;
456     deriv.dea[ic][2] += dedzic;
457    
458     deriv.dea[id][0] += dedxid;
459     deriv.dea[id][1] += dedyid;
460     deriv.dea[id][2] += dedzid;
461    
462     virial.virx += xad*dedxia + xbd*dedxib + xcd*dedxic;
463     virial.viry += yad*dedyia + ybd*dedyib + ycd*dedyic;
464     virial.virz += zad*dedzia + zbd*dedzib + zcd*dedzic;
465     }
466     }
467     }
468     }
469     // fixed angle code
470     for (i=0; i < fixangle.nfxangle; i++)
471     {
472     ia = fixangle.ifxangle[i][0];
473     ib = fixangle.ifxangle[i][1];
474     ic = fixangle.ifxangle[i][2];
475     id = 0;
476     if (atom[ia].use || atom[ib].use || atom[ic].use)
477     {
478     xia = atom[ia].x;
479     yia = atom[ia].y;
480     zia = atom[ia].z;
481     xib = atom[ib].x;
482     yib = atom[ib].y;
483     zib = atom[ib].z;
484     xic = atom[ic].x;
485     yic = atom[ic].y;
486     zic = atom[ic].z;
487     xab = xia - xib;
488     yab = yia - yib;
489     zab = zia - zib;
490     rab2 = xab*xab + yab*yab + zab*zab;
491     xcb = xic - xib;
492     ycb = yic - yib;
493     zcb = zic - zib;
494     rcb2 = xcb*xcb + ycb*ycb + zcb*zcb;
495     xp = ycb*zab - zcb*yab;
496     yp = zcb*xab - xcb*zab;
497     zp = xcb*yab - ycb*xab;
498     rp = sqrt(xp*xp + yp*yp + zp*zp);
499     if (rp != 0.0 )
500     {
501     dot = xab*xcb + yab*ycb + zab*zcb;
502     cosine = dot / sqrt(rab2*rcb2);
503     if (cosine > 1.0)
504     cosine = 1.0;
505     if (cosine < -1.0)
506     cosine = -1.0;
507     angle = radian*acos(cosine);
508     if (angle < 0.75*angles.anat[i])
509     if ( atom[ib].mmx_type != 80)
510     gmmx.hybrida = TRUE;
511    
512     dt = angle - fixangle.anat[i];
513     dt2 = dt*dt;
514     dt3 = dt2*dt;
515     dt4 = dt2*dt2;
516    
517     e = units.angunit*fixangle.acon[i]*dt2;
518     deddt = units.angunit*fixangle.acon[i]* dt*2.0;
519     deddt = deddt * radian;
520    
521     terma = sqrt(1.00-( (dot*dot)/(rab2*rcb2)));
522     termc = sqrt(rab2*rcb2);
523     dedxia = -deddt*( xcb/termc - (xab*cosine)/rab2)/terma;
524     dedyia = -deddt*( ycb/termc - (yab*cosine)/rab2)/terma;
525     dedzia = -deddt*( zcb/termc - (zab*cosine)/rab2)/terma;
526    
527     dedxic = -deddt*( xab/termc - (xcb*cosine)/rcb2)/terma;
528     dedyic = -deddt*( yab/termc - (ycb*cosine)/rcb2)/terma;
529     dedzic = -deddt*( zab/termc - (zcb*cosine)/rcb2)/terma;
530    
531     dedxib = -deddt*( (-xab-xcb)/termc + (dot*(xcb*rab2 + xab*rcb2))/(rab2*rcb2*termc))/terma;
532     dedyib = -deddt*( (-yab-ycb)/termc + (dot*(ycb*rab2 + yab*rcb2))/(rab2*rcb2*termc))/terma;
533     dedzib = -deddt*( (-zab-zcb)/termc + (dot*(zcb*rab2 + zab*rcb2))/(rab2*rcb2*termc))/terma;
534    
535    
536     deriv.dea[ia][0] += dedxia;
537     deriv.dea[ia][1] += dedyia;
538     deriv.dea[ia][2] += dedzia;
539    
540     deriv.dea[ib][0] += dedxib;
541     deriv.dea[ib][1] += dedyib;
542     deriv.dea[ib][2] += dedzib;
543    
544     deriv.dea[ic][0] += dedxic;
545     deriv.dea[ic][1] += dedyic;
546     deriv.dea[ic][2] += dedzic;
547    
548     energies.ebend += e;
549     }
550     }
551     }
552     // end of fixed angles
553     }
554    
555    
556     void eangle2(int i)
557     {
558     int j,k,ia,ib,ic,id;
559     double old,eps;
560     double **d0; // d0[MAXATOM][3];
561    
562     d0 = dmatrix(0,natom,0,3);
563    
564     eangle2a(i);
565    
566     eps = 1.0e-7;
567     for (k=0; k < angles.nang; k++)
568     {
569     if ( (angles.angin[k] == TRUE) && (field.type != MMFF94) )
570     {
571     ia = angles.i13[k][0];
572     ib = angles.i13[k][1];
573     ic = angles.i13[k][2];
574     id = angles.i13[k][3];
575     if (i == ia || i == ib || i == ic || i == id)
576     {
577     eangle2b(k);
578     for (j=0; j < 3; j++)
579     {
580     d0[ia][j] = deriv.dea[ia][j];
581     d0[ib][j] = deriv.dea[ib][j];
582     d0[ic][j] = deriv.dea[ic][j];
583     d0[id][j] = deriv.dea[id][j];
584     }
585     // numerical x derivative
586     old = atom[i].x;
587     atom[i].x += eps;
588     eangle2b(k);
589     atom[i].x = old;
590     for (j=0; j < 3; j++)
591     {
592     hess.hessx[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
593     hess.hessx[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
594     hess.hessx[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
595     hess.hessx[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
596     }
597     // numerical y derivative
598     old = atom[i].y;
599     atom[i].y += eps;
600     eangle2b(k);
601     atom[i].y = old;
602     for (j=0; j < 3; j++)
603     {
604     hess.hessy[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
605     hess.hessy[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
606     hess.hessy[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
607     hess.hessy[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
608     }
609     // numerical z derivative
610     old = atom[i].z;
611     atom[i].z += eps;
612     eangle2b(k);
613     atom[i].z = old;
614     for (j=0; j < 3; j++)
615     {
616     hess.hessz[ia][j] += (deriv.dea[ia][j]-d0[ia][j])/eps;
617     hess.hessz[ib][j] += (deriv.dea[ib][j]-d0[ib][j])/eps;
618     hess.hessz[ic][j] += (deriv.dea[ic][j]-d0[ic][j])/eps;
619     hess.hessz[id][j] += (deriv.dea[id][j]-d0[id][j])/eps;
620     }
621     }
622     }
623     }
624     // fixed angle code
625     if (fixangle.nfxangle > 0)
626     eangle2a_fixed(i);
627    
628     free_dmatrix(d0 ,0,natom,0,3);
629     }
630    
631     void eangle2a(int iatom)
632     {
633     int i,ia,ib,ic;
634     double angle,dot,cosine,factor,sine;
635     double dt,dt2,dt3,dt4;
636     double xia,yia,zia;
637     double xib,yib,zib;
638     double xic,yic,zic;
639     double xab,yab,zab,rab;
640     double xcb,ycb,zcb,rcb;
641     double xp,yp,zp,rp,rp2;
642     double xrab,yrab,zrab,rab2;
643     double xrcb,yrcb,zrcb,rcb2;
644     double xabp,yabp,zabp;
645     double xcbp,ycbp,zcbp;
646     double deddt,d2eddt2,terma,termc;
647     double ddtdxia,ddtdyia,ddtdzia;
648     double ddtdxib,ddtdyib,ddtdzib;
649     double ddtdxic,ddtdyic,ddtdzic;
650     double dxiaxia,dxiayia,dxiazia;
651     double dxibxib,dxibyib,dxibzib;
652     double dxicxic,dxicyic,dxiczic;
653     double dyiayia,dyiazia,dziazia;
654     double dyibyib,dyibzib,dzibzib;
655     double dyicyic,dyiczic,dziczic;
656     double dxibxia,dxibyia,dxibzia;
657     double dyibxia,dyibyia,dyibzia;
658     double dzibxia,dzibyia,dzibzia;
659     double dxibxic,dxibyic,dxibzic;
660     double dyibxic,dyibyic,dyibzic;
661     double dzibxic,dzibyic,dzibzic;
662     double dxiaxic,dxiayic,dxiazic;
663     double dyiaxic,dyiayic,dyiazic;
664     double dziaxic,dziayic,dziazic;
665    
666     for (i=0; i < angles.nang; i++)
667     {
668     ia = angles.i13[i][0];
669     ib = angles.i13[i][1];
670     ic = angles.i13[i][2];
671     if (iatom == ia || iatom == ib || iatom == ic)
672     {
673     xia = atom[ia].x;
674     yia = atom[ia].y;
675     zia = atom[ia].z;
676     xib = atom[ib].x;
677     yib = atom[ib].y;
678     zib = atom[ib].z;
679     xic = atom[ic].x;
680     yic = atom[ic].y;
681     zic = atom[ic].z;
682    
683     if ( (angles.angin[i] == FALSE) || field.type == MMFF94 )
684     {
685     xab = xia - xib;
686     yab = yia - yib;
687     zab = zia - zib;
688     rab = sqrt(xab*xab + yab*yab + zab*zab);
689     xcb = xic - xib;
690     ycb = yic - yib;
691     zcb = zic - zib;
692     rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
693     xp = ycb*zab - zcb*yab;
694     yp = zcb*xab - xcb*zab;
695     zp = xcb*yab - ycb*xab;
696     rp = sqrt(xp*xp + yp*yp + zp*zp);
697     if (rp != 0.0)
698     {
699     dot = xab*xcb + yab*ycb + zab*zcb;
700     cosine = dot / (rab*rcb);
701     if (cosine < -1.0)
702     cosine = -1.0;
703     if (cosine > 1.0)
704     cosine = 1.0;
705     angle = radian * acos(cosine);
706    
707     if (angles.angtype[i] == HARMONIC)
708     {
709     dt = angle - angles.anat[i];
710     dt2 = dt * dt;
711     dt3 = dt2 * dt;
712     dt4 = dt3 * dt;
713     deddt = units.angunit * angles.acon[i] * dt
714     * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
715     + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
716     d2eddt2 = units.angunit * angles.acon[i]
717     * (2.0 + 6.0*units.cang*dt + 12.0*units.qang*dt2
718     + 20.0*units.pang*dt3 + 30.0*units.sang*dt4);
719     terma = -radian/ (rab*rab*rp);
720     termc = radian / (rcb*rcb*rp);
721     } else
722     {
723     factor = angle/radian;
724     sine = sin(factor);
725     deddt = -units.angunit*angles.fcon[i]*
726     (angles.c1[i]*sine + 2*angles.c2[i]*sin(2*factor)
727     + 3*angles.c3[i]*sin(3*factor) + 4*angles.c4[i]*sin(4*factor)
728     + 5*angles.c5[i]*sin(5*factor) + 6*angles.c6[i]*sin(6*factor));
729     d2eddt2 = -units.angunit*angles.fcon[i]*
730     (angles.c1[i]*cosine + 4*angles.c2[i]*cos(2*factor)
731     + 9*angles.c3[i]*cos(3*factor) + 16*angles.c4[i]*cos(4*factor)
732     + 25*angles.c5[i]*cos(5*factor) + 36*angles.c6[i]*cos(6*factor));
733     terma = -1.0 / (rab*rab*rp);
734     termc = 1.0 / (rcb*rcb*rp);
735     }
736     ddtdxia = terma * (yab*zp-zab*yp);
737     ddtdyia = terma * (zab*xp-xab*zp);
738     ddtdzia = terma * (xab*yp-yab*xp);
739     ddtdxic = termc * (ycb*zp-zcb*yp);
740     ddtdyic = termc * (zcb*xp-xcb*zp);
741     ddtdzic = termc * (xcb*yp-ycb*xp);
742     ddtdxib = -ddtdxia - ddtdxic;
743     ddtdyib = -ddtdyia - ddtdyic;
744     ddtdzib = -ddtdzia - ddtdzic;
745     rab2 = 2.0 / (rab*rab);
746     xrab = xab * rab2;
747     yrab = yab * rab2;
748     zrab = zab * rab2;
749     rcb2 = 2.0 / (rcb*rcb);
750     xrcb = xcb * rcb2;
751     yrcb = ycb * rcb2;
752     zrcb = zcb * rcb2;
753     rp2 = 1.0 / (rp*rp);
754     xabp = (yab*zp-zab*yp) * rp2;
755     yabp = (zab*xp-xab*zp) * rp2;
756     zabp = (xab*yp-yab*xp) * rp2;
757     xcbp = (ycb*zp-zcb*yp) * rp2;
758     ycbp = (zcb*xp-xcb*zp) * rp2;
759     zcbp = (xcb*yp-ycb*xp) * rp2;
760    
761     dxiaxia = terma*(xab*xcb-dot) + ddtdxia*(xcbp-xrab);
762     dxiayia = terma*(zp+yab*xcb) + ddtdxia*(ycbp-yrab);
763     dxiazia = terma*(zab*xcb-yp) + ddtdxia*(zcbp-zrab);
764     dyiayia = terma*(yab*ycb-dot) + ddtdyia*(ycbp-yrab);
765     dyiazia = terma*(xp+zab*ycb) + ddtdyia*(zcbp-zrab);
766     dziazia = terma*(zab*zcb-dot) + ddtdzia*(zcbp-zrab);
767     dxicxic = termc*(dot-xab*xcb) - ddtdxic*(xabp+xrcb);
768     dxicyic = termc*(zp-ycb*xab) - ddtdxic*(yabp+yrcb);
769     dxiczic = -termc*(yp+zcb*xab) - ddtdxic*(zabp+zrcb);
770     dyicyic = termc*(dot-yab*ycb) - ddtdyic*(yabp+yrcb);
771     dyiczic = termc*(xp-zcb*yab) - ddtdyic*(zabp+zrcb);
772     dziczic = termc*(dot-zab*zcb) - ddtdzic*(zabp+zrcb);
773     dxiaxic = terma*(yab*yab+zab*zab) - ddtdxia*xabp;
774     dxiayic = -terma*xab*yab - ddtdxia*yabp;
775     dxiazic = -terma*xab*zab - ddtdxia*zabp;
776     dyiaxic = -terma*xab*yab - ddtdyia*xabp;
777     dyiayic = terma*(xab*xab+zab*zab) - ddtdyia*yabp;
778     dyiazic = -terma*yab*zab - ddtdyia*zabp;
779     dziaxic = -terma*xab*zab - ddtdzia*xabp;
780     dziayic = -terma*yab*zab - ddtdzia*yabp;
781     dziazic = terma*(xab*xab+yab*yab) - ddtdzia*zabp;
782    
783     dxibxia = -dxiaxia - dxiaxic;
784     dxibyia = -dxiayia - dyiaxic;
785     dxibzia = -dxiazia - dziaxic;
786     dyibxia = -dxiayia - dxiayic;
787     dyibyia = -dyiayia - dyiayic;
788     dyibzia = -dyiazia - dziayic;
789     dzibxia = -dxiazia - dxiazic;
790     dzibyia = -dyiazia - dyiazic;
791     dzibzia = -dziazia - dziazic;
792     dxibxic = -dxicxic - dxiaxic;
793     dxibyic = -dxicyic - dxiayic;
794     dxibzic = -dxiczic - dxiazic;
795     dyibxic = -dxicyic - dyiaxic;
796     dyibyic = -dyicyic - dyiayic;
797     dyibzic = -dyiczic - dyiazic;
798     dzibxic = -dxiczic - dziaxic;
799     dzibyic = -dyiczic - dziayic;
800     dzibzic = -dziczic - dziazic;
801     dxibxib = -dxibxia - dxibxic;
802     dxibyib = -dxibyia - dxibyic;
803     dxibzib = -dxibzia - dxibzic;
804     dyibyib = -dyibyia - dyibyic;
805     dyibzib = -dyibzia - dyibzic;
806     dzibzib = -dzibzia - dzibzic;
807    
808     if (ia == iatom)
809     {
810     hess.hessx[ia][0] += deddt*dxiaxia + d2eddt2*ddtdxia*ddtdxia;
811     hess.hessx[ia][1] += deddt*dxiayia + d2eddt2*ddtdxia*ddtdyia;
812     hess.hessx[ia][2] += deddt*dxiazia + d2eddt2*ddtdxia*ddtdzia;
813     hess.hessy[ia][0] += deddt*dxiayia + d2eddt2*ddtdyia*ddtdxia;
814     hess.hessy[ia][1] += deddt*dyiayia + d2eddt2*ddtdyia*ddtdyia;
815     hess.hessy[ia][2] += deddt*dyiazia + d2eddt2*ddtdyia*ddtdzia;
816     hess.hessz[ia][0] += deddt*dxiazia + d2eddt2*ddtdzia*ddtdxia;
817     hess.hessz[ia][1] += deddt*dyiazia + d2eddt2*ddtdzia*ddtdyia;
818     hess.hessz[ia][2] += deddt*dziazia + d2eddt2*ddtdzia*ddtdzia;
819     hess.hessx[ib][0] += deddt*dxibxia + d2eddt2*ddtdxia*ddtdxib;
820     hess.hessx[ib][1] += deddt*dyibxia + d2eddt2*ddtdxia*ddtdyib;
821     hess.hessx[ib][2] += deddt*dzibxia + d2eddt2*ddtdxia*ddtdzib;
822     hess.hessy[ib][0] += deddt*dxibyia + d2eddt2*ddtdyia*ddtdxib;
823     hess.hessy[ib][1] += deddt*dyibyia + d2eddt2*ddtdyia*ddtdyib;
824     hess.hessy[ib][2] += deddt*dzibyia + d2eddt2*ddtdyia*ddtdzib;
825     hess.hessz[ib][0] += deddt*dxibzia + d2eddt2*ddtdzia*ddtdxib;
826     hess.hessz[ib][1] += deddt*dyibzia + d2eddt2*ddtdzia*ddtdyib;
827     hess.hessz[ib][2] += deddt*dzibzia + d2eddt2*ddtdzia*ddtdzib;
828     hess.hessx[ic][0] += deddt*dxiaxic + d2eddt2*ddtdxia*ddtdxic;
829     hess.hessx[ic][1] += deddt*dxiayic + d2eddt2*ddtdxia*ddtdyic;
830     hess.hessx[ic][2] += deddt*dxiazic + d2eddt2*ddtdxia*ddtdzic;
831     hess.hessy[ic][0] += deddt*dyiaxic + d2eddt2*ddtdyia*ddtdxic;
832     hess.hessy[ic][1] += deddt*dyiayic + d2eddt2*ddtdyia*ddtdyic;
833     hess.hessy[ic][2] += deddt*dyiazic + d2eddt2*ddtdyia*ddtdzic;
834     hess.hessz[ic][0] += deddt*dziaxic + d2eddt2*ddtdzia*ddtdxic;
835     hess.hessz[ic][1] += deddt*dziayic + d2eddt2*ddtdzia*ddtdyic;
836     hess.hessz[ic][2] += deddt*dziazic + d2eddt2*ddtdzia*ddtdzic;
837     } else if (ib == iatom)
838     {
839     hess.hessx[ib][0] += deddt*dxibxib + d2eddt2*ddtdxib*ddtdxib;
840     hess.hessx[ib][1] += deddt*dxibyib + d2eddt2*ddtdxib*ddtdyib;
841     hess.hessx[ib][2] += deddt*dxibzib + d2eddt2*ddtdxib*ddtdzib;
842     hess.hessy[ib][0] += deddt*dxibyib + d2eddt2*ddtdyib*ddtdxib;
843     hess.hessy[ib][1] += deddt*dyibyib + d2eddt2*ddtdyib*ddtdyib;
844     hess.hessy[ib][2] += deddt*dyibzib + d2eddt2*ddtdyib*ddtdzib;
845     hess.hessz[ib][0] += deddt*dxibzib + d2eddt2*ddtdzib*ddtdxib;
846     hess.hessz[ib][1] += deddt*dyibzib + d2eddt2*ddtdzib*ddtdyib;
847     hess.hessz[ib][2] += deddt*dzibzib + d2eddt2*ddtdzib*ddtdzib;
848     hess.hessx[ia][0] += deddt*dxibxia + d2eddt2*ddtdxib*ddtdxia;
849     hess.hessx[ia][1] += deddt*dxibyia + d2eddt2*ddtdxib*ddtdyia;
850     hess.hessx[ia][2] += deddt*dxibzia + d2eddt2*ddtdxib*ddtdzia;
851     hess.hessy[ia][0] += deddt*dyibxia + d2eddt2*ddtdyib*ddtdxia;
852     hess.hessy[ia][1] += deddt*dyibyia + d2eddt2*ddtdyib*ddtdyia;
853     hess.hessy[ia][2] += deddt*dyibzia + d2eddt2*ddtdyib*ddtdzia;
854     hess.hessz[ia][0] += deddt*dzibxia + d2eddt2*ddtdzib*ddtdxia;
855     hess.hessz[ia][1] += deddt*dzibyia + d2eddt2*ddtdzib*ddtdyia;
856     hess.hessz[ia][2] += deddt*dzibzia + d2eddt2*ddtdzib*ddtdzia;
857     hess.hessx[ic][0] += deddt*dxibxic + d2eddt2*ddtdxib*ddtdxic;
858     hess.hessx[ic][1] += deddt*dxibyic + d2eddt2*ddtdxib*ddtdyic;
859     hess.hessx[ic][2] += deddt*dxibzic + d2eddt2*ddtdxib*ddtdzic;
860     hess.hessy[ic][0] += deddt*dyibxic + d2eddt2*ddtdyib*ddtdxic;
861     hess.hessy[ic][1] += deddt*dyibyic + d2eddt2*ddtdyib*ddtdyic;
862     hess.hessy[ic][2] += deddt*dyibzic + d2eddt2*ddtdyib*ddtdzic;
863     hess.hessz[ic][0] += deddt*dzibxic + d2eddt2*ddtdzib*ddtdxic;
864     hess.hessz[ic][1] += deddt*dzibyic + d2eddt2*ddtdzib*ddtdyic;
865     hess.hessz[ic][2] += deddt*dzibzic + d2eddt2*ddtdzib*ddtdzic;
866     }else if (ic == iatom)
867     {
868     hess.hessx[ic][0] += deddt*dxicxic + d2eddt2*ddtdxic*ddtdxic;
869     hess.hessx[ic][1] += deddt*dxicyic + d2eddt2*ddtdxic*ddtdyic;
870     hess.hessx[ic][2] += deddt*dxiczic + d2eddt2*ddtdxic*ddtdzic;
871     hess.hessy[ic][0] += deddt*dxicyic + d2eddt2*ddtdyic*ddtdxic;
872     hess.hessy[ic][1] += deddt*dyicyic + d2eddt2*ddtdyic*ddtdyic;
873     hess.hessy[ic][2] += deddt*dyiczic + d2eddt2*ddtdyic*ddtdzic;
874     hess.hessz[ic][0] += deddt*dxiczic + d2eddt2*ddtdzic*ddtdxic;
875     hess.hessz[ic][1] += deddt*dyiczic + d2eddt2*ddtdzic*ddtdyic;
876     hess.hessz[ic][2] += deddt*dziczic + d2eddt2*ddtdzic*ddtdzic;
877     hess.hessx[ib][0] += deddt*dxibxic + d2eddt2*ddtdxic*ddtdxib;
878     hess.hessx[ib][1] += deddt*dyibxic + d2eddt2*ddtdxic*ddtdyib;
879     hess.hessx[ib][2] += deddt*dzibxic + d2eddt2*ddtdxic*ddtdzib;
880     hess.hessy[ib][0] += deddt*dxibyic + d2eddt2*ddtdyic*ddtdxib;
881     hess.hessy[ib][1] += deddt*dyibyic + d2eddt2*ddtdyic*ddtdyib;
882     hess.hessy[ib][2] += deddt*dzibyic + d2eddt2*ddtdyic*ddtdzib;
883     hess.hessz[ib][0] += deddt*dxibzic + d2eddt2*ddtdzic*ddtdxib;
884     hess.hessz[ib][1] += deddt*dyibzic + d2eddt2*ddtdzic*ddtdyib;
885     hess.hessz[ib][2] += deddt*dzibzic + d2eddt2*ddtdzic*ddtdzib;
886     hess.hessx[ia][0] += deddt*dxiaxic + d2eddt2*ddtdxic*ddtdxia;
887     hess.hessx[ia][1] += deddt*dyiaxic + d2eddt2*ddtdxic*ddtdyia;
888     hess.hessx[ia][2] += deddt*dziaxic + d2eddt2*ddtdxic*ddtdzia;
889     hess.hessy[ia][0] += deddt*dxiayic + d2eddt2*ddtdyic*ddtdxia;
890     hess.hessy[ia][1] += deddt*dyiayic + d2eddt2*ddtdyic*ddtdyia;
891     hess.hessy[ia][2] += deddt*dziayic + d2eddt2*ddtdyic*ddtdzia;
892     hess.hessz[ia][0] += deddt*dxiazic + d2eddt2*ddtdzic*ddtdxia;
893     hess.hessz[ia][1] += deddt*dyiazic + d2eddt2*ddtdzic*ddtdyia;
894     hess.hessz[ia][2] += deddt*dziazic + d2eddt2*ddtdzic*ddtdzia;
895     }
896     }
897     }
898     }
899     }
900     }
901    
902     void eangle2b(int i)
903     {
904     int ia,ib,ic,id;
905     double angle,dot,cosine;
906     double dt,dt2,dt3,dt4;
907     double deddt,terma,termc,term;
908     double xia,yia,zia;
909     double xib,yib,zib;
910     double xic,yic,zic;
911     double xid,yid,zid;
912     double xad,yad,zad;
913     double xbd,ybd,zbd;
914     double xcd,ycd,zcd;
915     double xip,yip,zip;
916     double xap,yap,zap,rap2;
917     double xcp,ycp,zcp,rcp2;
918     double xt,yt,zt,rt2,ptrt2;
919     double xm,ym,zm,rm,delta,delta2;
920     double dedxia,dedyia,dedzia;
921     double dedxib,dedyib,dedzib;
922     double dedxic,dedyic,dedzic;
923     double dedxid,dedyid,dedzid;
924     double dedxip,dedyip,dedzip;
925     double dpdxia,dpdyia,dpdzia;
926     double dpdxic,dpdyic,dpdzic;
927    
928     ia = angles.i13[i][0];
929     ib = angles.i13[i][1];
930     ic = angles.i13[i][2];
931     id = angles.i13[i][3];
932     xia = atom[ia].x;
933     yia = atom[ia].y;
934     zia = atom[ia].z;
935     xib = atom[ib].x;
936     yib = atom[ib].y;
937     zib = atom[ib].z;
938     xic = atom[ic].x;
939     yic = atom[ic].y;
940     zic = atom[ic].z;
941     xid = atom[id].x;
942     yid = atom[id].y;
943     zid = atom[id].z;
944    
945     deriv.dea[ia][0] = 0.0;
946     deriv.dea[ia][1] = 0.0;
947     deriv.dea[ia][2] = 0.0;
948     deriv.dea[ib][0] = 0.0;
949     deriv.dea[ib][1] = 0.0;
950     deriv.dea[ib][2] = 0.0;
951     deriv.dea[ic][0] = 0.0;
952     deriv.dea[ic][1] = 0.0;
953     deriv.dea[ic][2] = 0.0;
954     deriv.dea[id][0] = 0.0;
955     deriv.dea[id][1] = 0.0;
956     deriv.dea[id][2] = 0.0;
957    
958     xad = xia - xid;
959     yad = yia - yid;
960     zad = zia - zid;
961     xbd = xib - xid;
962     ybd = yib - yid;
963     zbd = zib - zid;
964     xcd = xic - xid;
965     ycd = yic - yid;
966     zcd = zic - zid;
967     xt = yad*zcd - zad*ycd;
968     yt = zad*xcd - xad*zcd;
969     zt = xad*ycd - yad*xcd;
970     rt2 = xt*xt + yt*yt + zt*zt;
971     delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
972     xip = xib + xt*delta;
973     yip = yib + yt*delta;
974     zip = zib + zt*delta;
975     xap = xia - xip;
976     yap = yia - yip;
977     zap = zia - zip;
978     rap2 = xap*xap + yap*yap + zap*zap;
979     xcp = xic - xip;
980     ycp = yic - yip;
981     zcp = zic - zip;
982     rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
983     xm = ycp*zap - zcp*yap;
984     ym = zcp*xap - xcp*zap;
985     zm = xcp*yap - ycp*xap;
986     rm = sqrt(xm*xm + ym*ym + zm*zm);
987     if (rm != 0.0)
988     {
989     dot = xap*xcp + yap*ycp + zap*zcp;
990     cosine = dot / sqrt(rap2*rcp2);
991     if (cosine < -1.0)
992     cosine = -1.0;
993     if (cosine > 1.0)
994     cosine = 1.0;
995     angle = radian * acos(cosine);
996    
997     dt = angle - angles.anat[i];
998     dt2 = dt * dt;
999     dt3 = dt2 * dt;
1000     dt4 = dt2 * dt2;
1001     deddt = units.angunit * angles.acon[i] * dt
1002     * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
1003     + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
1004    
1005     deddt = deddt * radian;
1006     terma = -deddt / (rap2*rm);
1007     termc = deddt / (rcp2*rm);
1008     dedxia = terma * (yap*zm-zap*ym);
1009     dedyia = terma * (zap*xm-xap*zm);
1010     dedzia = terma * (xap*ym-yap*xm);
1011     dedxic = termc * (ycp*zm-zcp*ym);
1012     dedyic = termc * (zcp*xm-xcp*zm);
1013     dedzic = termc * (xcp*ym-ycp*xm);
1014     dedxip = -dedxia - dedxic;
1015     dedyip = -dedyia - dedyic;
1016     dedzip = -dedzia - dedzic;
1017    
1018     delta2 = 2.0 * delta;
1019     ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
1020     term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
1021     dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
1022     term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
1023     dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
1024     term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
1025     dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
1026     term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
1027     dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
1028     term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
1029     dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
1030     term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
1031     dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
1032     dedxia = dedxia + dpdxia;
1033     dedyia = dedyia + dpdyia;
1034     dedzia = dedzia + dpdzia;
1035     dedxib = dedxip;
1036     dedyib = dedyip;
1037     dedzib = dedzip;
1038     dedxic = dedxic + dpdxic;
1039     dedyic = dedyic + dpdyic;
1040     dedzic = dedzic + dpdzic;
1041     dedxid = -dedxia - dedxib - dedxic;
1042     dedyid = -dedyia - dedyib - dedyic;
1043     dedzid = -dedzia - dedzib - dedzic;
1044    
1045     deriv.dea[ia][0] = dedxia;
1046     deriv.dea[ia][1] = dedyia;
1047     deriv.dea[ia][2] = dedzia;
1048     deriv.dea[ib][0] = dedxib;
1049     deriv.dea[ib][1] = dedyib;
1050     deriv.dea[ib][2] = dedzib;
1051     deriv.dea[ic][0] = dedxic;
1052     deriv.dea[ic][1] = dedyic;
1053     deriv.dea[ic][2] = dedzic;
1054     deriv.dea[id][0] = dedxid;
1055     deriv.dea[id][1] = dedyid;
1056     deriv.dea[id][2] = dedzid;
1057     }
1058     }
1059    
1060     void eangle2a_fixed(int iatom)
1061     {
1062     int i,ia,ib,ic;
1063     double angle,dot,cosine;
1064     double dt,dt2,dt3,dt4;
1065     double xia,yia,zia;
1066     double xib,yib,zib;
1067     double xic,yic,zic;
1068     double xab,yab,zab,rab;
1069     double xcb,ycb,zcb,rcb;
1070     double xp,yp,zp,rp,rp2;
1071     double xrab,yrab,zrab,rab2;
1072     double xrcb,yrcb,zrcb,rcb2;
1073     double xabp,yabp,zabp;
1074     double xcbp,ycbp,zcbp;
1075     double deddt,d2eddt2,terma,termc;
1076     double ddtdxia,ddtdyia,ddtdzia;
1077     double ddtdxib,ddtdyib,ddtdzib;
1078     double ddtdxic,ddtdyic,ddtdzic;
1079     double dxiaxia,dxiayia,dxiazia;
1080     double dxibxib,dxibyib,dxibzib;
1081     double dxicxic,dxicyic,dxiczic;
1082     double dyiayia,dyiazia,dziazia;
1083     double dyibyib,dyibzib,dzibzib;
1084     double dyicyic,dyiczic,dziczic;
1085     double dxibxia,dxibyia,dxibzia;
1086     double dyibxia,dyibyia,dyibzia;
1087     double dzibxia,dzibyia,dzibzia;
1088     double dxibxic,dxibyic,dxibzic;
1089     double dyibxic,dyibyic,dyibzic;
1090     double dzibxic,dzibyic,dzibzic;
1091     double dxiaxic,dxiayic,dxiazic;
1092     double dyiaxic,dyiayic,dyiazic;
1093     double dziaxic,dziayic,dziazic;
1094    
1095     for (i=0; i < fixangle.nfxangle; i++)
1096     {
1097     ia = fixangle.ifxangle[i][0];
1098     ib = fixangle.ifxangle[i][1];
1099     ic = fixangle.ifxangle[i][2];
1100     if (iatom == ia || iatom == ib || iatom == ic)
1101     {
1102     xia = atom[ia].x;
1103     yia = atom[ia].y;
1104     zia = atom[ia].z;
1105     xib = atom[ib].x;
1106     yib = atom[ib].y;
1107     zib = atom[ib].z;
1108     xic = atom[ic].x;
1109     yic = atom[ic].y;
1110     zic = atom[ic].z;
1111    
1112     xab = xia - xib;
1113     yab = yia - yib;
1114     zab = zia - zib;
1115     rab = sqrt(xab*xab + yab*yab + zab*zab);
1116     xcb = xic - xib;
1117     ycb = yic - yib;
1118     zcb = zic - zib;
1119     rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
1120     xp = ycb*zab - zcb*yab;
1121     yp = zcb*xab - xcb*zab;
1122     zp = xcb*yab - ycb*xab;
1123     rp = sqrt(xp*xp + yp*yp + zp*zp);
1124     if (rp != 0.0)
1125     {
1126     dot = xab*xcb + yab*ycb + zab*zcb;
1127     cosine = dot / (rab*rcb);
1128     if (cosine < -1.0)
1129     cosine = -1.0;
1130     if (cosine > 1.0)
1131     cosine = 1.0;
1132     angle = radian * acos(cosine);
1133    
1134     dt = angle - fixangle.anat[i];
1135     dt2 = dt * dt;
1136     dt3 = dt2 * dt;
1137     dt4 = dt3 * dt;
1138     deddt = units.angunit * fixangle.acon[i] * dt*2.0;
1139     d2eddt2 = units.angunit * fixangle.acon[i]*2.0;
1140    
1141     terma = -radian / (rab*rab*rp);
1142     termc = radian / (rcb*rcb*rp);
1143     ddtdxia = terma * (yab*zp-zab*yp);
1144     ddtdyia = terma * (zab*xp-xab*zp);
1145     ddtdzia = terma * (xab*yp-yab*xp);
1146     ddtdxic = termc * (ycb*zp-zcb*yp);
1147     ddtdyic = termc * (zcb*xp-xcb*zp);
1148     ddtdzic = termc * (xcb*yp-ycb*xp);
1149     ddtdxib = -ddtdxia - ddtdxic;
1150     ddtdyib = -ddtdyia - ddtdyic;
1151     ddtdzib = -ddtdzia - ddtdzic;
1152     rab2 = 2.0 / (rab*rab);
1153     xrab = xab * rab2;
1154     yrab = yab * rab2;
1155     zrab = zab * rab2;
1156     rcb2 = 2.0 / (rcb*rcb);
1157     xrcb = xcb * rcb2;
1158     yrcb = ycb * rcb2;
1159     zrcb = zcb * rcb2;
1160     rp2 = 1.0 / (rp*rp);
1161     xabp = (yab*zp-zab*yp) * rp2;
1162     yabp = (zab*xp-xab*zp) * rp2;
1163     zabp = (xab*yp-yab*xp) * rp2;
1164     xcbp = (ycb*zp-zcb*yp) * rp2;
1165     ycbp = (zcb*xp-xcb*zp) * rp2;
1166     zcbp = (xcb*yp-ycb*xp) * rp2;
1167    
1168     dxiaxia = terma*(xab*xcb-dot) + ddtdxia*(xcbp-xrab);
1169     dxiayia = terma*(zp+yab*xcb) + ddtdxia*(ycbp-yrab);
1170     dxiazia = terma*(zab*xcb-yp) + ddtdxia*(zcbp-zrab);
1171     dyiayia = terma*(yab*ycb-dot) + ddtdyia*(ycbp-yrab);
1172     dyiazia = terma*(xp+zab*ycb) + ddtdyia*(zcbp-zrab);
1173     dziazia = terma*(zab*zcb-dot) + ddtdzia*(zcbp-zrab);
1174     dxicxic = termc*(dot-xab*xcb) - ddtdxic*(xabp+xrcb);
1175     dxicyic = termc*(zp-ycb*xab) - ddtdxic*(yabp+yrcb);
1176     dxiczic = -termc*(yp+zcb*xab) - ddtdxic*(zabp+zrcb);
1177     dyicyic = termc*(dot-yab*ycb) - ddtdyic*(yabp+yrcb);
1178     dyiczic = termc*(xp-zcb*yab) - ddtdyic*(zabp+zrcb);
1179     dziczic = termc*(dot-zab*zcb) - ddtdzic*(zabp+zrcb);
1180     dxiaxic = terma*(yab*yab+zab*zab) - ddtdxia*xabp;
1181     dxiayic = -terma*xab*yab - ddtdxia*yabp;
1182     dxiazic = -terma*xab*zab - ddtdxia*zabp;
1183     dyiaxic = -terma*xab*yab - ddtdyia*xabp;
1184     dyiayic = terma*(xab*xab+zab*zab) - ddtdyia*yabp;
1185     dyiazic = -terma*yab*zab - ddtdyia*zabp;
1186     dziaxic = -terma*xab*zab - ddtdzia*xabp;
1187     dziayic = -terma*yab*zab - ddtdzia*yabp;
1188     dziazic = terma*(xab*xab+yab*yab) - ddtdzia*zabp;
1189    
1190     dxibxia = -dxiaxia - dxiaxic;
1191     dxibyia = -dxiayia - dyiaxic;
1192     dxibzia = -dxiazia - dziaxic;
1193     dyibxia = -dxiayia - dxiayic;
1194     dyibyia = -dyiayia - dyiayic;
1195     dyibzia = -dyiazia - dziayic;
1196     dzibxia = -dxiazia - dxiazic;
1197     dzibyia = -dyiazia - dyiazic;
1198     dzibzia = -dziazia - dziazic;
1199     dxibxic = -dxicxic - dxiaxic;
1200     dxibyic = -dxicyic - dxiayic;
1201     dxibzic = -dxiczic - dxiazic;
1202     dyibxic = -dxicyic - dyiaxic;
1203     dyibyic = -dyicyic - dyiayic;
1204     dyibzic = -dyiczic - dyiazic;
1205     dzibxic = -dxiczic - dziaxic;
1206     dzibyic = -dyiczic - dziayic;
1207     dzibzic = -dziczic - dziazic;
1208     dxibxib = -dxibxia - dxibxic;
1209     dxibyib = -dxibyia - dxibyic;
1210     dxibzib = -dxibzia - dxibzic;
1211     dyibyib = -dyibyia - dyibyic;
1212     dyibzib = -dyibzia - dyibzic;
1213     dzibzib = -dzibzia - dzibzic;
1214    
1215     if (ia == iatom)
1216     {
1217     hess.hessx[ia][0] += deddt*dxiaxia + d2eddt2*ddtdxia*ddtdxia;
1218     hess.hessx[ia][1] += deddt*dxiayia + d2eddt2*ddtdxia*ddtdyia;
1219     hess.hessx[ia][2] += deddt*dxiazia + d2eddt2*ddtdxia*ddtdzia;
1220     hess.hessy[ia][0] += deddt*dxiayia + d2eddt2*ddtdyia*ddtdxia;
1221     hess.hessy[ia][1] += deddt*dyiayia + d2eddt2*ddtdyia*ddtdyia;
1222     hess.hessy[ia][2] += deddt*dyiazia + d2eddt2*ddtdyia*ddtdzia;
1223     hess.hessz[ia][0] += deddt*dxiazia + d2eddt2*ddtdzia*ddtdxia;
1224     hess.hessz[ia][1] += deddt*dyiazia + d2eddt2*ddtdzia*ddtdyia;
1225     hess.hessz[ia][2] += deddt*dziazia + d2eddt2*ddtdzia*ddtdzia;
1226     hess.hessx[ib][0] += deddt*dxibxia + d2eddt2*ddtdxia*ddtdxib;
1227     hess.hessx[ib][1] += deddt*dyibxia + d2eddt2*ddtdxia*ddtdyib;
1228     hess.hessx[ib][2] += deddt*dzibxia + d2eddt2*ddtdxia*ddtdzib;
1229     hess.hessy[ib][0] += deddt*dxibyia + d2eddt2*ddtdyia*ddtdxib;
1230     hess.hessy[ib][1] += deddt*dyibyia + d2eddt2*ddtdyia*ddtdyib;
1231     hess.hessy[ib][2] += deddt*dzibyia + d2eddt2*ddtdyia*ddtdzib;
1232     hess.hessz[ib][0] += deddt*dxibzia + d2eddt2*ddtdzia*ddtdxib;
1233     hess.hessz[ib][1] += deddt*dyibzia + d2eddt2*ddtdzia*ddtdyib;
1234     hess.hessz[ib][2] += deddt*dzibzia + d2eddt2*ddtdzia*ddtdzib;
1235     hess.hessx[ic][0] += deddt*dxiaxic + d2eddt2*ddtdxia*ddtdxic;
1236     hess.hessx[ic][1] += deddt*dxiayic + d2eddt2*ddtdxia*ddtdyic;
1237     hess.hessx[ic][2] += deddt*dxiazic + d2eddt2*ddtdxia*ddtdzic;
1238     hess.hessy[ic][0] += deddt*dyiaxic + d2eddt2*ddtdyia*ddtdxic;
1239     hess.hessy[ic][1] += deddt*dyiayic + d2eddt2*ddtdyia*ddtdyic;
1240     hess.hessy[ic][2] += deddt*dyiazic + d2eddt2*ddtdyia*ddtdzic;
1241     hess.hessz[ic][0] += deddt*dziaxic + d2eddt2*ddtdzia*ddtdxic;
1242     hess.hessz[ic][1] += deddt*dziayic + d2eddt2*ddtdzia*ddtdyic;
1243     hess.hessz[ic][2] += deddt*dziazic + d2eddt2*ddtdzia*ddtdzic;
1244     } else if (ib == iatom)
1245     {
1246     hess.hessx[ib][0] += deddt*dxibxib + d2eddt2*ddtdxib*ddtdxib;
1247     hess.hessx[ib][1] += deddt*dxibyib + d2eddt2*ddtdxib*ddtdyib;
1248     hess.hessx[ib][2] += deddt*dxibzib + d2eddt2*ddtdxib*ddtdzib;
1249     hess.hessy[ib][0] += deddt*dxibyib + d2eddt2*ddtdyib*ddtdxib;
1250     hess.hessy[ib][1] += deddt*dyibyib + d2eddt2*ddtdyib*ddtdyib;
1251     hess.hessy[ib][2] += deddt*dyibzib + d2eddt2*ddtdyib*ddtdzib;
1252     hess.hessz[ib][0] += deddt*dxibzib + d2eddt2*ddtdzib*ddtdxib;
1253     hess.hessz[ib][1] += deddt*dyibzib + d2eddt2*ddtdzib*ddtdyib;
1254     hess.hessz[ib][2] += deddt*dzibzib + d2eddt2*ddtdzib*ddtdzib;
1255     hess.hessx[ia][0] += deddt*dxibxia + d2eddt2*ddtdxib*ddtdxia;
1256     hess.hessx[ia][1] += deddt*dxibyia + d2eddt2*ddtdxib*ddtdyia;
1257     hess.hessx[ia][2] += deddt*dxibzia + d2eddt2*ddtdxib*ddtdzia;
1258     hess.hessy[ia][0] += deddt*dyibxia + d2eddt2*ddtdyib*ddtdxia;
1259     hess.hessy[ia][1] += deddt*dyibyia + d2eddt2*ddtdyib*ddtdyia;
1260     hess.hessy[ia][2] += deddt*dyibzia + d2eddt2*ddtdyib*ddtdzia;
1261     hess.hessz[ia][0] += deddt*dzibxia + d2eddt2*ddtdzib*ddtdxia;
1262     hess.hessz[ia][1] += deddt*dzibyia + d2eddt2*ddtdzib*ddtdyia;
1263     hess.hessz[ia][2] += deddt*dzibzia + d2eddt2*ddtdzib*ddtdzia;
1264     hess.hessx[ic][0] += deddt*dxibxic + d2eddt2*ddtdxib*ddtdxic;
1265     hess.hessx[ic][1] += deddt*dxibyic + d2eddt2*ddtdxib*ddtdyic;
1266     hess.hessx[ic][2] += deddt*dxibzic + d2eddt2*ddtdxib*ddtdzic;
1267     hess.hessy[ic][0] += deddt*dyibxic + d2eddt2*ddtdyib*ddtdxic;
1268     hess.hessy[ic][1] += deddt*dyibyic + d2eddt2*ddtdyib*ddtdyic;
1269     hess.hessy[ic][2] += deddt*dyibzic + d2eddt2*ddtdyib*ddtdzic;
1270     hess.hessz[ic][0] += deddt*dzibxic + d2eddt2*ddtdzib*ddtdxic;
1271     hess.hessz[ic][1] += deddt*dzibyic + d2eddt2*ddtdzib*ddtdyic;
1272     hess.hessz[ic][2] += deddt*dzibzic + d2eddt2*ddtdzib*ddtdzic;
1273     }else if (ic == iatom)
1274     {
1275     hess.hessx[ic][0] += deddt*dxicxic + d2eddt2*ddtdxic*ddtdxic;
1276     hess.hessx[ic][1] += deddt*dxicyic + d2eddt2*ddtdxic*ddtdyic;
1277     hess.hessx[ic][2] += deddt*dxiczic + d2eddt2*ddtdxic*ddtdzic;
1278     hess.hessy[ic][0] += deddt*dxicyic + d2eddt2*ddtdyic*ddtdxic;
1279     hess.hessy[ic][1] += deddt*dyicyic + d2eddt2*ddtdyic*ddtdyic;
1280     hess.hessy[ic][2] += deddt*dyiczic + d2eddt2*ddtdyic*ddtdzic;
1281     hess.hessz[ic][0] += deddt*dxiczic + d2eddt2*ddtdzic*ddtdxic;
1282     hess.hessz[ic][1] += deddt*dyiczic + d2eddt2*ddtdzic*ddtdyic;
1283     hess.hessz[ic][2] += deddt*dziczic + d2eddt2*ddtdzic*ddtdzic;
1284     hess.hessx[ib][0] += deddt*dxibxic + d2eddt2*ddtdxic*ddtdxib;
1285     hess.hessx[ib][1] += deddt*dyibxic + d2eddt2*ddtdxic*ddtdyib;
1286     hess.hessx[ib][2] += deddt*dzibxic + d2eddt2*ddtdxic*ddtdzib;
1287     hess.hessy[ib][0] += deddt*dxibyic + d2eddt2*ddtdyic*ddtdxib;
1288     hess.hessy[ib][1] += deddt*dyibyic + d2eddt2*ddtdyic*ddtdyib;
1289     hess.hessy[ib][2] += deddt*dzibyic + d2eddt2*ddtdyic*ddtdzib;
1290     hess.hessz[ib][0] += deddt*dxibzic + d2eddt2*ddtdzic*ddtdxib;
1291     hess.hessz[ib][1] += deddt*dyibzic + d2eddt2*ddtdzic*ddtdyib;
1292     hess.hessz[ib][2] += deddt*dzibzic + d2eddt2*ddtdzic*ddtdzib;
1293     hess.hessx[ia][0] += deddt*dxiaxic + d2eddt2*ddtdxic*ddtdxia;
1294     hess.hessx[ia][1] += deddt*dyiaxic + d2eddt2*ddtdxic*ddtdyia;
1295     hess.hessx[ia][2] += deddt*dziaxic + d2eddt2*ddtdxic*ddtdzia;
1296     hess.hessy[ia][0] += deddt*dxiayic + d2eddt2*ddtdyic*ddtdxia;
1297     hess.hessy[ia][1] += deddt*dyiayic + d2eddt2*ddtdyic*ddtdyia;
1298     hess.hessy[ia][2] += deddt*dziayic + d2eddt2*ddtdyic*ddtdzia;
1299     hess.hessz[ia][0] += deddt*dxiazic + d2eddt2*ddtdzic*ddtdxia;
1300     hess.hessz[ia][1] += deddt*dyiazic + d2eddt2*ddtdzic*ddtdyia;
1301     hess.hessz[ia][2] += deddt*dziazic + d2eddt2*ddtdzic*ddtdzia;
1302     }
1303     }
1304     }
1305     }
1306     }