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 561 by mmeineke, Fri Jun 20 20:29:36 2003 UTC vs.
Revision 782 by mmeineke, Tue Sep 23 20:34:31 2003 UTC

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

Diff Legend

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