mengine/src/kstrbnd.c

#define EXTERN extern
#include "pcwin.h"
#include "pcmod.h"

#include "angles.h"
#include "bonds_ff.h"
#include "field.h"

int find_bond(int, int);
void numeral(int, char *, int);

EXTERN struct t_crossterm_k {
        int nangang, nstrbnd, nstrtor;
        int ang_ang[MAXAA], stbnindex[MAXSTBN] ;
        char str_tor[MAXSTRTOR][7], stbn[MAXSTBN][10];
        float aacon[MAXAA][3], stbncon[MAXSTBN][3], str_torcon[MAXSTRTOR];
        } crossterm_k;

struct t_strbnd {
        int nstrbnd, isb[MAXANG][3];
        float ksb1[MAXANG], ksb2[MAXANG];
        } strbnd;
struct stbndat {
        int irow, krow, lrow;
        float con1, con2;
      } stbn[] =
      {
        0,    1,    0,      0.15,      0.15,
        0,    1,    1,      0.10,      0.30,
        0,    1,    2,      0.05,      0.35,
        0,    1,    3,      0.05,      0.35,
        0,    1,    4,      0.05,      0.35,
        0,    2,    0,      0.00,      0.00,
        0,    2,    1,      0.00,      0.15,
        0,    2,    2,      0.00,      0.15,
        0,    2,    3,      0.00,      0.15,
        0,    2,    4,      0.00,      0.15,
        1,    1,    1,      0.30,      0.30,
        1,    1,    2,      0.30,      0.50,
        1,    1,    3,      0.30,      0.50,
        1,    1,    4,      0.30,      0.50,
        2,    1,    2,      0.50,      0.50,
        2,    1,    3,      0.50,      0.50,
        2,    1,    4,      0.50,      0.50,
        3,    1,    3,      0.50,      0.50,
        3,    1,    4,      0.50,      0.50,
        4,    1,    4,      0.50,      0.50,
        1,    2,    1,      0.30,      0.30,
        1,    2,    2,      0.25,      0.25,
        1,    2,    3,      0.25,      0.25,
        1,    2,    4,      0.25,      0.25,
        2,    2,    2,      0.25,      0.25,
        2,    2,    3,      0.25,      0.25,
        2,    2,    4,      0.25,      0.25,
        3,    2,    3,      0.25,      0.25,
        3,    2,    4,      0.25,      0.25,
        4,    2,    4,      0.25,      0.25  };

EXTERN int Missing_constants;
//EXTERN FILE *errfile;

