ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/eopbend.c
Revision: 93
Committed: Sat Jan 17 23:24:55 2009 UTC (13 years, 4 months ago) by wdelano
File size: 16704 byte(s)
Log Message:
branch created
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
12 void eopbend2b(int);
13
14 EXTERN struct t_minim_values {
15 int iprint, ndc, nconst;
16 float dielc;
17 } minim_values;
18
19 void eopbend()
20 {
21 int i,ia,ib,ic,id;
22 double e,angle,dot,cosine;
23 double dt,dt2,dt3,dt4;
24 double xia,yia,zia;
25 double xib,yib,zib;
26 double xic,yic,zic;
27 double xid,yid,zid;
28 double xad,yad,zad;
29 double xbd,ybd,zbd,rbd2;
30 double xcd,ycd,zcd;
31 double xpd,ypd,zpd,rpd2;
32 double xt,yt,zt,rt2,delta;
33
34 energies.eopb = 0.0;
35
36 if (minim_values.iprint)
37 {
38 fprintf(pcmlogfile,"\nOut of Plane Bending Terms\n");
39 fprintf(pcmlogfile," At1 At2 At3 At4 Angle Oopconst Eoop\n");
40 }
41
42 for (i=0; i < angles.nang; i++)
43 {
44 if (angles.angin[i] == TRUE)
45 {
46 ia = angles.i13[i][0];
47 ib = angles.i13[i][1];
48 ic = angles.i13[i][2];
49 id = angles.i13[i][3];
50 if (atom[ia].use || atom[ib].use || atom[ic].use || atom[id].use)
51 {
52 xia = atom[ia].x;
53 yia = atom[ia].y;
54 zia = atom[ia].z;
55 xib = atom[ib].x;
56 yib = atom[ib].y;
57 zib = atom[ib].z;
58 xic = atom[ic].x;
59 yic = atom[ic].y;
60 zic = atom[ic].z;
61 xid = atom[id].x;
62 yid = atom[id].y;
63 zid = atom[id].z;
64 xad = xia - xid;
65 yad = yia - yid;
66 zad = zia - zid;
67 xbd = xib - xid;
68 ybd = yib - yid;
69 zbd = zib - zid;
70 xcd = xic - xid;
71 ycd = yic - yid;
72 zcd = zic - zid;
73 xt = yad*zcd - zad*ycd;
74 yt = zad*xcd - xad*zcd;
75 zt = xad*ycd - yad*xcd;
76 rt2 = xt*xt + yt*yt + zt*zt;
77 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
78 xpd = xbd + xt*delta;
79 ypd = ybd + yt*delta;
80 zpd = zbd + zt*delta;
81 rbd2 = xbd*xbd + ybd*ybd + zbd*zbd;
82 rpd2 = xpd*xpd + ypd*ypd + zpd*zpd;
83 if (rbd2 != 0.0 && rpd2 != 0.0)
84 {
85 dot = xbd*xpd + ybd*ypd + zbd*zpd;
86 cosine = dot / sqrt(rbd2*rpd2);
87 if (cosine < -1.0)
88 cosine = -1.0;
89 if (cosine > 1.0)
90 cosine = 1.0;
91 angle = radian * acos(cosine);
92 dt = angle;
93 dt2 = dt * dt;
94 dt3 = dt2 * dt;
95 dt4 = dt2 * dt2;
96 e = units.angunit * angles.copb[i] * dt2
97 * (1.0+units.cang*dt+units.qang*dt2+units.pang*dt3+units.sang*dt4);
98 energies.eopb += e;
99 atom[ia].energy += e;
100 atom[ib].energy += e;
101 atom[ic].energy += e;
102 atom[id].energy += e;
103 if (minim_values.iprint)
104 fprintf(pcmlogfile,"Oop: %-4d- %-4d- %-4d- %-4d %-8.3f %-8.3f = %-8.4f\n",ia,ib,ic,id,
105 angle, angles.copb[i],e);
106 }
107
108 }
109 }
110 }
111 }
112
113 void eopbend1()
114 {
115 int i,ia,ib,ic,id;
116 double e,angle,dot,cosine;
117 double dt,dt2,dt3,dt4;
118 double deddt,termb,termp,term;
119 double xia,yia,zia;
120 double xib,yib,zib;
121 double xic,yic,zic;
122 double xid,yid,zid;
123 double xad,yad,zad;
124 double xbd,ybd,zbd,rbd2;
125 double xcd,ycd,zcd;
126 double xpd,ypd,zpd,rpd2;
127 double xt,yt,zt,rt2,ptrt2;
128 double xm,ym,zm,rm,delta,delta2;
129 double dedxia,dedyia,dedzia;
130 double dedxib,dedyib,dedzib;
131 double dedxic,dedyic,dedzic;
132 double dedxid,dedyid,dedzid;
133 double dedxip,dedyip,dedzip;
134 double dpdxia,dpdyia,dpdzia;
135 double dpdxic,dpdyic,dpdzic;
136
137 energies.eopb = 0.0;
138 for (i=0; i <= natom; i++)
139 {
140 deriv.deopb[i][0] = 0.0;
141 deriv.deopb[i][1] = 0.0;
142 deriv.deopb[i][2] = 0.0;
143 }
144
145 for (i=0; i < angles.nang; i++)
146 {
147 if (angles.angin[i] == TRUE)
148 {
149 ia = angles.i13[i][0];
150 ib = angles.i13[i][1];
151 ic = angles.i13[i][2];
152 id = angles.i13[i][3];
153 if (atom[ia].use || atom[ib].use || atom[ic].use || atom[id].use)
154 {
155 xia = atom[ia].x;
156 yia = atom[ia].y;
157 zia = atom[ia].z;
158 xib = atom[ib].x;
159 yib = atom[ib].y;
160 zib = atom[ib].z;
161 xic = atom[ic].x;
162 yic = atom[ic].y;
163 zic = atom[ic].z;
164 xid = atom[id].x;
165 yid = atom[id].y;
166 zid = atom[id].z;
167 xad = xia - xid;
168 yad = yia - yid;
169 zad = zia - zid;
170 xbd = xib - xid;
171 ybd = yib - yid;
172 zbd = zib - zid;
173 xcd = xic - xid;
174 ycd = yic - yid;
175 zcd = zic - zid;
176 xt = yad*zcd - zad*ycd;
177 yt = zad*xcd - xad*zcd;
178 zt = xad*ycd - yad*xcd;
179 rt2 = xt*xt + yt*yt + zt*zt;
180 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
181 xpd = xbd + xt*delta;
182 ypd = ybd + yt*delta;
183 zpd = zbd + zt*delta;
184 rbd2 = xbd*xbd + ybd*ybd + zbd*zbd;
185 rpd2 = xpd*xpd + ypd*ypd + zpd*zpd;
186 xm = ypd*zbd - zpd*ybd;
187 ym = zpd*xbd - xpd*zbd;
188 zm = xpd*ybd - ypd*xbd;
189 rm = sqrt(xm*xm + ym*ym + zm*zm);
190 if (rm != 0.0)
191 {
192 dot = xbd*xpd + ybd*ypd + zbd*zpd;
193 cosine = dot / sqrt(rbd2*rpd2);
194 if (cosine < -1.0)
195 cosine = -1.0;
196 if (cosine > 1.0)
197 cosine = 1.0;
198 angle = radian * acos(cosine);
199 dt = angle;
200 dt2 = dt * dt;
201 dt3 = dt2 * dt;
202 dt4 = dt2 * dt2;
203 e = units.angunit * angles.copb[i] * dt2
204 * (1.0+units.cang*dt+units.qang*dt2+units.pang*dt3+units.sang*dt4);
205
206 deddt = units.angunit * angles.copb[i] * dt
207 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
208 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
209 deddt = deddt * radian;
210 termb = -deddt / (rbd2*rm);
211 termp = deddt / (rpd2*rm);
212 dedxib = termb * (ybd*zm-zbd*ym);
213 dedyib = termb * (zbd*xm-xbd*zm);
214 dedzib = termb * (xbd*ym-ybd*xm);
215 dedxip = termp * (ypd*zm-zpd*ym);
216 dedyip = termp * (zpd*xm-xpd*zm);
217 dedzip = termp * (xpd*ym-ypd*xm);
218
219 delta2 = 2.0 * delta;
220 ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
221 term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
222 dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
223 term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
224 dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
225 term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
226 dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
227 term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
228 dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
229 term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
230 dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
231 term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
232 dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
233
234 dedxia = dpdxia;
235 dedyia = dpdyia;
236 dedzia = dpdzia;
237 dedxic = dpdxic;
238 dedyic = dpdyic;
239 dedzic = dpdzic;
240 dedxid = -dedxia - dedxib - dedxic;
241 dedyid = -dedyia - dedyib - dedyic;
242 dedzid = -dedzia - dedzib - dedzic;
243
244 energies.eopb += e;
245 deriv.deopb[ia][0] += dedxia;
246 deriv.deopb[ia][1] += dedyia;
247 deriv.deopb[ia][2] += dedzia;
248
249 deriv.deopb[ib][0] += dedxib;
250 deriv.deopb[ib][1] += dedyib;
251 deriv.deopb[ib][2] += dedzib;
252
253 deriv.deopb[ic][0] += dedxic;
254 deriv.deopb[ic][1] += dedyic;
255 deriv.deopb[ic][2] += dedzic;
256
257 deriv.deopb[id][0] += dedxid;
258 deriv.deopb[id][1] += dedyid;
259 deriv.deopb[id][2] += dedzid;
260 virial.virx += xad*dedxia + xbd*dedxib + xcd*dedxic;
261 virial.viry += yad*dedyia + ybd*dedyib + ycd*dedyic;
262 virial.virz += zad*dedzia + zbd*dedzib + zcd*dedzic;
263 }
264
265 }
266 }
267 }
268 }
269
270 void eopbend2(int i)
271 {
272 int j,k;
273 int ia,ib,ic,id;
274 double old,eps;
275 double **d0; // d0[MAXATOM][3];
276
277 d0 = dmatrix(0,natom,0,3);
278 eps = 1.0e-7;
279
280 for (k=0; k < angles.nang; k++)
281 {
282 if (angles.angin[k] == TRUE)
283 {
284 ia = angles.i13[k][0];
285 ib = angles.i13[k][1];
286 ic = angles.i13[k][2];
287 id = angles.i13[k][3];
288 if (ia == i || ib == i || ic == i || id == i)
289 {
290 eopbend2b(k);
291 for (j = 0; j < 3; j++)
292 {
293 d0[ia][j] = deriv.deopb[ia][j];
294 d0[ib][j] = deriv.deopb[ib][j];
295 d0[ic][j] = deriv.deopb[ic][j];
296 d0[id][j] = deriv.deopb[id][j];
297 }
298 // numerical x
299 old = atom[i].x;
300 atom[i].x += eps;
301 eopbend2b(k);
302 atom[i].x = old;
303 for (j=0; j < 3; j++)
304 {
305 hess.hessx[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
306 hess.hessx[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
307 hess.hessx[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
308 hess.hessx[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
309 }
310 // numerical y
311 old = atom[i].y;
312 atom[i].y += eps;
313 eopbend2b(k);
314 atom[i].y = old;
315 for (j=0; j < 3; j++)
316 {
317 hess.hessy[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
318 hess.hessy[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
319 hess.hessy[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
320 hess.hessy[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
321 }
322 // numerical z
323 old = atom[i].z;
324 atom[i].z += eps;
325 eopbend2b(k);
326 atom[i].z = old;
327 for (j=0; j < 3; j++)
328 {
329 hess.hessz[ia][j] += (deriv.deopb[ia][j]-d0[ia][j])/eps;
330 hess.hessz[ib][j] += (deriv.deopb[ib][j]-d0[ib][j])/eps;
331 hess.hessz[ic][j] += (deriv.deopb[ic][j]-d0[ic][j])/eps;
332 hess.hessz[id][j] += (deriv.deopb[id][j]-d0[id][j])/eps;
333 }
334 }
335 }
336 }
337 free_dmatrix( d0 ,0,natom,0,3);
338 }
339
340 void eopbend2b(int k)
341 {
342 int j,ia,ib,ic,id;
343 double angle,dot,cosine;
344 double dt,dt2,dt3,dt4;
345 double deddt,termb,termp,term;
346 double xia,yia,zia;
347 double xib,yib,zib;
348 double xic,yic,zic;
349 double xid,yid,zid;
350 double xad,yad,zad;
351 double xbd,ybd,zbd,rbd2;
352 double xcd,ycd,zcd;
353 double xpd,ypd,zpd,rpd2;
354 double xt,yt,zt,rt2,ptrt2;
355 double xm,ym,zm,rm,delta,delta2;
356 double dedxia,dedyia,dedzia;
357 double dedxib,dedyib,dedzib;
358 double dedxic,dedyic,dedzic;
359 double dedxid,dedyid,dedzid;
360 double dedxip,dedyip,dedzip;
361 double dpdxia,dpdyia,dpdzia;
362 double dpdxic,dpdyic,dpdzic;
363
364 ia = angles.i13[k][0];
365 ib = angles.i13[k][1];
366 ic = angles.i13[k][2];
367 id = angles.i13[k][3];
368
369 xia = atom[ia].x;
370 yia = atom[ia].y;
371 zia = atom[ia].z;
372 xib = atom[ib].x;
373 yib = atom[ib].y;
374 zib = atom[ib].z;
375 xic = atom[ic].x;
376 yic = atom[ic].y;
377 zic = atom[ic].z;
378 xid = atom[id].x;
379 yid = atom[id].y;
380 zid = atom[id].z;
381
382 for (j=0; j < 3; j++)
383 {
384 deriv.deopb[ia][j] = 0.0;
385 deriv.deopb[ib][j] = 0.0;
386 deriv.deopb[ic][j] = 0.0;
387 deriv.deopb[id][j] = 0.0;
388 }
389
390 xad = xia - xid;
391 yad = yia - yid;
392 zad = zia - zid;
393 xbd = xib - xid;
394 ybd = yib - yid;
395 zbd = zib - zid;
396 xcd = xic - xid;
397 ycd = yic - yid;
398 zcd = zic - zid;
399 xt = yad*zcd - zad*ycd;
400 yt = zad*xcd - xad*zcd;
401 zt = xad*ycd - yad*xcd;
402 rt2 = xt*xt + yt*yt + zt*zt;
403 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
404 xpd = xbd + xt*delta;
405 ypd = ybd + yt*delta;
406 zpd = zbd + zt*delta;
407 rbd2 = xbd*xbd + ybd*ybd + zbd*zbd;
408 rpd2 = xpd*xpd + ypd*ypd + zpd*zpd;
409 xm = ypd*zbd - zpd*ybd;
410 ym = zpd*xbd - xpd*zbd;
411 zm = xpd*ybd - ypd*xbd;
412 rm = sqrt(xm*xm + ym*ym + zm*zm);
413 if (rm != 0.0)
414 {
415 dot = xbd*xpd + ybd*ypd + zbd*zpd;
416 cosine = dot / sqrt(rbd2*rpd2);
417 if (cosine < -1.0)
418 cosine = -1.0;
419 if (cosine > 1.0)
420 cosine = 1.0;
421 angle = radian * acos(cosine);
422 dt = angle;
423 dt2 = dt * dt;
424 dt3 = dt2 * dt;
425 dt4 = dt2 * dt2;
426 deddt = units.angunit * angles.copb[k] * dt
427 * (2.0 + 3.0*units.cang*dt + 4.0*units.qang*dt2
428 + 5.0*units.pang*dt3 + 6.0*units.sang*dt4);
429
430 deddt = deddt * radian;
431 termb = -deddt / (rbd2*rm);
432 termp = deddt / (rpd2*rm);
433 dedxib = termb * (ybd*zm-zbd*ym);
434 dedyib = termb * (zbd*xm-xbd*zm);
435 dedzib = termb * (xbd*ym-ybd*xm);
436 dedxip = termp * (ypd*zm-zpd*ym);
437 dedyip = termp * (zpd*xm-xpd*zm);
438 dedzip = termp * (xpd*ym-ypd*xm);
439
440 delta2 = 2.0 * delta;
441 ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
442 term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
443 dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
444 term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
445 dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
446 term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
447 dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
448 term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
449 dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
450 term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
451 dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
452 term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
453 dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
454
455 dedxia = dpdxia;
456 dedyia = dpdyia;
457 dedzia = dpdzia;
458 dedxic = dpdxic;
459 dedyic = dpdyic;
460 dedzic = dpdzic;
461 dedxid = -dedxia - dedxib - dedxic;
462 dedyid = -dedyia - dedyib - dedyic;
463 dedzid = -dedzia - dedzib - dedzic;
464
465 deriv.deopb[ia][0] = dedxia;
466 deriv.deopb[ia][1] = dedyia;
467 deriv.deopb[ia][2] = dedzia;
468 deriv.deopb[ib][0] = dedxib;
469 deriv.deopb[ib][1] = dedyib;
470 deriv.deopb[ib][2] = dedzib;
471 deriv.deopb[ic][0] = dedxic;
472 deriv.deopb[ic][1] = dedyic;
473 deriv.deopb[ic][2] = dedzic;
474 deriv.deopb[id][0] = dedxid;
475 deriv.deopb[id][1] = dedyid;
476 deriv.deopb[id][2] = dedzid;
477 }
478 }