mengine/src/pcm7.c

#define EXTERN extern

#include "pcwin.h"
#include "pcmod.h"

#include "attached.h"
#include "energies.h"
#include "bonds_ff.h"
#include "derivs.h"
#include "pot.h"
#include "utility.h"
#include "solv.h"   
#include "dipmom.h"
                
EXTERN struct t_minim_control {
        int type, method, field, added_const;
        char added_path[256],added_name[256];
        } minim_control;
        
EXTERN struct t_minim_values {
        int iprint, ndc, nconst;
        float dielc;
        } minim_values;
struct t_solvent {
    int type;
    double doffset, p1,p2,p3,p4,p5;
    double *shct,*asolv,*rsolv,*vsolv,*gpol,*rborn;
    } solvent;

int Missing_constants;
    
void message_alert(char *,char *);
void pcm7_min(void);
int isbond(int, int);
int initialise(void);
void get_bonds(void);
void get_torsions(void);
void get_rings(void);
void set_field(void);
double energy(void);
double print_energy(void);
void pcmfout(int);
void set_active(void);
int setup_calculation(void);
void end_calculation(void);

void egeom1(void);
void esolv1(void);
void egeom(void);
void esolv(void);

void ehal(void);
void ehal1(void);
void ebufcharge(void);
void ebufcharge1(void);
void eopbend_wilson(void);
void eopbend_wilson1(void);

void estrbnd(void);
void etorsion(void);
void eangang(void);
void eangle(void);
void ebond(void);

void estrbnd1(void);
void etorsion1(void);
void eangle1(void);
void ebond1(void);

void kangle(void);
void ktorsion(void);
void kstrbnd(void);
void kcharge(void);
void ksolv(void);
void kopbend(void);
void kvdw(void);
int  kbond(void);

void get_memory(void);
void free_memory(void);
void gradient(void);
void minimize(void);
void read_datafiles(char *);
void get_added_const(void);
void attach(void);
void hdel(int);
void set_atomtypes(int);
void type(void);
void zero_data(void);
void xlogp(float *);

