ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Integrator.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Integrator.cpp (file contents):
Revision 597 by mmeineke, Mon Jul 14 21:28:54 2003 UTC vs.
Revision 1284 by tim, Mon Jun 21 18:52:21 2004 UTC

# Line 1 | Line 1
1   #include <iostream>
2 < #include <cstdlib>
3 < #include <cmath>
4 <
2 > #include <stdlib.h>
3 > #include <math.h>
4 > #include "Rattle.hpp"
5 > #include "Roll.hpp"
6   #ifdef IS_MPI
7   #include "mpiSimulation.hpp"
8   #include <unistd.h>
9   #endif //is_mpi
10  
11 + #ifdef PROFILE
12 + #include "mdProfile.hpp"
13 + #endif // profile
14 +
15   #include "Integrator.hpp"
16   #include "simError.h"
17  
18  
19 < Integrator::Integrator( SimInfo *theInfo, ForceFields* the_ff ){
20 <  
19 > template<typename T> Integrator<T>::Integrator(SimInfo* theInfo,
20 >                                               ForceFields* the_ff){
21    info = theInfo;
22    myFF = the_ff;
23    isFirst = 1;
# Line 21 | Line 26 | Integrator::Integrator( SimInfo *theInfo, ForceFields*
26    nMols = info->n_mol;
27  
28    // give a little love back to the SimInfo object
24 –  
25 –  if( info->the_integrator != NULL ) delete info->the_integrator;
26 –  info->the_integrator = this;
29  
30 +  if (info->the_integrator != NULL){
31 +    delete info->the_integrator;
32 +  }
33 +
34    nAtoms = info->n_atoms;
35 +  integrableObjects = info->integrableObjects;
36  
37 <  std::cerr << "integ nAtoms = "  << nAtoms << "\n";
37 >  consFramework = new RattleFramework(info);
38  
39 <  // check for constraints
39 >  if(consFramework == NULL){
40 >    sprintf(painCave.errMsg,
41 >      "Integrator::Intergrator() Error: Memory allocation error for RattleFramework" );
42 >    painCave.isFatal = 1;
43 >    simError();
44 >  }
45    
46 <  constrainedA    = NULL;
47 <  constrainedB    = NULL;
46 > /*
47 >  // check for constraints
48 >
49 >  constrainedA = NULL;
50 >  constrainedB = NULL;
51    constrainedDsqr = NULL;
52 <  moving          = NULL;
53 <  moved           = NULL;
54 <  oldPos          = NULL;
55 <  
52 >  moving = NULL;
53 >  moved = NULL;
54 >  oldPos = NULL;
55 >
56    nConstrained = 0;
57  
58    checkConstraints();
59 + */
60   }
61  
62 < Integrator::~Integrator() {
63 <  
64 <  if( nConstrained ){
62 > template<typename T> Integrator<T>::~Integrator(){
63 >  if (consFramework != NULL)
64 >    delete consFramework;
65 > /*
66 >  if (nConstrained){
67      delete[] constrainedA;
68      delete[] constrainedB;
69      delete[] constrainedDsqr;
# Line 53 | Line 71 | Integrator::~Integrator() {
71      delete[] moved;
72      delete[] oldPos;
73    }
74 <  
74 > */
75   }
76  
77 < void Integrator::checkConstraints( void ){
78 <
61 <
77 > /*
78 > template<typename T> void Integrator<T>::checkConstraints(void){
79    isConstrained = 0;
80  
81 <  Constraint *temp_con;
82 <  Constraint *dummy_plug;
81 >  Constraint* temp_con;
82 >  Constraint* dummy_plug;
83    temp_con = new Constraint[info->n_SRI];
84    nConstrained = 0;
85    int constrained = 0;
86 <  
86 >
87    SRI** theArray;
88 <  for(int i = 0; i < nMols; i++){
89 <    
90 <    theArray = (SRI**) molecules[i].getMyBonds();
91 <    for(int j=0; j<molecules[i].getNBonds(); j++){
75 <      
88 >  for (int i = 0; i < nMols; i++){
89 >
90 >          theArray = (SRI * *) molecules[i].getMyBonds();
91 >    for (int j = 0; j < molecules[i].getNBonds(); j++){
92        constrained = theArray[j]->is_constrained();
93  
94 <      std::cerr << "Is the folowing bond constrained \n";
95 <      theArray[j]->printMe();
96 <      
97 <      if(constrained){
98 <        
83 <        std::cerr << "Yes\n";
94 >      if (constrained){
95 >        dummy_plug = theArray[j]->get_constraint();
96 >        temp_con[nConstrained].set_a(dummy_plug->get_a());
97 >        temp_con[nConstrained].set_b(dummy_plug->get_b());
98 >        temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());
99  
100 <        dummy_plug = theArray[j]->get_constraint();
101 <        temp_con[nConstrained].set_a( dummy_plug->get_a() );
102 <        temp_con[nConstrained].set_b( dummy_plug->get_b() );
88 <        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
89 <        
90 <        nConstrained++;
91 <        constrained = 0;
92 <      }
93 <      else std::cerr << "No.\n";
100 >        nConstrained++;
101 >        constrained = 0;
102 >      }
103      }
104  
105 <    theArray = (SRI**) molecules[i].getMyBends();
106 <    for(int j=0; j<molecules[i].getNBends(); j++){
98 <      
105 >    theArray = (SRI * *) molecules[i].getMyBends();
106 >    for (int j = 0; j < molecules[i].getNBends(); j++){
107        constrained = theArray[j]->is_constrained();
108 <      
109 <      if(constrained){
110 <        
111 <        dummy_plug = theArray[j]->get_constraint();
112 <        temp_con[nConstrained].set_a( dummy_plug->get_a() );
113 <        temp_con[nConstrained].set_b( dummy_plug->get_b() );
114 <        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
115 <        
116 <        nConstrained++;
109 <        constrained = 0;
108 >
109 >      if (constrained){
110 >        dummy_plug = theArray[j]->get_constraint();
111 >        temp_con[nConstrained].set_a(dummy_plug->get_a());
112 >        temp_con[nConstrained].set_b(Dummy_plug->get_b());
113 >        temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());
114 >
115 >        nConstrained++;
116 >        constrained = 0;
117        }
118      }
119  
120 <    theArray = (SRI**) molecules[i].getMyTorsions();
121 <    for(int j=0; j<molecules[i].getNTorsions(); j++){
115 <      
120 >    theArray = (SRI * *) molecules[i].getMyTorsions();
121 >    for (int j = 0; j < molecules[i].getNTorsions(); j++){
122        constrained = theArray[j]->is_constrained();
123 <      
124 <      if(constrained){
125 <        
126 <        dummy_plug = theArray[j]->get_constraint();
127 <        temp_con[nConstrained].set_a( dummy_plug->get_a() );
128 <        temp_con[nConstrained].set_b( dummy_plug->get_b() );
129 <        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
130 <        
131 <        nConstrained++;
126 <        constrained = 0;
123 >
124 >      if (constrained){
125 >        dummy_plug = theArray[j]->get_constraint();
126 >        temp_con[nConstrained].set_a(dummy_plug->get_a());
127 >        temp_con[nConstrained].set_b(dummy_plug->get_b());
128 >        temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());
129 >
130 >        nConstrained++;
131 >        constrained = 0;
132        }
133      }
134    }
135  
136 <  if(nConstrained > 0){
137 <    
136 >
137 >  if (nConstrained > 0){
138      isConstrained = 1;
139  
140 <    if(constrainedA != NULL )    delete[] constrainedA;
141 <    if(constrainedB != NULL )    delete[] constrainedB;
142 <    if(constrainedDsqr != NULL ) delete[] constrainedDsqr;
140 >    if (constrainedA != NULL)
141 >      delete[] constrainedA;
142 >    if (constrainedB != NULL)
143 >      delete[] constrainedB;
144 >    if (constrainedDsqr != NULL)
145 >      delete[] constrainedDsqr;
146  
147 <    constrainedA =    new int[nConstrained];
148 <    constrainedB =    new int[nConstrained];
147 >    constrainedA = new int[nConstrained];
148 >    constrainedB = new int[nConstrained];
149      constrainedDsqr = new double[nConstrained];
150 <    
151 <    for( int i = 0; i < nConstrained; i++){
144 <      
150 >
151 >    for (int i = 0; i < nConstrained; i++){
152        constrainedA[i] = temp_con[i].get_a();
153        constrainedB[i] = temp_con[i].get_b();
154        constrainedDsqr[i] = temp_con[i].get_dsqr();
148 –
155      }
156  
157 <    
158 <    // save oldAtoms to check for lode balanceing later on.
159 <    
157 >
158 >    // save oldAtoms to check for lode balancing later on.
159 >
160      oldAtoms = nAtoms;
161 <    
161 >
162      moving = new int[nAtoms];
163 <    moved  = new int[nAtoms];
163 >    moved = new int[nAtoms];
164  
165 <    oldPos = new double[nAtoms*3];
165 >    oldPos = new double[nAtoms * 3];
166    }
167 <  
167 >
168    delete[] temp_con;
169   }
170 + */
171  
172 + template<typename T> void Integrator<T>::integrate(void){
173  
174 < void Integrator::integrate( void ){
175 <
176 <  int i, j;                         // loop counters
169 <
170 <  double runTime     = info->run_time;
171 <  double sampleTime  = info->sampleTime;
172 <  double statusTime  = info->statusTime;
174 >  double runTime = info->run_time;
175 >  double sampleTime = info->sampleTime;
176 >  double statusTime = info->statusTime;
177    double thermalTime = info->thermalTime;
178 +  double resetTime = info->resetTime;
179  
180 +  double difference;
181    double currSample;
182    double currThermal;
183    double currStatus;
184 <  double currTime;
184 >  double currReset;
185  
186    int calcPot, calcStress;
181 –  int isError;
187  
188 +  tStats = new Thermo(info);
189 +  statOut = new StatWriter(info);
190 +  dumpOut = new DumpWriter(info);
191  
184 –
185 –  tStats   = new Thermo( info );
186 –  statOut  = new StatWriter( info );
187 –  dumpOut  = new DumpWriter( info );
188 –
192    atoms = info->atoms;
190 –  DirectionalAtom* dAtom;
193  
194    dt = info->dt;
195    dt2 = 0.5 * dt;
196  
197 +  readyCheck();
198 +
199 +  // remove center of mass drift velocity (in case we passed in a configuration
200 +  // that was drifting
201 +  tStats->removeCOMdrift();
202 +
203 +  // initialize the retraints if necessary
204 +  if (info->useSolidThermInt && !info->useLiquidThermInt) {
205 +    myFF->initRestraints();
206 +  }
207 +
208    // initialize the forces before the first step
209  
210 <  myFF->doForces(1,1);
210 >  calcForce(1, 1);
211 >
212 >  //execute constraint algorithm to make sure at the very beginning the system is constrained  
213 >  //consFramework->doPreConstraint();
214 >  //consFramework->doConstrainA();
215 >  //calcForce(1, 1);
216 >  //consFramework->doConstrainB();
217    
218 <  if( info->setTemp ){
219 <    
201 <    tStats->velocitize();
218 >  if (info->setTemp){
219 >    thermalize();
220    }
221 <  
204 <  dumpOut->writeDump( 0.0 );
205 <  statOut->writeStat( 0.0 );
206 <  
221 >
222    calcPot     = 0;
223    calcStress  = 0;
224 <  currSample  = sampleTime;
225 <  currThermal = thermalTime;
226 <  currStatus  = statusTime;
227 <  currTime    = 0.0;;
224 >  currSample  = sampleTime + info->getTime();
225 >  currThermal = thermalTime+ info->getTime();
226 >  currStatus  = statusTime + info->getTime();
227 >  currReset   = resetTime  + info->getTime();
228  
229 +  dumpOut->writeDump(info->getTime());
230 +  statOut->writeStat(info->getTime());
231  
215 –  readyCheck();
232  
233   #ifdef IS_MPI
234 <  strcpy( checkPointMsg,
219 <          "The integrator is ready to go." );
234 >  strcpy(checkPointMsg, "The integrator is ready to go.");
235    MPIcheckPoint();
236   #endif // is_mpi
237  
238 <
239 <  pos  = Atom::getPosArray();
240 <  vel  = Atom::getVelArray();
226 <  frc  = Atom::getFrcArray();
227 <
228 <  while( currTime < runTime ){
229 <
230 <    if( (currTime+dt) >= currStatus ){
238 >  while (info->getTime() < runTime && !stopIntegrator()){
239 >    difference = info->getTime() + dt - currStatus;
240 >    if (difference > 0 || fabs(difference) < 1e-4 ){
241        calcPot = 1;
242        calcStress = 1;
243      }
244  
245 <    std::cerr << currTime << "\n";
245 > #ifdef PROFILE
246 >    startProfile( pro1 );
247 > #endif
248 >    
249 >    integrateStep(calcPot, calcStress);
250  
251 <    integrateStep( calcPot, calcStress );
252 <      
239 <    currTime += dt;
251 > #ifdef PROFILE
252 >    endProfile( pro1 );
253  
254 <    if( info->setTemp ){
255 <      if( currTime >= currThermal ){
256 <        tStats->velocitize();
257 <        currThermal += thermalTime;
254 >    startProfile( pro2 );
255 > #endif // profile
256 >
257 >    info->incrTime(dt);
258 >
259 >    if (info->setTemp){
260 >      if (info->getTime() >= currThermal){
261 >        thermalize();
262 >        currThermal += thermalTime;
263        }
264      }
265  
266 <    if( currTime >= currSample ){
267 <      dumpOut->writeDump( currTime );
266 >    if (info->getTime() >= currSample){
267 >      dumpOut->writeDump(info->getTime());
268        currSample += sampleTime;
269      }
270  
271 <    if( currTime >= currStatus ){
272 <      statOut->writeStat( currTime );
273 <      calcPot = 0;
271 >    if (info->getTime() >= currStatus){
272 >      statOut->writeStat(info->getTime());
273 >      calcPot = 0;
274        calcStress = 0;
275        currStatus += statusTime;
276 <    }
276 >    }
277  
278 +    if (info->resetIntegrator){
279 +      if (info->getTime() >= currReset){
280 +        this->resetIntegrator();
281 +        currReset += resetTime;
282 +      }
283 +    }
284 +    
285 + #ifdef PROFILE
286 +    endProfile( pro2 );
287 + #endif //profile
288 +
289   #ifdef IS_MPI
290 <    strcpy( checkPointMsg,
262 <            "successfully took a time step." );
290 >    strcpy(checkPointMsg, "successfully took a time step.");
291      MPIcheckPoint();
292   #endif // is_mpi
265 –
293    }
294  
295 <  dumpOut->writeFinal(currTime);
295 >  // dump out a file containing the omega values for the final configuration
296 >  if (info->useSolidThermInt && !info->useLiquidThermInt)
297 >    myFF->dumpzAngle();
298 >  
299  
300    delete dumpOut;
301    delete statOut;
302   }
303  
304 < void Integrator::integrateStep( int calcPot, int calcStress ){
304 > template<typename T> void Integrator<T>::integrateStep(int calcPot,
305 >                                                       int calcStress){
306 >  // Position full step, and velocity half step
307  
308 + #ifdef PROFILE
309 +  startProfile(pro3);
310 + #endif //profile
311  
312 <      
313 <  // Position full step, and velocity half step
312 >  //save old state (position, velocity etc)
313 >  consFramework->doPreConstraint();
314  
315 <  preMove();
315 > #ifdef PROFILE
316 >  endProfile(pro3);
317 >
318 >  startProfile(pro4);
319 > #endif // profile
320 >
321    moveA();
282 –  //if( nConstrained ) constrainA();
322  
323 + #ifdef PROFILE
324 +  endProfile(pro4);
325 +  
326 +  startProfile(pro5);
327 + #endif//profile
328 +
329 +
330 + #ifdef IS_MPI
331 +  strcpy(checkPointMsg, "Succesful moveA\n");
332 +  MPIcheckPoint();
333 + #endif // is_mpi
334 +
335    // calc forces
336 +  calcForce(calcPot, calcStress);
337  
338 <  myFF->doForces(calcPot,calcStress);
338 > #ifdef IS_MPI
339 >  strcpy(checkPointMsg, "Succesful doForces\n");
340 >  MPIcheckPoint();
341 > #endif // is_mpi
342  
343 <  // finish the velocity  half step
344 <  
290 <  moveB();
291 <  if( nConstrained ) constrainB();
292 <  
293 < }
343 > #ifdef PROFILE
344 >  endProfile( pro5 );
345  
346 +  startProfile( pro6 );
347 + #endif //profile
348  
349 < void Integrator::moveA( void ){
297 <  
298 <  int i,j,k;
299 <  int atomIndex, aMatIndex;
300 <  DirectionalAtom* dAtom;
301 <  double Tb[3];
302 <  double ji[3];
303 <  double angle;
304 <  double A[3][3], At[3][3];
349 >  // finish the velocity  half step
350  
351 +  moveB();
352  
353 <  for( i=0; i<nAtoms; i++ ){
354 <    atomIndex = i * 3;
355 <    aMatIndex = i * 9;
353 > #ifdef PROFILE
354 >  endProfile(pro6);
355 > #endif // profile
356  
357 <    // velocity half step
358 <    for( j=atomIndex; j<(atomIndex+3); j++ )
359 <      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
357 > #ifdef IS_MPI
358 >  strcpy(checkPointMsg, "Succesful moveB\n");
359 >  MPIcheckPoint();
360 > #endif // is_mpi
361 > }
362  
363  
364 <    // position whole step    
365 <    for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j];
364 > template<typename T> void Integrator<T>::moveA(void){
365 >  size_t i, j;
366 >  DirectionalAtom* dAtom;
367 >  double Tb[3], ji[3];
368 >  double vel[3], pos[3], frc[3];
369 >  double mass;
370 >  double omega;
371 >
372 >  for (i = 0; i < integrableObjects.size() ; i++){
373 >    integrableObjects[i]->getVel(vel);
374 >    integrableObjects[i]->getPos(pos);
375 >    integrableObjects[i]->getFrc(frc);
376      
377 +    mass = integrableObjects[i]->getMass();
378  
379 <    if( atoms[i]->isDirectional() ){
379 >    for (j = 0; j < 3; j++){
380 >      // velocity half step
381 >      vel[j] += (dt2 * frc[j] / mass) * eConvert;
382 >      // position whole step
383 >      pos[j] += dt * vel[j];
384 >    }
385  
386 <      dAtom = (DirectionalAtom *)atoms[i];
387 <          
386 >    integrableObjects[i]->setVel(vel);
387 >    integrableObjects[i]->setPos(pos);
388 >
389 >    if (integrableObjects[i]->isDirectional()){
390 >
391        // get and convert the torque to body frame
325 –      
326 –      Tb[0] = dAtom->getTx();
327 –      Tb[1] = dAtom->getTy();
328 –      Tb[2] = dAtom->getTz();
392  
393 <      dAtom->lab2Body( Tb );
393 >      integrableObjects[i]->getTrq(Tb);
394 >      integrableObjects[i]->lab2Body(Tb);
395  
396        // get the angular momentum, and propagate a half step
333 –      
334 –      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
335 –      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
336 –      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
337 –      
338 –      // use the angular velocities to propagate the rotation matrix a
339 –      // full time step
340 –      
341 –      // rotate about the x-axis      
342 –      angle = dt2 * ji[0] / dAtom->getIxx();
343 –      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
397  
398 <      // rotate about the y-axis
346 <      angle = dt2 * ji[1] / dAtom->getIyy();
347 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
348 <      
349 <      // rotate about the z-axis
350 <      angle = dt * ji[2] / dAtom->getIzz();
351 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
352 <      
353 <      // rotate about the y-axis
354 <      angle = dt2 * ji[1] / dAtom->getIyy();
355 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
356 <      
357 <       // rotate about the x-axis
358 <      angle = dt2 * ji[0] / dAtom->getIxx();
359 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
360 <      
361 <      dAtom->setJx( ji[0] );
362 <      dAtom->setJy( ji[1] );
363 <      dAtom->setJz( ji[2] );
398 >      integrableObjects[i]->getJ(ji);
399  
400 <      std::cerr << "Amat[" << i << "]\n";
401 <      info->printMat9( &Amat[aMatIndex] );
402 <          
403 <      std::cerr << "ji[" << i << "]\t"
404 <                << ji[0] << "\t"
405 <                << ji[1] << "\t"
371 <                << ji[2] << "\n";
372 <          
400 >      for (j = 0; j < 3; j++)
401 >        ji[j] += (dt2 * Tb[j]) * eConvert;
402 >
403 >      this->rotationPropagation( integrableObjects[i], ji );
404 >
405 >      integrableObjects[i]->setJ(ji);
406      }
374 –    
407    }
408 +
409 +  consFramework->doConstrainA();
410   }
411  
412  
413 < void Integrator::moveB( void ){
414 <  int i,j,k;
415 <  int atomIndex, aMatIndex;
416 <  DirectionalAtom* dAtom;
417 <  double Tb[3];
384 <  double ji[3];
413 > template<typename T> void Integrator<T>::moveB(void){
414 >  int i, j;
415 >  double Tb[3], ji[3];
416 >  double vel[3], frc[3];
417 >  double mass;
418  
419 <  for( i=0; i<nAtoms; i++ ){
420 <    atomIndex = i * 3;
421 <    aMatIndex = i * 9;
419 >  for (i = 0; i < integrableObjects.size(); i++){
420 >    integrableObjects[i]->getVel(vel);
421 >    integrableObjects[i]->getFrc(frc);
422  
423 +    mass = integrableObjects[i]->getMass();
424 +
425      // velocity half step
426 <    for( j=atomIndex; j<(atomIndex+3); j++ )
427 <      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
426 >    for (j = 0; j < 3; j++)
427 >      vel[j] += (dt2 * frc[j] / mass) * eConvert;
428  
429 <
395 <    if( atoms[i]->isDirectional() ){
396 <      
397 <      dAtom = (DirectionalAtom *)atoms[i];
398 <      
399 <      // get and convert the torque to body frame
400 <      
401 <      Tb[0] = dAtom->getTx();
402 <      Tb[1] = dAtom->getTy();
403 <      Tb[2] = dAtom->getTz();
404 <      
405 <      std::cerr << "TrqB[" << i << "]\t"
406 <                << Tb[0] << "\t"
407 <                << Tb[1] << "\t"
408 <                << Tb[2] << "\n";
429 >    integrableObjects[i]->setVel(vel);
430  
431 <      dAtom->lab2Body( Tb );
411 <      
412 <      // get the angular momentum, and complete the angular momentum
413 <      // half step
414 <      
415 <      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
416 <      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
417 <      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
418 <      
419 <      dAtom->setJx( ji[0] );
420 <      dAtom->setJy( ji[1] );
421 <      dAtom->setJz( ji[2] );
431 >    if (integrableObjects[i]->isDirectional()){
432  
433 +      // get and convert the torque to body frame
434  
435 <      std::cerr << "Amat[" << i << "]\n";
436 <      info->printMat9( &Amat[aMatIndex] );
437 <          
438 <      std::cerr << "ji[" << i << "]\t"
439 <                << ji[0] << "\t"
440 <                << ji[1] << "\t"
441 <                << ji[2] << "\n";
435 >      integrableObjects[i]->getTrq(Tb);
436 >      integrableObjects[i]->lab2Body(Tb);
437 >
438 >      // get the angular momentum, and propagate a half step
439 >
440 >      integrableObjects[i]->getJ(ji);
441 >
442 >      for (j = 0; j < 3; j++)
443 >        ji[j] += (dt2 * Tb[j]) * eConvert;
444 >
445 >
446 >      integrableObjects[i]->setJ(ji);
447      }
448    }
449  
450 +  consFramework->doConstrainB();
451   }
452  
453 < void Integrator::preMove( void ){
454 <  int i;
453 > /*
454 > template<typename T> void Integrator<T>::preMove(void){
455 >  int i, j;
456 >  double pos[3];
457  
458 <  if( nConstrained ){
458 >  if (nConstrained){
459 >    for (i = 0; i < nAtoms; i++){
460 >      atoms[i]->getPos(pos);
461  
462 <    for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
462 >      for (j = 0; j < 3; j++){
463 >        oldPos[3 * i + j] = pos[j];
464 >      }
465 >    }
466    }
467 < }  
467 > }
468  
469 < void Integrator::constrainA(){
470 <
447 <  int i,j,k;
469 > template<typename T> void Integrator<T>::constrainA(){
470 >  int i, j;
471    int done;
472 +  double posA[3], posB[3];
473 +  double velA[3], velB[3];
474    double pab[3];
475    double rab[3];
476    int a, b, ax, ay, az, bx, by, bz;
# Line 457 | Line 482 | void Integrator::constrainA(){
482    double gab;
483    int iteration;
484  
485 <  for( i=0; i<nAtoms; i++){
461 <    
485 >  for (i = 0; i < nAtoms; i++){
486      moving[i] = 0;
487 <    moved[i]  = 1;
487 >    moved[i] = 1;
488    }
489  
490    iteration = 0;
491    done = 0;
492 <  while( !done && (iteration < maxIteration )){
469 <
492 >  while (!done && (iteration < maxIteration)){
493      done = 1;
494 <    for(i=0; i<nConstrained; i++){
472 <
494 >    for (i = 0; i < nConstrained; i++){
495        a = constrainedA[i];
496        b = constrainedB[i];
475 –      
476 –      ax = (a*3) + 0;
477 –      ay = (a*3) + 1;
478 –      az = (a*3) + 2;
497  
498 <      bx = (b*3) + 0;
499 <      by = (b*3) + 1;
500 <      bz = (b*3) + 2;
498 >      ax = (a * 3) + 0;
499 >      ay = (a * 3) + 1;
500 >      az = (a * 3) + 2;
501  
502 <      if( moved[a] || moved[b] ){
503 <        
504 <        pab[0] = pos[ax] - pos[bx];
487 <        pab[1] = pos[ay] - pos[by];
488 <        pab[2] = pos[az] - pos[bz];
502 >      bx = (b * 3) + 0;
503 >      by = (b * 3) + 1;
504 >      bz = (b * 3) + 2;
505  
506 <        //periodic boundary condition
506 >      if (moved[a] || moved[b]){
507 >        atoms[a]->getPos(posA);
508 >        atoms[b]->getPos(posB);
509  
510 <        info->wrapVector( pab );
510 >        for (j = 0; j < 3; j++)
511 >          pab[j] = posA[j] - posB[j];
512  
513 <        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
513 >        //periodic boundary condition
514  
515 <        rabsq = constrainedDsqr[i];
497 <        diffsq = rabsq - pabsq;
515 >        info->wrapVector(pab);
516  
517 <        // the original rattle code from alan tidesley
500 <        if (fabs(diffsq) > (tol*rabsq*2)) {
501 <          rab[0] = oldPos[ax] - oldPos[bx];
502 <          rab[1] = oldPos[ay] - oldPos[by];
503 <          rab[2] = oldPos[az] - oldPos[bz];
517 >        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
518  
519 <          info->wrapVector( rab );
519 >        rabsq = constrainedDsqr[i];
520 >        diffsq = rabsq - pabsq;
521  
522 <          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
522 >        // the original rattle code from alan tidesley
523 >        if (fabs(diffsq) > (tol * rabsq * 2)){
524 >          rab[0] = oldPos[ax] - oldPos[bx];
525 >          rab[1] = oldPos[ay] - oldPos[by];
526 >          rab[2] = oldPos[az] - oldPos[bz];
527  
528 <          rpabsq = rpab * rpab;
528 >          info->wrapVector(rab);
529  
530 +          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
531  
532 <          if (rpabsq < (rabsq * -diffsq)){
532 >          rpabsq = rpab * rpab;
533  
534 +
535 +          if (rpabsq < (rabsq * -diffsq)){
536   #ifdef IS_MPI
537 <            a = atoms[a]->getGlobalIndex();
538 <            b = atoms[b]->getGlobalIndex();
537 >            a = atoms[a]->getGlobalIndex();
538 >            b = atoms[b]->getGlobalIndex();
539   #endif //is_mpi
540 <            sprintf( painCave.errMsg,
541 <                     "Constraint failure in constrainA at atom %d and %d.\n",
542 <                     a, b );
543 <            painCave.isFatal = 1;
544 <            simError();
545 <          }
540 >            sprintf(painCave.errMsg,
541 >                    "Constraint failure in constrainA at atom %d and %d.\n", a,
542 >                    b);
543 >            painCave.isFatal = 1;
544 >            simError();
545 >          }
546  
547 <          rma = 1.0 / atoms[a]->getMass();
548 <          rmb = 1.0 / atoms[b]->getMass();
547 >          rma = 1.0 / atoms[a]->getMass();
548 >          rmb = 1.0 / atoms[b]->getMass();
549  
550 <          gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
550 >          gab = diffsq / (2.0 * (rma + rmb) * rpab);
551  
552            dx = rab[0] * gab;
553            dy = rab[1] * gab;
554            dz = rab[2] * gab;
555  
556 <          pos[ax] += rma * dx;
557 <          pos[ay] += rma * dy;
558 <          pos[az] += rma * dz;
556 >          posA[0] += rma * dx;
557 >          posA[1] += rma * dy;
558 >          posA[2] += rma * dz;
559  
560 <          pos[bx] -= rmb * dx;
539 <          pos[by] -= rmb * dy;
540 <          pos[bz] -= rmb * dz;
560 >          atoms[a]->setPos(posA);
561  
562 +          posB[0] -= rmb * dx;
563 +          posB[1] -= rmb * dy;
564 +          posB[2] -= rmb * dz;
565 +
566 +          atoms[b]->setPos(posB);
567 +
568            dx = dx / dt;
569            dy = dy / dt;
570            dz = dz / dt;
571  
572 <          vel[ax] += rma * dx;
547 <          vel[ay] += rma * dy;
548 <          vel[az] += rma * dz;
572 >          atoms[a]->getVel(velA);
573  
574 <          vel[bx] -= rmb * dx;
575 <          vel[by] -= rmb * dy;
576 <          vel[bz] -= rmb * dz;
574 >          velA[0] += rma * dx;
575 >          velA[1] += rma * dy;
576 >          velA[2] += rma * dz;
577  
578 <          moving[a] = 1;
579 <          moving[b] = 1;
580 <          done = 0;
581 <        }
578 >          atoms[a]->setVel(velA);
579 >
580 >          atoms[b]->getVel(velB);
581 >
582 >          velB[0] -= rmb * dx;
583 >          velB[1] -= rmb * dy;
584 >          velB[2] -= rmb * dz;
585 >
586 >          atoms[b]->setVel(velB);
587 >
588 >          moving[a] = 1;
589 >          moving[b] = 1;
590 >          done = 0;
591 >        }
592        }
593      }
594 <    
595 <    for(i=0; i<nAtoms; i++){
562 <      
594 >
595 >    for (i = 0; i < nAtoms; i++){
596        moved[i] = moving[i];
597        moving[i] = 0;
598      }
# Line 567 | Line 600 | void Integrator::constrainA(){
600      iteration++;
601    }
602  
603 <  if( !done ){
604 <
605 <    sprintf( painCave.errMsg,
606 <             "Constraint failure in constrainA, too many iterations: %d\n",
574 <             iteration );
603 >  if (!done){
604 >    sprintf(painCave.errMsg,
605 >            "Constraint failure in constrainA, too many iterations: %d\n",
606 >            iteration);
607      painCave.isFatal = 1;
608      simError();
609    }
610  
611   }
612  
613 < void Integrator::constrainB( void ){
614 <  
583 <  int i,j,k;
613 > template<typename T> void Integrator<T>::constrainB(void){
614 >  int i, j;
615    int done;
616 +  double posA[3], posB[3];
617 +  double velA[3], velB[3];
618    double vxab, vyab, vzab;
619    double rab[3];
620    int a, b, ax, ay, az, bx, by, bz;
621    double rma, rmb;
622    double dx, dy, dz;
623 <  double rabsq, pabsq, rvab;
591 <  double diffsq;
623 >  double rvab;
624    double gab;
625    int iteration;
626  
627 <  for(i=0; i<nAtoms; i++){
627 >  for (i = 0; i < nAtoms; i++){
628      moving[i] = 0;
629      moved[i] = 1;
630    }
631  
632    done = 0;
633    iteration = 0;
634 <  while( !done && (iteration < maxIteration ) ){
603 <
634 >  while (!done && (iteration < maxIteration)){
635      done = 1;
636  
637 <    for(i=0; i<nConstrained; i++){
607 <      
637 >    for (i = 0; i < nConstrained; i++){
638        a = constrainedA[i];
639        b = constrainedB[i];
640  
641 <      ax = (a*3) + 0;
642 <      ay = (a*3) + 1;
643 <      az = (a*3) + 2;
641 >      ax = (a * 3) + 0;
642 >      ay = (a * 3) + 1;
643 >      az = (a * 3) + 2;
644  
645 <      bx = (b*3) + 0;
646 <      by = (b*3) + 1;
647 <      bz = (b*3) + 2;
645 >      bx = (b * 3) + 0;
646 >      by = (b * 3) + 1;
647 >      bz = (b * 3) + 2;
648  
649 <      if( moved[a] || moved[b] ){
650 <        
651 <        vxab = vel[ax] - vel[bx];
622 <        vyab = vel[ay] - vel[by];
623 <        vzab = vel[az] - vel[bz];
649 >      if (moved[a] || moved[b]){
650 >        atoms[a]->getVel(velA);
651 >        atoms[b]->getVel(velB);
652  
653 <        rab[0] = pos[ax] - pos[bx];
654 <        rab[1] = pos[ay] - pos[by];
655 <        rab[2] = pos[az] - pos[bz];
628 <        
629 <        info->wrapVector( rab );
630 <        
631 <        rma = 1.0 / atoms[a]->getMass();
632 <        rmb = 1.0 / atoms[b]->getMass();
653 >        vxab = velA[0] - velB[0];
654 >        vyab = velA[1] - velB[1];
655 >        vzab = velA[2] - velB[2];
656  
657 <        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
658 <          
636 <        gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
657 >        atoms[a]->getPos(posA);
658 >        atoms[b]->getPos(posB);
659  
660 <        if (fabs(gab) > tol) {
661 <          
640 <          dx = rab[0] * gab;
641 <          dy = rab[1] * gab;
642 <          dz = rab[2] * gab;
643 <          
644 <          vel[ax] += rma * dx;
645 <          vel[ay] += rma * dy;
646 <          vel[az] += rma * dz;
660 >        for (j = 0; j < 3; j++)
661 >          rab[j] = posA[j] - posB[j];
662  
663 <          vel[bx] -= rmb * dx;
664 <          vel[by] -= rmb * dy;
665 <          vel[bz] -= rmb * dz;
666 <          
667 <          moving[a] = 1;
668 <          moving[b] = 1;
669 <          done = 0;
670 <        }
663 >        info->wrapVector(rab);
664 >
665 >        rma = 1.0 / atoms[a]->getMass();
666 >        rmb = 1.0 / atoms[b]->getMass();
667 >
668 >        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
669 >
670 >        gab = -rvab / ((rma + rmb) * constrainedDsqr[i]);
671 >
672 >        if (fabs(gab) > tol){
673 >          dx = rab[0] * gab;
674 >          dy = rab[1] * gab;
675 >          dz = rab[2] * gab;
676 >
677 >          velA[0] += rma * dx;
678 >          velA[1] += rma * dy;
679 >          velA[2] += rma * dz;
680 >
681 >          atoms[a]->setVel(velA);
682 >
683 >          velB[0] -= rmb * dx;
684 >          velB[1] -= rmb * dy;
685 >          velB[2] -= rmb * dz;
686 >
687 >          atoms[b]->setVel(velB);
688 >
689 >          moving[a] = 1;
690 >          moving[b] = 1;
691 >          done = 0;
692 >        }
693        }
694      }
695  
696 <    for(i=0; i<nAtoms; i++){
696 >    for (i = 0; i < nAtoms; i++){
697        moved[i] = moving[i];
698        moving[i] = 0;
699      }
700 <    
700 >
701      iteration++;
702    }
703  
704 <  if( !done ){
705 <
706 <  
707 <    sprintf( painCave.errMsg,
671 <             "Constraint failure in constrainB, too many iterations: %d\n",
672 <             iteration );
704 >  if (!done){
705 >    sprintf(painCave.errMsg,
706 >            "Constraint failure in constrainB, too many iterations: %d\n",
707 >            iteration);
708      painCave.isFatal = 1;
709      simError();
710 <  }
676 <
710 >  }
711   }
712 + */
713 + template<typename T> void Integrator<T>::rotationPropagation
714 + ( StuntDouble* sd, double ji[3] ){
715  
716 +  double angle;
717 +  double A[3][3], I[3][3];
718 +  int i, j, k;
719  
720 +  // use the angular velocities to propagate the rotation matrix a
721 +  // full time step
722  
723 +  sd->getA(A);
724 +  sd->getI(I);
725  
726 +  if (sd->isLinear()) {
727 +    i = sd->linearAxis();
728 +    j = (i+1)%3;
729 +    k = (i+2)%3;
730 +    
731 +    angle = dt2 * ji[j] / I[j][j];
732 +    this->rotate( k, i, angle, ji, A );
733  
734 +    angle = dt * ji[k] / I[k][k];
735 +    this->rotate( i, j, angle, ji, A);
736  
737 +    angle = dt2 * ji[j] / I[j][j];
738 +    this->rotate( k, i, angle, ji, A );
739  
740 < void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
741 <                         double A[9] ){
740 >  } else {
741 >    // rotate about the x-axis
742 >    angle = dt2 * ji[0] / I[0][0];
743 >    this->rotate( 1, 2, angle, ji, A );
744 >    
745 >    // rotate about the y-axis
746 >    angle = dt2 * ji[1] / I[1][1];
747 >    this->rotate( 2, 0, angle, ji, A );
748 >    
749 >    // rotate about the z-axis
750 >    angle = dt * ji[2] / I[2][2];
751 >    sd->addZangle(angle);
752 >    this->rotate( 0, 1, angle, ji, A);
753 >    
754 >    // rotate about the y-axis
755 >    angle = dt2 * ji[1] / I[1][1];
756 >    this->rotate( 2, 0, angle, ji, A );
757 >    
758 >    // rotate about the x-axis
759 >    angle = dt2 * ji[0] / I[0][0];
760 >    this->rotate( 1, 2, angle, ji, A );
761 >    
762 >  }
763 >  sd->setA( A  );
764 > }
765  
766 <  int i,j,k;
766 > template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
767 >                                                double angle, double ji[3],
768 >                                                double A[3][3]){
769 >  int i, j, k;
770    double sinAngle;
771    double cosAngle;
772    double angleSqr;
# Line 695 | Line 776 | void Integrator::rotate( int axes1, int axes2, double
776    double tempA[3][3];
777    double tempJ[3];
778  
698 –
779    // initialize the tempA
780  
781 <  for(i=0; i<3; i++){
782 <    for(j=0; j<3; j++){
783 <      tempA[j][i] = A[3*i+j];
781 >  for (i = 0; i < 3; i++){
782 >    for (j = 0; j < 3; j++){
783 >      tempA[j][i] = A[i][j];
784      }
785    }
786  
787    // initialize the tempJ
788  
789 <  for( i=0; i<3; i++) tempJ[i] = ji[i];
790 <  
789 >  for (i = 0; i < 3; i++)
790 >    tempJ[i] = ji[i];
791 >
792    // initalize rot as a unit matrix
793  
794    rot[0][0] = 1.0;
# Line 717 | Line 798 | void Integrator::rotate( int axes1, int axes2, double
798    rot[1][0] = 0.0;
799    rot[1][1] = 1.0;
800    rot[1][2] = 0.0;
801 <  
801 >
802    rot[2][0] = 0.0;
803    rot[2][1] = 0.0;
804    rot[2][2] = 1.0;
805 <  
805 >
806    // use a small angle aproximation for sin and cosine
807  
808 <  angleSqr  = angle * angle;
808 >  angleSqr = angle * angle;
809    angleSqrOver4 = angleSqr / 4.0;
810    top = 1.0 - angleSqrOver4;
811    bottom = 1.0 + angleSqrOver4;
# Line 737 | Line 818 | void Integrator::rotate( int axes1, int axes2, double
818  
819    rot[axes1][axes2] = sinAngle;
820    rot[axes2][axes1] = -sinAngle;
821 <  
821 >
822    // rotate the momentum acoording to: ji[] = rot[][] * ji[]
823 <  
824 <  for(i=0; i<3; i++){
823 >
824 >  for (i = 0; i < 3; i++){
825      ji[i] = 0.0;
826 <    for(k=0; k<3; k++){
826 >    for (k = 0; k < 3; k++){
827        ji[i] += rot[i][k] * tempJ[k];
828      }
829    }
830  
831 <  // rotate the Rotation matrix acording to:
831 >  // rotate the Rotation matrix acording to:
832    //            A[][] = A[][] * transpose(rot[][])
833  
834  
# Line 755 | Line 836 | void Integrator::rotate( int axes1, int axes2, double
836    // calculation as:
837    //                transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
838  
839 <  for(i=0; i<3; i++){
840 <    for(j=0; j<3; j++){
841 <      A[3*j+i] = 0.0;
842 <      for(k=0; k<3; k++){
843 <        A[3*j+i] += tempA[i][k] * rot[j][k];
839 >  for (i = 0; i < 3; i++){
840 >    for (j = 0; j < 3; j++){
841 >      A[j][i] = 0.0;
842 >      for (k = 0; k < 3; k++){
843 >        A[j][i] += tempA[i][k] * rot[j][k];
844        }
845      }
846    }
847   }
848 +
849 + template<typename T> void Integrator<T>::calcForce(int calcPot, int calcStress){
850 +  myFF->doForces(calcPot, calcStress);
851 + }
852 +
853 + template<typename T> void Integrator<T>::thermalize(){
854 +  tStats->velocitize();
855 + }
856 +
857 + template<typename T> double Integrator<T>::getConservedQuantity(void){
858 +  return tStats->getTotalE();
859 + }
860 + template<typename T> string Integrator<T>::getAdditionalParameters(void){
861 +  //By default, return a null string
862 +  //The reason we use string instead of char* is that if we use char*, we will
863 +  //return a pointer point to local variable which might cause problem
864 +  return string();
865 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines