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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines