OOPSE/libmdtools/Symplectic.cpp

#include <iostream>
#include <cstdlib>

#ifdef IS_MPI
#include "mpiSimulation.hpp"
#include <unistd.h>
#endif //is_mpi

#include "Integrator.hpp"
#include "simError.h"


Symplectic::Symplectic( SimInfo* theInfo, ForceFields* the_ff ){
  
  info = theInfo;
  myFF = the_ff;
  isFirst = 1;

  molecules = info->molecules;
  nMols = info->n_mol;

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

  // check for constraints
  
  constrainedI = NULL;
  constrainedJ = NULL;
  constrainedDsqr = NULL;
  nConstrained = 0;

  checkConstraints();
}

Symplectic::~Symplectic() {
  
  if( nConstrained ){
    delete[] constrainedI;
    delete[] constrainedJ;
    delete[] constrainedDsqr;
  }
  
}

void Symplectic::checkConstraints( void ){


  isConstrained = 0;

  Constraint *temp_con;
  Constraint *dummy_plug;
  temp_con = new Constraint[info->n_SRI];
  nConstrained = 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[nConstrained].set_a( dummy_plug->get_a() );
        temp_con[nConstrained].set_b( dummy_plug->get_b() );
        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
        
        nConstrained++;
        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[nConstrained].set_a( dummy_plug->get_a() );
        temp_con[nConstrained].set_b( dummy_plug->get_b() );
        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
        
        nConstrained++;
        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[nConstrained].set_a( dummy_plug->get_a() );
        temp_con[nConstrained].set_b( dummy_plug->get_b() );
        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
        
        nConstrained++;
        constrained = 0;
      }
    }
  }

  if(nConstrained > 0){
    
    isConstrained = 1;

    if(constrainedI != NULL )    delete[] constrainedI;
    if(constrainedJ != NULL )    delete[] constrainedJ;
    if(constrainedDsqr != NULL ) delete[] constrainedDsqr;

    constrainedI =    new int[nConstrained];
    constrainedJ =    new int[nConstrained];
    constrainedDsqr = new double[nConstrained];
    
    for( int i = 0; i < nConstrained; i++){
      
      constrainedI[i] = temp_con[i].get_a();
      constrainedJ[i] = temp_con[i].get_b();
      constrainedDsqr[i] = temp_con[i].get_dsqr();
    }
  }
  
  delete[] temp_con;
}


void Symplectic::integrate( void ){

  int i, j;                         // loop counters
  int nAtoms = info->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          = info->dt;
  double runTime     = info->run_time;
  double sampleTime  = info->sampleTime;
  double statusTime  = info->statusTime;
  double thermalTime = info->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;
  int isError;

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

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

  // initialize the forces before the first step

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

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

    integrateStep( calcPot, calcStress );
      
    time = tl + 1;
    
    if( info->setTemp ){
      if( !(time % vel_n) ) tStats->velocitize();
    }
    if( !(time % sample_n) ) dump_out->writeDump( time * dt );
    if( !((time+1) % status_n) ) {
      calcPot = 1;
      calcStress = 1;
    }
    if( !(time % status_n) ){ 
      e_out->writeStat( time * dt ); 
      calcPot = 0; 
      if (!strcasecmp(info->ensemble, "NPT")) calcStress = 1;
      else calcStress = 0;
    }    

   
  }

  dump_out->writeFinal();

  delete dump_out;
  delete e_out;
}


void Symplectic::moveA( void ){
  
  int i,j,k;
  int atomIndex, aMatIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];

  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;
    aMatIndex = i * 9;
    
    // velocity half step
    for( j=atomIndex; j<(atomIndex+3); j++ )
      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;

    // position whole step    
    for( j=atomIndex; j<(atomIndex+3); j++ )
      pos[j] += dt * vel[j];

   
    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] ) * eConvert;
      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
      
      // use the angular velocities to propagate the rotation matrix a
      // full time step
      
      // rotate about the x-axis      
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] ); 
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
      
      // rotate about the z-axis
      angle = dt * ji[2] / dAtom->getIzz();
      this->rotate( 0, 1, angle, ji, &aMat[aMatIndex] );
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
      
       // rotate about the x-axis
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
      
      dAtom->setJx( ji[0] );
      dAtom->setJy( ji[1] );
      dAtom->setJz( ji[2] );
    }
    
  }
}


void Integrator::moveB( void ){
  int i,j,k;
  int atomIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];

  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;

    // velocity half step
    for( j=atomIndex; j<(atomIndex+3); j++ )
      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;

    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] ) * eConvert;
      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
      
      jx2 = ji[0] * ji[0];
      jy2 = ji[1] * ji[1];
      jz2 = ji[2] * ji[2];
      
      dAtom->setJx( ji[0] );
      dAtom->setJy( ji[1] );
      dAtom->setJz( ji[2] );
    }
  }

}


void Integrator::constrainA(){

  
}


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:	542
Committed:	Fri May 30 21:31:48 2003 UTC (21 years, 11 months ago) by mmeineke
File size:	9903 byte(s)
Log Message:	currently modifiying Symplectic to become the basic integrator.
#	User	Rev	Content
1	mmeineke	377	#include <iostream>
2			#include <cstdlib>
3
4	mmeineke	542	#ifdef IS_MPI
5			#include "mpiSimulation.hpp"
6			#include <unistd.h>
7			#endif //is_mpi
8
9	mmeineke	377	#include "Integrator.hpp"
10			#include "simError.h"
11
12	mmeineke	542
13			Symplectic::Symplectic( SimInfo* theInfo, ForceFields* the_ff ){
14	mmeineke	377
15	mmeineke	542	info = theInfo;
16	mmeineke	377	myFF = the_ff;
17			isFirst = 1;
18
19	mmeineke	542	molecules = info->molecules;
20			nMols = info->n_mol;
21	mmeineke	377
22			// give a little love back to the SimInfo object
23
24	mmeineke	542	if( info->the_integrator != NULL ) delete info->the_integrator;
25			info->the_integrator = this;
26	mmeineke	377
27			// check for constraints
28	mmeineke	542
29			constrainedI = NULL;
30			constrainedJ = NULL;
31			constrainedDsqr = NULL;
32			nConstrained = 0;
33	mmeineke	377
34	mmeineke	542	checkConstraints();
35			}
36	mmeineke	377
37	mmeineke	542	Symplectic::~Symplectic() {
38
39			if( nConstrained ){
40			delete[] constrainedI;
41			delete[] constrainedJ;
42			delete[] constrainedDsqr;
43			}
44
45			}
46
47			void Symplectic::checkConstraints( void ){
48
49
50			isConstrained = 0;
51
52	mmeineke	377	Constraint *temp_con;
53			Constraint *dummy_plug;
54	mmeineke	542	temp_con = new Constraint[info->n_SRI];
55			nConstrained = 0;
56	mmeineke	377	int constrained = 0;
57
58	mmeineke	423	SRI** theArray;
59			for(int i = 0; i < nMols; i++){
60	mmeineke	377
61	mmeineke	428	theArray = (SRI**) molecules[i].getMyBonds();
62			for(int j=0; j<molecules[i].getNBonds(); j++){
63	mmeineke	377
64	mmeineke	423	constrained = theArray[j]->is_constrained();
65
66			if(constrained){
67
68			dummy_plug = theArray[j]->get_constraint();
69	mmeineke	542	temp_con[nConstrained].set_a( dummy_plug->get_a() );
70			temp_con[nConstrained].set_b( dummy_plug->get_b() );
71			temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
72	mmeineke	423
73	mmeineke	542	nConstrained++;
74	mmeineke	423	constrained = 0;
75			}
76			}
77	mmeineke	377
78	mmeineke	428	theArray = (SRI**) molecules[i].getMyBends();
79			for(int j=0; j<molecules[i].getNBends(); j++){
80	mmeineke	423
81			constrained = theArray[j]->is_constrained();
82
83			if(constrained){
84
85			dummy_plug = theArray[j]->get_constraint();
86	mmeineke	542	temp_con[nConstrained].set_a( dummy_plug->get_a() );
87			temp_con[nConstrained].set_b( dummy_plug->get_b() );
88			temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
89	mmeineke	423
90	mmeineke	542	nConstrained++;
91	mmeineke	423	constrained = 0;
92			}
93	mmeineke	377	}
94	mmeineke	423
95	mmeineke	428	theArray = (SRI**) molecules[i].getMyTorsions();
96			for(int j=0; j<molecules[i].getNTorsions(); j++){
97	mmeineke	423
98			constrained = theArray[j]->is_constrained();
99
100			if(constrained){
101
102			dummy_plug = theArray[j]->get_constraint();
103	mmeineke	542	temp_con[nConstrained].set_a( dummy_plug->get_a() );
104			temp_con[nConstrained].set_b( dummy_plug->get_b() );
105			temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
106	mmeineke	423
107	mmeineke	542	nConstrained++;
108	mmeineke	423	constrained = 0;
109			}
110			}
111	mmeineke	377	}
112
113	mmeineke	542	if(nConstrained > 0){
114	mmeineke	377
115	mmeineke	542	isConstrained = 1;
116
117			if(constrainedI != NULL ) delete[] constrainedI;
118			if(constrainedJ != NULL ) delete[] constrainedJ;
119			if(constrainedDsqr != NULL ) delete[] constrainedDsqr;
120
121			constrainedI = new int[nConstrained];
122			constrainedJ = new int[nConstrained];
123			constrainedDsqr = new double[nConstrained];
124	mmeineke	377
125	mmeineke	542	for( int i = 0; i < nConstrained; i++){
126	mmeineke	377
127	mmeineke	542	constrainedI[i] = temp_con[i].get_a();
128			constrainedJ[i] = temp_con[i].get_b();
129			constrainedDsqr[i] = temp_con[i].get_dsqr();
130	mmeineke	377	}
131			}
132
133			delete[] temp_con;
134			}
135
136
137			void Symplectic::integrate( void ){
138
139			int i, j; // loop counters
140	mmeineke	542	int nAtoms = info->n_atoms; // the number of atoms
141	mmeineke	377	double kE = 0.0; // the kinetic energy
142			double rot_kE;
143			double trans_kE;
144			int tl; // the time loop conter
145			double dt2; // half the dt
146
147			double vx, vy, vz; // the velocities
148	chuckv	482	double vx2, vy2, vz2; // the square of the velocities
149	mmeineke	377	double rx, ry, rz; // the postitions
150
151			double ji[3]; // the body frame angular momentum
152			double jx2, jy2, jz2; // the square of the angular momentums
153			double Tb[3]; // torque in the body frame
154			double angle; // the angle through which to rotate the rotation matrix
155			double A[3][3]; // the rotation matrix
156	gezelter	483	double press[9];
157	mmeineke	377
158			int time;
159
160	mmeineke	542	double dt = info->dt;
161			double runTime = info->run_time;
162			double sampleTime = info->sampleTime;
163			double statusTime = info->statusTime;
164			double thermalTime = info->thermalTime;
165	mmeineke	377
166			int n_loops = (int)( runTime / dt );
167			int sample_n = (int)( sampleTime / dt );
168			int status_n = (int)( statusTime / dt );
169			int vel_n = (int)( thermalTime / dt );
170
171	gezelter	468	int calcPot, calcStress;
172	mmeineke	542	int isError;
173	mmeineke	377
174	mmeineke	542	tStats = new Thermo( info );
175			e_out = new StatWriter( info );
176			dump_out = new DumpWriter( info );
177	mmeineke	377
178	mmeineke	542	Atom** atoms = info->atoms;
179	mmeineke	377	DirectionalAtom* dAtom;
180			dt2 = 0.5 * dt;
181
182	mmeineke	542	// initialize the forces before the first step
183	mmeineke	377
184	gezelter	468	myFF->doForces(1,1);
185	mmeineke	377
186	mmeineke	542	if( info->setTemp ){
187	mmeineke	377
188			tStats->velocitize();
189			}
190
191			dump_out->writeDump( 0.0 );
192			e_out->writeStat( 0.0 );
193
194			calcPot = 0;
195
196	mmeineke	542	for( tl=0; tl<nLoops; tl++){
197	gezelter	475
198	mmeineke	542	integrateStep( calcPot, calcStress );
199
200			time = tl + 1;
201	mmeineke	377
202	mmeineke	542	if( info->setTemp ){
203			if( !(time % vel_n) ) tStats->velocitize();
204			}
205			if( !(time % sample_n) ) dump_out->writeDump( time * dt );
206			if( !((time+1) % status_n) ) {
207			calcPot = 1;
208			calcStress = 1;
209			}
210			if( !(time % status_n) ){
211			e_out->writeStat( time * dt );
212			calcPot = 0;
213			if (!strcasecmp(info->ensemble, "NPT")) calcStress = 1;
214			else calcStress = 0;
215			}
216	mmeineke	377
217	mmeineke	542
218			}
219	gezelter	475
220	mmeineke	542	dump_out->writeFinal();
221	mmeineke	377
222	mmeineke	542	delete dump_out;
223			delete e_out;
224			}
225	mmeineke	377
226
227	mmeineke	542	void Symplectic::moveA( void ){
228
229			int i,j,k;
230			int atomIndex, aMatIndex;
231			DirectionalAtom* dAtom;
232			double Tb[3];
233			double ji[3];
234	mmeineke	377
235	mmeineke	542	for( i=0; i<nAtoms; i++ ){
236			atomIndex = i * 3;
237			aMatIndex = i * 9;
238
239			// velocity half step
240			for( j=atomIndex; j<(atomIndex+3); j++ )
241			vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
242	mmeineke	377
243	mmeineke	542	// position whole step
244			for( j=atomIndex; j<(atomIndex+3); j++ )
245			pos[j] += dt * vel[j];
246	mmeineke	377
247	mmeineke	542
248			if( atoms[i]->isDirectional() ){
249	mmeineke	377
250	mmeineke	542	dAtom = (DirectionalAtom *)atoms[i];
251	mmeineke	377
252	mmeineke	542	// get and convert the torque to body frame
253	mmeineke	377
254	mmeineke	542	Tb[0] = dAtom->getTx();
255			Tb[1] = dAtom->getTy();
256			Tb[2] = dAtom->getTz();
257	mmeineke	377
258	mmeineke	542	dAtom->lab2Body( Tb );
259	mmeineke	377
260	mmeineke	542	// get the angular momentum, and propagate a half step
261	mmeineke	377
262	mmeineke	542	ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
263			ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
264			ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
265	mmeineke	377
266	mmeineke	542	// use the angular velocities to propagate the rotation matrix a
267			// full time step
268	mmeineke	377
269	mmeineke	542	// rotate about the x-axis
270			angle = dt2 * ji[0] / dAtom->getIxx();
271			this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
272
273			// rotate about the y-axis
274			angle = dt2 * ji[1] / dAtom->getIyy();
275			this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
276
277			// rotate about the z-axis
278			angle = dt * ji[2] / dAtom->getIzz();
279			this->rotate( 0, 1, angle, ji, &aMat[aMatIndex] );
280
281			// rotate about the y-axis
282			angle = dt2 * ji[1] / dAtom->getIyy();
283			this->rotate( 2, 0, angle, ji, &aMat[aMatIndex] );
284
285			// rotate about the x-axis
286			angle = dt2 * ji[0] / dAtom->getIxx();
287			this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] );
288
289			dAtom->setJx( ji[0] );
290			dAtom->setJy( ji[1] );
291			dAtom->setJz( ji[2] );
292	mmeineke	377	}
293	mmeineke	542
294	mmeineke	377	}
295	mmeineke	542	}
296	mmeineke	377
297	gezelter	475
298	mmeineke	542	void Integrator::moveB( void ){
299			int i,j,k;
300			int atomIndex;
301			DirectionalAtom* dAtom;
302			double Tb[3];
303			double ji[3];
304	mmeineke	377
305	mmeineke	542	for( i=0; i<nAtoms; i++ ){
306			atomIndex = i * 3;
307	chuckv	497
308	mmeineke	542	// velocity half step
309			for( j=atomIndex; j<(atomIndex+3); j++ )
310			vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
311
312			if( atoms[i]->isDirectional() ){
313	mmeineke	377
314	mmeineke	542	dAtom = (DirectionalAtom *)atoms[i];
315	chuckv	497
316	mmeineke	542	// get and convert the torque to body frame
317	mmeineke	377
318	mmeineke	542	Tb[0] = dAtom->getTx();
319			Tb[1] = dAtom->getTy();
320			Tb[2] = dAtom->getTz();
321	mmeineke	377
322	mmeineke	542	dAtom->lab2Body( Tb );
323
324			// get the angular momentum, and complete the angular momentum
325			// half step
326
327			ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
328			ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
329			ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
330
331			jx2 = ji[0] * ji[0];
332			jy2 = ji[1] * ji[1];
333			jz2 = ji[2] * ji[2];
334
335			dAtom->setJx( ji[0] );
336			dAtom->setJy( ji[1] );
337			dAtom->setJz( ji[2] );
338			}
339			}
340	mmeineke	377
341	mmeineke	542	}
342	mmeineke	486
343	gezelter	475
344	mmeineke	542	void Integrator::constrainA(){
345
346	gezelter	475
347	mmeineke	377
348
349			}
350
351	mmeineke	542
352
353
354
355
356
357
358
359	mmeineke	377	void Symplectic::rotate( int axes1, int axes2, double angle, double ji[3],
360			double A[3][3] ){
361
362			int i,j,k;
363			double sinAngle;
364			double cosAngle;
365			double angleSqr;
366			double angleSqrOver4;
367			double top, bottom;
368			double rot[3][3];
369			double tempA[3][3];
370			double tempJ[3];
371
372			// initialize the tempA
373
374			for(i=0; i<3; i++){
375			for(j=0; j<3; j++){
376	mmeineke	443	tempA[j][i] = A[i][j];
377	mmeineke	377	}
378			}
379
380			// initialize the tempJ
381
382			for( i=0; i<3; i++) tempJ[i] = ji[i];
383
384			// initalize rot as a unit matrix
385
386			rot[0][0] = 1.0;
387			rot[0][1] = 0.0;
388			rot[0][2] = 0.0;
389
390			rot[1][0] = 0.0;
391			rot[1][1] = 1.0;
392			rot[1][2] = 0.0;
393
394			rot[2][0] = 0.0;
395			rot[2][1] = 0.0;
396			rot[2][2] = 1.0;
397
398			// use a small angle aproximation for sin and cosine
399
400			angleSqr = angle * angle;
401			angleSqrOver4 = angleSqr / 4.0;
402			top = 1.0 - angleSqrOver4;
403			bottom = 1.0 + angleSqrOver4;
404
405			cosAngle = top / bottom;
406			sinAngle = angle / bottom;
407
408			rot[axes1][axes1] = cosAngle;
409			rot[axes2][axes2] = cosAngle;
410
411			rot[axes1][axes2] = sinAngle;
412			rot[axes2][axes1] = -sinAngle;
413
414			// rotate the momentum acoording to: ji[] = rot[][] * ji[]
415
416			for(i=0; i<3; i++){
417			ji[i] = 0.0;
418			for(k=0; k<3; k++){
419			ji[i] += rot[i][k] * tempJ[k];
420			}
421			}
422
423			// rotate the Rotation matrix acording to:
424			// A[][] = A[][] * transpose(rot[][])
425
426
427			// NOte for as yet unknown reason, we are setting the performing the
428			// calculation as:
429			// transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
430
431			for(i=0; i<3; i++){
432			for(j=0; j<3; j++){
433			A[j][i] = 0.0;
434			for(k=0; k<3; k++){
435	mmeineke	443	A[j][i] += tempA[i][k] * rot[j][k];
436	mmeineke	377	}
437			}
438			}
439			}