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 677 by tim, Mon Aug 11 19:41:36 2003 UTC vs.
Revision 1125 by gezelter, Mon Apr 19 22:13:01 2004 UTC

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

Diff Legend

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