void kstrbnd()
{
   int i, j, iz, nb1, nb2, irow, krow, lrow;
   int ia, ib, ic, itb, ita, itc, ierr;
   int k, nstbn, index;
   char pa[4],pb[4],pc[4],pt[10],pt1[10];

   strbnd.nstrbnd = 0;
   nb1 = -1;
   nb2 = -1;
   
   for (i=0; i < angles.nang; i++)
   {
       ierr = FALSE;
       ia = angles.i13[i][0];
       ib = angles.i13[i][1];
       ic = angles.i13[i][2];
       ita = atom[ia].type;
       itb = atom[ib].type;
       itc = atom[ic].type;
       if (field.type == MMX && (atom[ia].type >= 300 || atom[ib].type >= 300 || atom[ic].type >= 300))
          goto L_10;
       if (field.type == MMX && (atom[ia].type >= 20 || atom[ic].type >= 20))
          goto L_10;
       if (field.type == MM3 && (atom[ia].type >= 300 || atom[ib].type >= 300 || atom[ic].type >= 300) )
          goto L_10;  
       nstbn = -1;
       numeral(ita,pa,3);
       numeral(itb,pb,3);
       numeral(itc,pc,3);
       if (ita <= itc)
       {
           strcpy(pt,pa);     strcat(pt,pb);  strcat(pt,pc);
       }else
       {
           strcpy(pt,pc);     strcat(pt,pb);  strcat(pt,pa);
       }
       itb = atom[ib].tclass;
       numeral(itb,pb,3);
       strcpy(pt1,"  0"); strcat(pt1,pb); strcat(pt1,"  0");
       
       for (j=0; j < crossterm_k.nstrbnd; j++)
       {
           if (strcmp(crossterm_k.stbn[j],pt) == 0 || strcmp(crossterm_k.stbn[j],pt1) == 0 )
           {
               nstbn = j;
               ierr = TRUE;
               break;
           }
       }
       for (j=0; j < bonds_ff.nbnd; j++)
       {
           iz = find_bond(ia, ib);
           if (iz != -1)
           {
              nb1 = iz;
              break;
           }
       }
      for (j=0; j < bonds_ff.nbnd; j++)
       {
           iz = find_bond(ib, ic);
           if (iz != -1)
           {
              nb2 = iz;
              break;
           }
       }
       if (angles.index[i] == 0)
          index = 0;
       else if (angles.index[i] == 1)
       {
          if (bonds_ff.index[nb1] == 1)
             index = 1;
          else
             index = 2;
       } else if (angles.index[i] == 2)
          index = 3;
       else if (angles.index[i] == 3)
          index = 5;
       else if (angles.index[i] == 4)
          index = 4;
       else if (angles.index[i] == 5)
       {
          if (bonds_ff.index[nb1] == 1)
             index = 6;
          else
             index = 7;
       } else if (angles.index[i] == 6)
          index = 8;
       else if (angles.index[i] == 7)
       {
          if (bonds_ff.index[nb1] == 1)
             index = 9;
          else
             index = 10;
       } else if (angles.index[i] == 8)
          index = 11;
       
       if (crossterm_k.stbnindex[nstbn] != index)
          nstbn = -1;
          
       if (field.type == MMFF94)
       {
           ierr = TRUE;
           if (angles.anat[i] < 175.0)
           {
              if (nstbn == -1)  // constant not found  assign constants base on periodic table
              {
                 if (atom[ia].atomnum <= 1)
                  irow = 0;
                 else if (atom[ia].atomnum <= 10)
                  irow = 1;
                 else if (atom[ia].atomnum <= 18)
                  irow = 2;
                 else if (atom[ia].atomnum <= 36)
                  irow = 3;
                 else if (atom[ia].atomnum <= 54)
                  irow = 4;
                  
                 if (atom[ib].atomnum <= 1)
                  krow = 0;
                 else if (atom[ib].atomnum <= 10)
                  krow = 1;
                 else if (atom[ib].atomnum <= 18)
                  krow = 2;
                 else if (atom[ib].atomnum <= 36)
                  krow = 3;
                 else if (atom[ib].atomnum <= 54)
                  krow = 4;
                  
                 if (atom[ic].atomnum <= 1)
                  lrow = 0;
                 else if (atom[ic].atomnum <= 10)
                  lrow = 1;
                 else if (atom[ic].atomnum <= 18)
                  lrow = 2;
                 else if (atom[ic].atomnum <= 36)
                  lrow = 3;
                 else if (atom[ic].atomnum <= 54)
                  lrow = 4;

                for (j=0; j < 30; j++)
                {
                  if ( irow == stbn[j].irow && krow == stbn[j].krow && lrow == stbn[j].lrow )
                  {
                      strbnd.isb[strbnd.nstrbnd][0] = i;
                      strbnd.isb[strbnd.nstrbnd][1] = nb1;
                      strbnd.isb[strbnd.nstrbnd][2] = nb2;
                      strbnd.ksb1[strbnd.nstrbnd] = stbn[j].con1;
                      strbnd.ksb2[strbnd.nstrbnd] = stbn[j].con2;
                      strbnd.nstrbnd++;
                      break;
                  } else if (lrow == stbn[j].irow && krow == stbn[j].krow && irow == stbn[j].lrow)
                  {
                      strbnd.isb[strbnd.nstrbnd][0] = i;
                      strbnd.isb[strbnd.nstrbnd][1] = nb1;
                      strbnd.isb[strbnd.nstrbnd][2] = nb2;
                      strbnd.ksb1[strbnd.nstrbnd] = stbn[j].con2;
                      strbnd.ksb2[strbnd.nstrbnd] = stbn[j].con1;
                      strbnd.nstrbnd++;
                      break;
                  }
                }   
             }else
             {
                strbnd.isb[strbnd.nstrbnd][0] = i;
                strbnd.isb[strbnd.nstrbnd][1] = nb1;
                strbnd.isb[strbnd.nstrbnd][2] = nb2;
                if (angles.anat[i] > 175.0)
                {
                   strbnd.ksb1[strbnd.nstrbnd] = 0.0;
                   strbnd.ksb2[strbnd.nstrbnd] = 0.0;
                } else
                { 
                 if (ita < itc)
                 {
                    strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][0];
                    strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][1];
                 } else
                 {
                    strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][1];
                    strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][0];
                 }
                }
                strbnd.nstrbnd++;
             }
           }
       } else
       {         
             k = 0;
             if ( atom[ia].atomnum <= 1)
             {
                 k++;
                 if (crossterm_k.stbncon[nstbn][2] == 0.0)
                   nb1 = -1;
             }
             if (atom[ic].atomnum <= 1)
             {
                 k++;
                 if (crossterm_k.stbncon[nstbn][2] == 0.0)
                    nb2 = -1;
             }
             if (crossterm_k.stbncon[nstbn][k] != 0.0)
             {
                 strbnd.isb[strbnd.nstrbnd][0] = i;
                 strbnd.isb[strbnd.nstrbnd][1] = nb1;
                 strbnd.isb[strbnd.nstrbnd][2] = nb2;
                 strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][k];           
                 strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][k];           
                 strbnd.nstrbnd++;
             }
       }
