OOPSE/libmdtools/Symplectic.cpp

#include <iostream>
#include <cstdlib>

#include "Integrator.hpp"
#include "Thermo.hpp"
#include "ReadWrite.hpp"
#include "ForceFields.hpp"
#include "ExtendedSystem.hpp"
#include "simError.h"

extern "C"{
  
  void v_constrain_a_( double &dt, int &n_atoms, double* mass, 
                       double* Rx, double* Ry, double* Rz, 
                       double* Vx, double* Vy, double* Vz, 
                       double* Fx, double* Fy, double* Fz, 
                       int &n_constrained, double *constr_sqr,
                       int* constr_i, int* constr_j, 
                       double &box_x, double &box_y, double &box_z );

  void v_constrain_b_( double &dt, int &n_atoms, double* mass, 
                       double* Rx, double* Ry, double* Rz, 
                       double* Vx, double* Vy, double* Vz, 
                       double* Fx, double* Fy, double* Fz, 
                       double &Kinetic, 
                       int &n_constrained, double *constr_sqr,
                       int* constr_i, int* constr_j, 
                       double &box_x, double &box_y, double &box_z );
}


Symplectic::Symplectic( SimInfo* the_entry_plug, ForceFields* the_ff,
                        ExtendedSystem* the_es ){
  entry_plug = the_entry_plug;
  myFF = the_ff;
  myES = the_es;
  isFirst = 1;

  std::cerr<< "calling symplectic constructor\n";

  molecules = entry_plug->molecules;
  nMols = entry_plug->n_mol;

  // give a little love back to the SimInfo object
  
  if( entry_plug->the_integrator != NULL ) delete entry_plug->the_integrator;
  entry_plug->the_integrator = this;

  // grab the masses

  mass = new double[entry_plug->n_atoms];
  for(int i = 0; i < entry_plug->n_atoms; i++){
    mass[i] = entry_plug->atoms[i]->getMass();
  }  
 
  // check for constraints

  is_constrained = 0;

  Constraint *temp_con;
  Constraint *dummy_plug;
  temp_con = new Constraint[entry_plug->n_SRI];
  n_constrained = 0;
  int constrained = 0;
  
  SRI** theArray;
  for(int i = 0; i < nMols; i++){
    
    theArray = (SRI**) molecules[i].getMyBonds();
    for(int j=0; j<molecules[i].getNBonds(); j++){
      
      constrained = theArray[j]->is_constrained();
      
      if(constrained){
        
        dummy_plug = theArray[j]->get_constraint();
        temp_con[n_constrained].set_a( dummy_plug->get_a() );
        temp_con[n_constrained].set_b( dummy_plug->get_b() );
        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
        
        n_constrained++;
        constrained = 0;
      }
    }

    theArray = (SRI**) molecules[i].getMyBends();
    for(int j=0; j<molecules[i].getNBends(); j++){
      
      constrained = theArray[j]->is_constrained();
      
      if(constrained){
        
        dummy_plug = theArray[j]->get_constraint();
        temp_con[n_constrained].set_a( dummy_plug->get_a() );
        temp_con[n_constrained].set_b( dummy_plug->get_b() );
        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
        
        n_constrained++;
        constrained = 0;
      }
    }

    theArray = (SRI**) molecules[i].getMyTorsions();
    for(int j=0; j<molecules[i].getNTorsions(); j++){
      
      constrained = theArray[j]->is_constrained();
      
      if(constrained){
        
        dummy_plug = theArray[j]->get_constraint();
        temp_con[n_constrained].set_a( dummy_plug->get_a() );
        temp_con[n_constrained].set_b( dummy_plug->get_b() );
        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
        
        n_constrained++;
        constrained = 0;
      }
    }
  }

  if(n_constrained > 0){
    
    is_constrained = 1;
    constrained_i = new int[n_constrained];
    constrained_j = new int[n_constrained];
    constrained_dsqr = new double[n_constrained];
    
    for( int i = 0; i < n_constrained; i++){
      
      /* add 1 to the index for the fortran arrays. */

      constrained_i[i] = temp_con[i].get_a() + 1;
      constrained_j[i] = temp_con[i].get_b() + 1;
      constrained_dsqr[i] = temp_con[i].get_dsqr();
    }
  }
  
  delete[] temp_con;
}

Symplectic::~Symplectic() {
  
  if( n_constrained ){
    delete[] constrained_i;
    delete[] constrained_j;
    delete[] constrained_dsqr;
  }
  
}


void Symplectic::integrate( void ){

  const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2

  int i, j;                         // loop counters
  int nAtoms = entry_plug->n_atoms; // the number of atoms
  double kE = 0.0;                  // the kinetic energy  
  double rot_kE;
  double trans_kE;
  int tl;                        // the time loop conter
  double dt2;                       // half the dt

  double vx, vy, vz;    // the velocities
//  double vx2, vy2, vz2; // the square of the velocities
  double rx, ry, rz;    // the postitions
  
  double ji[3];   // the body frame angular momentum
  double jx2, jy2, jz2; // the square of the angular momentums
  double Tb[3];   // torque in the body frame
  double angle;   // the angle through which to rotate the rotation matrix
  double A[3][3]; // the rotation matrix

  int time;

  double dt          = entry_plug->dt;
  double runTime     = entry_plug->run_time;
  double sampleTime  = entry_plug->sampleTime;
  double statusTime  = entry_plug->statusTime;
  double thermalTime = entry_plug->thermalTime;

  int n_loops  = (int)( runTime / dt );
  int sample_n = (int)( sampleTime / dt );
  int status_n = (int)( statusTime / dt );
  int vel_n    = (int)( thermalTime / dt );

  int calcPot, calcStress;

  Thermo *tStats;
  StatWriter*  e_out;
  DumpWriter*  dump_out;

  std::cerr << "about to call new thermo\n";

  tStats   = new Thermo( entry_plug );
  e_out    = new StatWriter( entry_plug );

  std::cerr << "calling dumpWriter \n";
  dump_out = new DumpWriter( entry_plug );
  std::cerr << "called dumpWriter \n";

  Atom** atoms = entry_plug->atoms;
  DirectionalAtom* dAtom;
  dt2 = 0.5 * dt;

  // initialize the forces the before the first step

  myFF->doForces(1,1);
  
  if( entry_plug->setTemp ){
    
    tStats->velocitize();
  }
  
  dump_out->writeDump( 0.0 );
  e_out->writeStat( 0.0 );
  
  calcPot = 0;

  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
    calcStress = 1;
  } else {
    calcStress = 0;
  }

  if( n_constrained ){

    double *Rx = new double[nAtoms];
    double *Ry = new double[nAtoms];
    double *Rz = new double[nAtoms];
    
    double *Vx = new double[nAtoms];
    double *Vy = new double[nAtoms];
    double *Vz = new double[nAtoms];
    
    double *Fx = new double[nAtoms];
    double *Fy = new double[nAtoms];
    double *Fz = new double[nAtoms];
    

    for( tl=0; tl < n_loops; tl++ ){

      if (!strcasecmp( entry_plug->ensemble, "NVT")) 
        myES->NoseHooverNVT( dt , tStats->getKinetic() );
      
      for( j=0; j<nAtoms; j++ ){

        Rx[j] = atoms[j]->getX();
        Ry[j] = atoms[j]->getY();
        Rz[j] = atoms[j]->getZ();

        Vx[j] = atoms[j]->get_vx();
        Vy[j] = atoms[j]->get_vy();
        Vz[j] = atoms[j]->get_vz();

        Fx[j] = atoms[j]->getFx();
        Fy[j] = atoms[j]->getFy();
        Fz[j] = atoms[j]->getFz();

      }
        
      v_constrain_a_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz, 
                      Fx, Fy, Fz,
                      n_constrained, constrained_dsqr, 
                      constrained_i, constrained_j,
                      entry_plug->box_x, 
                      entry_plug->box_y, 
                      entry_plug->box_z );
      
      for( j=0; j<nAtoms; j++ ){

        atoms[j]->setX(Rx[j]);
        atoms[j]->setY(Ry[j]);
        atoms[j]->setZ(Rz[j]);
        
        atoms[j]->set_vx(Vx[j]);
        atoms[j]->set_vy(Vy[j]);
        atoms[j]->set_vz(Vz[j]);
      }


//       for( i=0; i<nAtoms; i++ ){
// //   if( atoms[i]->isDirectional() ){
                  
// //     dAtom = (DirectionalAtom *)atoms[i];
          
// //     // get and convert the torque to body frame
          
// //     Tb[0] = dAtom->getTx();
// //     Tb[1] = dAtom->getTy();
// //     Tb[2] = dAtom->getTz();
          
// //     dAtom->lab2Body( Tb );
          
// //     // get the angular momentum, and propagate a half step
          
// //     ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
// //     ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
// //     ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
          
// //     // get the atom's rotation matrix
          
// //     A[0][0] = dAtom->getAxx();
// //     A[0][1] = dAtom->getAxy();
// //     A[0][2] = dAtom->getAxz();
          
// //     A[1][0] = dAtom->getAyx();
// //     A[1][1] = dAtom->getAyy();
// //     A[1][2] = dAtom->getAyz();
          
// //     A[2][0] = dAtom->getAzx();
// //     A[2][1] = dAtom->getAzy();
// //     A[2][2] = dAtom->getAzz();
          
          
// //     // use the angular velocities to propagate the rotation matrix a
// //     // full time step
          
          
// //     angle = dt2 * ji[0] / dAtom->getIxx();
// //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
          
// //     angle = dt2 * ji[1] / dAtom->getIyy();
// //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
          
// //     angle = dt * ji[2] / dAtom->getIzz();
// //     this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
          
// //     angle = dt2 * ji[1] / dAtom->getIyy();
// //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
          
// //     angle = dt2 * ji[0] / dAtom->getIxx();
// //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
          
          
// //     dAtom->setA( A );
// //     dAtom->setJx( ji[0] );
// //     dAtom->setJy( ji[1] );
// //     dAtom->setJz( ji[2] );
// //   }
//       }
      
      // calculate the forces
      
      myFF->doForces(calcPot, calcStress);
      
      // move b

      for( j=0; j<nAtoms; j++ ){

        Rx[j] = atoms[j]->getX();
        Ry[j] = atoms[j]->getY();
        Rz[j] = atoms[j]->getZ();

        Vx[j] = atoms[j]->get_vx();
        Vy[j] = atoms[j]->get_vy();
        Vz[j] = atoms[j]->get_vz();

        Fx[j] = atoms[j]->getFx();
        Fy[j] = atoms[j]->getFy();
        Fz[j] = atoms[j]->getFz();
      }
        
      v_constrain_b_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz, 
                      Fx, Fy, Fz,
                      kE, n_constrained, constrained_dsqr, 
                      constrained_i, constrained_j,
                      entry_plug->box_x,
                      entry_plug->box_y, 
                      entry_plug->box_z );
      
      for( j=0; j<nAtoms; j++ ){

        atoms[j]->setX(Rx[j]);
        atoms[j]->setY(Ry[j]);
        atoms[j]->setZ(Rz[j]);

        atoms[j]->set_vx(Vx[j]);
        atoms[j]->set_vy(Vy[j]);
        atoms[j]->set_vz(Vz[j]);
      }
      
//       for( i=0; i< nAtoms; i++ ){

//      if( atoms[i]->isDirectional() ){

//        dAtom = (DirectionalAtom *)atoms[i];
          
//        // get and convert the torque to body frame
          
//        Tb[0] = dAtom->getTx();
//        Tb[1] = dAtom->getTy();
//        Tb[2] = dAtom->getTz();
          
//        dAtom->lab2Body( Tb );
          
//        // get the angular momentum, and complete the angular momentum
//        // half step
          
//        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
//        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
//        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
          
//        dAtom->setJx( ji[0] );
//        dAtom->setJy( ji[1] );
//        dAtom->setJz( ji[2] );
//      }
//       }
     

      if (!strcasecmp( entry_plug->ensemble, "NVT")) 
        myES->NoseHooverNVT( dt , tStats->getKinetic() );

      if (!strcasecmp( entry_plug->ensemble, "NPT") )
        myES->NoseHooverAndersonNPT( dt, 
                                     tStats->getKinetic(), 
                                     tStats->getPressure());

      time = tl + 1;
      
      if( entry_plug->setTemp ){
        if( !(time % vel_n) ) tStats->velocitize();
      }
      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
      if( !((time+1) % status_n) ) {
        calcPot = 1;
        // bitwise masking in case we need it for NPT
        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
      }
      if( !(time % status_n) ){ 
        e_out->writeStat( time * dt ); 
        calcPot = 0; 
        // bitwise masking in case we need it for NPT
        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
      }
    }
  }
  else{

    for( tl=0; tl<n_loops; tl++ ){
      
      kE = 0.0;
      rot_kE= 0.0;
      trans_kE = 0.0;

      if (!strcasecmp( entry_plug->ensemble, "NVT")) 
        myES->NoseHooverNVT( dt , tStats->getKinetic() );
      
      for( i=0; i<nAtoms; i++ ){
        
        // velocity half step
        
        vx = atoms[i]->get_vx() + 
          ( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
        vy = atoms[i]->get_vy() + 
          ( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
        vz = atoms[i]->get_vz() + 
          ( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
        
        // position whole step
        
        rx = atoms[i]->getX() + dt * vx;
        ry = atoms[i]->getY() + dt * vy;
        rz = atoms[i]->getZ() + dt * vz;
        
        atoms[i]->setX( rx );
        atoms[i]->setY( ry );
        atoms[i]->setZ( rz );
        
        atoms[i]->set_vx( vx );
        atoms[i]->set_vy( vy );
        atoms[i]->set_vz( vz );
        
//      if( atoms[i]->isDirectional() ){

//        dAtom = (DirectionalAtom *)atoms[i];
          
//        // get and convert the torque to body frame
          
//        Tb[0] = dAtom->getTx();
//        Tb[1] = dAtom->getTy();
//        Tb[2] = dAtom->getTz();
          
//        dAtom->lab2Body( Tb );
          
//        // get the angular momentum, and propagate a half step
          
//        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
//        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
//        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
          
//        // get the atom's rotation matrix
          
//        A[0][0] = dAtom->getAxx();
//        A[0][1] = dAtom->getAxy();
//        A[0][2] = dAtom->getAxz();
          
//        A[1][0] = dAtom->getAyx();
//        A[1][1] = dAtom->getAyy();
//        A[1][2] = dAtom->getAyz();
          
//        A[2][0] = dAtom->getAzx();
//        A[2][1] = dAtom->getAzy();
//        A[2][2] = dAtom->getAzz();
          
          
//        // use the angular velocities to propagate the rotation matrix a
//        // full time step
          
          
//        angle = dt2 * ji[0] / dAtom->getIxx();
//        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
          
//        angle = dt2 * ji[1] / dAtom->getIyy();
//        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
          
//        angle = dt * ji[2] / dAtom->getIzz();
//        this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
          
//        angle = dt2 * ji[1] / dAtom->getIyy();
//        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
          
//        angle = dt2 * ji[0] / dAtom->getIxx();
//        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
          
          
//        dAtom->setA( A );
//        dAtom->setJx( ji[0] );
//        dAtom->setJy( ji[1] );
//        dAtom->setJz( ji[2] );
//      }
      }
      
      // calculate the forces
      
      myFF->doForces(calcPot,calcStress);
      
      for( i=0; i< nAtoms; i++ ){
        
        // complete the velocity half step
        
        vx = atoms[i]->get_vx() + 
          ( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
        vy = atoms[i]->get_vy() + 
          ( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
        vz = atoms[i]->get_vz() + 
          ( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
        
        atoms[i]->set_vx( vx );
        atoms[i]->set_vy( vy );
        atoms[i]->set_vz( vz );
        
//      vx2 = vx * vx;
//      vy2 = vy * vy;
//      vz2 = vz * vz;
        
//      if( atoms[i]->isDirectional() ){

//        dAtom = (DirectionalAtom *)atoms[i];
          
//        // get and convert the torque to body frame
          
//        Tb[0] = dAtom->getTx();
//        Tb[1] = dAtom->getTy();
//        Tb[2] = dAtom->getTz();
          
//        dAtom->lab2Body( Tb );
          
//        // get the angular momentum, and complete the angular momentum
//        // half step
          
//        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
//        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
//        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
          
//        jx2 = ji[0] * ji[0];
//        jy2 = ji[1] * ji[1];
//        jz2 = ji[2] * ji[2];
          
//        rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy()) 
//          + (jz2 / dAtom->getIzz());
          
//        dAtom->setJx( ji[0] );
//        dAtom->setJy( ji[1] );
//        dAtom->setJz( ji[2] );
//      }
      }
    
      if (!strcasecmp( entry_plug->ensemble, "NVT")) 
        myES->NoseHooverNVT( dt , tStats->getKinetic() );

      if (!strcasecmp( entry_plug->ensemble, "NPT") )
        myES->NoseHooverAndersonNPT( dt, 
                                     tStats->getKinetic(), 
                                     tStats->getPressure());
  
      time = tl + 1;
      
      if( entry_plug->setTemp ){
        if( !(time % vel_n) ) tStats->velocitize();
      }
      if( !(time % sample_n) ) dump_out->writeDump( time * dt );
      if( !((time+1) % status_n) ) {
        calcPot = 1;
        // bitwise masking in case we need it for NPT
        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
      }
      if( !(time % status_n) ){ 
        e_out->writeStat( time * dt ); 
        calcPot = 0; 
        // bitwise masking in case we need it for NPT
        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
      }      
    }
  }

  dump_out->writeFinal();

  delete dump_out;
  delete e_out;
}

void Symplectic::rotate( int axes1, int axes2, double angle, double ji[3], 
                         double A[3][3] ){

  int i,j,k;
  double sinAngle;
  double cosAngle;
  double angleSqr;
  double angleSqrOver4;
  double top, bottom;
  double rot[3][3];
  double tempA[3][3];
  double tempJ[3];

  // initialize the tempA

  for(i=0; i<3; i++){
    for(j=0; j<3; j++){
      tempA[j][i] = A[i][j];
    }
  }

  // initialize the tempJ

  for( i=0; i<3; i++) tempJ[i] = ji[i];
  
  // initalize rot as a unit matrix

  rot[0][0] = 1.0;
  rot[0][1] = 0.0;
  rot[0][2] = 0.0;

  rot[1][0] = 0.0;
  rot[1][1] = 1.0;
  rot[1][2] = 0.0;
  
  rot[2][0] = 0.0;
  rot[2][1] = 0.0;
  rot[2][2] = 1.0;
  
  // use a small angle aproximation for sin and cosine

  angleSqr  = angle * angle;
  angleSqrOver4 = angleSqr / 4.0;
  top = 1.0 - angleSqrOver4;
  bottom = 1.0 + angleSqrOver4;

  cosAngle = top / bottom;
  sinAngle = angle / bottom;

  rot[axes1][axes1] = cosAngle;
  rot[axes2][axes2] = cosAngle;

  rot[axes1][axes2] = sinAngle;
  rot[axes2][axes1] = -sinAngle;
  
  // rotate the momentum acoording to: ji[] = rot[][] * ji[]
  
  for(i=0; i<3; i++){
    ji[i] = 0.0;
    for(k=0; k<3; k++){
      ji[i] += rot[i][k] * tempJ[k];
    }
  }

  // rotate the Rotation matrix acording to: 
  //            A[][] = A[][] * transpose(rot[][])


  // NOte for as yet unknown reason, we are setting the performing the
  // calculation as:
  //                transpose(A[][]) = transpose(A[][]) * transpose(rot[][])

  for(i=0; i<3; i++){
    for(j=0; j<3; j++){
      A[j][i] = 0.0;
      for(k=0; k<3; k++){
        A[j][i] += tempA[i][k] * rot[j][k];
      }
    }
  }
}
Revision:	475
Committed:	Tue Apr 8 12:44:18 2003 UTC (21 years, 3 months ago) by gezelter
File size:	18147 byte(s)
Log Message:	Changes to integrate the NVT and NPT ensembles
#	User	Rev	Content
1	mmeineke	377	#include <iostream>
2			#include <cstdlib>
3
4			#include "Integrator.hpp"
5			#include "Thermo.hpp"
6			#include "ReadWrite.hpp"
7			#include "ForceFields.hpp"
8	gezelter	466	#include "ExtendedSystem.hpp"
9	mmeineke	377	#include "simError.h"
10
11			extern "C"{
12
13			void v_constrain_a_( double &dt, int &n_atoms, double* mass,
14			double* Rx, double* Ry, double* Rz,
15			double* Vx, double* Vy, double* Vz,
16			double* Fx, double* Fy, double* Fz,
17			int &n_constrained, double *constr_sqr,
18			int* constr_i, int* constr_j,
19			double &box_x, double &box_y, double &box_z );
20
21			void v_constrain_b_( double &dt, int &n_atoms, double* mass,
22			double* Rx, double* Ry, double* Rz,
23			double* Vx, double* Vy, double* Vz,
24			double* Fx, double* Fy, double* Fz,
25			double &Kinetic,
26			int &n_constrained, double *constr_sqr,
27			int* constr_i, int* constr_j,
28			double &box_x, double &box_y, double &box_z );
29			}
30
31
32
33
34	gezelter	466	Symplectic::Symplectic( SimInfo* the_entry_plug, ForceFields* the_ff,
35			ExtendedSystem* the_es ){
36	mmeineke	377	entry_plug = the_entry_plug;
37			myFF = the_ff;
38	gezelter	466	myES = the_es;
39	mmeineke	377	isFirst = 1;
40
41	mmeineke	469	std::cerr<< "calling symplectic constructor\n";
42
43	mmeineke	423	molecules = entry_plug->molecules;
44			nMols = entry_plug->n_mol;
45	mmeineke	377
46			// give a little love back to the SimInfo object
47
48			if( entry_plug->the_integrator != NULL ) delete entry_plug->the_integrator;
49			entry_plug->the_integrator = this;
50
51			// grab the masses
52
53			mass = new double[entry_plug->n_atoms];
54			for(int i = 0; i < entry_plug->n_atoms; i++){
55			mass[i] = entry_plug->atoms[i]->getMass();
56	gezelter	466	}
57	mmeineke	377
58			// check for constraints
59
60			is_constrained = 0;
61
62			Constraint *temp_con;
63			Constraint *dummy_plug;
64	mmeineke	423	temp_con = new Constraint[entry_plug->n_SRI];
65	mmeineke	377	n_constrained = 0;
66			int constrained = 0;
67
68	mmeineke	423	SRI** theArray;
69			for(int i = 0; i < nMols; i++){
70	mmeineke	377
71	mmeineke	428	theArray = (SRI**) molecules[i].getMyBonds();
72			for(int j=0; j<molecules[i].getNBonds(); j++){
73	mmeineke	377
74	mmeineke	423	constrained = theArray[j]->is_constrained();
75
76			if(constrained){
77
78			dummy_plug = theArray[j]->get_constraint();
79	mmeineke	428	temp_con[n_constrained].set_a( dummy_plug->get_a() );
80			temp_con[n_constrained].set_b( dummy_plug->get_b() );
81			temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
82	mmeineke	423
83	mmeineke	428	n_constrained++;
84	mmeineke	423	constrained = 0;
85			}
86			}
87	mmeineke	377
88	mmeineke	428	theArray = (SRI**) molecules[i].getMyBends();
89			for(int j=0; j<molecules[i].getNBends(); j++){
90	mmeineke	423
91			constrained = theArray[j]->is_constrained();
92
93			if(constrained){
94
95			dummy_plug = theArray[j]->get_constraint();
96	mmeineke	428	temp_con[n_constrained].set_a( dummy_plug->get_a() );
97			temp_con[n_constrained].set_b( dummy_plug->get_b() );
98			temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
99	mmeineke	423
100	mmeineke	428	n_constrained++;
101	mmeineke	423	constrained = 0;
102			}
103	mmeineke	377	}
104	mmeineke	423
105	mmeineke	428	theArray = (SRI**) molecules[i].getMyTorsions();
106			for(int j=0; j<molecules[i].getNTorsions(); j++){
107	mmeineke	423
108			constrained = theArray[j]->is_constrained();
109
110			if(constrained){
111
112			dummy_plug = theArray[j]->get_constraint();
113	mmeineke	428	temp_con[n_constrained].set_a( dummy_plug->get_a() );
114			temp_con[n_constrained].set_b( dummy_plug->get_b() );
115			temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
116	mmeineke	423
117	mmeineke	428	n_constrained++;
118	mmeineke	423	constrained = 0;
119			}
120			}
121	mmeineke	377	}
122
123			if(n_constrained > 0){
124
125			is_constrained = 1;
126			constrained_i = new int[n_constrained];
127			constrained_j = new int[n_constrained];
128			constrained_dsqr = new double[n_constrained];
129
130			for( int i = 0; i < n_constrained; i++){
131
132			/* add 1 to the index for the fortran arrays. */
133
134			constrained_i[i] = temp_con[i].get_a() + 1;
135			constrained_j[i] = temp_con[i].get_b() + 1;
136			constrained_dsqr[i] = temp_con[i].get_dsqr();
137			}
138			}
139
140			delete[] temp_con;
141			}
142
143			Symplectic::~Symplectic() {
144
145			if( n_constrained ){
146			delete[] constrained_i;
147			delete[] constrained_j;
148			delete[] constrained_dsqr;
149			}
150
151			}
152
153
154			void Symplectic::integrate( void ){
155
156			const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
157
158			int i, j; // loop counters
159			int nAtoms = entry_plug->n_atoms; // the number of atoms
160			double kE = 0.0; // the kinetic energy
161			double rot_kE;
162			double trans_kE;
163			int tl; // the time loop conter
164			double dt2; // half the dt
165
166			double vx, vy, vz; // the velocities
167			// double vx2, vy2, vz2; // the square of the velocities
168			double rx, ry, rz; // the postitions
169
170			double ji[3]; // the body frame angular momentum
171			double jx2, jy2, jz2; // the square of the angular momentums
172			double Tb[3]; // torque in the body frame
173			double angle; // the angle through which to rotate the rotation matrix
174			double A[3][3]; // the rotation matrix
175
176			int time;
177
178			double dt = entry_plug->dt;
179			double runTime = entry_plug->run_time;
180			double sampleTime = entry_plug->sampleTime;
181			double statusTime = entry_plug->statusTime;
182			double thermalTime = entry_plug->thermalTime;
183
184			int n_loops = (int)( runTime / dt );
185			int sample_n = (int)( sampleTime / dt );
186			int status_n = (int)( statusTime / dt );
187			int vel_n = (int)( thermalTime / dt );
188
189	gezelter	468	int calcPot, calcStress;
190	mmeineke	377
191	mmeineke	469	Thermo *tStats;
192			StatWriter* e_out;
193			DumpWriter* dump_out;
194	mmeineke	377
195	mmeineke	469	std::cerr << "about to call new thermo\n";
196	mmeineke	377
197	mmeineke	469	tStats = new Thermo( entry_plug );
198			e_out = new StatWriter( entry_plug );
199
200			std::cerr << "calling dumpWriter \n";
201			dump_out = new DumpWriter( entry_plug );
202			std::cerr << "called dumpWriter \n";
203
204	mmeineke	377	Atom** atoms = entry_plug->atoms;
205			DirectionalAtom* dAtom;
206			dt2 = 0.5 * dt;
207
208			// initialize the forces the before the first step
209
210	gezelter	468	myFF->doForces(1,1);
211	mmeineke	377
212			if( entry_plug->setTemp ){
213
214			tStats->velocitize();
215			}
216
217			dump_out->writeDump( 0.0 );
218			e_out->writeStat( 0.0 );
219
220			calcPot = 0;
221
222	gezelter	475	if (!strcasecmp( entry_plug->ensemble, "NPT")) {
223			calcStress = 1;
224			} else {
225			calcStress = 0;
226			}
227
228	mmeineke	377	if( n_constrained ){
229
230			double *Rx = new double[nAtoms];
231			double *Ry = new double[nAtoms];
232			double *Rz = new double[nAtoms];
233
234			double *Vx = new double[nAtoms];
235			double *Vy = new double[nAtoms];
236			double *Vz = new double[nAtoms];
237
238			double *Fx = new double[nAtoms];
239			double *Fy = new double[nAtoms];
240			double *Fz = new double[nAtoms];
241
242
243			for( tl=0; tl < n_loops; tl++ ){
244	gezelter	475
245			if (!strcasecmp( entry_plug->ensemble, "NVT"))
246			myES->NoseHooverNVT( dt , tStats->getKinetic() );
247	mmeineke	377
248			for( j=0; j<nAtoms; j++ ){
249
250			Rx[j] = atoms[j]->getX();
251			Ry[j] = atoms[j]->getY();
252			Rz[j] = atoms[j]->getZ();
253
254			Vx[j] = atoms[j]->get_vx();
255			Vy[j] = atoms[j]->get_vy();
256			Vz[j] = atoms[j]->get_vz();
257
258			Fx[j] = atoms[j]->getFx();
259			Fy[j] = atoms[j]->getFy();
260			Fz[j] = atoms[j]->getFz();
261
262			}
263
264			v_constrain_a_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz,
265			Fx, Fy, Fz,
266			n_constrained, constrained_dsqr,
267			constrained_i, constrained_j,
268			entry_plug->box_x,
269			entry_plug->box_y,
270			entry_plug->box_z );
271
272			for( j=0; j<nAtoms; j++ ){
273
274			atoms[j]->setX(Rx[j]);
275			atoms[j]->setY(Ry[j]);
276			atoms[j]->setZ(Rz[j]);
277
278			atoms[j]->set_vx(Vx[j]);
279			atoms[j]->set_vy(Vy[j]);
280			atoms[j]->set_vz(Vz[j]);
281			}
282
283
284	mmeineke	472	// for( i=0; i<nAtoms; i++ ){
285			// // if( atoms[i]->isDirectional() ){
286	mmeineke	377
287	mmeineke	472	// // dAtom = (DirectionalAtom *)atoms[i];
288	mmeineke	377
289	mmeineke	472	// // // get and convert the torque to body frame
290	mmeineke	377
291	mmeineke	472	// // Tb[0] = dAtom->getTx();
292			// // Tb[1] = dAtom->getTy();
293			// // Tb[2] = dAtom->getTz();
294	mmeineke	377
295	mmeineke	472	// // dAtom->lab2Body( Tb );
296	mmeineke	377
297	mmeineke	472	// // // get the angular momentum, and propagate a half step
298	mmeineke	377
299	mmeineke	472	// // ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
300			// // ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
301			// // ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
302	mmeineke	377
303	mmeineke	472	// // // get the atom's rotation matrix
304	mmeineke	377
305	mmeineke	472	// // A[0][0] = dAtom->getAxx();
306			// // A[0][1] = dAtom->getAxy();
307			// // A[0][2] = dAtom->getAxz();
308	mmeineke	377
309	mmeineke	472	// // A[1][0] = dAtom->getAyx();
310			// // A[1][1] = dAtom->getAyy();
311			// // A[1][2] = dAtom->getAyz();
312	mmeineke	377
313	mmeineke	472	// // A[2][0] = dAtom->getAzx();
314			// // A[2][1] = dAtom->getAzy();
315			// // A[2][2] = dAtom->getAzz();
316	mmeineke	377
317
318	mmeineke	472	// // // use the angular velocities to propagate the rotation matrix a
319			// // // full time step
320	mmeineke	377
321
322	mmeineke	472	// // angle = dt2 * ji[0] / dAtom->getIxx();
323			// // this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
324	mmeineke	377
325	mmeineke	472	// // angle = dt2 * ji[1] / dAtom->getIyy();
326			// // this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
327	mmeineke	377
328	mmeineke	472	// // angle = dt * ji[2] / dAtom->getIzz();
329			// // this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
330	mmeineke	377
331	mmeineke	472	// // angle = dt2 * ji[1] / dAtom->getIyy();
332			// // this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
333	mmeineke	377
334	mmeineke	472	// // angle = dt2 * ji[0] / dAtom->getIxx();
335			// // this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
336	mmeineke	377
337
338	mmeineke	472	// // dAtom->setA( A );
339			// // dAtom->setJx( ji[0] );
340			// // dAtom->setJy( ji[1] );
341			// // dAtom->setJz( ji[2] );
342			// // }
343			// }
344	mmeineke	377
345			// calculate the forces
346
347	gezelter	468	myFF->doForces(calcPot, calcStress);
348	mmeineke	377
349			// move b
350
351			for( j=0; j<nAtoms; j++ ){
352
353			Rx[j] = atoms[j]->getX();
354			Ry[j] = atoms[j]->getY();
355			Rz[j] = atoms[j]->getZ();
356
357			Vx[j] = atoms[j]->get_vx();
358			Vy[j] = atoms[j]->get_vy();
359			Vz[j] = atoms[j]->get_vz();
360
361			Fx[j] = atoms[j]->getFx();
362			Fy[j] = atoms[j]->getFy();
363			Fz[j] = atoms[j]->getFz();
364			}
365
366			v_constrain_b_( dt, nAtoms, mass, Rx, Ry, Rz, Vx, Vy, Vz,
367			Fx, Fy, Fz,
368			kE, n_constrained, constrained_dsqr,
369			constrained_i, constrained_j,
370			entry_plug->box_x,
371			entry_plug->box_y,
372			entry_plug->box_z );
373
374			for( j=0; j<nAtoms; j++ ){
375
376			atoms[j]->setX(Rx[j]);
377			atoms[j]->setY(Ry[j]);
378			atoms[j]->setZ(Rz[j]);
379
380			atoms[j]->set_vx(Vx[j]);
381			atoms[j]->set_vy(Vy[j]);
382			atoms[j]->set_vz(Vz[j]);
383			}
384
385	mmeineke	472	// for( i=0; i< nAtoms; i++ ){
386	mmeineke	377
387	mmeineke	472	// if( atoms[i]->isDirectional() ){
388	mmeineke	377
389	mmeineke	472	// dAtom = (DirectionalAtom *)atoms[i];
390	mmeineke	377
391	mmeineke	472	// // get and convert the torque to body frame
392	mmeineke	377
393	mmeineke	472	// Tb[0] = dAtom->getTx();
394			// Tb[1] = dAtom->getTy();
395			// Tb[2] = dAtom->getTz();
396	mmeineke	377
397	mmeineke	472	// dAtom->lab2Body( Tb );
398	mmeineke	377
399	mmeineke	472	// // get the angular momentum, and complete the angular momentum
400			// // half step
401	mmeineke	377
402	mmeineke	472	// ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
403			// ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
404			// ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
405	mmeineke	377
406	mmeineke	472	// dAtom->setJx( ji[0] );
407			// dAtom->setJy( ji[1] );
408			// dAtom->setJz( ji[2] );
409			// }
410			// }
411	mmeineke	377
412	gezelter	475
413			if (!strcasecmp( entry_plug->ensemble, "NVT"))
414			myES->NoseHooverNVT( dt , tStats->getKinetic() );
415
416			if (!strcasecmp( entry_plug->ensemble, "NPT") )
417			myES->NoseHooverAndersonNPT( dt,
418			tStats->getKinetic(),
419			tStats->getPressure());
420
421	mmeineke	377	time = tl + 1;
422
423			if( entry_plug->setTemp ){
424			if( !(time % vel_n) ) tStats->velocitize();
425			}
426			if( !(time % sample_n) ) dump_out->writeDump( time * dt );
427	gezelter	468	if( !((time+1) % status_n) ) {
428			calcPot = 1;
429	gezelter	475	// bitwise masking in case we need it for NPT
430			calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
431	gezelter	468	}
432			if( !(time % status_n) ){
433			e_out->writeStat( time * dt );
434			calcPot = 0;
435	gezelter	475	// bitwise masking in case we need it for NPT
436			calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
437	gezelter	468	}
438	mmeineke	377	}
439			}
440			else{
441
442			for( tl=0; tl<n_loops; tl++ ){
443
444			kE = 0.0;
445			rot_kE= 0.0;
446			trans_kE = 0.0;
447	gezelter	475
448			if (!strcasecmp( entry_plug->ensemble, "NVT"))
449			myES->NoseHooverNVT( dt , tStats->getKinetic() );
450	mmeineke	377
451			for( i=0; i<nAtoms; i++ ){
452
453			// velocity half step
454
455			vx = atoms[i]->get_vx() +
456			( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
457			vy = atoms[i]->get_vy() +
458			( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
459			vz = atoms[i]->get_vz() +
460			( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
461
462			// position whole step
463
464			rx = atoms[i]->getX() + dt * vx;
465			ry = atoms[i]->getY() + dt * vy;
466			rz = atoms[i]->getZ() + dt * vz;
467
468			atoms[i]->setX( rx );
469			atoms[i]->setY( ry );
470			atoms[i]->setZ( rz );
471
472			atoms[i]->set_vx( vx );
473			atoms[i]->set_vy( vy );
474			atoms[i]->set_vz( vz );
475
476	mmeineke	472	// if( atoms[i]->isDirectional() ){
477	mmeineke	377
478	mmeineke	472	// dAtom = (DirectionalAtom *)atoms[i];
479	mmeineke	377
480	mmeineke	472	// // get and convert the torque to body frame
481	mmeineke	377
482	mmeineke	472	// Tb[0] = dAtom->getTx();
483			// Tb[1] = dAtom->getTy();
484			// Tb[2] = dAtom->getTz();
485	mmeineke	377
486	mmeineke	472	// dAtom->lab2Body( Tb );
487	mmeineke	377
488	mmeineke	472	// // get the angular momentum, and propagate a half step
489	mmeineke	377
490	mmeineke	472	// ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
491			// ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
492			// ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
493	mmeineke	377
494	mmeineke	472	// // get the atom's rotation matrix
495	mmeineke	377
496	mmeineke	472	// A[0][0] = dAtom->getAxx();
497			// A[0][1] = dAtom->getAxy();
498			// A[0][2] = dAtom->getAxz();
499	mmeineke	377
500	mmeineke	472	// A[1][0] = dAtom->getAyx();
501			// A[1][1] = dAtom->getAyy();
502			// A[1][2] = dAtom->getAyz();
503	mmeineke	377
504	mmeineke	472	// A[2][0] = dAtom->getAzx();
505			// A[2][1] = dAtom->getAzy();
506			// A[2][2] = dAtom->getAzz();
507	mmeineke	377
508
509	mmeineke	472	// // use the angular velocities to propagate the rotation matrix a
510			// // full time step
511	mmeineke	377
512
513	mmeineke	472	// angle = dt2 * ji[0] / dAtom->getIxx();
514			// this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
515	mmeineke	377
516	mmeineke	472	// angle = dt2 * ji[1] / dAtom->getIyy();
517			// this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
518	mmeineke	377
519	mmeineke	472	// angle = dt * ji[2] / dAtom->getIzz();
520			// this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
521	mmeineke	377
522	mmeineke	472	// angle = dt2 * ji[1] / dAtom->getIyy();
523			// this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
524	mmeineke	377
525	mmeineke	472	// angle = dt2 * ji[0] / dAtom->getIxx();
526			// this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
527	mmeineke	377
528
529	mmeineke	472	// dAtom->setA( A );
530			// dAtom->setJx( ji[0] );
531			// dAtom->setJy( ji[1] );
532			// dAtom->setJz( ji[2] );
533			// }
534	mmeineke	377	}
535
536			// calculate the forces
537
538	gezelter	468	myFF->doForces(calcPot,calcStress);
539	mmeineke	377
540			for( i=0; i< nAtoms; i++ ){
541
542			// complete the velocity half step
543
544			vx = atoms[i]->get_vx() +
545			( dt2 * atoms[i]->getFx() / atoms[i]->getMass() ) * e_convert;
546			vy = atoms[i]->get_vy() +
547			( dt2 * atoms[i]->getFy() / atoms[i]->getMass() ) * e_convert;
548			vz = atoms[i]->get_vz() +
549			( dt2 * atoms[i]->getFz() / atoms[i]->getMass() ) * e_convert;
550
551			atoms[i]->set_vx( vx );
552			atoms[i]->set_vy( vy );
553			atoms[i]->set_vz( vz );
554
555			// vx2 = vx * vx;
556			// vy2 = vy * vy;
557			// vz2 = vz * vz;
558
559	mmeineke	472	// if( atoms[i]->isDirectional() ){
560	mmeineke	377
561	mmeineke	472	// dAtom = (DirectionalAtom *)atoms[i];
562	mmeineke	377
563	mmeineke	472	// // get and convert the torque to body frame
564	mmeineke	377
565	mmeineke	472	// Tb[0] = dAtom->getTx();
566			// Tb[1] = dAtom->getTy();
567			// Tb[2] = dAtom->getTz();
568	mmeineke	377
569	mmeineke	472	// dAtom->lab2Body( Tb );
570	mmeineke	377
571	mmeineke	472	// // get the angular momentum, and complete the angular momentum
572			// // half step
573	mmeineke	377
574	mmeineke	472	// ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
575			// ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
576			// ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
577	mmeineke	377
578	mmeineke	472	// jx2 = ji[0] * ji[0];
579			// jy2 = ji[1] * ji[1];
580			// jz2 = ji[2] * ji[2];
581	mmeineke	377
582	mmeineke	472	// rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
583			// + (jz2 / dAtom->getIzz());
584	mmeineke	377
585	mmeineke	472	// dAtom->setJx( ji[0] );
586			// dAtom->setJy( ji[1] );
587			// dAtom->setJz( ji[2] );
588			// }
589	mmeineke	377	}
590	gezelter	475
591			if (!strcasecmp( entry_plug->ensemble, "NVT"))
592			myES->NoseHooverNVT( dt , tStats->getKinetic() );
593
594			if (!strcasecmp( entry_plug->ensemble, "NPT") )
595			myES->NoseHooverAndersonNPT( dt,
596			tStats->getKinetic(),
597			tStats->getPressure());
598
599	mmeineke	377	time = tl + 1;
600
601			if( entry_plug->setTemp ){
602			if( !(time % vel_n) ) tStats->velocitize();
603			}
604			if( !(time % sample_n) ) dump_out->writeDump( time * dt );
605	gezelter	468	if( !((time+1) % status_n) ) {
606			calcPot = 1;
607	gezelter	475	// bitwise masking in case we need it for NPT
608			calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
609	gezelter	468	}
610			if( !(time % status_n) ){
611			e_out->writeStat( time * dt );
612			calcPot = 0;
613	gezelter	475	// bitwise masking in case we need it for NPT
614			calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
615	gezelter	468	}
616	mmeineke	377	}
617			}
618
619			dump_out->writeFinal();
620
621			delete dump_out;
622			delete e_out;
623			}
624
625			void Symplectic::rotate( int axes1, int axes2, double angle, double ji[3],
626			double A[3][3] ){
627
628			int i,j,k;
629			double sinAngle;
630			double cosAngle;
631			double angleSqr;
632			double angleSqrOver4;
633			double top, bottom;
634			double rot[3][3];
635			double tempA[3][3];
636			double tempJ[3];
637
638			// initialize the tempA
639
640			for(i=0; i<3; i++){
641			for(j=0; j<3; j++){
642	mmeineke	443	tempA[j][i] = A[i][j];
643	mmeineke	377	}
644			}
645
646			// initialize the tempJ
647
648			for( i=0; i<3; i++) tempJ[i] = ji[i];
649
650			// initalize rot as a unit matrix
651
652			rot[0][0] = 1.0;
653			rot[0][1] = 0.0;
654			rot[0][2] = 0.0;
655
656			rot[1][0] = 0.0;
657			rot[1][1] = 1.0;
658			rot[1][2] = 0.0;
659
660			rot[2][0] = 0.0;
661			rot[2][1] = 0.0;
662			rot[2][2] = 1.0;
663
664			// use a small angle aproximation for sin and cosine
665
666			angleSqr = angle * angle;
667			angleSqrOver4 = angleSqr / 4.0;
668			top = 1.0 - angleSqrOver4;
669			bottom = 1.0 + angleSqrOver4;
670
671			cosAngle = top / bottom;
672			sinAngle = angle / bottom;
673
674			rot[axes1][axes1] = cosAngle;
675			rot[axes2][axes2] = cosAngle;
676
677			rot[axes1][axes2] = sinAngle;
678			rot[axes2][axes1] = -sinAngle;
679
680			// rotate the momentum acoording to: ji[] = rot[][] * ji[]
681
682			for(i=0; i<3; i++){
683			ji[i] = 0.0;
684			for(k=0; k<3; k++){
685			ji[i] += rot[i][k] * tempJ[k];
686			}
687			}
688
689			// rotate the Rotation matrix acording to:
690			// A[][] = A[][] * transpose(rot[][])
691
692
693			// NOte for as yet unknown reason, we are setting the performing the
694			// calculation as:
695			// transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
696
697			for(i=0; i<3; i++){
698			for(j=0; j<3; j++){
699			A[j][i] = 0.0;
700			for(k=0; k<3; k++){
701	mmeineke	443	A[j][i] += tempA[i][k] * rot[j][k];
702	mmeineke	377	}
703			}
704			}
705			}