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;

  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
  double press[9];

  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;

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

  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 / 2.0 , 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 / 2.0, tStats->getKinetic() );

      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
        tStats->getPressureTensor(press);
        myES->NoseHooverAndersonNPT( dt, 
                                     tStats->getKinetic(), 
                                     press);
      }

      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 / 2.0, 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 / 2.0, tStats->getKinetic() );

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