//       if (ierr == FALSE)
//       {
//           Missing_constants = TRUE;
//           fprintf(errfile,"Missing strbnd constants for %d %d %d \n",ia,ib,ic);
//       }
L_10:
    continue;
   }
}
           
       
Revision:	27
Committed:	Tue Jul 8 19:10:25 2008 UTC (13 years, 11 months ago) by tjod
Original Path:	*trunk/src/mengine/src/kstrbnd.c*
File size:	9351 byte(s)
Log Message:	move mengine src to reflect new freemol directory structure
Line	User	Rev	File contents
1	tjod	3	#define EXTERN extern
2			#include "pcwin.h"
3			#include "pcmod.h"
4
5			#include "angles.h"
6			#include "bonds_ff.h"
7			#include "field.h"
8
9			int find_bond(int, int);
10			void numeral(int, char *, int);
11
12			EXTERN struct t_crossterm_k {
13			int nangang, nstrbnd, nstrtor;
14			int ang_ang[MAXAA], stbnindex[MAXSTBN] ;
15			char str_tor[MAXSTRTOR][7], stbn[MAXSTBN][10];
16			float aacon[MAXAA][3], stbncon[MAXSTBN][3], str_torcon[MAXSTRTOR];
17			} crossterm_k;
18
19			struct t_strbnd {
20			int nstrbnd, isb[MAXANG][3];
21			float ksb1[MAXANG], ksb2[MAXANG];
22			} strbnd;
23			struct stbndat {
24			int irow, krow, lrow;
25			float con1, con2;
26			} stbn[] =
27			{
28			0, 1, 0, 0.15, 0.15,
29			0, 1, 1, 0.10, 0.30,
30			0, 1, 2, 0.05, 0.35,
31			0, 1, 3, 0.05, 0.35,
32			0, 1, 4, 0.05, 0.35,
33			0, 2, 0, 0.00, 0.00,
34			0, 2, 1, 0.00, 0.15,
35			0, 2, 2, 0.00, 0.15,
36			0, 2, 3, 0.00, 0.15,
37			0, 2, 4, 0.00, 0.15,
38			1, 1, 1, 0.30, 0.30,
39			1, 1, 2, 0.30, 0.50,
40			1, 1, 3, 0.30, 0.50,
41			1, 1, 4, 0.30, 0.50,
42			2, 1, 2, 0.50, 0.50,
43			2, 1, 3, 0.50, 0.50,
44			2, 1, 4, 0.50, 0.50,
45			3, 1, 3, 0.50, 0.50,
46			3, 1, 4, 0.50, 0.50,
47			4, 1, 4, 0.50, 0.50,
48			1, 2, 1, 0.30, 0.30,
49			1, 2, 2, 0.25, 0.25,
50			1, 2, 3, 0.25, 0.25,
51			1, 2, 4, 0.25, 0.25,
52			2, 2, 2, 0.25, 0.25,
53			2, 2, 3, 0.25, 0.25,
54			2, 2, 4, 0.25, 0.25,
55			3, 2, 3, 0.25, 0.25,
56			3, 2, 4, 0.25, 0.25,
57			4, 2, 4, 0.25, 0.25 };
58
59			EXTERN int Missing_constants;
60			//EXTERN FILE *errfile;
61
62			void kstrbnd()
63			{
64			int i, j, iz, nb1, nb2, irow, krow, lrow;
65			int ia, ib, ic, itb, ita, itc, ierr;
66			int k, nstbn, index;
67			char pa[4],pb[4],pc[4],pt[10],pt1[10];
68
69			strbnd.nstrbnd = 0;
70			nb1 = -1;
71			nb2 = -1;
72
73			for (i=0; i < angles.nang; i++)
74			{
75			ierr = FALSE;
76			ia = angles.i13[i][0];
77			ib = angles.i13[i][1];
78			ic = angles.i13[i][2];
79			ita = atom[ia].type;
80			itb = atom[ib].type;
81			itc = atom[ic].type;
82			if (field.type == MMX && (atom[ia].type >= 300 \|\| atom[ib].type >= 300 \|\| atom[ic].type >= 300))
83			goto L_10;
84			if (field.type == MMX && (atom[ia].type >= 20 \|\| atom[ic].type >= 20))
85			goto L_10;
86			if (field.type == MM3 && (atom[ia].type >= 300 \|\| atom[ib].type >= 300 \|\| atom[ic].type >= 300) )
87			goto L_10;
88			nstbn = -1;
89			numeral(ita,pa,3);
90			numeral(itb,pb,3);
91			numeral(itc,pc,3);
92			if (ita <= itc)
93			{
94			strcpy(pt,pa); strcat(pt,pb); strcat(pt,pc);
95			}else
96			{
97			strcpy(pt,pc); strcat(pt,pb); strcat(pt,pa);
98			}
99			itb = atom[ib].tclass;
100			numeral(itb,pb,3);
101			strcpy(pt1," 0"); strcat(pt1,pb); strcat(pt1," 0");
102
103			for (j=0; j < crossterm_k.nstrbnd; j++)
104			{
105			if (strcmp(crossterm_k.stbn[j],pt) == 0 \|\| strcmp(crossterm_k.stbn[j],pt1) == 0 )
106			{
107			nstbn = j;
108			ierr = TRUE;
109			break;
110			}
111			}
112			for (j=0; j < bonds_ff.nbnd; j++)
113			{
114			iz = find_bond(ia, ib);
115			if (iz != -1)
116			{
117			nb1 = iz;
118			break;
119			}
120			}
121			for (j=0; j < bonds_ff.nbnd; j++)
122			{
123			iz = find_bond(ib, ic);
124			if (iz != -1)
125			{
126			nb2 = iz;
127			break;
128			}
129			}
130			if (angles.index[i] == 0)
131			index = 0;
132			else if (angles.index[i] == 1)
133			{
134			if (bonds_ff.index[nb1] == 1)
135			index = 1;
136			else
137			index = 2;
138			} else if (angles.index[i] == 2)
139			index = 3;
140			else if (angles.index[i] == 3)
141			index = 5;
142			else if (angles.index[i] == 4)
143			index = 4;
144			else if (angles.index[i] == 5)
145			{
146			if (bonds_ff.index[nb1] == 1)
147			index = 6;
148			else
149			index = 7;
150			} else if (angles.index[i] == 6)
151			index = 8;
152			else if (angles.index[i] == 7)
153			{
154			if (bonds_ff.index[nb1] == 1)
155			index = 9;
156			else
157			index = 10;
158			} else if (angles.index[i] == 8)
159			index = 11;
160
161			if (crossterm_k.stbnindex[nstbn] != index)
162			nstbn = -1;
163
164			if (field.type == MMFF94)
165			{
166			ierr = TRUE;
167			if (angles.anat[i] < 175.0)
168			{
169			if (nstbn == -1) // constant not found assign constants base on periodic table
170			{
171			if (atom[ia].atomnum <= 1)
172			irow = 0;
173			else if (atom[ia].atomnum <= 10)
174			irow = 1;
175			else if (atom[ia].atomnum <= 18)
176			irow = 2;
177			else if (atom[ia].atomnum <= 36)
178			irow = 3;
179			else if (atom[ia].atomnum <= 54)
180			irow = 4;
181
182			if (atom[ib].atomnum <= 1)
183			krow = 0;
184			else if (atom[ib].atomnum <= 10)
185			krow = 1;
186			else if (atom[ib].atomnum <= 18)
187			krow = 2;
188			else if (atom[ib].atomnum <= 36)
189			krow = 3;
190			else if (atom[ib].atomnum <= 54)
191			krow = 4;
192
193			if (atom[ic].atomnum <= 1)
194			lrow = 0;
195			else if (atom[ic].atomnum <= 10)
196			lrow = 1;
197			else if (atom[ic].atomnum <= 18)
198			lrow = 2;
199			else if (atom[ic].atomnum <= 36)
200			lrow = 3;
201			else if (atom[ic].atomnum <= 54)
202			lrow = 4;
203
204			for (j=0; j < 30; j++)
205			{
206			if ( irow == stbn[j].irow && krow == stbn[j].krow && lrow == stbn[j].lrow )
207			{
208			strbnd.isb[strbnd.nstrbnd][0] = i;
209			strbnd.isb[strbnd.nstrbnd][1] = nb1;
210			strbnd.isb[strbnd.nstrbnd][2] = nb2;
211			strbnd.ksb1[strbnd.nstrbnd] = stbn[j].con1;
212			strbnd.ksb2[strbnd.nstrbnd] = stbn[j].con2;
213			strbnd.nstrbnd++;
214			break;
215			} else if (lrow == stbn[j].irow && krow == stbn[j].krow && irow == stbn[j].lrow)
216			{
217			strbnd.isb[strbnd.nstrbnd][0] = i;
218			strbnd.isb[strbnd.nstrbnd][1] = nb1;
219			strbnd.isb[strbnd.nstrbnd][2] = nb2;
220			strbnd.ksb1[strbnd.nstrbnd] = stbn[j].con2;
221			strbnd.ksb2[strbnd.nstrbnd] = stbn[j].con1;
222			strbnd.nstrbnd++;
223			break;
224			}
225			}
226			}else
227			{
228			strbnd.isb[strbnd.nstrbnd][0] = i;
229			strbnd.isb[strbnd.nstrbnd][1] = nb1;
230			strbnd.isb[strbnd.nstrbnd][2] = nb2;
231			if (angles.anat[i] > 175.0)
232			{
233			strbnd.ksb1[strbnd.nstrbnd] = 0.0;
234			strbnd.ksb2[strbnd.nstrbnd] = 0.0;
235			} else
236			{
237			if (ita < itc)
238			{
239			strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][0];
240			strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][1];
241			} else
242			{
243			strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][1];
244			strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][0];
245			}
246			}
247			strbnd.nstrbnd++;
248			}
249			}
250			} else
251			{
252			k = 0;
253			if ( atom[ia].atomnum <= 1)
254			{
255			k++;
256			if (crossterm_k.stbncon[nstbn][2] == 0.0)
257			nb1 = -1;
258			}
259			if (atom[ic].atomnum <= 1)
260			{
261			k++;
262			if (crossterm_k.stbncon[nstbn][2] == 0.0)
263			nb2 = -1;
264			}
265			if (crossterm_k.stbncon[nstbn][k] != 0.0)
266			{
267			strbnd.isb[strbnd.nstrbnd][0] = i;
268			strbnd.isb[strbnd.nstrbnd][1] = nb1;
269			strbnd.isb[strbnd.nstrbnd][2] = nb2;
270			strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][k];
271			strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][k];
272			strbnd.nstrbnd++;
273			}
274			}
275			// if (ierr == FALSE)
276			// {
277			// Missing_constants = TRUE;
278			// fprintf(errfile,"Missing strbnd constants for %d %d %d \n",ia,ib,ic);
279			// }
280			L_10:
281			continue;
282			}
283			}
284
285