mengine/src/pcm7.c

#define EXTERN extern

#include "pcwin.h"
#include "pcmod.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;
FILE *errfile;
    
FILE * fopen_path ( char * , char * , char * ) ;
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 generate_bonds(void);
void zero_data(void);
void xlogp(float *);

// ==================================================================
int setup_calculation()
{
    int nRet;
    char string[30];
    double etot;
    if (minim_control.field == MM3)
    {
        // copy mm3 types into type file and retype
        pot.use_hbond = FALSE;
        hbondreset();
        hdel(1);
        strcpy(string,"mm3.prm");
        zero_data();
        read_datafiles(string);
        set_field();
        default_intype = MM3;
        default_outtype = MM3;
        set_atomtypes(MM3);
        generate_bonds();
        if (minim_control.added_const)
           get_added_const();
    } else if (minim_control.field == MMFF94)
    {
        // copy mm3 types into type file and retype
        pot.use_hbond = FALSE;
        pot.use_picalc = FALSE;
        hbondreset();
        pireset();
        hdel(1);
        set_field();
        default_intype = MMFF94;
        default_outtype = MMFF94;
        set_atomtypes(MMFF94);
        generate_bonds();
        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();
        generate_bonds();
    }

            
  // 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();

/*  assign local geometry potential function parameters  */
    Missing_constants = FALSE;
    //     errfile = fopen_path(pcwindir,"pcmod.err","w");
     if (pot.use_bond || pot.use_strbnd)  nRet = kbond();
     if (nRet == FALSE) // use end_calc to free memory
     {
       //         message_alert("Parameters missing. See pcmod.err for information","Error");
         //         fclose(errfile);
         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)
       {
         //           message_alert("Parameters missing. See pcmod.err for information","Error");
         //           fclose(errfile);
           energies.total = 9999.0;
           return (FALSE);
       } else
       {
         //           fclose(errfile);
         //           remove("pcmod.err");
       }    
       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;

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