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