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, 9 months ago) by wdelano
File size: 8024 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
Line File contents
1 #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