ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/ebond.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 11 months ago) by tjod
File size: 12371 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 "bonds_ff.h"
8     #include "derivs.h"
9     #include "hess.h"
10     #include "field.h"
11     #include "fix.h"
12    
13     EXTERN struct t_units {
14     double bndunit, cbnd, qbnd;
15     double angunit, cang, qang, pang, sang, aaunit;
16     double stbnunit, ureyunit, torsunit, storunit, v14scale;
17     double aterm, bterm, cterm, dielec, chgscale;
18     } units;
19    
20     EXTERN struct t_minim_values {
21     int iprint, ndc, nconst;
22     float dielc;
23     } minim_values;
24    
25     int find_bond(int,int);
26     int find_fixed_bond(int,int);
27     void bnd_corr(int,int, int,int,double *);
28     double dihdrl(int,int,int,int);
29     double deltaks(double,double);
30    
31     void ebond()
32     {
33     /* compute stretching energy */
34    
35     int i, it, kt, itype, ktype, classi,classk;
36     double xr, yr, zr, rik, rik2, bcorr;
37     double dt, dt2, e;
38    
39     energies.estr = 0.0F;
40    
41     if (minim_values.iprint)
42     {
43     fprintf(pcmoutfile,"\nBond Terms \n");
44     fprintf(pcmoutfile," At1 At2 R BLen Bconst Eb\n");
45     }
46    
47     for (i=0; i < bonds_ff.nbnd; i++)
48     {
49     it = bonds_ff.i12[i][0];
50     kt = bonds_ff.i12[i][1];
51     itype = atom[it].type;
52     ktype = atom[kt].type;
53     classi = atom[it].tclass;
54     classk = atom[kt].tclass;
55     bcorr = 0.0;
56     if ( atom[it].use || atom[kt].use )
57     {
58     xr = atom[it].x - atom[kt].x;
59     yr = atom[it].y - atom[kt].y;
60     zr = atom[it].z - atom[kt].z;
61     rik2 = xr*xr + yr*yr + zr*zr;
62     rik = sqrt(rik2);
63     if (field.type == MM3 && ( (itype ==1 && ktype == 6) || (itype == 6 && ktype ==1) ||
64     (classi == 1 && classk == 6) || (classi == 6 && classk == 1) ) )
65     bnd_corr(it,kt, itype, ktype, &bcorr);
66    
67     dt = rik - (bonds_ff.bl[i] - bcorr);
68     dt2 = dt*dt;
69     e = units.bndunit *bonds_ff.bk[i]*dt2*
70     (1.0 + units.cbnd*dt + units.qbnd*dt2);
71     energies.estr += e;
72     atom[it].energy += e;
73     atom[kt].energy += e;
74     if (minim_values.iprint)
75     fprintf(pcmoutfile,"Bond: %2s(%-3d) - %2s(%-3d) %-8.3f %-8.3f %-8.3f = %-8.4f\n",atom[it].name,it,atom[kt].name
76     ,kt, rik, bonds_ff.bl[i]-bcorr, bonds_ff.bk[i],e);
77     }
78     }
79     // fixed distances
80     for (i=0; i < fixdis.nfxstr; i++)
81     {
82     it = fixdis.ifxstr[i][0];
83     kt = fixdis.ifxstr[i][1];
84     itype = atom[it].type;
85     ktype = atom[kt].type;
86     bcorr = 0.0;
87     if ( atom[it].use || atom[kt].use )
88     {
89     xr = atom[it].x - atom[kt].x;
90     yr = atom[it].y - atom[kt].y;
91     zr = atom[it].z - atom[kt].z;
92     rik2 = xr*xr + yr*yr + zr*zr;
93     rik = sqrt(rik2);
94     dt = rik - fixdis.fxstrd[i];
95    
96     dt2 = dt*dt;
97     e = units.bndunit *fixdis.fxstrc[i]*dt2;
98    
99     energies.estr += e;
100     atom[it].energy += e;
101     atom[kt].energy += e;
102     if (minim_values.iprint)
103     fprintf(pcmoutfile,"Fixed Dist: %2s(%-3d) - %2s(%-3d) %-8.3f %-8.3f %-8.3f = %-8.4f\n",atom[it].name,it,atom[kt].name
104     ,kt, rik, fixdis.fxstrd[i], fixdis.fxstrc[i],e);
105     }
106     }
107     }
108    
109     void ebond1()
110     {
111     /* compute stretching energy and first derivatives */
112     int i, it, kt, itype, ktype, classi,classk;
113     double xr, yr, zr, rik, rik2, bcorr;
114     double dt, dt2, e, deddt;
115     double de,dedx,dedy,dedz;
116    
117     energies.estr = 0.0F;
118     for (i=0; i <= natom; i++)
119     {
120     deriv.deb[i][0] = 0.0;
121     deriv.deb[i][1] = 0.0;
122     deriv.deb[i][2] = 0.0;
123     }
124    
125     for (i=0; i < bonds_ff.nbnd; i++)
126     {
127     it = bonds_ff.i12[i][0];
128     kt = bonds_ff.i12[i][1];
129     itype = atom[it].type;
130     ktype = atom[kt].type;
131     classi = atom[it].tclass;
132     classk = atom[kt].tclass;
133     bcorr = 0.0;
134     if ( atom[it].use || atom[kt].use )
135     {
136     xr = atom[it].x - atom[kt].x;
137     yr = atom[it].y - atom[kt].y;
138     zr = atom[it].z - atom[kt].z;
139     rik2 = xr*xr + yr*yr + zr*zr;
140     rik = sqrt(rik2);
141    
142     if (field.type == MM3 && ( (itype ==1 && ktype == 6) || (itype == 6 && ktype ==1) ||
143     (classi == 1 && classk == 6) || (classi == 6 && classk == 1) ) )
144     bnd_corr(it,kt, itype, ktype, &bcorr);
145     dt = rik - (bonds_ff.bl[i] - bcorr);
146    
147     dt2 = dt*dt;
148     e = units.bndunit *bonds_ff.bk[i]*dt2*
149     (1.0 + units.cbnd*dt + units.qbnd*dt2);
150     deddt = 2.0 * units.bndunit * bonds_ff.bk[i] * dt
151     * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
152     if (rik == 0.0)
153     de = 0.0;
154     else
155     de = deddt/rik;
156     dedx = de*xr;
157     dedy = de*yr;
158     dedz = de*zr;
159     energies.estr += e;
160     deriv.deb[it][0] += dedx;
161     deriv.deb[it][1] += dedy;
162     deriv.deb[it][2] += dedz;
163     deriv.deb[kt][0] -= dedx;
164     deriv.deb[kt][1] -= dedy;
165     deriv.deb[kt][2] -= dedz;
166    
167     virial.virx += xr*dedx;
168     virial.viry += yr*dedy;
169     virial.virz += zr*dedz;
170     }
171     }
172     // fixed distances
173     for (i=0; i < fixdis.nfxstr; i++)
174     {
175     it = fixdis.ifxstr[i][0];
176     kt = fixdis.ifxstr[i][1];
177     if ( atom[it].use || atom[kt].use )
178     {
179     xr = atom[it].x - atom[kt].x;
180     yr = atom[it].y - atom[kt].y;
181     zr = atom[it].z - atom[kt].z;
182     rik2 = xr*xr + yr*yr + zr*zr;
183     rik = sqrt(rik2);
184     dt = rik - fixdis.fxstrd[i];
185    
186     dt2 = dt*dt;
187     e = units.bndunit *fixdis.fxstrc[i]*dt2;
188     deddt = 2.0 * units.bndunit * fixdis.fxstrc[i] * dt;
189     if (rik == 0.0)
190     de = 0.0;
191     else
192     de = deddt/rik;
193     dedx = de*xr;
194     dedy = de*yr;
195     dedz = de*zr;
196     energies.estr += e;
197     deriv.deb[it][0] += dedx;
198     deriv.deb[it][1] += dedy;
199     deriv.deb[it][2] += dedz;
200     deriv.deb[kt][0] -= dedx;
201     deriv.deb[kt][1] -= dedy;
202     deriv.deb[kt][2] -= dedz;
203    
204     virial.virx += xr*dedx;
205     virial.viry += yr*dedy;
206     virial.virz += zr*dedz;
207     }
208     }
209     }
210    
211    
212     void ebond2(int ia)
213     {
214     int j,k,m,ibond;
215     double dt,dt2;
216     double xr,yr,zr,rik,rik2,deddt,d2eddt2;
217     double de,term,termx,termy,termz,d2e[3][3];
218    
219     for (m=0; m < MAXIAT; m++)
220     {
221     if (atom[ia].iat[m]!= 0 && atom[ia].bo[m] != 9)
222     {
223     ibond = find_bond(ia, atom[ia].iat[m]);
224     if (bonds_ff.i12[ibond][0] == ia)
225     k = bonds_ff.i12[ibond][1];
226     else
227     k = bonds_ff.i12[ibond][0];
228    
229     xr = atom[ia].x - atom[k].x;
230     yr = atom[ia].y - atom[k].y;
231     zr = atom[ia].z - atom[k].z;
232     rik2 = xr*xr + yr*yr + zr*zr;
233     rik = sqrt(rik2);
234     dt = rik - bonds_ff.bl[ibond];
235     dt2 = dt * dt;
236    
237     deddt = 2.0 * units.bndunit * bonds_ff.bk[ibond] * dt
238     * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
239     d2eddt2 = 2.0 * units.bndunit * bonds_ff.bk[ibond]
240     * (1.0+3.0*units.cbnd*dt+6.0*units.qbnd*dt2);
241    
242     if (rik2 == 0.0)
243     {
244     de = 0.0;
245     term = 0.0;
246     }else
247     {
248     de = deddt / rik;
249     term = (d2eddt2-de) / rik2;
250     }
251    
252     termx = term * xr;
253     termy = term * yr;
254     termz = term * zr;
255     d2e[0][0] = termx*xr + de;
256     d2e[1][0] = termx*yr;
257     d2e[2][0] = termx*zr;
258     d2e[0][1] = d2e[1][0];
259     d2e[1][1] = termy*yr + de;
260     d2e[2][1] = termy*zr;
261     d2e[0][2] = d2e[2][0];
262     d2e[1][2] = d2e[2][1];
263     d2e[2][2] = termz*zr + de;
264    
265     for (j=0; j < 3; j++)
266     {
267     hess.hessx[ia][j] += d2e[j][0];
268     hess.hessy[ia][j] += d2e[j][1];
269     hess.hessz[ia][j] += d2e[j][2];
270     hess.hessx[k][j] -= d2e[j][0];
271     hess.hessy[k][j] -= d2e[j][1];
272     hess.hessz[k][j] -= d2e[j][2];
273     }
274     }
275     }
276     // fixed distances
277     for (m=0; m < MAXIAT; m++)
278     {
279     if (atom[ia].iat[m]!= 0 && atom[ia].bo[m] != 9)
280     {
281     ibond = find_fixed_bond(ia, atom[ia].iat[m]);
282     if (ibond > -1)
283     {
284     if (fixdis.ifxstr[ibond][0] == ia)
285     k = fixdis.ifxstr[ibond][1];
286     else
287     k = fixdis.ifxstr[ibond][0];
288    
289     xr = atom[ia].x - atom[k].x;
290     yr = atom[ia].y - atom[k].y;
291     zr = atom[ia].z - atom[k].z;
292     rik2 = xr*xr + yr*yr + zr*zr;
293     rik = sqrt(rik2);
294     dt = rik - fixdis.fxstrd[ibond];
295     dt2 = dt * dt;
296    
297     deddt = 2.0 * units.bndunit * fixdis.fxstrc[ibond] * dt;
298     d2eddt2 = 2.0 * units.bndunit * fixdis.fxstrc[ibond];
299    
300     if (rik2 == 0.0)
301     {
302     de = 0.0;
303     term = 0.0;
304     }else
305     {
306     de = deddt / rik;
307     term = (d2eddt2-de) / rik2;
308     }
309    
310     termx = term * xr;
311     termy = term * yr;
312     termz = term * zr;
313     d2e[0][0] = termx*xr + de;
314     d2e[1][0] = termx*yr;
315     d2e[2][0] = termx*zr;
316     d2e[0][1] = d2e[1][0];
317     d2e[1][1] = termy*yr + de;
318     d2e[2][1] = termy*zr;
319     d2e[0][2] = d2e[2][0];
320     d2e[1][2] = d2e[2][1];
321     d2e[2][2] = termz*zr + de;
322    
323     for (j=0; j < 3; j++)
324     {
325     hess.hessx[ia][j] += d2e[j][0];
326     hess.hessy[ia][j] += d2e[j][1];
327     hess.hessz[ia][j] += d2e[j][2];
328     hess.hessx[k][j] -= d2e[j][0];
329     hess.hessy[k][j] -= d2e[j][1];
330     hess.hessz[k][j] -= d2e[j][2];
331     }
332     }
333     }
334     }
335    
336     }
337    
338     /* ------------------------------------------ */
339     int find_fixed_bond(int ia, int ib)
340     {
341     int i;
342     int iz;
343    
344     iz = -1;
345     if (fixdis.nfxstr <= 0)
346     return (iz);
347     for (i=0; i < fixdis.nfxstr; i++)
348     {
349     if (ia == fixdis.ifxstr[i][0] && ib == fixdis.ifxstr[i][1])
350     return(i);
351     if (ib == fixdis.ifxstr[i][0] && ia == fixdis.ifxstr[i][1])
352     return(i);
353     }
354     return(iz);
355     }
356     /* ------------------------------------------- */
357     void bnd_corr(int ia, int ib, int itype, int ktype, double *bcorr)
358     {
359     // two terms - first O-ia-ib-R
360     // then X-O-ia-ib
361     // it = carbon, kt = oxygen
362    
363     int i, j,k, iatt, jatt, latt,iatype, jatype, iox;
364     int it, kt;
365     double angle;
366    
367     *bcorr = 0.0;
368    
369     if (itype == 1)
370     {
371     it = ia;
372     kt = ib;
373     } else
374     {
375     it = ib;
376     kt = ia;
377     }
378    
379     for (i=0; i < MAXIAT; i++)
380     {
381     if (atom[kt].iat[i] != 0 && atom[kt].iat[i] != it)
382     {
383     iatt = atom[kt].iat[i];
384     iatype = atom[iatt].type;
385     for(j=0; j < MAXIAT; j++)
386     {
387     iox = 0;
388     if (atom[it].iat[j] != 0 && atom[it].iat[j] != kt)
389     {
390     jatt = atom[it].iat[j];
391     jatype = atom[jatt].type;
392    
393     if (itype == 6 && jatype == 6) iox = 1;
394     if (ktype == 6 && iatype == 6) iox = 2;
395     if (iox != 0)
396     {
397     angle = dihdrl(iatt,kt,it,jatt)/radian;
398     *bcorr += 0.5* 0.01*(1.0-cos(2.0*angle));
399     if (iox == 1)
400     {
401     for (k = 0 ; k < MAXIAT; k++)
402     {
403     if (atom[jatt].iat[k] != 0 && atom[jatt].iat[k] != it)
404     {
405     latt = atom[jatt].iat[k];
406     angle = dihdrl(kt,it,jatt,latt)/radian;
407     *bcorr += 0.5*0.01*(1.0 - cos(2.0*angle))*(-1.30)+0.005;
408     break;
409     }
410     }
411     }
412     }
413     }
414     }
415     }
416     }
417     }