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 601 by gezelter, Mon Jul 14 23:06:09 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 < 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  
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 <  
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 | Integrator::~Integrator() {
55      delete[] moved;
56      delete[] oldPos;
57    }
54 –  
58   }
59  
60 < void Integrator::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 <      std::cerr << "Is the folowing bond constrained \n";
77 <      theArray[j]->printMe();
78 <      
79 <      if(constrained){
80 <        
81 <        std::cerr << "Yes\n";
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() );
86 <        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
87 <        
88 <        nConstrained++;
89 <        constrained = 0;
90 <      }
91 <      else std::cerr << "No.\n";
82 >        nConstrained++;
83 >        constrained = 0;
84 >      }
85      }
86  
87 <    theArray = (SRI**) molecules[i].getMyBends();
88 <    for(int j=0; j<molecules[i].getNBends(); j++){
96 <      
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++;
107 <        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++){
113 <      
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++;
124 <        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++){
142 <      
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();
146 –
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 < void Integrator::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;
169 <  double sampleTime  = info->sampleTime;
170 <  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 currTime;
166 >  double currReset;
167  
168    int calcPot, calcStress;
179 –  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;
186 –  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 <  myFF->doForces(1,1);
183 >  calcForce(1, 1);
184 >
185 >  //temp test
186 >  tStats->getPotential();
187    
188 <  if( info->setTemp ){
189 <    
190 <    tStats->velocitize();
188 >  if (nConstrained){
189 >    preMove();
190 >    constrainA();
191 >    calcForce(1, 1);
192 >    constrainB();
193    }
194    
195 <  dumpOut->writeDump( 0.0 );
196 <  statOut->writeStat( 0.0 );
197 <  
195 >  if (info->setTemp){
196 >    thermalize();
197 >  }
198 >
199    calcPot     = 0;
200    calcStress  = 0;
201 <  currSample  = sampleTime;
202 <  currThermal = thermalTime;
203 <  currStatus  = statusTime;
204 <  currTime    = 0.0;;
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());
208  
211 –  readyCheck();
209  
210   #ifdef IS_MPI
211 <  strcpy( checkPointMsg,
215 <          "The integrator is ready to go." );
211 >  strcpy(checkPointMsg, "The integrator is ready to go.");
212    MPIcheckPoint();
213   #endif // is_mpi
214  
215 <  while( currTime < runTime ){
216 <
221 <    if( (currTime+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 <      
223 <    currTime += dt;
221 > #ifdef PROFILE
222 >    startProfile( pro1 );
223 > #endif
224 >    
225 >    integrateStep(calcPot, calcStress);
226  
227 <    if( info->setTemp ){
228 <      if( currTime >= currThermal ){
229 <        tStats->velocitize();
230 <        currThermal += thermalTime;
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;
239        }
240      }
241  
242 <    if( currTime >= currSample ){
243 <      dumpOut->writeDump( currTime );
242 >    if (info->getTime() >= currSample){
243 >      dumpOut->writeDump(info->getTime());
244        currSample += sampleTime;
245      }
246  
247 <    if( currTime >= currStatus ){
248 <      statOut->writeStat( currTime );
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,
251 <            "successfully took a time step." );
266 >    strcpy(checkPointMsg, "successfully took a time step.");
267      MPIcheckPoint();
268   #endif // is_mpi
254 –
269    }
270  
257 –  dumpOut->writeFinal(currTime);
258 –
271    delete dumpOut;
272    delete statOut;
273   }
274  
275 < void Integrator::integrateStep( int calcPot, int calcStress ){
276 <
265 <
266 <      
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();
271 –  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");
302 +  MPIcheckPoint();
303 + #endif // is_mpi
304 +
305 +
306    // calc forces
307  
308 <  myFF->doForces(calcPot,calcStress);
308 >  calcForce(calcPot, calcStress);
309  
310 + #ifdef IS_MPI
311 +  strcpy(checkPointMsg, "Succesful doForces\n");
312 +  MPIcheckPoint();
313 + #endif // is_mpi
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");
331 >  MPIcheckPoint();
332 > #endif // is_mpi
333   }
334  
335  
336 < void Integrator::moveA( void ){
286 <  
336 > template<typename T> void Integrator<T>::moveA(void){
337    int i, j;
338    DirectionalAtom* dAtom;
339    double Tb[3], ji[3];
290 –  double A[3][3], I[3][3];
291 –  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  
297 –    atoms[i]->getVel( vel );
298 –    atoms[i]->getPos( pos );
299 –    atoms[i]->getFrc( frc );
300 –
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  
315 –      dAtom = (DirectionalAtom *)atoms[i];
316 –          
363        // get and convert the torque to body frame
318 –      
319 –      dAtom->getTrq( Tb );
320 –      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;
328 –      
329 –      // use the angular velocities to propagate the rotation matrix a
330 –      // full time step
374  
375 <      dAtom->getA(A);
333 <      dAtom->getI(I);
334 <    
335 <      // rotate about the x-axis      
336 <      angle = dt2 * ji[0] / I[0][0];
337 <      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 );
342 <      
343 <      // rotate about the z-axis
344 <      angle = dt * ji[2] / I[2][2];
345 <      this->rotate( 0, 1, angle, ji, A);
346 <      
347 <      // rotate about the y-axis
348 <      angle = dt2 * ji[1] / I[1][1];
349 <      this->rotate( 2, 0, angle, ji, A );
350 <      
351 <       // rotate about the x-axis
352 <      angle = dt2 * ji[0] / I[0][0];
353 <      this->rotate( 1, 2, angle, ji, A );
354 <      
377 >      dAtom->setJ(ji);
378 >    }
379 >  }
380  
381 <      dAtom->setJ( ji );
382 <      dAtom->setA( A  );
358 <          
359 <    }    
381 >  if (nConstrained){
382 >    constrainA();
383    }
384   }
385  
386  
387 < void Integrator::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 );
374 <    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;
381 <    
382 <    atoms[i]->setVel( vel );
383 <
384 <    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 );
391 <      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;
399 –      
420  
421 <      dAtom->setJ( ji );
421 >
422 >      dAtom->setJ(ji);
423      }
424    }
425 +
426 +  if (nConstrained){
427 +    constrainB();
428 +  }
429   }
430  
431 < void Integrator::preMove( void ){
431 > template<typename T> void Integrator<T>::preMove(void){
432    int i, j;
433    double pos[3];
434  
435 <  if( nConstrained ){
436 <
437 <    for(i=0; i < nAtoms; i++) {
413 <
414 <      atoms[i]->getPos( pos );
415 <
416 <      for (j = 0; j < 3; j++) {        
417 <        oldPos[3*i + j] = pos[j];
418 <      }
435 >  if (nConstrained){
436 >    for (i = 0; i < nAtoms; i++){
437 >      atoms[i]->getPos(pos);
438  
439 +      for (j = 0; j < 3; j++){
440 +        oldPos[3 * i + j] = pos[j];
441 +      }
442      }
443 <  }  
443 >  }
444   }
445  
446 < void Integrator::constrainA(){
447 <
426 <  int i,j,k;
446 > template<typename T> void Integrator<T>::constrainA(){
447 >  int i, j;
448    int done;
449    double posA[3], posB[3];
450    double velA[3], velB[3];
# Line 438 | Line 459 | void Integrator::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 )){
449 <
469 >  while (!done && (iteration < maxIteration)){
470      done = 1;
471 <    for(i=0; i<nConstrained; i++){
452 <
471 >    for (i = 0; i < nConstrained; i++){
472        a = constrainedA[i];
473        b = constrainedB[i];
455 –      
456 –      ax = (a*3) + 0;
457 –      ay = (a*3) + 1;
458 –      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];
471 –        
472 –        //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];
479 <        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)) {
483 <          rab[0] = oldPos[ax] - oldPos[bx];
484 <          rab[1] = oldPos[ay] - oldPos[by];
485 <          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  
494 –          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++){
556 <      
571 >
572 >    for (i = 0; i < nAtoms; i++){
573        moved[i] = moving[i];
574        moving[i] = 0;
575      }
# Line 561 | Line 577 | void Integrator::constrainA(){
577      iteration++;
578    }
579  
580 <  if( !done ){
581 <
582 <    sprintf( painCave.errMsg,
583 <             "Constraint failure in constrainA, too many iterations: %d\n",
568 <             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 < void Integrator::constrainB( void ){
591 <  
577 <  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 583 | Line 597 | void Integrator::constrainB( void ){
597    int a, b, ax, ay, az, bx, by, bz;
598    double rma, rmb;
599    double dx, dy, dz;
600 <  double rabsq, pabsq, rvab;
587 <  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 ) ){
599 <
611 >  while (!done && (iteration < maxIteration)){
612      done = 1;
613  
614 <    for(i=0; i<nConstrained; i++){
603 <      
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 <          
620 <        vxab = velA[0] - velB[0];
621 <        vyab = velA[1] - velB[1];
622 <        vzab = velA[2] - velB[2];
623 <
624 <        atoms[a]->getPos( posA );
625 <        atoms[b]->getPos( posB );
630 >        vxab = velA[0] - velB[0];
631 >        vyab = velA[1] - velB[1];
632 >        vzab = velA[2] - velB[2];
633  
634 <        for (j = 0; j < 3; j++)
634 >        atoms[a]->getPos(posA);
635 >        atoms[b]->getPos(posB);
636 >
637 >        for (j = 0; j < 3; j++)
638            rab[j] = posA[j] - posB[j];
629 –          
630 –        info->wrapVector( rab );
631 –        
632 –        rma = 1.0 / atoms[a]->getMass();
633 –        rmb = 1.0 / atoms[b]->getMass();
639  
640 <        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
636 <          
637 <        gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
640 >        info->wrapVector(rab);
641  
642 <        if (fabs(gab) > tol) {
643 <          
641 <          dx = rab[0] * gab;
642 <          dy = rab[1] * gab;
643 <          dz = rab[2] * gab;
644 <        
645 <          velA[0] += rma * dx;
646 <          velA[1] += rma * dy;
647 <          velA[2] += rma * dz;
642 >        rma = 1.0 / atoms[a]->getMass();
643 >        rmb = 1.0 / atoms[b]->getMass();
644  
645 <          atoms[a]->setVel( velA );
645 >        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
646  
647 <          velB[0] -= rmb * dx;
652 <          velB[1] -= rmb * dy;
653 <          velB[2] -= rmb * dz;
647 >        gab = -rvab / ((rma + rmb) * constrainedDsqr[i]);
648  
649 <          atoms[b]->setVel( velB );
650 <          
651 <          moving[a] = 1;
652 <          moving[b] = 1;
653 <          done = 0;
654 <        }
649 >        if (fabs(gab) > tol){
650 >          dx = rab[0] * gab;
651 >          dy = rab[1] * gab;
652 >          dz = rab[2] * gab;
653 >
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    }
671 –  
672 –  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 <  }
681 <
687 >  }
688   }
689  
690 < void Integrator::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 696 | Line 737 | void Integrator::rotate( int axes1, int axes2, double
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 715 | Line 757 | void Integrator::rotate( int axes1, int axes2, double
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 735 | Line 777 | void Integrator::rotate( int axes1, int axes2, double
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 753 | Line 795 | void Integrator::rotate( int axes1, int axes2, double
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);
810 + }
811 +
812 + template<typename T> void Integrator<T>::thermalize(){
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