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

Diff Legend

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