NamdOneTools.C File Reference

#include "InfoStream.h"
#include "common.h"
#include "NamdTypes.h"
#include "NamdOneTools.h"
#include "Vector.h"
#include "PDB.h"
#include "Molecule.h"
#include "Debug.h"

Go to the source code of this file.

Macros
#define	MIN_DEBUG_LEVEL   4

Functions
void	read_binary_coors (char *fname, PDB *pdbobj)
void	read_binary_file (const char *fname, Vector *data, int n)
void	vec_rotation_matrix (BigReal angle, Vector v, BigReal m[])
Vector	mat_multiply_vec (const Vector &v, BigReal m[])

Macro Definition Documentation

MIN_DEBUG_LEVEL

#define MIN_DEBUG_LEVEL   4

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 18 of file NamdOneTools.C.

Function Documentation

mat_multiply_vec()

Vector mat_multiply_vec	(	const Vector &	v,
		BigReal	m[]
	)

Definition at line 232 of file NamdOneTools.C.

References Vector::x, Vector::y, and Vector::z.

Referenced by ComputeRestraints::doForce().

                                                      {
      return Vector( m[0]*v.x + m[1]*v.y + m[2]*v.z,
                     m[3]*v.x + m[4]*v.y + m[5]*v.z,
                     m[6]*v.x + m[7]*v.y + m[8]*v.z );
}

read_binary_coors()

void read_binary_coors	(	char *	fname,
		PDB *	pdbobj
	)

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 34 of file NamdOneTools.C.

References PDB::num_atoms(), read_binary_file(), and PDB::set_all_positions().

Referenced by NamdState::loadStructure().

                                                 {
  Vector *newcoords;    //  Array of vectors to hold coordinates from file

  //  Allocate an array to hold the new coordinates
  newcoords = new Vector[pdbobj->num_atoms()];

  //  Read the coordinate from the file
  read_binary_file(fname,newcoords,pdbobj->num_atoms());

  //  Set the coordinates in the PDB object to the new coordinates
  pdbobj->set_all_positions(newcoords);

  //  Clean up
  delete [] newcoords;

} // END OF FUNCTION read_binary_coors()

read_binary_file()

void read_binary_file	(	const char *	fname,
		Vector *	data,
		int	n
	)

Definition at line 52 of file NamdOneTools.C.

References endi(), Fclose(), Fopen(), iINFO(), iout, iWARN(), and NAMD_die().

Referenced by WorkDistrib::createAtomLists(), and read_binary_coors().

{
  int32 filen;          //  Number of atoms read from file
  FILE *fp;             //  File descriptor
  int needToFlip = 0;

  iout << iINFO << "Reading from binary file " << fname << "\n" << endi;

  //  Open the file and die if the open fails
  if ( (fp = Fopen(fname, "rb")) == NULL)
  {
    char errmsg[256];
    sprintf(errmsg, "Unable to open binary file %s", fname);
    NAMD_die(errmsg);
  }

  //  read the number of coordinates in this file
  if (fread(&filen, sizeof(int32), 1, fp) != (size_t)1)
  {
    char errmsg[256];
    sprintf(errmsg, "Error reading binary file %s", fname);
    NAMD_die(errmsg);
  }

  //  read the number of coordinates in this file
  //  check for palindromic number of atoms
  char lenbuf[4];
  memcpy(lenbuf, (const char *)&filen, 4);
  char tmpc;
  tmpc = lenbuf[0]; lenbuf[0] = lenbuf[3]; lenbuf[3] = tmpc;
  tmpc = lenbuf[1]; lenbuf[1] = lenbuf[2]; lenbuf[2] = tmpc;
  if ( ! memcmp((const char *)&filen, lenbuf, 4) ) {
    iout << iWARN << "Number of atoms in binary file " << fname <<
                " is palindromic, assuming same endian.\n" << endi;
  }

  //  Die if this doesn't match the number in our system
  if (filen != n)
  {
    needToFlip = 1;
    memcpy((char *)&filen, lenbuf, 4);
  }
  if (filen != n)
  {
    char errmsg[256];
    sprintf(errmsg, "Incorrect atom count in binary file %s", fname);
    NAMD_die(errmsg);
  }

  if (fread(data, sizeof(Vector), n, fp) != (size_t)n)
  {
    char errmsg[256];
    sprintf(errmsg, "Error reading binary file %s", fname);
    NAMD_die(errmsg);
  }

  Fclose(fp);

  if (needToFlip) { 
    iout << iWARN << "Converting binary file " << fname << "\n" << endi;
    int i;
    char *cdata = (char *) data;
    for ( i=0; i<3*n; ++i, cdata+=8 ) {
      char tmp0, tmp1, tmp2, tmp3;
      tmp0 = cdata[0]; tmp1 = cdata[1];
      tmp2 = cdata[2]; tmp3 = cdata[3];
      cdata[0] = cdata[7]; cdata[1] = cdata[6];
      cdata[2] = cdata[5]; cdata[3] = cdata[4];
      cdata[7] = tmp0; cdata[6] = tmp1;
      cdata[5] = tmp2; cdata[4] = tmp3;
    }
  }

}

