OOPSE/libmdtools/Symplectic.cpp

#include <iostream>
#include <cstdlib>

#include "Integrator.hpp"
#include "Thermo.hpp"
#include "ReadWrite.hpp"
#include "ForceFields.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 ){
  entry_plug = the_entry_plug;
  myFF = the_ff;
  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

  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;

   Thermo *tStats = new Thermo( entry_plug );

  StatWriter*  e_out    = new StatWriter( entry_plug );
  DumpWriter*  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,0);
  
  if( entry_plug->setTemp ){
    
    tStats->velocitize();
  }
  
  dump_out->writeDump( 0.0 );
  e_out->writeStat( 0.0 );
  
  calcPot = 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++ ){
      
      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, 0);
      
      // 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] );
        }
      }
     
      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;
      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
    }
  }
  else{

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