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 "energies.h"
9 #include "rings.h"
10 #include "torsions.h"
7   #include "bonds_ff.h"
8   #include "derivs.h"
13 #include "hess.h"
9   #include "pot.h"
15 #include "gmmx.h"
10   #include "utility.h"
11   #include "solv.h"  
18 #include "field.h"
12   #include "dipmom.h"
13                  
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
14   EXTERN struct t_minim_control {
15          int type, method, field, added_const;
16          char added_path[256],added_name[256];
# Line 46 | Line 31
31      
32   FILE * fopen_path ( char * , char * , char * ) ;
33   void message_alert(char *,char *);
49 void get_added_const(void);
50 void assign_gast_charge();
51 void lattice(void);
52 void bounds(void);
34   void pcm7_min(void);
54 int iscoord_bond(int, int);
55 void get_coordbonds(void);
35   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);
36   int initialise(void);
37   void get_bonds(void);
38   void get_torsions(void);
# Line 67 | Line 41
41   double energy(void);
42   double print_energy(void);
43   void pcmfout(int);
70 void orbital(void);
71 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);
80 void edipole(void);
61   void etorsion(void);
62   void eangang(void);
63   void eangle(void);
64   void ebond(void);
65 < 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);
65 >
66   void estrbnd1(void);
100 void edipole1(void);
67   void etorsion1(void);
102 void eangang1(void);
68   void eangle1(void);
69   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);
70  
71   void kangle(void);
118 void korbit(void);
72   void ktorsion(void);
120 void kdipole(void);
73   void kstrbnd(void);
74   void kcharge(void);
75   void ksolv(void);
124 void kangang(void);
76   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);
77   void kvdw(void);
78   int  kbond(void);
134 void kewald(void);
135 void kurey(void);
79  
80   void get_memory(void);
81   void free_memory(void);
82   void gradient(void);
83   void minimize(void);
141 void dynamics(void);
84   void read_datafiles(char *);
85 + void get_added_const(void);
86   void attach(void);
144 void eheat(void);
145 void pirite(void);
146 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);
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);
92   void xlogp(float *);
93  
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 }
94   // ==================================================================
95   int setup_calculation()
96   {
# Line 256 | Line 128
128          generate_bonds();
129          if (minim_control.added_const)
130             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();
131      } else
132      {
133          if (minim_control.added_const)
# Line 300 | Line 150
150    get_bonds();
151    get_angles();
152    get_torsions();
303  get_coordbonds();
304
305  if (high_coord.ncoord > 0)
306      pot.use_highcoord = TRUE;
153  
154    attach();
309
310  
155    get_rings();
156    // need allene
157  
# Line 526 | Line 370
370   }    
371  
372   /* ================================================================== */
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 }
373   int isbond(int i, int j)
374   {
375      int k;
# Line 621 | Line 380
380      }
381      return FALSE;
382   }
383 <
383 > // ==========================
384   void get_bonds()
385   {
386      int i, j;
628    long int mask;
387  
388 <    bonds_ff.nbnd = 0;
631 <    mask = 1L << 0;     // flag for pi atoms
632 <    
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