vec_rotation_matrix()

void vec_rotation_matrix	(	BigReal	angle,
		Vector	v,
		BigReal	m[]
	)

Definition at line 132 of file NamdOneTools.C.

References Vector::length(), PI, Vector::x, xx, xy, Vector::y, yy, yz, Vector::z, zx, and zz.

Referenced by ComputeRestraints::doForce().

                                                                 {
   /* This function contributed by Erich Boleyn (erich@uruk.org) */
   BigReal mag, s, c;
   BigReal xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c;

   s = sin(angle * PI/180.0);
   c = cos(angle * PI/180.0);

   mag = v.length();

   if (mag == 0.0) {
      /* generate an identity matrix and return */
      for ( int i = 0; i < 9; ++i ) m[i] = 0.0;
      m[0] = m[4] = m[8] = 1.0;
      return;
   }

   // normalize the vector 
   v /= mag;

   /*
    *     Arbitrary axis rotation matrix.
    *
    *  This is composed of 5 matrices, Rz, Ry, T, Ry', Rz', multiplied
    *  like so:  Rz * Ry * T * Ry' * Rz'.  T is the final rotation
    *  (which is about the X-axis), and the two composite transforms
    *  Ry' * Rz' and Rz * Ry are (respectively) the rotations necessary
    *  from the arbitrary axis to the X-axis then back.  They are
    *  all elementary rotations.
    *
    *  Rz' is a rotation about the Z-axis, to bring the axis vector
    *  into the x-z plane.  Then Ry' is applied, rotating about the
    *  Y-axis to bring the axis vector parallel with the X-axis.  The
    *  rotation about the X-axis is then performed.  Ry and Rz are
    *  simply the respective inverse transforms to bring the arbitrary
    *  axis back to it's original orientation.  The first transforms
    *  Rz' and Ry' are considered inverses, since the data from the
    *  arbitrary axis gives you info on how to get to it, not how
    *  to get away from it, and an inverse must be applied.
    *
    *  The basic calculation used is to recognize that the arbitrary
    *  axis vector (x, y, z), since it is of unit length, actually
    *  represents the sines and cosines of the angles to rotate the
    *  X-axis to the same orientation, with theta being the angle about
    *  Z and phi the angle about Y (in the order described above)
    *  as follows:
    *
    *  cos ( theta ) = x / sqrt ( 1 - z^2 )
    *  sin ( theta ) = y / sqrt ( 1 - z^2 )
    *
    *  cos ( phi ) = sqrt ( 1 - z^2 )
    *  sin ( phi ) = z
    *
    *  Note that cos ( phi ) can further be inserted to the above
    *  formulas:
    *
    *  cos ( theta ) = x / cos ( phi )
    *  sin ( theta ) = y / sin ( phi )
    *
    *  ...etc.  Because of those relations and the standard trigonometric
    *  relations, it is pssible to reduce the transforms down to what
    *  is used below.  It may be that any primary axis chosen will give the
    *  same results (modulo a sign convention) using thie method.
    *
    *  Particularly nice is to notice that all divisions that might
    *  have caused trouble when parallel to certain planes or
    *  axis go away with care paid to reducing the expressions.
    *  After checking, it does perform correctly under all cases, since
    *  in all the cases of division where the denominator would have
    *  been zero, the numerator would have been zero as well, giving
    *  the expected result.
    */

   // store matrix in (row, col) form, i.e., m(row,col) = m[row*3+col]

   xx = v.x * v.x;
   yy = v.y * v.y;
   zz = v.z * v.z;
   xy = v.x * v.y;
   yz = v.y * v.z;
   zx = v.z * v.x;
   xs = v.x * s;
   ys = v.y * s;
   zs = v.z * s;
   one_c = 1.0 - c;

   m[0] = (one_c * xx) + c;
   m[1] = (one_c * xy) - zs;
   m[2] = (one_c * zx) + ys;
   
   m[3] = (one_c * xy) + zs;
   m[4] = (one_c * yy) + c;
   m[5] = (one_c * yz) - xs;
   
   m[6] = (one_c * zx) - ys;
   m[7] = (one_c * yz) + xs;
   m[8] = (one_c * zz) + c;
}