3
\brief function which places one fragment into the coordinate system of another
5
int natoms_A, natoms_B - number of atoms in each fragment
6
int P_A, P_B - number of ref pts to worry about, less than 3 for linear fragments
7
double **geom_A, double **geom_B - geometry of fragments A and B - geom B is changed
8
double **ref_coeff_A, double **ref_coeff_B - linear combinations which specify reference atoms
9
double R_AB - distance between reference atoms #1 on each fragment
10
theta_A, theta_B, tau, chi-A, chi-B - interfragment angles
14
#include <libciomr/libciomr.h>
17
#include <physconst.h>
22
double dot_prod(double *v1, double *v2);
23
void cross_prod(double *v1, double *v2, double *out);
24
void unit_vec(double *B, double *A, double *AB);
26
void orient_fragment(int natom_A, int natom_B, int P_A, int P_B, double **geom_A, double **geom_B,
27
double **ref_coeff_A, double **ref_coeff_B, double R_AB, double theta_A, double theta_B,
28
double tau, double phi_A, double phi_B, FILE *outfile)
30
int i, j, errcod, pts, xyz;
31
double tval, norm, B_angle, R_B1B2, R_B2B3, e12[3], e12b[3], e12c[3], e12d[3], erot[3];
32
double **ref_A, **ref_B, **ref_B_final;
33
double sign, cross1[3], cross2[3], cross3[3], phi2, phi3;
35
ref_A = block_matrix(3,3);
36
ref_B = block_matrix(3,3);
37
ref_B_final = block_matrix(3,3);
39
/* stick SOMETHING in for non-specified reference atoms - necessary to prevent complaints
40
about collinear reference atoms and to make zmat_point() work in such cases */
42
for (xyz=0; xyz<3; ++xyz)
43
ref_A[2][xyz] = (xyz+1)/_pi;
46
for (xyz=0; xyz<3; ++xyz)
47
ref_A[1][xyz] = (xyz+1)/(2*_pi);
50
for (pts=0; pts<P_A; ++pts)
51
for (xyz=0; xyz<3; ++xyz)
52
for (i=0; i<natom_A; ++i)
53
ref_A[pts][xyz] += ref_coeff_A[pts][i] * geom_A[i][xyz];
55
for (pts=0; pts<P_B; ++pts)
56
for (xyz=0; xyz<3; ++xyz)
57
for (i=0; i<natom_B;++i)
58
ref_B[pts][xyz] += ref_coeff_B[pts][i] * geom_B[i][xyz];
60
fprintf(outfile,"\tCoordinates for reference points on fragment A\n");
61
print_mat(ref_A,P_A,3,outfile);
62
fprintf(outfile,"\tCoordinates for reference points on fragment B (initial) \n");
63
print_mat(ref_B,P_B,3,outfile);
64
fprintf(outfile,"\tInterfragment coordinates:\n");
65
fprintf(outfile,"\t(1/)R_AB:%10.5f, theta_A:%10.5f, theta_B:%10.5f\n", R_AB, theta_A, theta_B);
66
fprintf(outfile,"\t tau:%10.5f, phi_A:%10.5f, phi_B:%10.5f\n", tau, phi_A, phi_B);
68
/* compute B1-B2 distance, B2-B3 distance, and B1-B2-B3 angle */
71
for (xyz=0; xyz<3; ++xyz)
72
R_B1B2 += (ref_B[1][xyz]-ref_B[0][xyz])*(ref_B[1][xyz]-ref_B[0][xyz]);
73
R_B1B2 = sqrt(R_B1B2);
78
for (xyz=0; xyz<3; ++xyz)
79
R_B2B3 += (ref_B[2][xyz]-ref_B[1][xyz])*(ref_B[2][xyz]-ref_B[1][xyz]);
80
R_B2B3 = sqrt(R_B2B3);
81
unit_vec(ref_B[1],ref_B[0],e12);
82
unit_vec(ref_B[1],ref_B[2],e12b);
83
B_angle = acos(dot_prod(e12,e12b))*180.0/_pi;
86
/* determine location of reference pts for B in coordinate system of A */
87
zmat_point(ref_A[2], ref_A[1], ref_A[0], R_AB, theta_A, phi_A, ref_B_final[0]);
89
zmat_point(ref_A[1], ref_A[0], ref_B_final[0], R_B1B2, theta_B, tau, ref_B_final[1]);
91
zmat_point(ref_A[0], ref_B_final[0], ref_B_final[1], R_B2B3, B_angle, phi_B, ref_B_final[2]);
93
fprintf(outfile,"Coordinates for references points for fragment B in coordinate system of A\n");
94
print_mat(ref_B_final,P_B,3,outfile);
96
/* translate geom_B to place B1 in correct location */
97
for (xyz=0; xyz<3; ++xyz) {
98
tval = ref_B_final[0][xyz] - ref_B[0][xyz];
99
for (i=0; i<natom_B; ++i)
100
geom_B[i][xyz] += tval;
103
for (pts=0; pts<P_B; ++pts)
104
for (xyz=0; xyz<3; ++xyz) {
105
ref_B[pts][xyz] = 0.0;
106
for (i=0; i<natom_B;++i)
107
ref_B[pts][xyz] += ref_coeff_B[pts][i] * geom_B[i][xyz];
109
//fprintf(outfile,"Reference points after translation (to fix point B1):\n");
110
//print_mat(ref_B,P_B,3,outfile);
112
if (P_B>1) { /* move fragment B to place reference point B2 in correct location */
113
/* Determine rotational angle and axis */
114
unit_vec(ref_B[1], ref_B[0], e12); /* v B1->B2 */
115
unit_vec(ref_B_final[1], ref_B[0], e12b); /* v B1->B2_final */
116
B_angle = acos(dot_prod(e12b,e12));
117
fprintf(outfile,"Rotation by %f degrees (to fix point B2)\n", 180.0*B_angle/_pi);
118
if (fabs(B_angle) > 1.0e-7) {
119
cross_prod(e12,e12b,erot);
121
/* Move B to put B1 at origin */
122
for (xyz=0; xyz<3; ++xyz)
123
for (i=0; i<natom_B;++i)
124
geom_B[i][xyz] -= ref_B[0][xyz];
127
rotate_vecs(erot, B_angle, geom_B, natom_B);
129
/* Move B back to coordinate system of A */
130
for (xyz=0; xyz<3; ++xyz)
131
for (i=0; i<natom_B;++i)
132
geom_B[i][xyz] += ref_B[0][xyz];
134
for (pts=0; pts<P_B; ++pts)
135
for (xyz=0; xyz<3; ++xyz) {
136
ref_B[pts][xyz] = 0.0;
137
for (i=0; i<natom_B;++i)
138
ref_B[pts][xyz] += ref_coeff_B[pts][i] * geom_B[i][xyz];
140
//fprintf(outfile,"Reference points after rotation (to fix point B2) \n");
141
//print_mat(ref_B,P_B,3,outfile);
145
if (P_B==3) { /* move fragment B to place reference point B3 in correct location */
146
/* Determine rotational angle and axis */
147
unit_vec(ref_B[1], ref_B[0], erot); /* B1 -> B2 is rotation axis */
149
/* Calculate B3-B1-B2-B3' torsion angle */
150
unit_vec(ref_B[2], ref_B[0], e12); /* v B1->B3 */
151
unit_vec(ref_B[1], ref_B[0], e12b); /* v B1->B2 */
152
phi2 = acos(dot_prod(e12,e12b));
153
unit_vec(ref_B[0], ref_B[1], e12c); /* v B2->B1 */
154
unit_vec(ref_B_final[2], ref_B[1], e12d); /* v B2->B3' */
155
phi3 = acos(dot_prod(e12c,e12d));
157
cross_prod(e12 , e12b, cross1) ; /* B3->B1 x B1->B2 */
158
cross_prod(e12c, e12d, cross2) ; /* B1->B2 x B2->B3 */
159
tval = dot_prod(cross1, cross2) ;
161
if ((sin(phi2) > 0.00001) && (sin(phi3) > 0.00001)) {
167
if (tval > 0.99999) B_angle = 0.0000;
168
else if (tval < -0.99999) B_angle = _pi;
169
else B_angle = acos(tval) ;
171
sign = 1.0; /* check sign */
172
cross_prod(cross1, cross2, cross3);
173
norm = sqrt(dot_prod(cross3, cross3));
174
if (fabs(norm) > 0.00001) {
175
for (xyz=0; xyz<3; ++xyz)
176
cross3[xyz] *= 1.0/norm;
177
tval = dot_prod(cross3, e12b);
178
if (tval < 0.0) sign = -1.0;
182
if (fabs(B_angle) > 1.0e-7) {
183
fprintf(outfile,"Rotation by %f degrees (to fix point B3)\n", 180.0*B_angle/_pi);
185
/* Move B to put B2 at origin */
186
for (xyz=0; xyz<3; ++xyz)
187
for (i=0; i<natom_B;++i)
188
geom_B[i][xyz] -= ref_B[1][xyz];
190
rotate_vecs(erot, B_angle, geom_B, natom_B);
192
/* Translate B1 back to coordinate system of A */
193
for (xyz=0; xyz<3; ++xyz)
194
for (i=0; i<natom_B;++i)
195
geom_B[i][xyz] += ref_B[1][xyz];
197
for (pts=0; pts<P_B; ++pts)
198
for (xyz=0; xyz<3; ++xyz) {
199
ref_B[pts][xyz] = 0.0;
200
for (i=0; i<natom_B;++i)
201
ref_B[pts][xyz] += ref_coeff_B[pts][i] * geom_B[i][xyz];
203
//fprintf(outfile,"Reference points on B after rotation for B2 \n");
204
//print_mat(ref_B,P_B,3,outfile);
208
/* check to see if desired reference points were obtained */
210
for (i=0; i<P_B; ++i)
211
for (xyz=0; xyz<3; ++xyz)
212
tval += fabs(ref_B[i][xyz] - ref_B_final[i][xyz]);
214
fprintf(outfile, "Unable to construct multi-fragment geometry.");
215
exit(PSI_RETURN_FAILURE);
218
fprintf(outfile,"Successfully constructed multifragment geometry.\n");
222
free_block(ref_B_final);
226
double dot_prod(double *v1, double *v2) {
227
return v1[0]*v2[0]+ v1[1]*v2[1]+ v1[2]*v2[2];
229
void cross_prod(double *v1, double *v2, double *out) {
230
out[0] = v1[1]*v2[2]-v1[2]*v2[1];
231
out[1] = -v1[0]*v2[2]+v1[2]*v2[0];
232
out[2] = v1[0]*v2[1]-v1[1]*v2[0];
235
void unit_vec(double *B, double *A, double *AB) {
240
norm += (A[i]-B[i])*(A[i]-B[i]);
243
AB[i] = (B[i] - A[i]) / norm;
added
removed
src/lib/libqt/orient_fragment.cc