// ==================================================================
int setup_calculation()
{
  int i,j,nRet;
    char string[30];
    double etot;
    if (minim_control.field == MMFF94)
    {
        // copy mm3 types into type file and retype
        pot.use_hbond = FALSE;
        pot.use_picalc = FALSE;
        hdel(1);
        set_field();
        set_atomtypes(MMFF94);
        if (minim_control.added_const)
           get_added_const();
    } else
    {
        if (minim_control.added_const)
        {
           strcpy(string,"mmxconst.prm");
           zero_data();
           read_datafiles(string);
           get_added_const();
        }
        type();
        set_field();
    }

            
  // allocate memeory for derivatives
  get_memory();
 
  // setup bonds, angles, torsions, improper torsions and angles, allenes
  get_bonds();
  get_angles();
  get_torsions();

  attach();
  get_rings();
  // need allene

  // set active atoms
  set_active();
  // setup non_bonded list of atoms to skip
  for (i=1; i <= natom; i++)
    {
         for (j=0; j < MAXIAT; j++)
             if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
               {
                    skip[i][atom[i].iat[j]] = i;
                    skip[atom[i].iat[j]][i] = atom[i].iat[j];
               }
         for (j=0; j < attached.n13[i]; j++)
           {
               skip[i][attached.i13[j][i]] = i;
               skip[attached.i13[j][i]][i] = attached.i13[j][i];
           }
         for(j=0; j < attached.n14[i]; j++)
           {
            skip[i][attached.i14[j][i]] = -i;
            skip[attached.i14[j][i]][i] = -attached.i14[j][i];
           }
    }
/*  assign local geometry potential function parameters  */
    Missing_constants = FALSE;
     if (pot.use_bond || pot.use_strbnd)  nRet = kbond();
     if (nRet == FALSE) // use end_calc to free memory
     {
         energies.total = 9999.;
         return FALSE;
     }
     if (pot.use_angle || pot.use_strbnd || pot.use_angang) kangle();
     if (pot.use_angle || pot.use_opbend || pot.use_opbend_wilson) kopbend();
     if (pot.use_tors  || pot.use_strtor || pot.use_tortor) ktorsion();
     if (pot.use_strbnd) kstrbnd();
      
      if (pot.use_buck || pot.use_hal || pot.use_lj) kvdw();
      if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();

       if (Missing_constants == TRUE)
       {
           energies.total = 9999.0;
           return (FALSE);
       }    
       if (minim_values.iprint == TRUE)
          etot = print_energy();
       else
          etot = energy();
       return TRUE;
}
// =================================================================
void end_calculation()
{
 
    free_memory();
}
// ==================================================================
void gradient()
{
    int i, j;

      energies.estr = 0.0;
      energies.ebend = 0.0;
      energies.estrbnd = 0.0;
      energies.e14 = 0.0;
      energies.evdw = 0.0;
      energies.etor = 0.0;
      energies.eu = 0.0;
      energies.eopb = 0.0;
      energies.eangang = 0.0;
      energies.estrtor = 0.0;
      energies.ehbond = 0.0;
      energies.efix = 0.0;
      energies.eimprop = 0.0;
      energies.eimptors = 0.0;
      energies.eurey = 0.0;
      energies.esolv = 0.0;
      energies.egeom =0.0;
      
      for (i=1; i <= natom; i++)
           atom[i].energy = 0.0;
      
      for (i=1; i <= natom; i++)
      {
        for (j=0; j < 3; j++)
        {
            deriv.deb[i][j] = 0.0;
            deriv.dea[i][j] = 0.0;
            deriv.destb[i][j] = 0.0;
            deriv.deopb[i][j] = 0.0;
            deriv.detor[i][j] = 0.0;
            deriv.de14[i][j] = 0.0;
            deriv.devdw[i][j]= 0.0;
            deriv.deqq[i][j] = 0.0;
            deriv.deaa[i][j] = 0.0;
            deriv.destor[i][j] = 0.0;
            deriv.deimprop[i][j] = 0.0;
            deriv.dehb[i][j] = 0.0;
            deriv.desolv[i][j] = 0.0;
            deriv.degeom[i][j] = 0.0;
        }
      }
   
      if (pot.use_bond)ebond1();
      if (pot.use_angle)eangle1();
      if (pot.use_opbend_wilson) eopbend_wilson1();
      if (pot.use_tors)etorsion1();
      if (pot.use_strbnd)estrbnd1();
     
      // mmff 
      if (pot.use_hal) ehal1();
      if (pot.use_bufcharge) ebufcharge1();

      if (pot.use_geom) egeom1();
     
      energies.total =  energies.estr + energies.ebend + energies.etor + energies.estrbnd + energies.e14+
                        energies.evdw + energies.eu + energies.ehbond + energies.eangang + energies.estrtor +
                        energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
      for (i=1; i <= natom; i++)
      {
          for (j=0; j < 3; j++)
          {
              deriv.d1[i][j] = deriv.deb[i][j] + deriv.dea[i][j] + deriv.deaa[i][j] +
                               deriv.destb[i][j] + deriv.detor[i][j] + deriv.deopb[i][j] + deriv.dehb[i][j] +
                               deriv.destor[i][j] + deriv.deqq[i][j] + deriv.devdw[i][j] + deriv.de14[i][j] +
                               deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];

         }
      }
}
// ==================================================================
double energy()
{
    double etot;
    int i;
    
      energies.total = 0.0;
      energies.estr = 0.0;
      energies.ebend = 0.0;
      energies.etor = 0.0;
      energies.estrbnd = 0.0;
      energies.evdw = 0.0;
      energies.e14 = 0.0;
      energies.ehbond = 0.0;
      energies.eu = 0.0;
      energies.eangang = 0.0;
      energies.estrtor = 0.0;
      energies.eimprop = 0.0;
      energies.eimptors = 0.0;
      energies.eopb = 0.0;
      energies.eurey = 0.0;
      energies.esolv = 0.0;
      energies.egeom = 0.0;

      for (i=1; i <= natom; i++)
           atom[i].energy = 0.0;

      if (pot.use_bond)ebond();
      if (pot.use_angle)eangle();
      if (pot.use_opbend_wilson) eopbend_wilson();
      if (pot.use_tors)etorsion();
      if (pot.use_strbnd)estrbnd();
      // mmff 
      if (pot.use_hal) ehal();
      if (pot.use_bufcharge) ebufcharge();

      if (pot.use_geom) egeom();
                              
      energies.total =  energies.estr + energies.ebend + energies.etor + energies.estrbnd
             + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
               energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
      etot = energies.total;
      return (etot);
}
/* =========================================================================== */ 
double print_energy()
{
    double etot;
    int i;
    
      energies.total = 0.0;
      energies.estr = 0.0;
      energies.ebend = 0.0;
      energies.etor = 0.0;
      energies.estrbnd = 0.0;
      energies.evdw = 0.0;
      energies.e14 = 0.0;
      energies.ehbond = 0.0;
      energies.eu = 0.0;
      energies.eangang = 0.0;
      energies.estrtor = 0.0;
      energies.eimprop = 0.0;
      energies.eimptors = 0.0;
      energies.eopb = 0.0;
      energies.eurey = 0.0;
      energies.esolv = 0.0;
      energies.egeom = 0.0;

      for (i=1; i <= natom; i++)
           atom[i].energy = 0.0;

      if (pot.use_bond)ebond();
      if (pot.use_angle)eangle();
      if (pot.use_opbend_wilson) eopbend_wilson();
      if (pot.use_tors)etorsion();
      if (pot.use_strbnd)estrbnd();
      // mmff 
      if (pot.use_hal) ehal();
      if (pot.use_bufcharge) ebufcharge();
      
      if (pot.use_geom) egeom();

      energies.total =  energies.estr + energies.ebend + energies.etor + energies.estrbnd
             + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
               energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
      etot = energies.total;
      
      return (etot);
}    

