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

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