ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcm7.c
Revision: 90
Committed: Fri Jan 16 21:55:34 2009 UTC (13 years ago) by gilbertke
File size: 10584 byte(s)
Log Message:
code cleanup
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    
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 gilbertke 85 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 tjod 3 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 gilbertke 85 void get_added_const(void);
85 tjod 3 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 gilbertke 89 int i,j,nRet;
96 tjod 3 char string[30];
97     double etot;
98 gilbertke 90 if (minim_control.field == MMFF94)
99 tjod 3 {
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 gilbertke 89 // 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 tjod 3 /* 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 wdelano 58 if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();
171 tjod 3
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 wdelano 58 energies.egeom =0.0;
211 tjod 3
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 wdelano 58 deriv.degeom[i][j] = 0.0;
233 tjod 3 }
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 wdelano 58
246     if (pot.use_geom) egeom1();
247 tjod 3
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 wdelano 58 energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
251 tjod 3 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 wdelano 58 deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
259 tjod 3
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 wdelano 58 energies.egeom = 0.0;
286 tjod 3
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 wdelano 58
299     if (pot.use_geom) egeom();
300 tjod 3
301     energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
302     + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
303 wdelano 58 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
304 tjod 3 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 wdelano 58 energies.egeom = 0.0;
330 tjod 3
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 wdelano 58 if (pot.use_geom) egeom();
344    
345 tjod 3 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
346     + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
347 wdelano 58 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
348 tjod 3 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 gilbertke 85 // ==========================
365 tjod 3 void get_bonds()
366     {
367     int i, j;
368    
369 gilbertke 85 bonds_ff.nbnd = 0;
370 tjod 3 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 wdelano 58 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
382 tjod 3 exit(0);
383     }
384     }
385     }
386     }
387     }