ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcm7.c
(Generate patch)
# Line 3 | Line 3
3   #include "pcwin.h"
4   #include "pcmod.h"
5  
6 + #include "attached.h"
7   #include "energies.h"
7 #include "angles.h"
8 #include "rings.h"
9 #include "torsions.h"
8   #include "bonds_ff.h"
9   #include "derivs.h"
12 #include "hess.h"
10   #include "pot.h"
14 #include "gmmx.h"
11   #include "utility.h"
12   #include "solv.h"  
13 < #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;
13 > #include "dipmom.h"
14                  
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
15   EXTERN struct t_minim_control {
16          int type, method, field, added_const;
17          char added_path[256],added_name[256];
18          } minim_control;
42
43 struct t_dipolemom {
44        double total, xdipole, ydipole, zdipole;
45       }  dipolemom;
19          
20   EXTERN struct t_minim_values {
21          int iprint, ndc, nconst;
# Line 59 | Line 32
32      
33   FILE * fopen_path ( char * , char * , char * ) ;
34   void message_alert(char *,char *);
62 void get_added_const(void);
63 void assign_gast_charge();
64 void lattice(void);
65 void bounds(void);
35   void pcm7_min(void);
67 int iscoord_bond(int, int);
68 void get_coordbonds(void);
36   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);
37   int initialise(void);
38   void get_bonds(void);
78 void get_angles(void);
39   void get_torsions(void);
40   void get_rings(void);
41   void set_field(void);
42   double energy(void);
43   double print_energy(void);
44   void pcmfout(int);
85 void orbital(void);
86 void free_orbit(void);
45   void set_active(void);
46 < void estrtor(void);
47 < void echarge(void);
48 < void ewald(void);
49 < void eurey(void);
50 < void ebuck(void);
51 < void ebuckmm3(void);
46 > int setup_calculation(void);
47 > void end_calculation(void);
48 >
49 > void egeom1(void);
50 > void esolv1(void);
51 > void egeom(void);
52 > void esolv(void);
53 >
54 > void ehal(void);
55 > void ehal1(void);
56 > void ebufcharge(void);
57 > void ebufcharge1(void);
58 > void eopbend_wilson(void);
59 > void eopbend_wilson1(void);
60 >
61   void estrbnd(void);
95 void edipole(void);
62   void etorsion(void);
63   void eangang(void);
64   void eangle(void);
65   void ebond(void);
66 < 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);
66 >
67   void estrbnd1(void);
115 void edipole1(void);
68   void etorsion1(void);
117 void eangang1(void);
69   void eangle1(void);
70   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);
71  
72   void kangle(void);
133 void korbit(void);
73   void ktorsion(void);
135 void kdipole(void);
74   void kstrbnd(void);
75   void kcharge(void);
76   void ksolv(void);
139 void kangang(void);
77   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);
78   void kvdw(void);
79   int  kbond(void);
149 void kewald(void);
150 void kurey(void);
80  
81   void get_memory(void);
82   void free_memory(void);
83   void gradient(void);
84   void minimize(void);
156 void dynamics(void);
85   void read_datafiles(char *);
86 + void get_added_const(void);
87   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);
88   void hdel(int);
89   void set_atomtypes(int);
90   void type(void);
91   void generate_bonds(void);
92   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);
93   void xlogp(float *);
94  
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 }
95   // ==================================================================
96   int setup_calculation()
97   {
98 <    int nRet;
98 >  int i,j,nRet;
99      char string[30];
100      double etot;
101      if (minim_control.field == MM3)
# Line 277 | Line 129
129          generate_bonds();
130          if (minim_control.added_const)
131             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();
132      } else
133      {
134          if (minim_control.added_const)
# Line 321 | Line 151
151    get_bonds();
152    get_angles();
153    get_torsions();
324  get_coordbonds();
325
326  if (high_coord.ncoord > 0)
327      pot.use_highcoord = TRUE;
154  
155    attach();
330
331  
156    get_rings();
157    // need allene
158  
159    // set active atoms
160    set_active();
161 <
161 >  // setup non_bonded list of atoms to skip
162 >  for (i=1; i <= natom; i++)
163 >    {
164 >         for (j=0; j < MAXIAT; j++)
165 >             if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
166 >               {
167 >                    skip[i][atom[i].iat[j]] = i;
168 >                    skip[atom[i].iat[j]][i] = atom[i].iat[j];
169 >               }
170 >         for (j=0; j < attached.n13[i]; j++)
171 >           {
172 >               skip[i][attached.i13[j][i]] = i;
173 >               skip[attached.i13[j][i]][i] = attached.i13[j][i];
174 >           }
175 >         for(j=0; j < attached.n14[i]; j++)
176 >           {
177 >            skip[i][attached.i14[j][i]] = -i;
178 >            skip[attached.i14[j][i]][i] = -attached.i14[j][i];
179 >           }
180 >    }
181   /*  assign local geometry potential function parameters  */
182      Missing_constants = FALSE;
183      //     errfile = fopen_path(pcwindir,"pcmod.err","w");
# Line 342 | Line 185
185       if (nRet == FALSE) // use end_calc to free memory
186       {
187         //         message_alert("Parameters missing. See pcmod.err for information","Error");
345         optimize_data.param_avail = 1;
188           //         fclose(errfile);
189           energies.total = 9999.;
190           return FALSE;
# Line 353 | Line 195
195       if (pot.use_strbnd) kstrbnd();
196        
197        if (pot.use_buck || pot.use_hal || pot.use_lj) kvdw();
198 <      if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald) && !use_external_chrg && !use_gast_chrg) kcharge();
198 >      if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();
199  
200         if (Missing_constants == TRUE)
201         {
202           //           message_alert("Parameters missing. See pcmod.err for information","Error");
203           //           fclose(errfile);
362           optimize_data.param_avail = 1;
204             energies.total = 9999.0;
205             return (FALSE);
206         } else
# Line 400 | Line 241
241        energies.eimptors = 0.0;
242        energies.eurey = 0.0;
243        energies.esolv = 0.0;
244 <
404 <      virial.virx = 0.0;
405 <      virial.viry = 0.0;
406 <      virial.virz = 0.0;
244 >      energies.egeom =0.0;
245        
246        for (i=1; i <= natom; i++)
247             atom[i].energy = 0.0;
# Line 425 | Line 263
263              deriv.deimprop[i][j] = 0.0;
264              deriv.dehb[i][j] = 0.0;
265              deriv.desolv[i][j] = 0.0;
266 +            deriv.degeom[i][j] = 0.0;
267          }
268        }
269    
# Line 437 | Line 276
276        // mmff
277        if (pot.use_hal) ehal1();
278        if (pot.use_bufcharge) ebufcharge1();
279 +
280 +      if (pot.use_geom) egeom1();
281      
282        energies.total =  energies.estr + energies.ebend + energies.etor + energies.estrbnd + energies.e14+
283                          energies.evdw + energies.eu + energies.ehbond + energies.eangang + energies.estrtor +
284 <                        energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv;
284 >                        energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
285        for (i=1; i <= natom; i++)
286        {
287            for (j=0; j < 3; j++)
# Line 448 | Line 289
289                deriv.d1[i][j] = deriv.deb[i][j] + deriv.dea[i][j] + deriv.deaa[i][j] +
290                                 deriv.destb[i][j] + deriv.detor[i][j] + deriv.deopb[i][j] + deriv.dehb[i][j] +
291                                 deriv.destor[i][j] + deriv.deqq[i][j] + deriv.devdw[i][j] + deriv.de14[i][j] +
292 <                               deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j];
292 >                               deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
293  
294           }
295        }
# Line 458 | Line 299
299   {
300      double etot;
301      int i;
461 //    struct timeb start,end;
302      
303        energies.total = 0.0;
304        energies.estr = 0.0;
# Line 476 | Line 316
316        energies.eopb = 0.0;
317        energies.eurey = 0.0;
318        energies.esolv = 0.0;
319 +      energies.egeom = 0.0;
320  
321        for (i=1; i <= natom; i++)
322             atom[i].energy = 0.0;
# Line 488 | Line 329
329        // mmff
330        if (pot.use_hal) ehal();
331        if (pot.use_bufcharge) ebufcharge();
332 +
333 +      if (pot.use_geom) egeom();
334                                
335        energies.total =  energies.estr + energies.ebend + energies.etor + energies.estrbnd
336               + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
337 <               energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv;
337 >               energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
338        etot = energies.total;
339        return (etot);
340   }
# Line 517 | Line 360
360        energies.eopb = 0.0;
361        energies.eurey = 0.0;
362        energies.esolv = 0.0;
363 +      energies.egeom = 0.0;
364  
365        for (i=1; i <= natom; i++)
366             atom[i].energy = 0.0;
# Line 530 | Line 374
374        if (pot.use_hal) ehal();
375        if (pot.use_bufcharge) ebufcharge();
376        
377 +      if (pot.use_geom) egeom();
378 +
379        energies.total =  energies.estr + energies.ebend + energies.etor + energies.estrbnd
380               + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
381 <               energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv;
381 >               energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
382        etot = energies.total;
383        
384        return (etot);
385   }    
386  
387   /* ================================================================== */
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 }
388   int isbond(int i, int j)
389   {
390      int k;
# Line 634 | Line 395
395      }
396      return FALSE;
397   }
398 <
398 > // ==========================
399   void get_bonds()
400   {
401      int i, j;
641    long int mask;
402  
403 <    bonds_ff.nbnd = 0;
644 <    mask = 1L << 0;     // flag for pi atoms
645 <    
403 >    bonds_ff.nbnd = 0;    
404      for (i=1; i <= natom; i++)
405      {
406          for (j=i+1; j <= natom; j++)
# Line 654 | Line 412
412                bonds_ff.nbnd++;
413                if (bonds_ff.nbnd > MAXBND)
414                {
415 <                 fprintf(pcmoutfile,"Error - Too many bonds!\nProgram will now quit.\n");
415 >                 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
416                   exit(0);
417                 }
418             }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines