mdtools/md_code/Symplectic.cpp

#include <iostream>
#include <cstdlib>

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

  srInteractions = entry_plug->sr_interactions;
  longRange      =       entry_plug->longRange;
  nSRI           =           entry_plug->n_SRI;

  // 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[nSRI];
  n_constrained = 0;
  int constrained = 0;
  
  for(int i = 0; i < nSRI; i++){
    
    constrained = srInteractions[i]->is_constrained();
    
    if(constrained){
      
      dummy_plug = srInteractions[i]->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 );

  Thermo *tStats = new Thermo( entry_plug );

  StatWriter*  e_out    = new StatWriter( entry_plug );
  DumpWriter*  dump_out = new DumpWriter( entry_plug );

  strcpy( checkPointMsg, "Fast foward!" );
  MPIcheckPoint();

  strcpy( checkPointMsg, "Fast fowarding!" );
  MPIcheckPoint();

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

  // initialize the forces the before the first step

  
  for(i = 0; i < nAtoms; i++){
    atoms[i]->zeroForces();
  }
  
  for(i = 0; i < nSRI; i++){
    
    srInteractions[i]->calc_forces();
  }
  
  longRange->calc_forces();

  if( entry_plug->setTemp ){
    
    tStats->velocitize();
  }
  
  dump_out->writeDump( 0.0 );
  e_out->writeStat( 0.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
      
      for(j = 0; j < nAtoms; j++){
        atoms[j]->zeroForces();
      }
      
      for(j = 0; j < nSRI; j++){
        srInteractions[j]->calc_forces();
      }
      
      longRange->calc_forces();
      
      // 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 % status_n) ) e_out->writeStat( time * dt );
    }
  }
  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
      
      for(j = 0; j < nAtoms; j++){
        atoms[j]->zeroForces();
      }
      
      for(j = 0; j < nSRI; j++){
        srInteractions[j]->calc_forces();
      }
      
      longRange->calc_forces();
      
      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 % status_n) ) e_out->writeStat( time * dt );
    }
  }

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