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 46 | Line 32
32      
33   FILE * fopen_path ( char * , char * , char * ) ;
34   void message_alert(char *,char *);
49 void get_added_const(void);
50 void assign_gast_charge();
51 void lattice(void);
52 void bounds(void);
35   void pcm7_min(void);
54 int iscoord_bond(int, int);
55 void get_coordbonds(void);
36   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);
37   int initialise(void);
38   void get_bonds(void);
39   void get_torsions(void);
# Line 67 | Line 42
42   double energy(void);
43   double print_energy(void);
44   void pcmfout(int);
70 void orbital(void);
71 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);
80 void edipole(void);
62   void etorsion(void);
63   void eangang(void);
64   void eangle(void);
65   void ebond(void);
66 < 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);
66 >
67   void estrbnd1(void);
100 void edipole1(void);
68   void etorsion1(void);
102 void eangang1(void);
69   void eangle1(void);
70   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);
71  
72   void kangle(void);
118 void korbit(void);
73   void ktorsion(void);
120 void kdipole(void);
74   void kstrbnd(void);
75   void kcharge(void);
76   void ksolv(void);
124 void kangang(void);
77   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);
78   void kvdw(void);
79   int  kbond(void);
134 void kewald(void);
135 void kurey(void);
80  
81   void get_memory(void);
82   void free_memory(void);
83   void gradient(void);
84   void minimize(void);
141 void dynamics(void);
85   void read_datafiles(char *);
86 + void get_added_const(void);
87   void attach(void);
144 void eheat(void);
145 void pirite(void);
146 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);
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);
93   void xlogp(float *);
94  
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 }
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 256 | Line 129
129          generate_bonds();
130          if (minim_control.added_const)
131             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();
132      } else
133      {
134          if (minim_control.added_const)
# Line 300 | Line 151
151    get_bonds();
152    get_angles();
153    get_torsions();
303  get_coordbonds();
304
305  if (high_coord.ncoord > 0)
306      pot.use_highcoord = TRUE;
154  
155    attach();
309
310  
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 378 | Line 242
242        energies.eurey = 0.0;
243        energies.esolv = 0.0;
244        energies.egeom =0.0;
381
382      virial.virx = 0.0;
383      virial.viry = 0.0;
384      virial.virz = 0.0;
245        
246        for (i=1; i <= natom; i++)
247             atom[i].energy = 0.0;
# Line 439 | Line 299
299   {
300      double etot;
301      int i;
442 //    struct timeb start,end;
302      
303        energies.total = 0.0;
304        energies.estr = 0.0;
# Line 526 | Line 385
385   }    
386  
387   /* ================================================================== */
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 }
388   int isbond(int i, int j)
389   {
390      int k;
# Line 621 | Line 395
395      }
396      return FALSE;
397   }
398 <
398 > // ==========================
399   void get_bonds()
400   {
401      int i, j;
628    long int mask;
402  
403 <    bonds_ff.nbnd = 0;
631 <    mask = 1L << 0;     // flag for pi atoms
632 <    
403 >    bonds_ff.nbnd = 0;    
404      for (i=1; i <= natom; i++)
405      {
406          for (j=i+1; j <= natom; j++)

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines