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:	64
Committed:	Thu Dec 4 02:03:23 2008 UTC (13 years, 6 months ago) by gilbertke
Original Path:	*trunk/src/mengine/src/kstrbnd.c*
File size:	9335 byte(s)
Log Message:	fixed bug in routine to recognize delocalized bonds in kbond.c
Line	File contents
1	#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	for (j=0; j < crossterm_k.nstrbnd; j++)
103	{
104	if (strcmp(crossterm_k.stbn[j],pt) == 0 \|\| strcmp(crossterm_k.stbn[j],pt1) == 0 )
105	{
106	nstbn = j;
107	ierr = TRUE;
108	break;
109	}
110	}
111	for (j=0; j < bonds_ff.nbnd; j++)
112	{
113	iz = find_bond(ia, ib);
114	if (iz != -1)
115	{
116	nb1 = iz;
117	break;
118	}
119	}
120	for (j=0; j < bonds_ff.nbnd; j++)
121	{
122	iz = find_bond(ib, ic);
123	if (iz != -1)
124	{
125	nb2 = iz;
126	break;
127	}
128	}
129	if (angles.index[i] == 0)
130	index = 0;
131	else if (angles.index[i] == 1)
132	{
133	if (bonds_ff.index[nb1] == 1)
134	index = 1;
135	else
136	index = 2;
137	} else if (angles.index[i] == 2)
138	index = 3;
139	else if (angles.index[i] == 3)
140	index = 5;
141	else if (angles.index[i] == 4)
142	index = 4;
143	else if (angles.index[i] == 5)
144	{
145	if (bonds_ff.index[nb1] == 1)
146	index = 6;
147	else
148	index = 7;
149	} else if (angles.index[i] == 6)
150	index = 8;
151	else if (angles.index[i] == 7)
152	{
153	if (bonds_ff.index[nb1] == 1)
154	index = 9;
155	else
156	index = 10;
157	} else if (angles.index[i] == 8)
158	index = 11;
159	if (crossterm_k.stbnindex[nstbn] != index)
160	nstbn = -1;
161
162	if (field.type == MMFF94)
163	{
164	ierr = TRUE;
165	if (angles.anat[i] < 175.0)
166	{
167	if (nstbn == -1) // constant not found assign constants base on periodic table
168	{
169	if (atom[ia].atomnum <= 1)
170	irow = 0;
171	else if (atom[ia].atomnum <= 10)
172	irow = 1;
173	else if (atom[ia].atomnum <= 18)
174	irow = 2;
175	else if (atom[ia].atomnum <= 36)
176	irow = 3;
177	else if (atom[ia].atomnum <= 54)
178	irow = 4;
179
180	if (atom[ib].atomnum <= 1)
181	krow = 0;
182	else if (atom[ib].atomnum <= 10)
183	krow = 1;
184	else if (atom[ib].atomnum <= 18)
185	krow = 2;
186	else if (atom[ib].atomnum <= 36)
187	krow = 3;
188	else if (atom[ib].atomnum <= 54)
189	krow = 4;
190
191	if (atom[ic].atomnum <= 1)
192	lrow = 0;
193	else if (atom[ic].atomnum <= 10)
194	lrow = 1;
195	else if (atom[ic].atomnum <= 18)
196	lrow = 2;
197	else if (atom[ic].atomnum <= 36)
198	lrow = 3;
199	else if (atom[ic].atomnum <= 54)
200	lrow = 4;
201
202	for (j=0; j < 30; j++)
203	{
204	if ( irow == stbn[j].irow && krow == stbn[j].krow && lrow == stbn[j].lrow )
205	{
206	strbnd.isb[strbnd.nstrbnd][0] = i;
207	strbnd.isb[strbnd.nstrbnd][1] = nb1;
208	strbnd.isb[strbnd.nstrbnd][2] = nb2;
209	strbnd.ksb1[strbnd.nstrbnd] = stbn[j].con1;
210	strbnd.ksb2[strbnd.nstrbnd] = stbn[j].con2;
211	strbnd.nstrbnd++;
212	break;
213	} else if (lrow == stbn[j].irow && krow == stbn[j].krow && irow == stbn[j].lrow)
214	{
215	strbnd.isb[strbnd.nstrbnd][0] = i;
216	strbnd.isb[strbnd.nstrbnd][1] = nb1;
217	strbnd.isb[strbnd.nstrbnd][2] = nb2;
218	strbnd.ksb1[strbnd.nstrbnd] = stbn[j].con2;
219	strbnd.ksb2[strbnd.nstrbnd] = stbn[j].con1;
220	strbnd.nstrbnd++;
221	break;
222	}
223	}
224	}else
225	{
226	strbnd.isb[strbnd.nstrbnd][0] = i;
227	strbnd.isb[strbnd.nstrbnd][1] = nb1;
228	strbnd.isb[strbnd.nstrbnd][2] = nb2;
229	if (angles.anat[i] > 175.0)
230	{
231	strbnd.ksb1[strbnd.nstrbnd] = 0.0;
232	strbnd.ksb2[strbnd.nstrbnd] = 0.0;
233	} else
234	{
235	if (ita < itc)
236	{
237	strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][0];
238	strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][1];
239	} else
240	{
241	strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][1];
242	strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][0];
243	}
244	}
245	strbnd.nstrbnd++;
246	}
247	}
248	} else
249	{
250	k = 0;
251	if ( atom[ia].atomnum <= 1)
252	{
253	k++;
254	if (crossterm_k.stbncon[nstbn][2] == 0.0)
255	nb1 = -1;
256	}
257	if (atom[ic].atomnum <= 1)
258	{
259	k++;
260	if (crossterm_k.stbncon[nstbn][2] == 0.0)
261	nb2 = -1;
262	}
263	if (crossterm_k.stbncon[nstbn][k] != 0.0)
264	{
265	strbnd.isb[strbnd.nstrbnd][0] = i;
266	strbnd.isb[strbnd.nstrbnd][1] = nb1;
267	strbnd.isb[strbnd.nstrbnd][2] = nb2;
268	strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][k];
269	strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][k];
270	strbnd.nstrbnd++;
271	}
272	}
273	// if (ierr == FALSE)
274	// {
275	// Missing_constants = TRUE;
276	// fprintf(errfile,"Missing strbnd constants for %d %d %d \n",ia,ib,ic);
277	// }
278	L_10:
279	continue;
280	}
281	}
282
283