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 (10 years, 9 months ago) by gilbertke
File size: 11131 byte(s)
Log Message:
more cleanup and code removal
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5
6 #include "energies.h"
7 #include "bonds_ff.h"
8 #include "derivs.h"
9 #include "pot.h"
10 #include "utility.h"
11 #include "solv.h"
12 #include "dipmom.h"
13
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 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 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 void get_added_const(void);
86 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 if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();
179
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 energies.egeom =0.0;
225
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 deriv.degeom[i][j] = 0.0;
251 }
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
264 if (pot.use_geom) egeom1();
265
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 energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
269 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 deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
277
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 energies.egeom = 0.0;
305
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
318 if (pot.use_geom) egeom();
319
320 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
321 + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
322 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
323 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 energies.egeom = 0.0;
349
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 if (pot.use_geom) egeom();
363
364 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
365 + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
366 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
367 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 // ==========================
384 void get_bonds()
385 {
386 int i, j;
387
388 bonds_ff.nbnd = 0;
389 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 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
401 exit(0);
402 }
403 }
404 }
405 }
406 }