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