ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcm7.c
Revision: 63
Committed: Wed Dec 3 03:46:32 2008 UTC (13 years ago) by gilbertke
File size: 17549 byte(s)
Log Message:
updated read_sdf type_mmx and first pass at using best practices
Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5    
6     #include "angles.h"
7 gilbertke 63
8     #include "energies.h"
9 tjod 3 #include "rings.h"
10     #include "torsions.h"
11     #include "bonds_ff.h"
12     #include "derivs.h"
13     #include "hess.h"
14     #include "pot.h"
15     #include "gmmx.h"
16     #include "utility.h"
17     #include "solv.h"
18     #include "field.h"
19 gilbertke 63 #include "dipmom.h"
20 tjod 3
21     struct t_cbonds {
22     int nbnd, icbonds[MAXCBND][3], lpd[MAXCBND], ia1[MAXCBND],ia2[MAXCBND];
23     double vrad[MAXCBND],veps[MAXCBND];
24     } cbonds;
25     EXTERN struct t_high_coord {
26     int ncoord, i13[400][3];
27     } high_coord;
28    
29     EXTERN struct t_minim_control {
30     int type, method, field, added_const;
31     char added_path[256],added_name[256];
32     } minim_control;
33    
34     EXTERN struct t_minim_values {
35     int iprint, ndc, nconst;
36     float dielc;
37     } minim_values;
38     struct t_solvent {
39     int type;
40     double doffset, p1,p2,p3,p4,p5;
41     double *shct,*asolv,*rsolv,*vsolv,*gpol,*rborn;
42     } solvent;
43    
44     int Missing_constants;
45     FILE *errfile;
46    
47     FILE * fopen_path ( char * , char * , char * ) ;
48     void message_alert(char *,char *);
49     void get_added_const(void);
50     void assign_gast_charge();
51     void lattice(void);
52     void bounds(void);
53     void pcm7_min(void);
54     int iscoord_bond(int, int);
55     void get_coordbonds(void);
56     int isbond(int, int);
57     int ishbond(int, int);
58     int istorsion(int, int);
59     void get_hbond(void);
60     void mark_hbond(void);
61     void reset_atom_data(void);
62     int initialise(void);
63     void get_bonds(void);
64     void get_torsions(void);
65     void get_rings(void);
66     void set_field(void);
67     double energy(void);
68     double print_energy(void);
69     void pcmfout(int);
70     void orbital(void);
71     void free_orbit(void);
72     void set_active(void);
73     void estrtor(void);
74     void echarge(void);
75     void ewald(void);
76     void eurey(void);
77     void ebuck(void);
78     void ebuckmm3(void);
79     void estrbnd(void);
80     void edipole(void);
81     void etorsion(void);
82     void eangang(void);
83     void eangle(void);
84     void ebond(void);
85     void eopbend(void);
86     void ehbond(void);
87     void eimprop(void);
88     void eimptors(void);
89     void elj(void);
90     void ecoord_bond(void);
91     void elj_qq(void);
92     void ebuck_charge(void);
93     void ehigh_coord(void);
94    
95     void estrtor1(void);
96     void echarge1(void);
97     void ebuck1(void);
98     void ebuckmm31(void);
99     void estrbnd1(void);
100     void edipole1(void);
101     void etorsion1(void);
102     void eangang1(void);
103     void eangle1(void);
104     void ebond1(void);
105     void eopbend1(void);
106     void ehbond1(void);
107     void eimprop1(void);
108     void eimptors1(void);
109     void elj1(void);
110     void ecoord_bond1(void);
111     void ewald1(void);
112     void eurey1(void);
113     void elj_qq1(void);
114     void ebuck_charge1(void);
115     void ehigh_coord1(void);
116    
117     void kangle(void);
118     void korbit(void);
119     void ktorsion(void);
120     void kdipole(void);
121     void kstrbnd(void);
122     void kcharge(void);
123     void ksolv(void);
124     void kangang(void);
125     void kopbend(void);
126     void kstrtor(void);
127     void khbond(void);
128     void piseq(int, int);
129     void kimprop(void);
130     void kimptors(void);
131     void kcoord_bonds(void);
132     void kvdw(void);
133     int kbond(void);
134     void kewald(void);
135     void kurey(void);
136    
137     void get_memory(void);
138     void free_memory(void);
139     void gradient(void);
140     void minimize(void);
141     void dynamics(void);
142     void read_datafiles(char *);
143     void attach(void);
144     void eheat(void);
145     void pirite(void);
146     void piden(void);
147     void hdel(int);
148     void set_atomtypes(int);
149     void type(void);
150     void generate_bonds(void);
151     void zero_data(void);
152     void ehal(void);
153     void ehal1(void);
154     void ebufcharge(void);
155     void ebufcharge1(void);
156     void eopbend_wilson(void);
157     void eopbend_wilson1(void);
158     void hbondreset(void);
159     void vibrate(void);
160     void pireset(void);
161     void adjust_mmfftypes(void);
162     void build_neighbor_list(int);
163     int setup_calculation(void);
164     void end_calculation(void);
165     void egeom1(void);
166     void esolv1(void);
167     void egeom(void);
168     void esolv(void);
169     void xlogp(float *);
170    
171     void mmxsub(int ia)
172     {
173     int nret;
174    
175     minim_control.type = ia; // 0 to minimize, 1 for single point
176     nret = setup_calculation();
177     if (nret == FALSE)
178     {
179     end_calculation();
180     return;
181     }
182    
183     pcm7_min();
184     end_calculation();
185     }
186     // ================================================================
187     void pcm7_min()
188     {
189     int i, print;
190     double etot;
191    
192     print = FALSE;
193     if (minim_values.iprint == TRUE)
194     {
195     print = TRUE;
196     minim_values.iprint = FALSE;
197     }
198    
199     if (minim_control.type == 1) // single point calculation
200     return;
201    
202     minimize();
203    
204     if (print == TRUE)
205     minim_values.iprint = TRUE;
206    
207     if (minim_values.iprint == TRUE)
208     etot = print_energy();
209     else
210     etot = energy();
211    
212     // compute dipole moment here
213     dipolemom.xdipole = 0.0;
214     dipolemom.ydipole = 0.0;
215     dipolemom.zdipole = 0.0;
216 wdelano 58 if (pot.use_charge || pot.use_bufcharge) charge_dipole();
217     if (pot.use_dipole) dipole_dipole();
218 tjod 3 dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
219     dipolemom.ydipole*dipolemom.ydipole +
220     dipolemom.zdipole*dipolemom.zdipole);
221     }
222     // ==================================================================
223     int setup_calculation()
224     {
225     int nRet;
226     char string[30];
227     double etot;
228     if (minim_control.field == MM3)
229     {
230     // copy mm3 types into type file and retype
231     pot.use_hbond = FALSE;
232     hbondreset();
233     hdel(1);
234     strcpy(string,"mm3.prm");
235     zero_data();
236     read_datafiles(string);
237     set_field();
238     default_intype = MM3;
239     default_outtype = MM3;
240     set_atomtypes(MM3);
241     generate_bonds();
242     if (minim_control.added_const)
243     get_added_const();
244     } else if (minim_control.field == MMFF94)
245     {
246     // copy mm3 types into type file and retype
247     pot.use_hbond = FALSE;
248     pot.use_picalc = FALSE;
249     hbondreset();
250     pireset();
251     hdel(1);
252     set_field();
253     default_intype = MMFF94;
254     default_outtype = MMFF94;
255     set_atomtypes(MMFF94);
256     generate_bonds();
257     if (minim_control.added_const)
258     get_added_const();
259     } else if (minim_control.field == AMBER)
260     {
261     pot.use_hbond = FALSE;
262     pot.use_picalc = FALSE;
263     hbondreset();
264     pireset();
265     hdel(2);
266     set_field();
267     generate_bonds();
268     if (minim_control.added_const)
269     get_added_const();
270     } else if (minim_control.field == OPLSAA)
271     {
272     pot.use_hbond = FALSE;
273     pot.use_picalc = FALSE;
274     hbondreset();
275     pireset();
276     hdel(1);
277     set_field();
278     generate_bonds();
279     if (minim_control.added_const)
280     get_added_const();
281     } else
282     {
283     if (minim_control.added_const)
284     {
285     strcpy(string,"mmxconst.prm");
286     zero_data();
287     read_datafiles(string);
288     get_added_const();
289     }
290     type();
291     set_field();
292     generate_bonds();
293     }
294    
295    
296     // allocate memeory for derivatives
297     get_memory();
298    
299     // setup bonds, angles, torsions, improper torsions and angles, allenes
300     get_bonds();
301     get_angles();
302     get_torsions();
303     get_coordbonds();
304    
305     if (high_coord.ncoord > 0)
306     pot.use_highcoord = TRUE;
307    
308     attach();
309    
310    
311     get_rings();
312     // need allene
313    
314     // set active atoms
315     set_active();
316    
317     /* assign local geometry potential function parameters */
318     Missing_constants = FALSE;
319     // errfile = fopen_path(pcwindir,"pcmod.err","w");
320     if (pot.use_bond || pot.use_strbnd) nRet = kbond();
321     if (nRet == FALSE) // use end_calc to free memory
322     {
323     // message_alert("Parameters missing. See pcmod.err for information","Error");
324     // fclose(errfile);
325     energies.total = 9999.;
326     return FALSE;
327     }
328     if (pot.use_angle || pot.use_strbnd || pot.use_angang) kangle();
329     if (pot.use_angle || pot.use_opbend || pot.use_opbend_wilson) kopbend();
330     if (pot.use_tors || pot.use_strtor || pot.use_tortor) ktorsion();
331     if (pot.use_strbnd) kstrbnd();
332    
333     if (pot.use_buck || pot.use_hal || pot.use_lj) kvdw();
334 wdelano 58 if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();
335 tjod 3
336     if (Missing_constants == TRUE)
337     {
338     // message_alert("Parameters missing. See pcmod.err for information","Error");
339     // fclose(errfile);
340     energies.total = 9999.0;
341     return (FALSE);
342     } else
343     {
344     // fclose(errfile);
345     // remove("pcmod.err");
346     }
347     if (minim_values.iprint == TRUE)
348     etot = print_energy();
349     else
350     etot = energy();
351     return TRUE;
352     }
353     // =================================================================
354     void end_calculation()
355     {
356    
357     free_memory();
358     }
359     // ==================================================================
360     void gradient()
361     {
362     int i, j;
363    
364     energies.estr = 0.0;
365     energies.ebend = 0.0;
366     energies.estrbnd = 0.0;
367     energies.e14 = 0.0;
368     energies.evdw = 0.0;
369     energies.etor = 0.0;
370     energies.eu = 0.0;
371     energies.eopb = 0.0;
372     energies.eangang = 0.0;
373     energies.estrtor = 0.0;
374     energies.ehbond = 0.0;
375     energies.efix = 0.0;
376     energies.eimprop = 0.0;
377     energies.eimptors = 0.0;
378     energies.eurey = 0.0;
379     energies.esolv = 0.0;
380 wdelano 58 energies.egeom =0.0;
381 tjod 3
382     virial.virx = 0.0;
383     virial.viry = 0.0;
384     virial.virz = 0.0;
385    
386     for (i=1; i <= natom; i++)
387     atom[i].energy = 0.0;
388    
389     for (i=1; i <= natom; i++)
390     {
391     for (j=0; j < 3; j++)
392     {
393     deriv.deb[i][j] = 0.0;
394     deriv.dea[i][j] = 0.0;
395     deriv.destb[i][j] = 0.0;
396     deriv.deopb[i][j] = 0.0;
397     deriv.detor[i][j] = 0.0;
398     deriv.de14[i][j] = 0.0;
399     deriv.devdw[i][j]= 0.0;
400     deriv.deqq[i][j] = 0.0;
401     deriv.deaa[i][j] = 0.0;
402     deriv.destor[i][j] = 0.0;
403     deriv.deimprop[i][j] = 0.0;
404     deriv.dehb[i][j] = 0.0;
405     deriv.desolv[i][j] = 0.0;
406 wdelano 58 deriv.degeom[i][j] = 0.0;
407 tjod 3 }
408     }
409    
410     if (pot.use_bond)ebond1();
411     if (pot.use_angle)eangle1();
412     if (pot.use_opbend_wilson) eopbend_wilson1();
413     if (pot.use_tors)etorsion1();
414     if (pot.use_strbnd)estrbnd1();
415    
416     // mmff
417     if (pot.use_hal) ehal1();
418     if (pot.use_bufcharge) ebufcharge1();
419 wdelano 58
420     if (pot.use_geom) egeom1();
421 tjod 3
422     energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd + energies.e14+
423     energies.evdw + energies.eu + energies.ehbond + energies.eangang + energies.estrtor +
424 wdelano 58 energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
425 tjod 3 for (i=1; i <= natom; i++)
426     {
427     for (j=0; j < 3; j++)
428     {
429     deriv.d1[i][j] = deriv.deb[i][j] + deriv.dea[i][j] + deriv.deaa[i][j] +
430     deriv.destb[i][j] + deriv.detor[i][j] + deriv.deopb[i][j] + deriv.dehb[i][j] +
431     deriv.destor[i][j] + deriv.deqq[i][j] + deriv.devdw[i][j] + deriv.de14[i][j] +
432 wdelano 58 deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
433 tjod 3
434     }
435     }
436     }
437     // ==================================================================
438     double energy()
439     {
440     double etot;
441     int i;
442     // struct timeb start,end;
443    
444     energies.total = 0.0;
445     energies.estr = 0.0;
446     energies.ebend = 0.0;
447     energies.etor = 0.0;
448     energies.estrbnd = 0.0;
449     energies.evdw = 0.0;
450     energies.e14 = 0.0;
451     energies.ehbond = 0.0;
452     energies.eu = 0.0;
453     energies.eangang = 0.0;
454     energies.estrtor = 0.0;
455     energies.eimprop = 0.0;
456     energies.eimptors = 0.0;
457     energies.eopb = 0.0;
458     energies.eurey = 0.0;
459     energies.esolv = 0.0;
460 wdelano 58 energies.egeom = 0.0;
461 tjod 3
462     for (i=1; i <= natom; i++)
463     atom[i].energy = 0.0;
464    
465     if (pot.use_bond)ebond();
466     if (pot.use_angle)eangle();
467     if (pot.use_opbend_wilson) eopbend_wilson();
468     if (pot.use_tors)etorsion();
469     if (pot.use_strbnd)estrbnd();
470     // mmff
471     if (pot.use_hal) ehal();
472     if (pot.use_bufcharge) ebufcharge();
473 wdelano 58
474     if (pot.use_geom) egeom();
475 tjod 3
476     energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
477     + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
478 wdelano 58 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
479 tjod 3 etot = energies.total;
480     return (etot);
481     }
482     /* =========================================================================== */
483     double print_energy()
484     {
485     double etot;
486     int i;
487    
488     energies.total = 0.0;
489     energies.estr = 0.0;
490     energies.ebend = 0.0;
491     energies.etor = 0.0;
492     energies.estrbnd = 0.0;
493     energies.evdw = 0.0;
494     energies.e14 = 0.0;
495     energies.ehbond = 0.0;
496     energies.eu = 0.0;
497     energies.eangang = 0.0;
498     energies.estrtor = 0.0;
499     energies.eimprop = 0.0;
500     energies.eimptors = 0.0;
501     energies.eopb = 0.0;
502     energies.eurey = 0.0;
503     energies.esolv = 0.0;
504 wdelano 58 energies.egeom = 0.0;
505 tjod 3
506     for (i=1; i <= natom; i++)
507     atom[i].energy = 0.0;
508    
509     if (pot.use_bond)ebond();
510     if (pot.use_angle)eangle();
511     if (pot.use_opbend_wilson) eopbend_wilson();
512     if (pot.use_tors)etorsion();
513     if (pot.use_strbnd)estrbnd();
514     // mmff
515     if (pot.use_hal) ehal();
516     if (pot.use_bufcharge) ebufcharge();
517    
518 wdelano 58 if (pot.use_geom) egeom();
519    
520 tjod 3 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
521     + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
522 wdelano 58 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
523 tjod 3 etot = energies.total;
524    
525     return (etot);
526     }
527    
528     /* ================================================================== */
529    
530     int iscoord_bond(int i, int j)
531     {
532     int k;
533     long int mask;
534    
535     mask = 1 << METCOORD_MASK ;
536    
537     for (k=0; k < MAXIAT; k++)
538     {
539     if (atom[i].iat[k] == j && atom[i].bo[k] == 9)
540     {
541     if (atom[i].type >= 300)
542     {
543     if (atom[j].type == 2 || atom[j].type == 4 || atom[j].type == 29 ||
544     atom[j].type == 30 || atom[j].type == 40 || atom[j].type == 48 )
545     return TRUE;
546     }
547     if (atom[j].type >= 300)
548     {
549     if (atom[i].type == 2 || atom[i].type == 4 || atom[i].type == 29 ||
550     atom[i].type == 30 || atom[i].type == 40 || atom[i].type == 48 )
551     return TRUE;
552     }
553     }
554     }
555     // throughout Metal-lone pairs if coordinated bond
556     // metal to donor atom
557     if ( (atom[i].type >= 300) && atom[j].flags & mask)
558     return TRUE;
559     if ( (atom[j].type >= 300) && atom[i].flags & mask)
560     return TRUE;
561    
562     return FALSE;
563     }
564    
565     void get_coordbonds()
566     {
567     int i, j, k, iatype, jatm;
568    
569     cbonds.nbnd = 0;
570     pot.use_coordb = FALSE;
571    
572     for (i=1; i <= natom; i++)
573     {
574     if (atom[i].mmx_type >= 300)
575     {
576     for (j=0; j <= MAXIAT; j++)
577     {
578     if (atom[i].bo[j] == 9) // coord bond
579     {
580     iatype = atom[atom[i].iat[j]].mmx_type;
581     if (iatype == 2 || iatype == 4 || iatype == 29 ||
582     iatype == 30 || iatype == 40 || iatype == 48 )
583     {
584     cbonds.icbonds[cbonds.nbnd][0] = i;
585     cbonds.icbonds[cbonds.nbnd][1] = atom[i].iat[j];
586     cbonds.icbonds[cbonds.nbnd][2] = 0;
587     cbonds.nbnd++;
588     } else
589     {
590     jatm = atom[i].iat[j];
591     for (k = 0; k < MAXIAT; k++)
592     {
593     if (atom[atom[jatm].iat[k]].mmx_type == 20) // lp
594     {
595     cbonds.icbonds[cbonds.nbnd][0] = i;
596     cbonds.icbonds[cbonds.nbnd][1] = jatm;
597     cbonds.icbonds[cbonds.nbnd][2] = atom[jatm].iat[k];
598     cbonds.nbnd++;
599     }
600     }
601     }
602     if (cbonds.nbnd > MAXCBND)
603     {
604 wdelano 58 fprintf(pcmlogfile,"Error - Too many coordinated bonds!\nProgram will now quit.\n");
605 tjod 3 exit(0);
606     }
607     }
608     }
609     }
610     }
611     if (cbonds.nbnd > 0)
612     pot.use_coordb = TRUE;
613     }
614     int isbond(int i, int j)
615     {
616     int k;
617     for (k=0; k < MAXIAT; k++)
618     {
619     if (atom[i].iat[k] == j && atom[i].bo[k] != 9)
620     return TRUE;
621     }
622     return FALSE;
623     }
624    
625     void get_bonds()
626     {
627     int i, j;
628     long int mask;
629    
630     bonds_ff.nbnd = 0;
631     mask = 1L << 0; // flag for pi atoms
632    
633     for (i=1; i <= natom; i++)
634     {
635     for (j=i+1; j <= natom; j++)
636     {
637     if (isbond(i,j))
638     {
639     bonds_ff.i12[bonds_ff.nbnd][0] = i;
640     bonds_ff.i12[bonds_ff.nbnd][1] = j;
641     bonds_ff.nbnd++;
642     if (bonds_ff.nbnd > MAXBND)
643     {
644 wdelano 58 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
645 tjod 3 exit(0);
646     }
647     }
648     }
649     }
650     }