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