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

Diff Legend

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