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

Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5
6 #include "energies.h"
7 #include "angles.h"
8 #include "derivs.h"
9 #include "hess.h"
10 #include "utility.h"
11 #include "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 }