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 "angles.h"
7 <
6 > #include "attached.h"
7   #include "energies.h"
9 #include "rings.h"
10 #include "torsions.h"
8   #include "bonds_ff.h"
9   #include "derivs.h"
13 #include "hess.h"
10   #include "pot.h"
15 #include "gmmx.h"
11   #include "utility.h"
12   #include "solv.h"  
18 #include "field.h"
13   #include "dipmom.h"
14                  
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
15   EXTERN struct t_minim_control {
16          int type, method, field, added_const;
17          char added_path[256],added_name[256];
# Line 42 | Line 28
28      } solvent;
29  
30   int Missing_constants;
45 FILE *errfile;
31      
47 FILE * fopen_path ( char * , char * , char * ) ;
32   void message_alert(char *,char *);
49 void get_added_const(void);
50 void assign_gast_charge();
51 void lattice(void);
52 void bounds(void);
33   void pcm7_min(void);
54 int iscoord_bond(int, int);
55 void get_coordbonds(void);
34   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);
35   int initialise(void);
36   void get_bonds(void);
37   void get_torsions(void);
# Line 67 | Line 40
40   double energy(void);
41   double print_energy(void);
42   void pcmfout(int);
70 void orbital(void);
71 void free_orbit(void);
43   void set_active(void);
44 < void estrtor(void);
45 < void echarge(void);
46 < void ewald(void);
47 < void eurey(void);
48 < void ebuck(void);
49 < void ebuckmm3(void);
44 > int setup_calculation(void);
45 > void end_calculation(void);
46 >
47 > void egeom1(void);
48 > void esolv1(void);
49 > void egeom(void);
50 > void esolv(void);
51 >
52 > void ehal(void);
53 > void ehal1(void);
54 > void ebufcharge(void);
55 > void ebufcharge1(void);
56 > void eopbend_wilson(void);
57 > void eopbend_wilson1(void);
58 >
59   void estrbnd(void);
80 void edipole(void);
60   void etorsion(void);
61   void eangang(void);
62   void eangle(void);
63   void ebond(void);
64 < 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);
64 >
65   void estrbnd1(void);
100 void edipole1(void);
66   void etorsion1(void);
102 void eangang1(void);
67   void eangle1(void);
68   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);
69  
70   void kangle(void);
118 void korbit(void);
71   void ktorsion(void);
120 void kdipole(void);
72   void kstrbnd(void);
73   void kcharge(void);
74   void ksolv(void);
124 void kangang(void);
75   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);
76   void kvdw(void);
77   int  kbond(void);
134 void kewald(void);
135 void kurey(void);
78  
79   void get_memory(void);
80   void free_memory(void);
81   void gradient(void);
82   void minimize(void);
141 void dynamics(void);
83   void read_datafiles(char *);
84 + void get_added_const(void);
85   void attach(void);
144 void eheat(void);
145 void pirite(void);
146 void piden(void);
86   void hdel(int);
87   void set_atomtypes(int);
88   void type(void);
150 void generate_bonds(void);
89   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);
90   void xlogp(float *);
91  
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      if (pot.use_charge || pot.use_bufcharge) charge_dipole();
217      if (pot.use_dipole) dipole_dipole();
218      dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
219                             dipolemom.ydipole*dipolemom.ydipole +
220                             dipolemom.zdipole*dipolemom.zdipole);
221 }
92   // ==================================================================
93   int setup_calculation()
94   {
95 <    int nRet;
95 >  int i,j,nRet;
96      char string[30];
97      double etot;
98 <    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)
98 >    if (minim_control.field == MMFF94)
99      {
100          // copy mm3 types into type file and retype
101          pot.use_hbond = FALSE;
102          pot.use_picalc = FALSE;
249        hbondreset();
250        pireset();
103          hdel(1);
104          set_field();
253        default_intype = MMFF94;
254        default_outtype = MMFF94;
105          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();
106          if (minim_control.added_const)
107             get_added_const();
108      } else
# Line 289 | Line 116
116          }
117          type();
118          set_field();
292        generate_bonds();
119      }
120  
121              
# Line 300 | Line 126
126    get_bonds();
127    get_angles();
128    get_torsions();
303  get_coordbonds();
304
305  if (high_coord.ncoord > 0)
306      pot.use_highcoord = TRUE;
129  
130    attach();
309
310  
131    get_rings();
132    // need allene
133  
134    // set active atoms
135    set_active();
136 <
136 >  // setup non_bonded list of atoms to skip
137 >  for (i=1; i <= natom; i++)
138 >    {
139 >         for (j=0; j < MAXIAT; j++)
140 >             if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
141 >               {
142 >                    skip[i][atom[i].iat[j]] = i;
143 >                    skip[atom[i].iat[j]][i] = atom[i].iat[j];
144 >               }
145 >         for (j=0; j < attached.n13[i]; j++)
146 >           {
147 >               skip[i][attached.i13[j][i]] = i;
148 >               skip[attached.i13[j][i]][i] = attached.i13[j][i];
149 >           }
150 >         for(j=0; j < attached.n14[i]; j++)
151 >           {
152 >            skip[i][attached.i14[j][i]] = -i;
153 >            skip[attached.i14[j][i]][i] = -attached.i14[j][i];
154 >           }
155 >    }
156   /*  assign local geometry potential function parameters  */
157      Missing_constants = FALSE;
319    //     errfile = fopen_path(pcwindir,"pcmod.err","w");
158       if (pot.use_bond || pot.use_strbnd)  nRet = kbond();
159       if (nRet == FALSE) // use end_calc to free memory
160       {
323       //         message_alert("Parameters missing. See pcmod.err for information","Error");
324         //         fclose(errfile);
161           energies.total = 9999.;
162           return FALSE;
163       }
# Line 335 | Line 171
171  
172         if (Missing_constants == TRUE)
173         {
338         //           message_alert("Parameters missing. See pcmod.err for information","Error");
339         //           fclose(errfile);
174             energies.total = 9999.0;
175             return (FALSE);
342       } else
343       {
344         //           fclose(errfile);
345         //           remove("pcmod.err");
176         }    
177         if (minim_values.iprint == TRUE)
178            etot = print_energy();
# Line 378 | Line 208
208        energies.eurey = 0.0;
209        energies.esolv = 0.0;
210        energies.egeom =0.0;
381
382      virial.virx = 0.0;
383      virial.viry = 0.0;
384      virial.virz = 0.0;
211        
212        for (i=1; i <= natom; i++)
213             atom[i].energy = 0.0;
# Line 439 | Line 265
265   {
266      double etot;
267      int i;
442 //    struct timeb start,end;
268      
269        energies.total = 0.0;
270        energies.estr = 0.0;
# Line 526 | Line 351
351   }    
352  
353   /* ================================================================== */
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                      fprintf(pcmlogfile,"Error - Too many coordinated bonds!\nProgram will now quit.\n");
605                       exit(0);
606                   }
607               }
608           }
609        }
610    }
611    if (cbonds.nbnd > 0)
612       pot.use_coordb = TRUE;
613 }
354   int isbond(int i, int j)
355   {
356      int k;
# Line 621 | Line 361
361      }
362      return FALSE;
363   }
364 <
364 > // ==========================
365   void get_bonds()
366   {
367      int i, j;
628    long int mask;
368  
369 <    bonds_ff.nbnd = 0;
631 <    mask = 1L << 0;     // flag for pi atoms
632 <    
369 >    bonds_ff.nbnd = 0;    
370      for (i=1; i <= natom; i++)
371      {
372          for (j=i+1; j <= natom; j++)

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines