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

Diff Legend

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