ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcm7.c
(Generate patch)
# Line 4 | Line 4
4   #include "pcmod.h"
5  
6   #include "energies.h"
7 #include "angles.h"
8 #include "rings.h"
9 #include "torsions.h"
7   #include "bonds_ff.h"
8   #include "derivs.h"
12 #include "hess.h"
9   #include "pot.h"
14 #include "gmmx.h"
10   #include "utility.h"
11   #include "solv.h"  
12 < #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;
12 > #include "dipmom.h"
13                  
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
14   EXTERN struct t_minim_control {
15          int type, method, field, added_const;
16          char added_path[256],added_name[256];
17          } minim_control;
42
43 struct t_dipolemom {
44        double total, xdipole, ydipole, zdipole;
45       }  dipolemom;
18          
19   EXTERN struct t_minim_values {
20          int iprint, ndc, nconst;
# Line 59 | Line 31
31      
32   FILE * fopen_path ( char * , char * , char * ) ;
33   void message_alert(char *,char *);
62 void get_added_const(void);
63 void assign_gast_charge();
64 void lattice(void);
65 void bounds(void);
34   void pcm7_min(void);
67 int iscoord_bond(int, int);
68 void get_coordbonds(void);
35   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);
36   int initialise(void);
37   void get_bonds(void);
78 void get_angles(void);
38   void get_torsions(void);
39   void get_rings(void);
40   void set_field(void);
41   double energy(void);
42   double print_energy(void);
43   void pcmfout(int);
85 void orbital(void);
86 void free_orbit(void);
44   void set_active(void);
45 < void estrtor(void);
46 < void echarge(void);
47 < void ewald(void);
48 < void eurey(void);
49 < void ebuck(void);
50 < void ebuckmm3(void);
45 > int setup_calculation(void);
46 > void end_calculation(void);
47 >
48 > void egeom1(void);
49 > void esolv1(void);
50 > void egeom(void);
51 > void esolv(void);
52 >
53 > void ehal(void);
54 > void ehal1(void);
55 > void ebufcharge(void);
56 > void ebufcharge1(void);
57 > void eopbend_wilson(void);
58 > void eopbend_wilson1(void);
59 >
60   void estrbnd(void);
95 void edipole(void);
61   void etorsion(void);
62   void eangang(void);
63   void eangle(void);
64   void ebond(void);
65 < 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);
65 >
66   void estrbnd1(void);
115 void edipole1(void);
67   void etorsion1(void);
117 void eangang1(void);
68   void eangle1(void);
69   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);
70  
71   void kangle(void);
133 void korbit(void);
72   void ktorsion(void);
135 void kdipole(void);
73   void kstrbnd(void);
74   void kcharge(void);
75   void ksolv(void);
139 void kangang(void);
76   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);
77   void kvdw(void);
78   int  kbond(void);
149 void kewald(void);
150 void kurey(void);
79  
80   void get_memory(void);
81   void free_memory(void);
82   void gradient(void);
83   void minimize(void);
156 void dynamics(void);
84   void read_datafiles(char *);
85 + void get_added_const(void);
86   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);
87   void hdel(int);
88   void set_atomtypes(int);
89   void type(void);
90   void generate_bonds(void);
91   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);
92   void xlogp(float *);
93  
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) charge_dipole();
238      if (pot.use_dipole) dipole_dipole();
239      dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
240                             dipolemom.ydipole*dipolemom.ydipole +
241                             dipolemom.zdipole*dipolemom.zdipole);
242 }
94   // ==================================================================
95   int setup_calculation()
96   {
# Line 277 | Line 128
128          generate_bonds();
129          if (minim_control.added_const)
130             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();
131      } else
132      {
133          if (minim_control.added_const)
# Line 321 | Line 150
150    get_bonds();
151    get_angles();
152    get_torsions();
324  get_coordbonds();
325
326  if (high_coord.ncoord > 0)
327      pot.use_highcoord = TRUE;
153  
154    attach();
330
331  
155    get_rings();
156    // need allene
157  
# Line 342 | Line 165
165       if (nRet == FALSE) // use end_calc to free memory
166       {
167         //         message_alert("Parameters missing. See pcmod.err for information","Error");
345         optimize_data.param_avail = 1;
168           //         fclose(errfile);
169           energies.total = 9999.;
170           return FALSE;
# Line 359 | Line 181
181         {
182           //           message_alert("Parameters missing. See pcmod.err for information","Error");
183           //           fclose(errfile);
362           optimize_data.param_avail = 1;
184             energies.total = 9999.0;
185             return (FALSE);
186         } else
# Line 549 | Line 370
370   }    
371  
372   /* ================================================================== */
552
553 int iscoord_bond(int i, int j)
554 {
555    int k;
556    long int mask;
557
558    mask = 1 << METCOORD_MASK ;
559
560    for (k=0; k < MAXIAT; k++)
561    {
562        if (atom[i].iat[k] == j && atom[i].bo[k] == 9)
563        {
564            if (atom[i].type >= 300)
565            {
566                if (atom[j].type == 2 || atom[j].type == 4 || atom[j].type == 29 ||
567                    atom[j].type == 30 || atom[j].type == 40 || atom[j].type == 48 )
568                      return TRUE;
569            }
570            if (atom[j].type >= 300)
571            {
572                if (atom[i].type == 2 || atom[i].type == 4 || atom[i].type == 29 ||
573                    atom[i].type == 30 || atom[i].type == 40 || atom[i].type == 48 )
574                      return TRUE;
575            }
576        }
577    }
578 // throughout Metal-lone pairs if coordinated bond
579 //  metal to donor atom
580    if ( (atom[i].type >= 300) && atom[j].flags & mask)
581      return TRUE;
582    if ( (atom[j].type >= 300) && atom[i].flags & mask)
583      return TRUE;
584  
585    return FALSE;
586 }
587
588 void get_coordbonds()
589 {
590    int i, j, k, iatype, jatm;
591
592    cbonds.nbnd = 0;
593    pot.use_coordb = FALSE;
594        
595    for (i=1; i <= natom; i++)
596    {
597        if (atom[i].mmx_type >= 300)
598        {
599           for (j=0; j <= MAXIAT; j++)
600           {
601               if (atom[i].bo[j] == 9)  // coord bond
602               {
603                   iatype = atom[atom[i].iat[j]].mmx_type;
604                   if (iatype == 2 || iatype == 4 || iatype == 29 ||
605                       iatype == 30 || iatype == 40 || iatype == 48 )
606                   {
607                      cbonds.icbonds[cbonds.nbnd][0] = i;
608                      cbonds.icbonds[cbonds.nbnd][1] = atom[i].iat[j];
609                      cbonds.icbonds[cbonds.nbnd][2] = 0;
610                      cbonds.nbnd++;
611                   } else
612                   {
613                      jatm = atom[i].iat[j];
614                      for (k = 0; k < MAXIAT; k++)
615                      {
616                          if (atom[atom[jatm].iat[k]].mmx_type == 20) // lp
617                          {
618                              cbonds.icbonds[cbonds.nbnd][0] = i;
619                              cbonds.icbonds[cbonds.nbnd][1] = jatm;
620                              cbonds.icbonds[cbonds.nbnd][2] = atom[jatm].iat[k];
621                              cbonds.nbnd++;
622                          }
623                      }
624                   }
625                   if (cbonds.nbnd > MAXCBND)
626                   {
627                      fprintf(pcmlogfile,"Error - Too many coordinated bonds!\nProgram will now quit.\n");
628                       exit(0);
629                   }
630               }
631           }
632        }
633    }
634    if (cbonds.nbnd > 0)
635       pot.use_coordb = TRUE;
636 }
373   int isbond(int i, int j)
374   {
375      int k;
# Line 644 | Line 380
380      }
381      return FALSE;
382   }
383 <
383 > // ==========================
384   void get_bonds()
385   {
386      int i, j;
651    long int mask;
387  
388 <    bonds_ff.nbnd = 0;
654 <    mask = 1L << 0;     // flag for pi atoms
655 <    
388 >    bonds_ff.nbnd = 0;    
389      for (i=1; i <= natom; i++)
390      {
391          for (j=i+1; j <= natom; j++)

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines