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