NAMD
NamdOneTools.C
Go to the documentation of this file.
1 
7 /*
8  NAMD 1.X functions copied without change for use primarily in WorkDistrib
9 */
10 
11 #include "InfoStream.h"
12 #include "common.h"
13 #include "NamdTypes.h"
14 #include "NamdOneTools.h"
15 #include "Vector.h"
16 #include "PDB.h"
17 #include "Molecule.h"
18 #define MIN_DEBUG_LEVEL 4
19 //#define DEBUGM
20 #include "Debug.h"
21 
22 //#define DEBUGM
23 #include "Debug.h"
24 
25 /************************************************************************/
26 /* FUNCTION read_binary_coors */
27 /* INPUTS: */
28 /* fname - Filename to read coordinates from */
29 /* pdbobj - PDB object to place coordinates into */
30 /* This function reads initial coordinates from a binary */
31 /* restart file */
32 /************************************************************************/
33 
34 void read_binary_coors(char *fname, PDB *pdbobj) {
35  Vector *newcoords; // Array of vectors to hold coordinates from file
36 
37  // Allocate an array to hold the new coordinates
38  newcoords = new Vector[pdbobj->num_atoms()];
39 
40  // Read the coordinate from the file
41  read_binary_file(fname,newcoords,pdbobj->num_atoms());
42 
43  // Set the coordinates in the PDB object to the new coordinates
44  pdbobj->set_all_positions(newcoords);
45 
46  // Clean up
47  delete [] newcoords;
48 
49 } // END OF FUNCTION read_binary_coors()
50 
51 
52 void read_binary_file(const char *fname, Vector *data, int n)
53 {
54  int32 filen; // Number of atoms read from file
55  FILE *fp; // File descriptor
56  int needToFlip = 0;
57 
58  iout << iINFO << "Reading from binary file " << fname << "\n" << endi;
59 
60  // Open the file and die if the open fails
61  if ( (fp = Fopen(fname, "rb")) == NULL)
62  {
63  char errmsg[256];
64  sprintf(errmsg, "Unable to open binary file %s", fname);
65  NAMD_die(errmsg);
66  }
67 
68  // read the number of coordinates in this file
69  if (fread(&filen, sizeof(int32), 1, fp) != (size_t)1)
70  {
71  char errmsg[256];
72  sprintf(errmsg, "Error reading binary file %s", fname);
73  NAMD_die(errmsg);
74  }
75 
76  // read the number of coordinates in this file
77  // check for palindromic number of atoms
78  char lenbuf[4];
79  memcpy(lenbuf, (const char *)&filen, 4);
80  char tmpc;
81  tmpc = lenbuf[0]; lenbuf[0] = lenbuf[3]; lenbuf[3] = tmpc;
82  tmpc = lenbuf[1]; lenbuf[1] = lenbuf[2]; lenbuf[2] = tmpc;
83  if ( ! memcmp((const char *)&filen, lenbuf, 4) ) {
84  iout << iWARN << "Number of atoms in binary file " << fname <<
85  " is palindromic, assuming same endian.\n" << endi;
86  }
87 
88  // Die if this doesn't match the number in our system
89  if (filen != n)
90  {
91  needToFlip = 1;
92  memcpy((char *)&filen, lenbuf, 4);
93  }
94  if (filen != n)
95  {
96  char errmsg[256];
97  sprintf(errmsg, "Incorrect atom count in binary file %s", fname);
98  NAMD_die(errmsg);
99  }
100 
101  if (fread(data, sizeof(Vector), n, fp) != (size_t)n)
102  {
103  char errmsg[256];
104  sprintf(errmsg, "Error reading binary file %s", fname);
105  NAMD_die(errmsg);
106  }
107 
108  Fclose(fp);
109 
110  if (needToFlip) {
111  iout << iWARN << "Converting binary file " << fname << "\n" << endi;
112  int i;
113  char *cdata = (char *) data;
114  for ( i=0; i<3*n; ++i, cdata+=8 ) {
115  char tmp0, tmp1, tmp2, tmp3;
116  tmp0 = cdata[0]; tmp1 = cdata[1];
117  tmp2 = cdata[2]; tmp3 = cdata[3];
118  cdata[0] = cdata[7]; cdata[1] = cdata[6];
119  cdata[2] = cdata[5]; cdata[3] = cdata[4];
120  cdata[7] = tmp0; cdata[6] = tmp1;
121  cdata[5] = tmp2; cdata[4] = tmp3;
122  }
123  }
124 
125 }
126 
127 
128 /*
129  * Generate a 3x3 transformation matrix for rotation about a given
130  * vector v by an angle given in degrees
131  */
133  /* This function contributed by Erich Boleyn (erich@uruk.org) */
134  BigReal mag, s, c;
135  BigReal xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c;
136 
137  s = sin(angle * PI/180.0);
138  c = cos(angle * PI/180.0);
139 
140  mag = v.length();
141 
142  if (mag == 0.0) {
143  /* generate an identity matrix and return */
144  for ( int i = 0; i < 9; ++i ) m[i] = 0.0;
145  m[0] = m[4] = m[8] = 1.0;
146  return;
147  }
148 
149  // normalize the vector
150  v /= mag;
151 
152  /*
153  * Arbitrary axis rotation matrix.
154  *
155  * This is composed of 5 matrices, Rz, Ry, T, Ry', Rz', multiplied
156  * like so: Rz * Ry * T * Ry' * Rz'. T is the final rotation
157  * (which is about the X-axis), and the two composite transforms
158  * Ry' * Rz' and Rz * Ry are (respectively) the rotations necessary
159  * from the arbitrary axis to the X-axis then back. They are
160  * all elementary rotations.
161  *
162  * Rz' is a rotation about the Z-axis, to bring the axis vector
163  * into the x-z plane. Then Ry' is applied, rotating about the
164  * Y-axis to bring the axis vector parallel with the X-axis. The
165  * rotation about the X-axis is then performed. Ry and Rz are
166  * simply the respective inverse transforms to bring the arbitrary
167  * axis back to it's original orientation. The first transforms
168  * Rz' and Ry' are considered inverses, since the data from the
169  * arbitrary axis gives you info on how to get to it, not how
170  * to get away from it, and an inverse must be applied.
171  *
172  * The basic calculation used is to recognize that the arbitrary
173  * axis vector (x, y, z), since it is of unit length, actually
174  * represents the sines and cosines of the angles to rotate the
175  * X-axis to the same orientation, with theta being the angle about
176  * Z and phi the angle about Y (in the order described above)
177  * as follows:
178  *
179  * cos ( theta ) = x / sqrt ( 1 - z^2 )
180  * sin ( theta ) = y / sqrt ( 1 - z^2 )
181  *
182  * cos ( phi ) = sqrt ( 1 - z^2 )
183  * sin ( phi ) = z
184  *
185  * Note that cos ( phi ) can further be inserted to the above
186  * formulas:
187  *
188  * cos ( theta ) = x / cos ( phi )
189  * sin ( theta ) = y / sin ( phi )
190  *
191  * ...etc. Because of those relations and the standard trigonometric
192  * relations, it is pssible to reduce the transforms down to what
193  * is used below. It may be that any primary axis chosen will give the
194  * same results (modulo a sign convention) using thie method.
195  *
196  * Particularly nice is to notice that all divisions that might
197  * have caused trouble when parallel to certain planes or
198  * axis go away with care paid to reducing the expressions.
199  * After checking, it does perform correctly under all cases, since
200  * in all the cases of division where the denominator would have
201  * been zero, the numerator would have been zero as well, giving
202  * the expected result.
203  */
204 
205  // store matrix in (row, col) form, i.e., m(row,col) = m[row*3+col]
206 
207  xx = v.x * v.x;
208  yy = v.y * v.y;
209  zz = v.z * v.z;
210  xy = v.x * v.y;
211  yz = v.y * v.z;
212  zx = v.z * v.x;
213  xs = v.x * s;
214  ys = v.y * s;
215  zs = v.z * s;
216  one_c = 1.0 - c;
217 
218  m[0] = (one_c * xx) + c;
219  m[1] = (one_c * xy) - zs;
220  m[2] = (one_c * zx) + ys;
221 
222  m[3] = (one_c * xy) + zs;
223  m[4] = (one_c * yy) + c;
224  m[5] = (one_c * yz) - xs;
225 
226  m[6] = (one_c * zx) - ys;
227  m[7] = (one_c * yz) + xs;
228  m[8] = (one_c * zz) + c;
229 }
230 
231 // multiply vector v by a 3x3 matrix m stored so that m(row,col) = m[row*3+col]
233  return Vector( m[0]*v.x + m[1]*v.y + m[2]*v.z,
234  m[3]*v.x + m[4]*v.y + m[5]*v.z,
235  m[6]*v.x + m[7]*v.y + m[8]*v.z );
236 }
237 
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
virial xy
Definition: PDB.h:35
short int32
Definition: dumpdcd.c:24
Definition: Vector.h:64
BigReal z
Definition: Vector.h:66
void vec_rotation_matrix(BigReal angle, Vector v, BigReal m[])
Definition: NamdOneTools.C:132
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
void set_all_positions(Vector *)
Definition: PDB.C:331
#define iout
Definition: InfoStream.h:87
int num_atoms(void)
Definition: PDB.C:323
BigReal length(void) const
Definition: Vector.h:169
virial zz
#define PI
Definition: common.h:81
void read_binary_file(const char *fname, Vector *data, int n)
Definition: NamdOneTools.C:52
void read_binary_coors(char *fname, PDB *pdbobj)
Definition: NamdOneTools.C:34
Vector mat_multiply_vec(const Vector &v, BigReal m[])
Definition: NamdOneTools.C:232
FILE * Fopen(const char *filename, const char *mode)
Definition: common.C:265
virial xx
BigReal x
Definition: Vector.h:66
void NAMD_die(const char *err_msg)
Definition: common.C:83
virial yz
virial yy
int Fclose(FILE *fout)
Definition: common.C:359
virial zx
BigReal y
Definition: Vector.h:66
infostream & endi(infostream &s)
Definition: InfoStream.C:38
double BigReal
Definition: common.h:112