ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/pcm7.c
Revision: 99
Committed: Mon Jan 19 04:31:29 2009 UTC (10 years, 8 months ago) by wdelano
File size: 10584 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 "attached.h"
7 #include "energies.h"
8 #include "bonds_ff.h"
9 #include "derivs.h"
10 #include "pot.h"
11 #include "utility.h"
12 #include "solv.h"
13 #include "dipmom.h"
14
15 EXTERN struct t_minim_control {
16 int type, method, field, added_const;
17 char added_path[256],added_name[256];
18 } minim_control;
19
20 EXTERN struct t_minim_values {
21 int iprint, ndc, nconst;
22 float dielc;
23 } minim_values;
24 struct t_solvent {
25 int type;
26 double doffset, p1,p2,p3,p4,p5;
27 double *shct,*asolv,*rsolv,*vsolv,*gpol,*rborn;
28 } solvent;
29
30 int Missing_constants;
31
32 void message_alert(char *,char *);
33 void pcm7_min(void);
34 int isbond(int, int);
35 int initialise(void);
36 void get_bonds(void);
37 void get_torsions(void);
38 void get_rings(void);
39 void set_field(void);
40 double energy(void);
41 double print_energy(void);
42 void pcmfout(int);
43 void set_active(void);
44 int setup_calculation(void);
45 void end_calculation(void);
46
47 void egeom1(void);
48 void esolv1(void);
49 void egeom(void);
50 void esolv(void);
51
52 void ehal(void);
53 void ehal1(void);
54 void ebufcharge(void);
55 void ebufcharge1(void);
56 void eopbend_wilson(void);
57 void eopbend_wilson1(void);
58
59 void estrbnd(void);
60 void etorsion(void);
61 void eangang(void);
62 void eangle(void);
63 void ebond(void);
64
65 void estrbnd1(void);
66 void etorsion1(void);
67 void eangle1(void);
68 void ebond1(void);
69
70 void kangle(void);
71 void ktorsion(void);
72 void kstrbnd(void);
73 void kcharge(void);
74 void ksolv(void);
75 void kopbend(void);
76 void kvdw(void);
77 int kbond(void);
78
79 void get_memory(void);
80 void free_memory(void);
81 void gradient(void);
82 void minimize(void);
83 void read_datafiles(char *);
84 void get_added_const(void);
85 void attach(void);
86 void hdel(int);
87 void set_atomtypes(int);
88 void type(void);
89 void zero_data(void);
90 void xlogp(float *);
91
92 // ==================================================================
93 int setup_calculation()
94 {
95 int i,j,nRet;
96 char string[30];
97 double etot;
98 if (minim_control.field == MMFF94)
99 {
100 // copy mm3 types into type file and retype
101 pot.use_hbond = FALSE;
102 pot.use_picalc = FALSE;
103 hdel(1);
104 set_field();
105 set_atomtypes(MMFF94);
106 if (minim_control.added_const)
107 get_added_const();
108 } else
109 {
110 if (minim_control.added_const)
111 {
112 strcpy(string,"mmxconst.prm");
113 zero_data();
114 read_datafiles(string);
115 get_added_const();
116 }
117 type();
118 set_field();
119 }
120
121
122 // allocate memeory for derivatives
123 get_memory();
124
125 // setup bonds, angles, torsions, improper torsions and angles, allenes
126 get_bonds();
127 get_angles();
128 get_torsions();
129
130 attach();
131 get_rings();
132 // need allene
133
134 // set active atoms
135 set_active();
136 // setup non_bonded list of atoms to skip
137 for (i=1; i <= natom; i++)
138 {
139 for (j=0; j < MAXIAT; j++)
140 if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
141 {
142 skip[i][atom[i].iat[j]] = i;
143 skip[atom[i].iat[j]][i] = atom[i].iat[j];
144 }
145 for (j=0; j < attached.n13[i]; j++)
146 {
147 skip[i][attached.i13[j][i]] = i;
148 skip[attached.i13[j][i]][i] = attached.i13[j][i];
149 }
150 for(j=0; j < attached.n14[i]; j++)
151 {
152 skip[i][attached.i14[j][i]] = -i;
153 skip[attached.i14[j][i]][i] = -attached.i14[j][i];
154 }
155 }
156 /* assign local geometry potential function parameters */
157 Missing_constants = FALSE;
158 if (pot.use_bond || pot.use_strbnd) nRet = kbond();
159 if (nRet == FALSE) // use end_calc to free memory
160 {
161 energies.total = 9999.;
162 return FALSE;
163 }
164 if (pot.use_angle || pot.use_strbnd || pot.use_angang) kangle();
165 if (pot.use_angle || pot.use_opbend || pot.use_opbend_wilson) kopbend();
166 if (pot.use_tors || pot.use_strtor || pot.use_tortor) ktorsion();
167 if (pot.use_strbnd) kstrbnd();
168
169 if (pot.use_buck || pot.use_hal || pot.use_lj) kvdw();
170 if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();
171
172 if (Missing_constants == TRUE)
173 {
174 energies.total = 9999.0;
175 return (FALSE);
176 }
177 if (minim_values.iprint == TRUE)
178 etot = print_energy();
179 else
180 etot = energy();
181 return TRUE;
182 }
183 // =================================================================
184 void end_calculation()
185 {
186
187 free_memory();
188 }
189 // ==================================================================
190 void gradient()
191 {
192 int i, j;
193
194 energies.estr = 0.0;
195 energies.ebend = 0.0;
196 energies.estrbnd = 0.0;
197 energies.e14 = 0.0;
198 energies.evdw = 0.0;
199 energies.etor = 0.0;
200 energies.eu = 0.0;
201 energies.eopb = 0.0;
202 energies.eangang = 0.0;
203 energies.estrtor = 0.0;
204 energies.ehbond = 0.0;
205 energies.efix = 0.0;
206 energies.eimprop = 0.0;
207 energies.eimptors = 0.0;
208 energies.eurey = 0.0;
209 energies.esolv = 0.0;
210 energies.egeom =0.0;
211
212 for (i=1; i <= natom; i++)
213 atom[i].energy = 0.0;
214
215 for (i=1; i <= natom; i++)
216 {
217 for (j=0; j < 3; j++)
218 {
219 deriv.deb[i][j] = 0.0;
220 deriv.dea[i][j] = 0.0;
221 deriv.destb[i][j] = 0.0;
222 deriv.deopb[i][j] = 0.0;
223 deriv.detor[i][j] = 0.0;
224 deriv.de14[i][j] = 0.0;
225 deriv.devdw[i][j]= 0.0;
226 deriv.deqq[i][j] = 0.0;
227 deriv.deaa[i][j] = 0.0;
228 deriv.destor[i][j] = 0.0;
229 deriv.deimprop[i][j] = 0.0;
230 deriv.dehb[i][j] = 0.0;
231 deriv.desolv[i][j] = 0.0;
232 deriv.degeom[i][j] = 0.0;
233 }
234 }
235
236 if (pot.use_bond)ebond1();
237 if (pot.use_angle)eangle1();
238 if (pot.use_opbend_wilson) eopbend_wilson1();
239 if (pot.use_tors)etorsion1();
240 if (pot.use_strbnd)estrbnd1();
241
242 // mmff
243 if (pot.use_hal) ehal1();
244 if (pot.use_bufcharge) ebufcharge1();
245
246 if (pot.use_geom) egeom1();
247
248 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd + energies.e14+
249 energies.evdw + energies.eu + energies.ehbond + energies.eangang + energies.estrtor +
250 energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
251 for (i=1; i <= natom; i++)
252 {
253 for (j=0; j < 3; j++)
254 {
255 deriv.d1[i][j] = deriv.deb[i][j] + deriv.dea[i][j] + deriv.deaa[i][j] +
256 deriv.destb[i][j] + deriv.detor[i][j] + deriv.deopb[i][j] + deriv.dehb[i][j] +
257 deriv.destor[i][j] + deriv.deqq[i][j] + deriv.devdw[i][j] + deriv.de14[i][j] +
258 deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
259
260 }
261 }
262 }
263 // ==================================================================
264 double energy()
265 {
266 double etot;
267 int i;
268
269 energies.total = 0.0;
270 energies.estr = 0.0;
271 energies.ebend = 0.0;
272 energies.etor = 0.0;
273 energies.estrbnd = 0.0;
274 energies.evdw = 0.0;
275 energies.e14 = 0.0;
276 energies.ehbond = 0.0;
277 energies.eu = 0.0;
278 energies.eangang = 0.0;
279 energies.estrtor = 0.0;
280 energies.eimprop = 0.0;
281 energies.eimptors = 0.0;
282 energies.eopb = 0.0;
283 energies.eurey = 0.0;
284 energies.esolv = 0.0;
285 energies.egeom = 0.0;
286
287 for (i=1; i <= natom; i++)
288 atom[i].energy = 0.0;
289
290 if (pot.use_bond)ebond();
291 if (pot.use_angle)eangle();
292 if (pot.use_opbend_wilson) eopbend_wilson();
293 if (pot.use_tors)etorsion();
294 if (pot.use_strbnd)estrbnd();
295 // mmff
296 if (pot.use_hal) ehal();
297 if (pot.use_bufcharge) ebufcharge();
298
299 if (pot.use_geom) egeom();
300
301 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
302 + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
303 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
304 etot = energies.total;
305 return (etot);
306 }
307 /* =========================================================================== */
308 double print_energy()
309 {
310 double etot;
311 int i;
312
313 energies.total = 0.0;
314 energies.estr = 0.0;
315 energies.ebend = 0.0;
316 energies.etor = 0.0;
317 energies.estrbnd = 0.0;
318 energies.evdw = 0.0;
319 energies.e14 = 0.0;
320 energies.ehbond = 0.0;
321 energies.eu = 0.0;
322 energies.eangang = 0.0;
323 energies.estrtor = 0.0;
324 energies.eimprop = 0.0;
325 energies.eimptors = 0.0;
326 energies.eopb = 0.0;
327 energies.eurey = 0.0;
328 energies.esolv = 0.0;
329 energies.egeom = 0.0;
330
331 for (i=1; i <= natom; i++)
332 atom[i].energy = 0.0;
333
334 if (pot.use_bond)ebond();
335 if (pot.use_angle)eangle();
336 if (pot.use_opbend_wilson) eopbend_wilson();
337 if (pot.use_tors)etorsion();
338 if (pot.use_strbnd)estrbnd();
339 // mmff
340 if (pot.use_hal) ehal();
341 if (pot.use_bufcharge) ebufcharge();
342
343 if (pot.use_geom) egeom();
344
345 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
346 + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
347 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
348 etot = energies.total;
349
350 return (etot);
351 }
352
353 /* ================================================================== */
354 int isbond(int i, int j)
355 {
356 int k;
357 for (k=0; k < MAXIAT; k++)
358 {
359 if (atom[i].iat[k] == j && atom[i].bo[k] != 9)
360 return TRUE;
361 }
362 return FALSE;
363 }
364 // ==========================
365 void get_bonds()
366 {
367 int i, j;
368
369 bonds_ff.nbnd = 0;
370 for (i=1; i <= natom; i++)
371 {
372 for (j=i+1; j <= natom; j++)
373 {
374 if (isbond(i,j))
375 {
376 bonds_ff.i12[bonds_ff.nbnd][0] = i;
377 bonds_ff.i12[bonds_ff.nbnd][1] = j;
378 bonds_ff.nbnd++;
379 if (bonds_ff.nbnd > MAXBND)
380 {
381 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
382 exit(0);
383 }
384 }
385 }
386 }
387 }