NamdOneTools.C

Go to the documentation of this file.
00001 
00007 /*
00008    NAMD 1.X functions copied without change for use primarily in WorkDistrib
00009 */
00010 
00011 #include "InfoStream.h"
00012 #include "common.h"
00013 #include "NamdTypes.h"
00014 #include "NamdOneTools.h"
00015 #include "Vector.h"
00016 #include "PDB.h"
00017 #include "Molecule.h"
00018 #define MIN_DEBUG_LEVEL 4
00019 //#define DEBUGM
00020 #include "Debug.h"
00021 
00022 //#define DEBUGM
00023 #include "Debug.h"
00024 
00025 /************************************************************************/
00026 /*   FUNCTION read_binary_coors                                         */
00027 /*   INPUTS:                                                            */
00028 /*      fname - Filename to read coordinates from                       */
00029 /*      pdbobj - PDB object to place coordinates into                   */
00030 /*      This function reads initial coordinates from a binary           */
00031 /*      restart file                                                    */
00032 /************************************************************************/
00033 
00034 void read_binary_coors(char *fname, PDB *pdbobj) {
00035   Vector *newcoords;    //  Array of vectors to hold coordinates from file
00036 
00037   //  Allocate an array to hold the new coordinates
00038   newcoords = new Vector[pdbobj->num_atoms()];
00039 
00040   //  Read the coordinate from the file
00041   read_binary_file(fname,newcoords,pdbobj->num_atoms());
00042 
00043   //  Set the coordinates in the PDB object to the new coordinates
00044   pdbobj->set_all_positions(newcoords);
00045 
00046   //  Clean up
00047   delete [] newcoords;
00048 
00049 } // END OF FUNCTION read_binary_coors()
00050 
00051 
00052 void read_binary_file(const char *fname, Vector *data, int n)
00053 {
00054   int32 filen;          //  Number of atoms read from file
00055   FILE *fp;             //  File descriptor
00056   int needToFlip = 0;
00057 
00058   iout << iINFO << "Reading from binary file " << fname << "\n" << endi;
00059 
00060   //  Open the file and die if the open fails
00061   if ( (fp = Fopen(fname, "rb")) == NULL)
00062   {
00063     char errmsg[256];
00064     sprintf(errmsg, "Unable to open binary file %s", fname);
00065     NAMD_die(errmsg);
00066   }
00067 
00068   //  read the number of coordinates in this file
00069   if (fread(&filen, sizeof(int32), 1, fp) != (size_t)1)
00070   {
00071     char errmsg[256];
00072     sprintf(errmsg, "Error reading binary file %s", fname);
00073     NAMD_die(errmsg);
00074   }
00075 
00076   //  read the number of coordinates in this file
00077   //  check for palindromic number of atoms
00078   char lenbuf[4];
00079   memcpy(lenbuf, (const char *)&filen, 4);
00080   char tmpc;
00081   tmpc = lenbuf[0]; lenbuf[0] = lenbuf[3]; lenbuf[3] = tmpc;
00082   tmpc = lenbuf[1]; lenbuf[1] = lenbuf[2]; lenbuf[2] = tmpc;
00083   if ( ! memcmp((const char *)&filen, lenbuf, 4) ) {
00084     iout << iWARN << "Number of atoms in binary file " << fname <<
00085                 " is palindromic, assuming same endian.\n" << endi;
00086   }
00087 
00088   //  Die if this doesn't match the number in our system
00089   if (filen != n)
00090   {
00091     needToFlip = 1;
00092     memcpy((char *)&filen, lenbuf, 4);
00093   }
00094   if (filen != n)
00095   {
00096     char errmsg[256];
00097     sprintf(errmsg, "Incorrect atom count in binary file %s", fname);
00098     NAMD_die(errmsg);
00099   }
00100 
00101   if (fread(data, sizeof(Vector), n, fp) != (size_t)n)
00102   {
00103     char errmsg[256];
00104     sprintf(errmsg, "Error reading binary file %s", fname);
00105     NAMD_die(errmsg);
00106   }
00107 
00108   Fclose(fp);
00109 
00110   if (needToFlip) { 
00111     iout << iWARN << "Converting binary file " << fname << "\n" << endi;
00112     int i;
00113     char *cdata = (char *) data;
00114     for ( i=0; i<3*n; ++i, cdata+=8 ) {
00115       char tmp0, tmp1, tmp2, tmp3;
00116       tmp0 = cdata[0]; tmp1 = cdata[1];
00117       tmp2 = cdata[2]; tmp3 = cdata[3];
00118       cdata[0] = cdata[7]; cdata[1] = cdata[6];
00119       cdata[2] = cdata[5]; cdata[3] = cdata[4];
00120       cdata[7] = tmp0; cdata[6] = tmp1;
00121       cdata[5] = tmp2; cdata[4] = tmp3;
00122     }
00123   }
00124 
00125 }
00126 
00127 
00128 /* 
00129  * Generate a 3x3 transformation matrix for rotation about a given
00130  * vector v by an angle given in degrees
00131  */
00132 void vec_rotation_matrix( BigReal angle, Vector v, BigReal m[] ) {
00133    /* This function contributed by Erich Boleyn (erich@uruk.org) */
00134    BigReal mag, s, c;
00135    BigReal xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c;
00136 
00137    s = sin(angle * PI/180.0);
00138    c = cos(angle * PI/180.0);
00139 
00140    mag = v.length();
00141 
00142    if (mag == 0.0) {
00143       /* generate an identity matrix and return */
00144       for ( int i = 0; i < 9; ++i ) m[i] = 0.0;
00145       m[0] = m[4] = m[8] = 1.0;
00146       return;
00147    }
00148 
00149    // normalize the vector 
00150    v /= mag;
00151 
00152    /*
00153     *     Arbitrary axis rotation matrix.
00154     *
00155     *  This is composed of 5 matrices, Rz, Ry, T, Ry', Rz', multiplied
00156     *  like so:  Rz * Ry * T * Ry' * Rz'.  T is the final rotation
00157     *  (which is about the X-axis), and the two composite transforms
00158     *  Ry' * Rz' and Rz * Ry are (respectively) the rotations necessary
00159     *  from the arbitrary axis to the X-axis then back.  They are
00160     *  all elementary rotations.
00161     *
00162     *  Rz' is a rotation about the Z-axis, to bring the axis vector
00163     *  into the x-z plane.  Then Ry' is applied, rotating about the
00164     *  Y-axis to bring the axis vector parallel with the X-axis.  The
00165     *  rotation about the X-axis is then performed.  Ry and Rz are
00166     *  simply the respective inverse transforms to bring the arbitrary
00167     *  axis back to it's original orientation.  The first transforms
00168     *  Rz' and Ry' are considered inverses, since the data from the
00169     *  arbitrary axis gives you info on how to get to it, not how
00170     *  to get away from it, and an inverse must be applied.
00171     *
00172     *  The basic calculation used is to recognize that the arbitrary
00173     *  axis vector (x, y, z), since it is of unit length, actually
00174     *  represents the sines and cosines of the angles to rotate the
00175     *  X-axis to the same orientation, with theta being the angle about
00176     *  Z and phi the angle about Y (in the order described above)
00177     *  as follows:
00178     *
00179     *  cos ( theta ) = x / sqrt ( 1 - z^2 )
00180     *  sin ( theta ) = y / sqrt ( 1 - z^2 )
00181     *
00182     *  cos ( phi ) = sqrt ( 1 - z^2 )
00183     *  sin ( phi ) = z
00184     *
00185     *  Note that cos ( phi ) can further be inserted to the above
00186     *  formulas:
00187     *
00188     *  cos ( theta ) = x / cos ( phi )
00189     *  sin ( theta ) = y / sin ( phi )
00190     *
00191     *  ...etc.  Because of those relations and the standard trigonometric
00192     *  relations, it is pssible to reduce the transforms down to what
00193     *  is used below.  It may be that any primary axis chosen will give the
00194     *  same results (modulo a sign convention) using thie method.
00195     *
00196     *  Particularly nice is to notice that all divisions that might
00197     *  have caused trouble when parallel to certain planes or
00198     *  axis go away with care paid to reducing the expressions.
00199     *  After checking, it does perform correctly under all cases, since
00200     *  in all the cases of division where the denominator would have
00201     *  been zero, the numerator would have been zero as well, giving
00202     *  the expected result.
00203     */
00204 
00205    // store matrix in (row, col) form, i.e., m(row,col) = m[row*3+col]
00206 
00207    xx = v.x * v.x;
00208    yy = v.y * v.y;
00209    zz = v.z * v.z;
00210    xy = v.x * v.y;
00211    yz = v.y * v.z;
00212    zx = v.z * v.x;
00213    xs = v.x * s;
00214    ys = v.y * s;
00215    zs = v.z * s;
00216    one_c = 1.0 - c;
00217 
00218    m[0] = (one_c * xx) + c;
00219    m[1] = (one_c * xy) - zs;
00220    m[2] = (one_c * zx) + ys;
00221    
00222    m[3] = (one_c * xy) + zs;
00223    m[4] = (one_c * yy) + c;
00224    m[5] = (one_c * yz) - xs;
00225    
00226    m[6] = (one_c * zx) - ys;
00227    m[7] = (one_c * yz) + xs;
00228    m[8] = (one_c * zz) + c;
00229 }
00230 
00231 // multiply vector v by a 3x3 matrix m stored so that m(row,col) = m[row*3+col]
00232 Vector mat_multiply_vec(const Vector &v, BigReal m[]) {
00233       return Vector( m[0]*v.x + m[1]*v.y + m[2]*v.z,
00234                      m[3]*v.x + m[4]*v.y + m[5]*v.z,
00235                      m[6]*v.x + m[7]*v.y + m[8]*v.z );
00236 }
00237 

Generated on Sat Nov 18 01:17:15 2017 for NAMD by  doxygen 1.4.7