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, 5 months ago) by wdelano
File size: 16704 byte(s)
Log Message:
branch created
Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5    
6     #include "energies.h"
7     #include "angles.h"
8     #include "derivs.h"
9     #include "hess.h"
10     #include "utility.h"
11    
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 wdelano 58 fprintf(pcmlogfile,"\nOut of Plane Bending Terms\n");
39     fprintf(pcmlogfile," At1 At2 At3 At4 Angle Oopconst Eoop\n");
40 tjod 3 }
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 wdelano 58 fprintf(pcmlogfile,"Oop: %-4d- %-4d- %-4d- %-4d %-8.3f %-8.3f = %-8.4f\n",ia,ib,ic,id,
105 tjod 3 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     }