ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/ebond.c
Revision: 99
Committed: Mon Jan 19 04:31:29 2009 UTC (10 years, 7 months ago) by wdelano
File size: 8024 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
Line User Rev File contents
1 wdelano 58 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5    
6     #include "energies.h"
7     #include "bonds_ff.h"
8     #include "derivs.h"
9     #include "hess.h"
10     #include "field.h"
11    
12    
13     EXTERN struct t_minim_values {
14     int iprint, ndc, nconst;
15     float dielc;
16     } minim_values;
17    
18     int find_bond(int,int);
19     void bnd_corr(int,int, int,int,double *);
20     double dihdrl(int,int,int,int);
21     double deltaks(double,double);
22    
23     void ebond()
24     {
25     /* compute stretching energy */
26    
27     int i, it, kt, itype, ktype, classi,classk;
28     double xr, yr, zr, rik, rik2, bcorr;
29     double dt, dt2, e;
30    
31     energies.estr = 0.0F;
32    
33     if (minim_values.iprint)
34     {
35     fprintf(pcmlogfile,"\nBond Terms \n");
36     fprintf(pcmlogfile," At1 At2 R BLen Bconst Eb\n");
37     }
38    
39     for (i=0; i < bonds_ff.nbnd; i++)
40     {
41     it = bonds_ff.i12[i][0];
42     kt = bonds_ff.i12[i][1];
43     itype = atom[it].type;
44     ktype = atom[kt].type;
45     classi = atom[it].tclass;
46     classk = atom[kt].tclass;
47     bcorr = 0.0;
48     if ( atom[it].use || atom[kt].use )
49     {
50     xr = atom[it].x - atom[kt].x;
51     yr = atom[it].y - atom[kt].y;
52     zr = atom[it].z - atom[kt].z;
53     rik2 = xr*xr + yr*yr + zr*zr;
54     rik = sqrt(rik2);
55     if (field.type == MM3 && ( (itype ==1 && ktype == 6) || (itype == 6 && ktype ==1) ||
56     (classi == 1 && classk == 6) || (classi == 6 && classk == 1) ) )
57     bnd_corr(it,kt, itype, ktype, &bcorr);
58    
59     dt = rik - (bonds_ff.bl[i] - bcorr);
60     dt2 = dt*dt;
61     e = units.bndunit *bonds_ff.bk[i]*dt2*
62     (1.0 + units.cbnd*dt + units.qbnd*dt2);
63     energies.estr += e;
64     atom[it].energy += e;
65     atom[kt].energy += e;
66     if (minim_values.iprint)
67     fprintf(pcmlogfile,"Bond: %2s(%-3d) - %2s(%-3d) %-8.3f %-8.3f %-8.3f = %-8.4f\n",atom[it].name,it,atom[kt].name
68     ,kt, rik, bonds_ff.bl[i]-bcorr, bonds_ff.bk[i],e);
69     }
70     }
71     }
72    
73     void ebond1()
74     {
75     /* compute stretching energy and first derivatives */
76     int i, it, kt, itype, ktype, classi,classk;
77     double xr, yr, zr, rik, rik2, bcorr;
78     double dt, dt2, e, deddt;
79     double de,dedx,dedy,dedz;
80    
81     energies.estr = 0.0F;
82     for (i=0; i <= natom; i++)
83     {
84     deriv.deb[i][0] = 0.0;
85     deriv.deb[i][1] = 0.0;
86     deriv.deb[i][2] = 0.0;
87     }
88    
89     for (i=0; i < bonds_ff.nbnd; i++)
90     {
91     it = bonds_ff.i12[i][0];
92     kt = bonds_ff.i12[i][1];
93     itype = atom[it].type;
94     ktype = atom[kt].type;
95     classi = atom[it].tclass;
96     classk = atom[kt].tclass;
97     bcorr = 0.0;
98     if ( atom[it].use || atom[kt].use )
99     {
100     xr = atom[it].x - atom[kt].x;
101     yr = atom[it].y - atom[kt].y;
102     zr = atom[it].z - atom[kt].z;
103     rik2 = xr*xr + yr*yr + zr*zr;
104     rik = sqrt(rik2);
105    
106     if (field.type == MM3 && ( (itype ==1 && ktype == 6) || (itype == 6 && ktype ==1) ||
107     (classi == 1 && classk == 6) || (classi == 6 && classk == 1) ) )
108     bnd_corr(it,kt, itype, ktype, &bcorr);
109     dt = rik - (bonds_ff.bl[i] - bcorr);
110    
111     dt2 = dt*dt;
112     e = units.bndunit *bonds_ff.bk[i]*dt2*
113     (1.0 + units.cbnd*dt + units.qbnd*dt2);
114     deddt = 2.0 * units.bndunit * bonds_ff.bk[i] * dt
115     * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
116     if (rik == 0.0)
117     de = 0.0;
118     else
119     de = deddt/rik;
120     dedx = de*xr;
121     dedy = de*yr;
122     dedz = de*zr;
123     energies.estr += e;
124     deriv.deb[it][0] += dedx;
125     deriv.deb[it][1] += dedy;
126     deriv.deb[it][2] += dedz;
127     deriv.deb[kt][0] -= dedx;
128     deriv.deb[kt][1] -= dedy;
129     deriv.deb[kt][2] -= dedz;
130    
131     }
132     }
133     }
134     // =============================================
135     void ebond2(int ia)
136     {
137     int j,k,m,ibond;
138     double dt,dt2;
139     double xr,yr,zr,rik,rik2,deddt,d2eddt2;
140     double de,term,termx,termy,termz,d2e[3][3];
141    
142     for (m=0; m < MAXIAT; m++)
143     {
144     if (atom[ia].iat[m]!= 0 && atom[ia].bo[m] != 9)
145     {
146     ibond = find_bond(ia, atom[ia].iat[m]);
147     if (bonds_ff.i12[ibond][0] == ia)
148     k = bonds_ff.i12[ibond][1];
149     else
150     k = bonds_ff.i12[ibond][0];
151    
152     xr = atom[ia].x - atom[k].x;
153     yr = atom[ia].y - atom[k].y;
154     zr = atom[ia].z - atom[k].z;
155     rik2 = xr*xr + yr*yr + zr*zr;
156     rik = sqrt(rik2);
157     dt = rik - bonds_ff.bl[ibond];
158     dt2 = dt * dt;
159    
160     deddt = 2.0 * units.bndunit * bonds_ff.bk[ibond] * dt
161     * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
162     d2eddt2 = 2.0 * units.bndunit * bonds_ff.bk[ibond]
163     * (1.0+3.0*units.cbnd*dt+6.0*units.qbnd*dt2);
164    
165     if (rik2 == 0.0)
166     {
167     de = 0.0;
168     term = 0.0;
169     }else
170     {
171     de = deddt / rik;
172     term = (d2eddt2-de) / rik2;
173     }
174    
175     termx = term * xr;
176     termy = term * yr;
177     termz = term * zr;
178     d2e[0][0] = termx*xr + de;
179     d2e[1][0] = termx*yr;
180     d2e[2][0] = termx*zr;
181     d2e[0][1] = d2e[1][0];
182     d2e[1][1] = termy*yr + de;
183     d2e[2][1] = termy*zr;
184     d2e[0][2] = d2e[2][0];
185     d2e[1][2] = d2e[2][1];
186     d2e[2][2] = termz*zr + de;
187    
188     for (j=0; j < 3; j++)
189     {
190     hess.hessx[ia][j] += d2e[j][0];
191     hess.hessy[ia][j] += d2e[j][1];
192     hess.hessz[ia][j] += d2e[j][2];
193     hess.hessx[k][j] -= d2e[j][0];
194     hess.hessy[k][j] -= d2e[j][1];
195     hess.hessz[k][j] -= d2e[j][2];
196     }
197     }
198     }
199     }
200     /* ------------------------------------------- */
201     void bnd_corr(int ia, int ib, int itype, int ktype, double *bcorr)
202     {
203     // two terms - first O-ia-ib-R
204     // then X-O-ia-ib
205     // it = carbon, kt = oxygen
206    
207     int i, j,k, iatt, jatt, latt,iatype, jatype, iox;
208     int it, kt;
209     double angle;
210    
211     *bcorr = 0.0;
212    
213     if (itype == 1)
214     {
215     it = ia;
216     kt = ib;
217     } else
218     {
219     it = ib;
220     kt = ia;
221     }
222    
223     for (i=0; i < MAXIAT; i++)
224     {
225     if (atom[kt].iat[i] != 0 && atom[kt].iat[i] != it)
226     {
227     iatt = atom[kt].iat[i];
228     iatype = atom[iatt].type;
229     for(j=0; j < MAXIAT; j++)
230     {
231     iox = 0;
232     if (atom[it].iat[j] != 0 && atom[it].iat[j] != kt)
233     {
234     jatt = atom[it].iat[j];
235     jatype = atom[jatt].type;
236    
237     if (itype == 6 && jatype == 6) iox = 1;
238     if (ktype == 6 && iatype == 6) iox = 2;
239     if (iox != 0)
240     {
241     angle = dihdrl(iatt,kt,it,jatt)/radian;
242     *bcorr += 0.5* 0.01*(1.0-cos(2.0*angle));
243     if (iox == 1)
244     {
245     for (k = 0 ; k < MAXIAT; k++)
246     {
247     if (atom[jatt].iat[k] != 0 && atom[jatt].iat[k] != it)
248     {
249     latt = atom[jatt].iat[k];
250     angle = dihdrl(kt,it,jatt,latt)/radian;
251     *bcorr += 0.5*0.01*(1.0 - cos(2.0*angle))*(-1.30)+0.005;
252     break;
253     }
254     }
255     }
256     }
257     }
258     }
259     }
260     }
261     }
262