ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/ebond.c
Revision: 63
Committed: Wed Dec 3 03:46:32 2008 UTC (11 years, 1 month ago) by gilbertke
File size: 8126 byte(s)
Log Message:
updated read_sdf type_mmx and first pass at using best practices
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     virial.virx += xr*dedx;
132     virial.viry += yr*dedy;
133     virial.virz += zr*dedz;
134     }
135     }
136     }
137     // =============================================
138     void ebond2(int ia)
139     {
140     int j,k,m,ibond;
141     double dt,dt2;
142     double xr,yr,zr,rik,rik2,deddt,d2eddt2;
143     double de,term,termx,termy,termz,d2e[3][3];
144    
145     for (m=0; m < MAXIAT; m++)
146     {
147     if (atom[ia].iat[m]!= 0 && atom[ia].bo[m] != 9)
148     {
149     ibond = find_bond(ia, atom[ia].iat[m]);
150     if (bonds_ff.i12[ibond][0] == ia)
151     k = bonds_ff.i12[ibond][1];
152     else
153     k = bonds_ff.i12[ibond][0];
154    
155     xr = atom[ia].x - atom[k].x;
156     yr = atom[ia].y - atom[k].y;
157     zr = atom[ia].z - atom[k].z;
158     rik2 = xr*xr + yr*yr + zr*zr;
159     rik = sqrt(rik2);
160     dt = rik - bonds_ff.bl[ibond];
161     dt2 = dt * dt;
162    
163     deddt = 2.0 * units.bndunit * bonds_ff.bk[ibond] * dt
164     * (1.0+1.5*units.cbnd*dt+2.0*units.qbnd*dt2);
165     d2eddt2 = 2.0 * units.bndunit * bonds_ff.bk[ibond]
166     * (1.0+3.0*units.cbnd*dt+6.0*units.qbnd*dt2);
167    
168     if (rik2 == 0.0)
169     {
170     de = 0.0;
171     term = 0.0;
172     }else
173     {
174     de = deddt / rik;
175     term = (d2eddt2-de) / rik2;
176     }
177    
178     termx = term * xr;
179     termy = term * yr;
180     termz = term * zr;
181     d2e[0][0] = termx*xr + de;
182     d2e[1][0] = termx*yr;
183     d2e[2][0] = termx*zr;
184     d2e[0][1] = d2e[1][0];
185     d2e[1][1] = termy*yr + de;
186     d2e[2][1] = termy*zr;
187     d2e[0][2] = d2e[2][0];
188     d2e[1][2] = d2e[2][1];
189     d2e[2][2] = termz*zr + de;
190    
191     for (j=0; j < 3; j++)
192     {
193     hess.hessx[ia][j] += d2e[j][0];
194     hess.hessy[ia][j] += d2e[j][1];
195     hess.hessz[ia][j] += d2e[j][2];
196     hess.hessx[k][j] -= d2e[j][0];
197     hess.hessy[k][j] -= d2e[j][1];
198     hess.hessz[k][j] -= d2e[j][2];
199     }
200     }
201     }
202     }
203     /* ------------------------------------------- */
204     void bnd_corr(int ia, int ib, int itype, int ktype, double *bcorr)
205     {
206     // two terms - first O-ia-ib-R
207     // then X-O-ia-ib
208     // it = carbon, kt = oxygen
209    
210     int i, j,k, iatt, jatt, latt,iatype, jatype, iox;
211     int it, kt;
212     double angle;
213    
214     *bcorr = 0.0;
215    
216     if (itype == 1)
217     {
218     it = ia;
219     kt = ib;
220     } else
221     {
222     it = ib;
223     kt = ia;
224     }
225    
226     for (i=0; i < MAXIAT; i++)
227     {
228     if (atom[kt].iat[i] != 0 && atom[kt].iat[i] != it)
229     {
230     iatt = atom[kt].iat[i];
231     iatype = atom[iatt].type;
232     for(j=0; j < MAXIAT; j++)
233     {
234     iox = 0;
235     if (atom[it].iat[j] != 0 && atom[it].iat[j] != kt)
236     {
237     jatt = atom[it].iat[j];
238     jatype = atom[jatt].type;
239    
240     if (itype == 6 && jatype == 6) iox = 1;
241     if (ktype == 6 && iatype == 6) iox = 2;
242     if (iox != 0)
243     {
244     angle = dihdrl(iatt,kt,it,jatt)/radian;
245     *bcorr += 0.5* 0.01*(1.0-cos(2.0*angle));
246     if (iox == 1)
247     {
248     for (k = 0 ; k < MAXIAT; k++)
249     {
250     if (atom[jatt].iat[k] != 0 && atom[jatt].iat[k] != it)
251     {
252     latt = atom[jatt].iat[k];
253     angle = dihdrl(kt,it,jatt,latt)/radian;
254     *bcorr += 0.5*0.01*(1.0 - cos(2.0*angle))*(-1.30)+0.005;
255     break;
256     }
257     }
258     }
259     }
260     }
261     }
262     }
263     }
264     }
265