ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/rotbond.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 11 months ago) by tjod
File size: 5963 byte(s)
Log Message:
test

Line User Rev File contents
1 tjod 3 #define EXTERN extern
2     #include "pcwin.h"
3     #include "pcmod.h"
4     #include "utility.h"
5     #include "substr.h"
6    
7     void matident(float (*)[4]);
8     void matrotsc(float (*)[4], int, float, float);
9     void matrotat(float (*)[4], int, float);
10     void mattrans(float (*)[4], float, float, float);
11     void matriply(float (*)[4], float (*)[4], float (*)[4]);
12     void matxform(float (*)[4], int);
13     void rotb(int, int, int, int, float);
14     void chain(int, int, int);
15     double dihdrl(int,int,int,int);
16     void RefreshScreen(void);
17     void message_alert(char *, char *);
18    
19    
20     void rotb(int ia1,int n1,int n2,int ia4, float degree)
21     {
22     int issnum;
23     long int mask;
24     int i;
25     float aa, cos1, cos2, denom, rot1[4][4], sin1,
26     sin2, trans1[4][4], x1, x2, xform1[4][4], xform2[4][4], y1, y2,
27     z1, z2, xtmp, ytmp, ztmp;
28    
29     /* subroutine to rotate about a bond. the atoms of the bond are n1 and n2,
30     * the degree to rotate is idegree, and issnum is the substructure number
31     * we get from chain to tell matxform which atoms to rotate
32     * we assume chain has been called to set up the atoms to rotate
33     * we use xaxis and xyplan to generate the rotation matrices to put
34     * n1 at the origin and n2 on the x axis
35     * */
36    
37     for( i = 0; i < MAXSS; i++ )
38     {
39     if(!substr.istract[i])
40     goto L_1;
41     }
42     fprintf( pcmoutfile, "Error in rotbond- no substructs left\n" );
43     goto L_2;
44    
45     L_1:
46     issnum = i;
47     mask = (1L << issnum);
48     chain( n1, n2, issnum );
49    
50     x1 = atom[n1].x;
51     y1 = atom[n1].y;
52     z1 = atom[n1].z;
53     x2 = atom[n2].x - x1;
54     y2 = atom[n2].y - y1;
55     z2 = atom[n2].z - z1;
56     denom = sqrt( x2*x2 + y2*y2 );
57     if( fabs( denom ) < .00001 )
58     {
59     cos1 = 1.0;
60     sin1 = 0.0;
61     }
62     else
63     {
64     sin1 = -y2/denom;
65     cos1 = x2/denom;
66     }
67     x2 = x2*cos1 - y2*sin1;
68     denom = sqrt( x2*x2 + z2*z2 );
69     if( fabs( denom ) < .00001 )
70     {
71     cos2 = 1.0;
72     sin2 = 0.0;
73     }
74     else
75     {
76     sin2 = z2/denom;
77     cos2 = x2/denom;
78     }
79     aa = degree - dihdrl( ia1, n1, n2, ia4 );
80    
81     i = 3;
82     xtmp = -x1;
83     ytmp = -y1;
84     ztmp = -z1;
85     mattrans( trans1, xtmp, ytmp, ztmp );
86     matrotsc( rot1, i, sin1, cos1 );
87     matriply( trans1, rot1, xform1 );
88    
89     i = 2;
90     matrotsc( rot1, i, sin2, cos2 );
91     matriply( xform1, rot1, xform2 );
92    
93     i = 1;
94     matrotat( rot1, i, aa );
95     matriply( xform2, rot1, xform1 );
96    
97     i = 2;
98     xtmp = -sin2;
99     matrotsc( rot1, i, xtmp, cos2 );
100     matriply( xform1, rot1, xform2 );
101    
102     i = 3;
103     xtmp = -sin1;
104     matrotsc( rot1, i, xtmp, cos1 );
105     matriply( xform2, rot1, xform1 );
106    
107     mattrans( trans1, x1, y1, z1 );
108     matriply( xform1, trans1, xform2 );
109    
110     matxform( xform2, issnum );
111    
112     for( i = 1; i <= natom; i++ )
113     {
114     atom[i].substr[0] &= ~mask ;
115     }
116     L_2:
117     return;
118     }
119     /* ------------------------------------------ */
120     void chain(int ia,int ib,int issnum)
121     {
122     int newmem;
123     int i, ij, j;
124     long int mask;
125    
126     /**** this routine starts at atom ia and fills
127     * iarray with the atoms attached directly
128     * or indirectly to ia but not attached either
129     * directly or indirectly to ib. iarray is
130     * filled such that a non-zero value for the
131     * i-th member means that atom i is either
132     * directly or indirectly attached to ia. *** */
133     /**** zero the output array *** */
134    
135     mask = 1L << issnum ;
136    
137     /**** find the atoms adjacent to ia other than atom ib *** */
138    
139     for( j = 0; j < MAXIAT; j++ )
140     {
141     if( atom[ia].iat[j])
142     {
143     if( atom[ia].iat[j] != ib )
144     {
145     atom[atom[ia].iat[j]].substr[0] = 0 ;
146     atom[atom[ia].iat[j]].substr[0] |= mask ;
147     }
148     }
149     }
150    
151     /**** repetitively loop through the output array
152     * looking for adjacent atoms until no new
153     * attachments are found *** */
154     while( TRUE )
155     {
156     newmem = FALSE;
157     for( i = 1; i <= natom; i++ )
158     {
159     if( atom[i].substr[0] & mask) // atom not previously marked
160     {
161     for( j = 0; j < MAXIAT; j++ )
162     {
163     if( atom[i].iat[j] != 0 && atom[i].iat[j] != ia && !(atom[atom[i].iat[j]].substr[0] & mask))
164     {
165     /**** a new atom has been found *** */
166     atom[atom[i].iat[j]].substr[0] = 0 ;
167     atom[atom[i].iat[j]].substr[0] |= mask ;
168     newmem = TRUE;
169     }
170     }
171     }
172     }
173     /**** check for newly added atoms *** */
174     if( !newmem )
175     break;
176     }
177     if( ib )
178     {
179     if( atom[ib].substr[0] & mask)
180     {
181     atom[ib].substr[0] &= ~mask;
182     for( ij = 0; ij < MAXIAT; ij++ )
183     {
184     if( atom[ib].iat[ij] )
185     {
186     // if( atom[ib].bo[ij] != 9 )
187     {
188     if( atom[ib].iat[ij] != ia )
189     atom[atom[ib].iat[ij]].substr[0] &= ~mask;
190     }
191     }
192     }
193     }
194     }
195     return;
196     }