ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcm7.c
Revision: 27
Committed: Tue Jul 8 19:10:25 2008 UTC (13 years, 5 months ago) by tjod
File size: 18103 byte(s)
Log Message:
move mengine src to reflect new freemol directory structure

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