/* ================================================================== */
int isbond(int i, int j)
{
    int k;
    for (k=0; k < MAXIAT; k++)
    {
        if (atom[i].iat[k] == j && atom[i].bo[k] != 9)
            return TRUE;
    }
    return FALSE;
}
// ==========================
void get_bonds()
{
    int i, j;

    bonds_ff.nbnd = 0;    
    for (i=1; i <= natom; i++)
    {
        for (j=i+1; j <= natom; j++)
        {
           if (isbond(i,j))
           {
              bonds_ff.i12[bonds_ff.nbnd][0] = i;
              bonds_ff.i12[bonds_ff.nbnd][1] = j;
              bonds_ff.nbnd++;
              if (bonds_ff.nbnd > MAXBND)
              {
                 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
                 exit(0);
               }
           }
        }
    }
}
Revision:	93
Committed:	Sat Jan 17 23:24:55 2009 UTC (13 years, 4 months ago) by wdelano
File size:	10584 byte(s)
Log Message:	branch created
Line	File contents
1	#define EXTERN extern
2
3	#include "pcwin.h"
4	#include "pcmod.h"
5
6	#include "attached.h"
7	#include "energies.h"
8	#include "bonds_ff.h"
9	#include "derivs.h"
10	#include "pot.h"
11	#include "utility.h"
12	#include "solv.h"
13	#include "dipmom.h"
14
15	EXTERN struct t_minim_control {
16	int type, method, field, added_const;
17	char added_path[256],added_name[256];
18	} minim_control;
19
20	EXTERN struct t_minim_values {
21	int iprint, ndc, nconst;
22	float dielc;
23	} minim_values;
24	struct t_solvent {
25	int type;
26	double doffset, p1,p2,p3,p4,p5;
27	double shct,asolv,rsolv,vsolv,gpol,rborn;
28	} solvent;
29
30	int Missing_constants;
31
32	void message_alert(char ,char );
33	void pcm7_min(void);
34	int isbond(int, int);
35	int initialise(void);
36	void get_bonds(void);
37	void get_torsions(void);
38	void get_rings(void);
39	void set_field(void);
40	double energy(void);
41	double print_energy(void);
42	void pcmfout(int);
43	void set_active(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);
60	void etorsion(void);
61	void eangang(void);
62	void eangle(void);
63	void ebond(void);
64
65	void estrbnd1(void);
66	void etorsion1(void);
67	void eangle1(void);
68	void ebond1(void);
69
70	void kangle(void);
71	void ktorsion(void);
72	void kstrbnd(void);
73	void kcharge(void);
74	void ksolv(void);
75	void kopbend(void);
76	void kvdw(void);
77	int kbond(void);
78
79	void get_memory(void);
80	void free_memory(void);
81	void gradient(void);
82	void minimize(void);
83	void read_datafiles(char *);
84	void get_added_const(void);
85	void attach(void);
86	void hdel(int);
87	void set_atomtypes(int);
88	void type(void);
89	void zero_data(void);
90	void xlogp(float *);
91
92	// ==================================================================
93	int setup_calculation()
94	{
95	int i,j,nRet;
96	char string[30];
97	double etot;
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;
103	hdel(1);
104	set_field();
105	set_atomtypes(MMFF94);
106	if (minim_control.added_const)
107	get_added_const();
108	} else
109	{
110	if (minim_control.added_const)
111	{
112	strcpy(string,"mmxconst.prm");
113	zero_data();
114	read_datafiles(string);
115	get_added_const();
116	}
117	type();
118	set_field();
119	}
120
121
122	// allocate memeory for derivatives
123	get_memory();
124
125	// setup bonds, angles, torsions, improper torsions and angles, allenes
126	get_bonds();
127	get_angles();
128	get_torsions();
129
130	attach();
131	get_rings();
132	// need allene
133
134	// set active atoms
135	set_active();
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;
158	if (pot.use_bond \|\| pot.use_strbnd) nRet = kbond();
159	if (nRet == FALSE) // use end_calc to free memory
160	{
161	energies.total = 9999.;
162	return FALSE;
163	}
164	if (pot.use_angle \|\| pot.use_strbnd \|\| pot.use_angang) kangle();
165	if (pot.use_angle \|\| pot.use_opbend \|\| pot.use_opbend_wilson) kopbend();
166	if (pot.use_tors \|\| pot.use_strtor \|\| pot.use_tortor) ktorsion();
167	if (pot.use_strbnd) kstrbnd();
168
169	if (pot.use_buck \|\| pot.use_hal \|\| pot.use_lj) kvdw();
170	if ((pot.use_charge \|\| pot.use_bufcharge \|\| pot.use_ewald)) kcharge();
171
172	if (Missing_constants == TRUE)
173	{
174	energies.total = 9999.0;
175	return (FALSE);
176	}
177	if (minim_values.iprint == TRUE)
178	etot = print_energy();
179	else
180	etot = energy();
181	return TRUE;
182	}
183	// =================================================================
184	void end_calculation()
185	{
186
187	free_memory();
188	}
189	// ==================================================================
190	void gradient()
191	{
192	int i, j;
193
194	energies.estr = 0.0;
195	energies.ebend = 0.0;
196	energies.estrbnd = 0.0;
197	energies.e14 = 0.0;
198	energies.evdw = 0.0;
199	energies.etor = 0.0;
200	energies.eu = 0.0;
201	energies.eopb = 0.0;
202	energies.eangang = 0.0;
203	energies.estrtor = 0.0;
204	energies.ehbond = 0.0;
205	energies.efix = 0.0;
206	energies.eimprop = 0.0;
207	energies.eimptors = 0.0;
208	energies.eurey = 0.0;
209	energies.esolv = 0.0;
210	energies.egeom =0.0;
211
212	for (i=1; i <= natom; i++)
213	atom[i].energy = 0.0;
214
215	for (i=1; i <= natom; i++)
216	{
217	for (j=0; j < 3; j++)
218	{
219	deriv.deb[i][j] = 0.0;
220	deriv.dea[i][j] = 0.0;
221	deriv.destb[i][j] = 0.0;
222	deriv.deopb[i][j] = 0.0;
223	deriv.detor[i][j] = 0.0;
224	deriv.de14[i][j] = 0.0;
225	deriv.devdw[i][j]= 0.0;
226	deriv.deqq[i][j] = 0.0;
227	deriv.deaa[i][j] = 0.0;
228	deriv.destor[i][j] = 0.0;
229	deriv.deimprop[i][j] = 0.0;
230	deriv.dehb[i][j] = 0.0;
231	deriv.desolv[i][j] = 0.0;
232	deriv.degeom[i][j] = 0.0;
233	}
234	}
235
236	if (pot.use_bond)ebond1();
237	if (pot.use_angle)eangle1();
238	if (pot.use_opbend_wilson) eopbend_wilson1();
239	if (pot.use_tors)etorsion1();
240	if (pot.use_strbnd)estrbnd1();
241
242	// mmff
243	if (pot.use_hal) ehal1();
244	if (pot.use_bufcharge) ebufcharge1();
245
246	if (pot.use_geom) egeom1();
247
248	energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd + energies.e14+
249	energies.evdw + energies.eu + energies.ehbond + energies.eangang + energies.estrtor +
250	energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
251	for (i=1; i <= natom; i++)
252	{
253	for (j=0; j < 3; j++)
254	{
255	deriv.d1[i][j] = deriv.deb[i][j] + deriv.dea[i][j] + deriv.deaa[i][j] +
256	deriv.destb[i][j] + deriv.detor[i][j] + deriv.deopb[i][j] + deriv.dehb[i][j] +
257	deriv.destor[i][j] + deriv.deqq[i][j] + deriv.devdw[i][j] + deriv.de14[i][j] +
258	deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
259
260	}
261	}
262	}
263	// ==================================================================
264	double energy()
265	{
266	double etot;
267	int i;
268
269	energies.total = 0.0;
270	energies.estr = 0.0;
271	energies.ebend = 0.0;
272	energies.etor = 0.0;
273	energies.estrbnd = 0.0;
274	energies.evdw = 0.0;
275	energies.e14 = 0.0;
276	energies.ehbond = 0.0;
277	energies.eu = 0.0;
278	energies.eangang = 0.0;
279	energies.estrtor = 0.0;
280	energies.eimprop = 0.0;
281	energies.eimptors = 0.0;
282	energies.eopb = 0.0;
283	energies.eurey = 0.0;
284	energies.esolv = 0.0;
285	energies.egeom = 0.0;
286
287	for (i=1; i <= natom; i++)
288	atom[i].energy = 0.0;
289
290	if (pot.use_bond)ebond();
291	if (pot.use_angle)eangle();
292	if (pot.use_opbend_wilson) eopbend_wilson();
293	if (pot.use_tors)etorsion();
294	if (pot.use_strbnd)estrbnd();
295	// mmff
296	if (pot.use_hal) ehal();
297	if (pot.use_bufcharge) ebufcharge();
298
299	if (pot.use_geom) egeom();
300
301	energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
302	+ energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
303	energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
304	etot = energies.total;
305	return (etot);
306	}
307	/* =========================================================================== */
308	double print_energy()
309	{
310	double etot;
311	int i;
312
313	energies.total = 0.0;
314	energies.estr = 0.0;
315	energies.ebend = 0.0;
316	energies.etor = 0.0;
317	energies.estrbnd = 0.0;
318	energies.evdw = 0.0;
319	energies.e14 = 0.0;
320	energies.ehbond = 0.0;
321	energies.eu = 0.0;
322	energies.eangang = 0.0;
323	energies.estrtor = 0.0;
324	energies.eimprop = 0.0;
325	energies.eimptors = 0.0;
326	energies.eopb = 0.0;
327	energies.eurey = 0.0;
328	energies.esolv = 0.0;
329	energies.egeom = 0.0;
330
331	for (i=1; i <= natom; i++)
332	atom[i].energy = 0.0;
333
334	if (pot.use_bond)ebond();
335	if (pot.use_angle)eangle();
336	if (pot.use_opbend_wilson) eopbend_wilson();
337	if (pot.use_tors)etorsion();
338	if (pot.use_strbnd)estrbnd();
339	// mmff
340	if (pot.use_hal) ehal();
341	if (pot.use_bufcharge) ebufcharge();
342
343	if (pot.use_geom) egeom();
344
345	energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
346	+ energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
347	energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
348	etot = energies.total;
349
350	return (etot);
351	}
352
353	/* ================================================================== */
354	int isbond(int i, int j)
355	{
356	int k;
357	for (k=0; k < MAXIAT; k++)
358	{
359	if (atom[i].iat[k] == j && atom[i].bo[k] != 9)
360	return TRUE;
361	}
362	return FALSE;
363	}
364	// ==========================
365	void get_bonds()
366	{
367	int i, j;
368
369	bonds_ff.nbnd = 0;
370	for (i=1; i <= natom; i++)
371	{
372	for (j=i+1; j <= natom; j++)
373	{
374	if (isbond(i,j))
375	{
376	bonds_ff.i12[bonds_ff.nbnd][0] = i;
377	bonds_ff.i12[bonds_ff.nbnd][1] = j;
378	bonds_ff.nbnd++;
379	if (bonds_ff.nbnd > MAXBND)
380	{
381	fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
382	exit(0);
383	}
384	}
385	}
386	}
387	}