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