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 User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5    
6 gilbertke 89 #include "attached.h"
7 gilbertke 63 #include "energies.h"
8 tjod 3 #include "bonds_ff.h"
9     #include "derivs.h"
10     #include "pot.h"
11     #include "utility.h"
12     #include "solv.h"
13 gilbertke 63 #include "dipmom.h"
14 tjod 3
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 gilbertke 85 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 tjod 3 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 gilbertke 85 void get_added_const(void);
87 tjod 3 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 gilbertke 89 int i,j,nRet;
99 tjod 3 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 gilbertke 89 // 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 tjod 3 /* 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 wdelano 58 if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();
199 tjod 3
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 wdelano 58 energies.egeom =0.0;
245 tjod 3
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 wdelano 58 deriv.degeom[i][j] = 0.0;
267 tjod 3 }
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 wdelano 58
280     if (pot.use_geom) egeom1();
281 tjod 3
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 wdelano 58 energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
285 tjod 3 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 wdelano 58 deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
293 tjod 3
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 wdelano 58 energies.egeom = 0.0;
320 tjod 3
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 wdelano 58
333     if (pot.use_geom) egeom();
334 tjod 3
335     energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
336     + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
337 wdelano 58 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
338 tjod 3 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 wdelano 58 energies.egeom = 0.0;
364 tjod 3
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 wdelano 58 if (pot.use_geom) egeom();
378    
379 tjod 3 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
380     + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
381 wdelano 58 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
382 tjod 3 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 gilbertke 85 // ==========================
399 tjod 3 void get_bonds()
400     {
401     int i, j;
402    
403 gilbertke 85 bonds_ff.nbnd = 0;
404 tjod 3 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 wdelano 58 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
416 tjod 3 exit(0);
417     }
418     }
419     }
420     }
421     }