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, 9 months ago) by gilbertke
File size: 10547 byte(s)
Log Message:
combine rings and srings to eliminate duplication
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5
6 #include "attached.h"
7 #include "energies.h"
8 #include "bonds_ff.h"
9 #include "derivs.h"
10 #include "pot.h"
11 #include "utility.h"
12 #include "solv.h"
13 #include "dipmom.h"
14
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 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 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 void get_added_const(void);
84 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 int i,j,nRet;
95 char string[30];
96 double etot;
97 if (minim_control.field == MMFF94)
98 {
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 // 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 /* 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 if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();
169
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 energies.egeom =0.0;
209
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 deriv.degeom[i][j] = 0.0;
231 }
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
244 if (pot.use_geom) egeom1();
245
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 energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
249 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 deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
257
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 energies.egeom = 0.0;
284
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
297 if (pot.use_geom) egeom();
298
299 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
300 + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
301 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
302 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 energies.egeom = 0.0;
328
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 if (pot.use_geom) egeom();
342
343 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
344 + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
345 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
346 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 // ==========================
363 void get_bonds()
364 {
365 int i, j;
366
367 bonds_ff.nbnd = 0;
368 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 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
380 exit(0);
381 }
382 }
383 }
384 }
385 }