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, 4 months ago) by tjod
File size: 5963 byte(s)
Log Message:
test

Line File contents
1 #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 }