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 567 by mmeineke, Wed Jun 25 21:12:14 2003 UTC vs.
Revision 784 by mmeineke, Wed Sep 24 19:34:39 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();
140 –
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 ){
159 <
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  
176 –
177 –  tStats   = new Thermo( info );
178 –  statOut  = new StatWriter( info );
179 –  dumpOut  = new DumpWriter( info );
180 –
170    atoms = info->atoms;
171    DirectionalAtom* dAtom;
172  
173    dt = info->dt;
174    dt2 = 0.5 * dt;
175  
176 +  readyCheck();
177 +
178    // initialize the forces before the first step
179  
180 <  myFF->doForces(1,1);
181 <  
182 <  if( info->setTemp ){
183 <    
184 <    tStats->velocitize();
180 >  calcForce(1, 1);
181 >
182 >  if (nConstrained){
183 >    preMove();
184 >    constrainA();
185 >    calcForce(1, 1);    
186 >    constrainB();
187    }
188    
189 <  dumpOut->writeDump( 0.0 );
190 <  statOut->writeStat( 0.0 );
191 <  
189 >  if (info->setTemp){
190 >    thermalize();
191 >  }
192 >
193    calcPot     = 0;
194    calcStress  = 0;
195 <  currSample  = sampleTime;
196 <  currThermal = thermalTime;
197 <  currStatus  = statusTime;
198 <  currTime    = 0.0;;
195 >  currSample  = sampleTime + info->getTime();
196 >  currThermal = thermalTime+ info->getTime();
197 >  currStatus  = statusTime + info->getTime();
198 >  currReset   = resetTime  + info->getTime();
199  
200 +  dumpOut->writeDump(info->getTime());
201 +  statOut->writeStat(info->getTime());
202  
207 –  readyCheck();
203  
204 +
205   #ifdef IS_MPI
206 <  strcpy( checkPointMsg,
211 <          "The integrator is ready to go." );
206 >  strcpy(checkPointMsg, "The integrator is ready to go.");
207    MPIcheckPoint();
208   #endif // is_mpi
209  
210 <
211 <  pos  = Atom::getPosArray();
217 <  vel  = Atom::getVelArray();
218 <  frc  = Atom::getFrcArray();
219 <  trq  = Atom::getTrqArray();
220 <  Amat = Atom::getAmatArray();
221 <
222 <  while( currTime < runTime ){
223 <
224 <    if( (currTime+dt) >= currStatus ){
210 >  while (info->getTime() < runTime){
211 >    if ((info->getTime() + dt) >= currStatus){
212        calcPot = 1;
213        calcStress = 1;
214      }
215  
216 <    integrateStep( calcPot, calcStress );
230 <      
231 <    currTime += dt;
216 >    integrateStep(calcPot, calcStress);
217  
218 <    if( info->setTemp ){
219 <      if( currTime >= currThermal ){
220 <        tStats->velocitize();
221 <        currThermal += thermalTime;
218 >    info->incrTime(dt);
219 >
220 >    if (info->setTemp){
221 >      if (info->getTime() >= currThermal){
222 >        thermalize();
223 >        currThermal += thermalTime;
224        }
225      }
226  
227 <    if( currTime >= currSample ){
228 <      dumpOut->writeDump( currTime );
227 >    if (info->getTime() >= currSample){
228 >      dumpOut->writeDump(info->getTime());
229        currSample += sampleTime;
230      }
231  
232 <    if( currTime >= currStatus ){
233 <      statOut->writeStat( currTime );
232 >    if (info->getTime() >= currStatus){
233 >      statOut->writeStat(info->getTime());
234        calcPot = 0;
235        calcStress = 0;
236        currStatus += statusTime;
237      }
238  
239 +    if (info->resetIntegrator){
240 +      if (info->getTime() >= currReset){
241 +        this->resetIntegrator();
242 +        currReset += resetTime;
243 +      }
244 +    }
245 +
246   #ifdef IS_MPI
247 <    strcpy( checkPointMsg,
254 <            "successfully took a time step." );
247 >    strcpy(checkPointMsg, "successfully took a time step.");
248      MPIcheckPoint();
249   #endif // is_mpi
257 –
250    }
251  
252 <  dumpOut->writeFinal();
252 >  dumpOut->writeFinal(info->getTime());
253  
254    delete dumpOut;
255    delete statOut;
256   }
257  
258 < void Integrator::integrateStep( int calcPot, int calcStress ){
259 <
268 <
269 <      
258 > template<typename T> void Integrator<T>::integrateStep(int calcPot,
259 >                                                       int calcStress){
260    // Position full step, and velocity half step
271 –
261    preMove();
262 +
263    moveA();
274 –  if( nConstrained ) constrainA();
264  
265 +
266 +
267 +
268 + #ifdef IS_MPI
269 +  strcpy(checkPointMsg, "Succesful moveA\n");
270 +  MPIcheckPoint();
271 + #endif // is_mpi
272 +
273 +
274    // calc forces
275  
276 <  myFF->doForces(calcPot,calcStress);
276 >  calcForce(calcPot, calcStress);
277  
278 + #ifdef IS_MPI
279 +  strcpy(checkPointMsg, "Succesful doForces\n");
280 +  MPIcheckPoint();
281 + #endif // is_mpi
282 +
283 +
284    // finish the velocity  half step
285 <  
285 >
286    moveB();
287 <  if( nConstrained ) constrainB();
288 <  
287 >
288 >
289 >
290 > #ifdef IS_MPI
291 >  strcpy(checkPointMsg, "Succesful moveB\n");
292 >  MPIcheckPoint();
293 > #endif // is_mpi
294   }
295  
296  
297 < void Integrator::moveA( void ){
298 <  
290 <  int i,j,k;
291 <  int atomIndex, aMatIndex;
297 > template<typename T> void Integrator<T>::moveA(void){
298 >  int i, j;
299    DirectionalAtom* dAtom;
300 <  double Tb[3];
301 <  double ji[3];
302 <  double angle;
300 >  double Tb[3], ji[3];
301 >  double vel[3], pos[3], frc[3];
302 >  double mass;
303  
304 +  for (i = 0; i < nAtoms; i++){
305 +    atoms[i]->getVel(vel);
306 +    atoms[i]->getPos(pos);
307 +    atoms[i]->getFrc(frc);
308  
309 +    mass = atoms[i]->getMass();
310  
311 <  for( i=0; i<nAtoms; i++ ){
312 <    atomIndex = i * 3;
313 <    aMatIndex = i * 9;
311 >    for (j = 0; j < 3; j++){
312 >      // velocity half step
313 >      vel[j] += (dt2 * frc[j] / mass) * eConvert;
314 >      // position whole step
315 >      pos[j] += dt * vel[j];
316 >    }
317  
318 <    // velocity half step
319 <    for( j=atomIndex; j<(atomIndex+3); j++ )
305 <      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
318 >    atoms[i]->setVel(vel);
319 >    atoms[i]->setPos(pos);
320  
321 <    // position whole step    
322 <    for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j];
309 <    
310 <    if( atoms[i]->isDirectional() ){
321 >    if (atoms[i]->isDirectional()){
322 >      dAtom = (DirectionalAtom *) atoms[i];
323  
312 –      dAtom = (DirectionalAtom *)atoms[i];
313 –          
324        // get and convert the torque to body frame
325 <      
326 <      Tb[0] = dAtom->getTx();
327 <      Tb[1] = dAtom->getTy();
328 <      Tb[2] = dAtom->getTz();
319 <      
320 <      dAtom->lab2Body( Tb );
321 <      
325 >
326 >      dAtom->getTrq(Tb);
327 >      dAtom->lab2Body(Tb);
328 >
329        // get the angular momentum, and propagate a half step
330 <      
331 <      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
332 <      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
333 <      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
334 <      
335 <      // use the angular velocities to propagate the rotation matrix a
336 <      // full time step
337 <      
338 <      // rotate about the x-axis      
332 <      angle = dt2 * ji[0] / dAtom->getIxx();
333 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
334 <      
335 <      // rotate about the y-axis
336 <      angle = dt2 * ji[1] / dAtom->getIyy();
337 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
338 <      
339 <      // rotate about the z-axis
340 <      angle = dt * ji[2] / dAtom->getIzz();
341 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
342 <      
343 <      // rotate about the y-axis
344 <      angle = dt2 * ji[1] / dAtom->getIyy();
345 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
346 <      
347 <       // rotate about the x-axis
348 <      angle = dt2 * ji[0] / dAtom->getIxx();
349 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
350 <      
351 <      dAtom->setJx( ji[0] );
352 <      dAtom->setJy( ji[1] );
353 <      dAtom->setJz( ji[2] );
330 >
331 >      dAtom->getJ(ji);
332 >
333 >      for (j = 0; j < 3; j++)
334 >        ji[j] += (dt2 * Tb[j]) * eConvert;
335 >
336 >      this->rotationPropagation( dAtom, ji );
337 >
338 >      dAtom->setJ(ji);
339      }
355 –    
340    }
341 +
342 +  if (nConstrained){
343 +    constrainA();
344 +  }
345   }
346  
347  
348 < void Integrator::moveB( void ){
349 <  int i,j,k;
362 <  int atomIndex;
348 > template<typename T> void Integrator<T>::moveB(void){
349 >  int i, j;
350    DirectionalAtom* dAtom;
351 <  double Tb[3];
352 <  double ji[3];
351 >  double Tb[3], ji[3];
352 >  double vel[3], frc[3];
353 >  double mass;
354  
355 <  for( i=0; i<nAtoms; i++ ){
356 <    atomIndex = i * 3;
355 >  for (i = 0; i < nAtoms; i++){
356 >    atoms[i]->getVel(vel);
357 >    atoms[i]->getFrc(frc);
358  
359 +    mass = atoms[i]->getMass();
360 +
361      // velocity half step
362 <    for( j=atomIndex; j<(atomIndex+3); j++ )
363 <      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
362 >    for (j = 0; j < 3; j++)
363 >      vel[j] += (dt2 * frc[j] / mass) * eConvert;
364  
365 <    if( atoms[i]->isDirectional() ){
375 <      
376 <      dAtom = (DirectionalAtom *)atoms[i];
377 <      
378 <      // get and convert the torque to body frame
379 <      
380 <      Tb[0] = dAtom->getTx();
381 <      Tb[1] = dAtom->getTy();
382 <      Tb[2] = dAtom->getTz();
383 <      
384 <      dAtom->lab2Body( Tb );
385 <      
386 <      // get the angular momentum, and complete the angular momentum
387 <      // half step
388 <      
389 <      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
390 <      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
391 <      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
392 <      
393 <      dAtom->setJx( ji[0] );
394 <      dAtom->setJy( ji[1] );
395 <      dAtom->setJz( ji[2] );
396 <    }
397 <  }
365 >    atoms[i]->setVel(vel);
366  
367 < }
367 >    if (atoms[i]->isDirectional()){
368 >      dAtom = (DirectionalAtom *) atoms[i];
369  
370 < void Integrator::preMove( void ){
402 <  int i;
370 >      // get and convert the torque to body frame      
371  
372 <  if( nConstrained ){
372 >      dAtom->getTrq(Tb);
373 >      dAtom->lab2Body(Tb);
374  
375 <    for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
375 >      // get the angular momentum, and propagate a half step
376 >
377 >      dAtom->getJ(ji);
378 >
379 >      for (j = 0; j < 3; j++)
380 >        ji[j] += (dt2 * Tb[j]) * eConvert;
381 >
382 >
383 >      dAtom->setJ(ji);
384 >    }
385    }
408 – }  
386  
387 < void Integrator::constrainA(){
387 >  if (nConstrained){
388 >    constrainB();
389 >  }
390 > }
391  
392 <  int i,j,k;
392 > template<typename T> void Integrator<T>::preMove(void){
393 >  int i, j;
394 >  double pos[3];
395 >
396 >  if (nConstrained){
397 >    for (i = 0; i < nAtoms; i++){
398 >      atoms[i]->getPos(pos);
399 >
400 >      for (j = 0; j < 3; j++){
401 >        oldPos[3 * i + j] = pos[j];
402 >      }
403 >    }
404 >  }
405 > }
406 >
407 > template<typename T> void Integrator<T>::constrainA(){
408 >  int i, j, k;
409    int done;
410 <  double pxab, pyab, pzab;
411 <  double rxab, ryab, rzab;
410 >  double posA[3], posB[3];
411 >  double velA[3], velB[3];
412 >  double pab[3];
413 >  double rab[3];
414    int a, b, ax, ay, az, bx, by, bz;
415    double rma, rmb;
416    double dx, dy, dz;
# Line 422 | Line 420 | void Integrator::constrainA(){
420    double gab;
421    int iteration;
422  
423 <
426 <  
427 <  for( i=0; i<nAtoms; i++){
428 <    
423 >  for (i = 0; i < nAtoms; i++){
424      moving[i] = 0;
425 <    moved[i]  = 1;
425 >    moved[i] = 1;
426    }
427  
428    iteration = 0;
429    done = 0;
430 <  while( !done && (iteration < maxIteration )){
436 <
430 >  while (!done && (iteration < maxIteration)){
431      done = 1;
432 <    for(i=0; i<nConstrained; i++){
439 <
432 >    for (i = 0; i < nConstrained; i++){
433        a = constrainedA[i];
434        b = constrainedB[i];
442 –      
443 –      ax = (a*3) + 0;
444 –      ay = (a*3) + 1;
445 –      az = (a*3) + 2;
435  
436 <      bx = (b*3) + 0;
437 <      by = (b*3) + 1;
438 <      bz = (b*3) + 2;
436 >      ax = (a * 3) + 0;
437 >      ay = (a * 3) + 1;
438 >      az = (a * 3) + 2;
439  
440 <      if( moved[a] || moved[b] ){
441 <        
442 <        pxab = pos[ax] - pos[bx];
454 <        pyab = pos[ay] - pos[by];
455 <        pzab = pos[az] - pos[bz];
440 >      bx = (b * 3) + 0;
441 >      by = (b * 3) + 1;
442 >      bz = (b * 3) + 2;
443  
444 <        //periodic boundary condition
445 <        pxab = pxab - info->box_x * copysign(1, pxab)
446 <          * (int)( fabs(pxab / info->box_x) + 0.5);
460 <        pyab = pyab - info->box_y * copysign(1, pyab)
461 <          * (int)( fabs(pyab / info->box_y) + 0.5);
462 <        pzab = pzab - info->box_z * copysign(1, pzab)
463 <          * (int)( fabs(pzab / info->box_z) + 0.5);
444 >      if (moved[a] || moved[b]){
445 >        atoms[a]->getPos(posA);
446 >        atoms[b]->getPos(posB);
447  
448 <        pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
448 >        for (j = 0; j < 3; j++)
449 >          pab[j] = posA[j] - posB[j];
450  
451 <        rabsq = constrainedDsqr[i];
468 <        diffsq = rabsq - pabsq;
451 >        //periodic boundary condition
452  
453 <        // the original rattle code from alan tidesley
471 <        if (fabs(diffsq) > (tol*rabsq*2)) {
472 <          rxab = oldPos[ax] - oldPos[bx];
473 <          ryab = oldPos[ay] - oldPos[by];
474 <          rzab = oldPos[az] - oldPos[bz];
453 >        info->wrapVector(pab);
454  
455 <          rxab = rxab - info->box_x * copysign(1, rxab)
477 <            * (int)( fabs(rxab / info->box_x) + 0.5);
478 <          ryab = ryab - info->box_y * copysign(1, ryab)
479 <            * (int)( fabs(ryab / info->box_y) + 0.5);
480 <          rzab = rzab - info->box_z * copysign(1, rzab)
481 <            * (int)( fabs(rzab / info->box_z) + 0.5);
455 >        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
456  
457 <          rpab = rxab * pxab + ryab * pyab + rzab * pzab;
457 >        rabsq = constrainedDsqr[i];
458 >        diffsq = rabsq - pabsq;
459  
460 <          rpabsq = rpab * rpab;
460 >        // the original rattle code from alan tidesley
461 >        if (fabs(diffsq) > (tol * rabsq * 2)){
462 >          rab[0] = oldPos[ax] - oldPos[bx];
463 >          rab[1] = oldPos[ay] - oldPos[by];
464 >          rab[2] = oldPos[az] - oldPos[bz];
465  
466 +          info->wrapVector(rab);
467  
468 <          if (rpabsq < (rabsq * -diffsq)){
468 >          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
469  
470 +          rpabsq = rpab * rpab;
471 +
472 +
473 +          if (rpabsq < (rabsq * -diffsq)){
474   #ifdef IS_MPI
475 <            a = atoms[a]->getGlobalIndex();
476 <            b = atoms[b]->getGlobalIndex();
475 >            a = atoms[a]->getGlobalIndex();
476 >            b = atoms[b]->getGlobalIndex();
477   #endif //is_mpi
478 <            sprintf( painCave.errMsg,
479 <                     "Constraint failure in constrainA at atom %d and %d.\n",
480 <                     a, b );
481 <            painCave.isFatal = 1;
482 <            simError();
483 <          }
478 >            sprintf(painCave.errMsg,
479 >                    "Constraint failure in constrainA at atom %d and %d.\n", a,
480 >                    b);
481 >            painCave.isFatal = 1;
482 >            simError();
483 >          }
484  
485 <          rma = 1.0 / atoms[a]->getMass();
486 <          rmb = 1.0 / atoms[b]->getMass();
485 >          rma = 1.0 / atoms[a]->getMass();
486 >          rmb = 1.0 / atoms[b]->getMass();
487  
488 <          gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
488 >          gab = diffsq / (2.0 * (rma + rmb) * rpab);
489  
490 <          dx = rxab * gab;
491 <          dy = ryab * gab;
492 <          dz = rzab * gab;
490 >          dx = rab[0] * gab;
491 >          dy = rab[1] * gab;
492 >          dz = rab[2] * gab;
493  
494 <          pos[ax] += rma * dx;
495 <          pos[ay] += rma * dy;
496 <          pos[az] += rma * dz;
494 >          posA[0] += rma * dx;
495 >          posA[1] += rma * dy;
496 >          posA[2] += rma * dz;
497  
498 <          pos[bx] -= rmb * dx;
515 <          pos[by] -= rmb * dy;
516 <          pos[bz] -= rmb * dz;
498 >          atoms[a]->setPos(posA);
499  
500 +          posB[0] -= rmb * dx;
501 +          posB[1] -= rmb * dy;
502 +          posB[2] -= rmb * dz;
503 +
504 +          atoms[b]->setPos(posB);
505 +
506            dx = dx / dt;
507            dy = dy / dt;
508            dz = dz / dt;
509  
510 <          vel[ax] += rma * dx;
523 <          vel[ay] += rma * dy;
524 <          vel[az] += rma * dz;
510 >          atoms[a]->getVel(velA);
511  
512 <          vel[bx] -= rmb * dx;
513 <          vel[by] -= rmb * dy;
514 <          vel[bz] -= rmb * dz;
512 >          velA[0] += rma * dx;
513 >          velA[1] += rma * dy;
514 >          velA[2] += rma * dz;
515  
516 <          moving[a] = 1;
517 <          moving[b] = 1;
518 <          done = 0;
519 <        }
516 >          atoms[a]->setVel(velA);
517 >
518 >          atoms[b]->getVel(velB);
519 >
520 >          velB[0] -= rmb * dx;
521 >          velB[1] -= rmb * dy;
522 >          velB[2] -= rmb * dz;
523 >
524 >          atoms[b]->setVel(velB);
525 >
526 >          moving[a] = 1;
527 >          moving[b] = 1;
528 >          done = 0;
529 >        }
530        }
531      }
532 <    
533 <    for(i=0; i<nAtoms; i++){
538 <      
532 >
533 >    for (i = 0; i < nAtoms; i++){
534        moved[i] = moving[i];
535        moving[i] = 0;
536      }
# Line 543 | Line 538 | void Integrator::constrainA(){
538      iteration++;
539    }
540  
541 <  if( !done ){
542 <
543 <    sprintf( painCave.errMsg,
544 <             "Constraint failure in constrainA, too many iterations: %d\n",
550 <             iteration );
541 >  if (!done){
542 >    sprintf(painCave.errMsg,
543 >            "Constraint failure in constrainA, too many iterations: %d\n",
544 >            iteration);
545      painCave.isFatal = 1;
546      simError();
547    }
548  
549   }
550  
551 < void Integrator::constrainB( void ){
552 <  
559 <  int i,j,k;
551 > template<typename T> void Integrator<T>::constrainB(void){
552 >  int i, j, k;
553    int done;
554 +  double posA[3], posB[3];
555 +  double velA[3], velB[3];
556    double vxab, vyab, vzab;
557 <  double rxab, ryab, rzab;
557 >  double rab[3];
558    int a, b, ax, ay, az, bx, by, bz;
559    double rma, rmb;
560    double dx, dy, dz;
# Line 568 | Line 563 | void Integrator::constrainB( void ){
563    double gab;
564    int iteration;
565  
566 <  for(i=0; i<nAtoms; i++){
566 >  for (i = 0; i < nAtoms; i++){
567      moving[i] = 0;
568      moved[i] = 1;
569    }
570  
571    done = 0;
572    iteration = 0;
573 <  while( !done && (iteration < maxIteration ) ){
579 <
573 >  while (!done && (iteration < maxIteration)){
574      done = 1;
575  
576 <    for(i=0; i<nConstrained; i++){
583 <      
576 >    for (i = 0; i < nConstrained; i++){
577        a = constrainedA[i];
578        b = constrainedB[i];
579  
580 <      ax = (a*3) + 0;
581 <      ay = (a*3) + 1;
582 <      az = (a*3) + 2;
580 >      ax = (a * 3) + 0;
581 >      ay = (a * 3) + 1;
582 >      az = (a * 3) + 2;
583 >
584 >      bx = (b * 3) + 0;
585 >      by = (b * 3) + 1;
586 >      bz = (b * 3) + 2;
587  
588 <      bx = (b*3) + 0;
589 <      by = (b*3) + 1;
590 <      bz = (b*3) + 2;
588 >      if (moved[a] || moved[b]){
589 >        atoms[a]->getVel(velA);
590 >        atoms[b]->getVel(velB);
591  
592 <      if( moved[a] || moved[b] ){
593 <        
594 <        vxab = vel[ax] - vel[bx];
598 <        vyab = vel[ay] - vel[by];
599 <        vzab = vel[az] - vel[bz];
592 >        vxab = velA[0] - velB[0];
593 >        vyab = velA[1] - velB[1];
594 >        vzab = velA[2] - velB[2];
595  
596 <        rxab = pos[ax] - pos[bx];
597 <        ryab = pos[ay] - pos[by];
603 <        rzab = pos[az] - pos[bz];
604 <        
596 >        atoms[a]->getPos(posA);
597 >        atoms[b]->getPos(posB);
598  
599 <        rxab = rxab - info->box_x * copysign(1, rxab)
600 <          * (int)( fabs(rxab / info->box_x) + 0.5);
608 <        ryab = ryab - info->box_y * copysign(1, ryab)
609 <          * (int)( fabs(ryab / info->box_y) + 0.5);
610 <        rzab = rzab - info->box_z * copysign(1, rzab)
611 <          * (int)( fabs(rzab / info->box_z) + 0.5);
612 <        
613 <        rma = 1.0 / atoms[a]->getMass();
614 <        rmb = 1.0 / atoms[b]->getMass();
599 >        for (j = 0; j < 3; j++)
600 >          rab[j] = posA[j] - posB[j];
601  
602 <        rvab = rxab * vxab + ryab * vyab + rzab * vzab;
617 <          
618 <        gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
602 >        info->wrapVector(rab);
603  
604 <        if (fabs(gab) > tol) {
605 <          
622 <          dx = rxab * gab;
623 <          dy = ryab * gab;
624 <          dz = rzab * gab;
625 <          
626 <          vel[ax] += rma * dx;
627 <          vel[ay] += rma * dy;
628 <          vel[az] += rma * dz;
604 >        rma = 1.0 / atoms[a]->getMass();
605 >        rmb = 1.0 / atoms[b]->getMass();
606  
607 <          vel[bx] -= rmb * dx;
608 <          vel[by] -= rmb * dy;
609 <          vel[bz] -= rmb * dz;
610 <          
611 <          moving[a] = 1;
612 <          moving[b] = 1;
613 <          done = 0;
614 <        }
607 >        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
608 >
609 >        gab = -rvab / ((rma + rmb) * constrainedDsqr[i]);
610 >
611 >        if (fabs(gab) > tol){
612 >          dx = rab[0] * gab;
613 >          dy = rab[1] * gab;
614 >          dz = rab[2] * gab;
615 >
616 >          velA[0] += rma * dx;
617 >          velA[1] += rma * dy;
618 >          velA[2] += rma * dz;
619 >
620 >          atoms[a]->setVel(velA);
621 >
622 >          velB[0] -= rmb * dx;
623 >          velB[1] -= rmb * dy;
624 >          velB[2] -= rmb * dz;
625 >
626 >          atoms[b]->setVel(velB);
627 >
628 >          moving[a] = 1;
629 >          moving[b] = 1;
630 >          done = 0;
631 >        }
632        }
633      }
634  
635 <    for(i=0; i<nAtoms; i++){
635 >    for (i = 0; i < nAtoms; i++){
636        moved[i] = moving[i];
637        moving[i] = 0;
638      }
639 <    
639 >
640      iteration++;
641    }
642  
643 <  if( !done ){
644 <
645 <  
646 <    sprintf( painCave.errMsg,
653 <             "Constraint failure in constrainB, too many iterations: %d\n",
654 <             iteration );
643 >  if (!done){
644 >    sprintf(painCave.errMsg,
645 >            "Constraint failure in constrainB, too many iterations: %d\n",
646 >            iteration);
647      painCave.isFatal = 1;
648      simError();
649 <  }
658 <
649 >  }
650   }
651  
652 + template<typename T> void Integrator<T>::rotationPropagation
653 + ( DirectionalAtom* dAtom, double ji[3] ){
654  
655 +  double angle;
656 +  double A[3][3], I[3][3];
657  
658 +  // use the angular velocities to propagate the rotation matrix a
659 +  // full time step
660  
661 +  dAtom->getA(A);
662 +  dAtom->getI(I);
663 +  
664 +  // rotate about the x-axis      
665 +  angle = dt2 * ji[0] / I[0][0];
666 +  this->rotate( 1, 2, angle, ji, A );
667 +  
668 +  // rotate about the y-axis
669 +  angle = dt2 * ji[1] / I[1][1];
670 +  this->rotate( 2, 0, angle, ji, A );
671 +  
672 +  // rotate about the z-axis
673 +  angle = dt * ji[2] / I[2][2];
674 +  this->rotate( 0, 1, angle, ji, A);
675 +  
676 +  // rotate about the y-axis
677 +  angle = dt2 * ji[1] / I[1][1];
678 +  this->rotate( 2, 0, angle, ji, A );
679 +  
680 +  // rotate about the x-axis
681 +  angle = dt2 * ji[0] / I[0][0];
682 +  this->rotate( 1, 2, angle, ji, A );
683 +  
684 +  dAtom->setA( A  );    
685 + }
686  
687 <
688 <
689 < void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
690 <                         double A[9] ){
669 <
670 <  int i,j,k;
687 > template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
688 >                                                double angle, double ji[3],
689 >                                                double A[3][3]){
690 >  int i, j, k;
691    double sinAngle;
692    double cosAngle;
693    double angleSqr;
# Line 679 | Line 699 | void Integrator::rotate( int axes1, int axes2, double
699  
700    // initialize the tempA
701  
702 <  for(i=0; i<3; i++){
703 <    for(j=0; j<3; j++){
704 <      tempA[j][i] = A[3*i + j];
702 >  for (i = 0; i < 3; i++){
703 >    for (j = 0; j < 3; j++){
704 >      tempA[j][i] = A[i][j];
705      }
706    }
707  
708    // initialize the tempJ
709  
710 <  for( i=0; i<3; i++) tempJ[i] = ji[i];
711 <  
710 >  for (i = 0; i < 3; i++)
711 >    tempJ[i] = ji[i];
712 >
713    // initalize rot as a unit matrix
714  
715    rot[0][0] = 1.0;
# Line 698 | Line 719 | void Integrator::rotate( int axes1, int axes2, double
719    rot[1][0] = 0.0;
720    rot[1][1] = 1.0;
721    rot[1][2] = 0.0;
722 <  
722 >
723    rot[2][0] = 0.0;
724    rot[2][1] = 0.0;
725    rot[2][2] = 1.0;
726 <  
726 >
727    // use a small angle aproximation for sin and cosine
728  
729 <  angleSqr  = angle * angle;
729 >  angleSqr = angle * angle;
730    angleSqrOver4 = angleSqr / 4.0;
731    top = 1.0 - angleSqrOver4;
732    bottom = 1.0 + angleSqrOver4;
# Line 718 | Line 739 | void Integrator::rotate( int axes1, int axes2, double
739  
740    rot[axes1][axes2] = sinAngle;
741    rot[axes2][axes1] = -sinAngle;
742 <  
742 >
743    // rotate the momentum acoording to: ji[] = rot[][] * ji[]
744 <  
745 <  for(i=0; i<3; i++){
744 >
745 >  for (i = 0; i < 3; i++){
746      ji[i] = 0.0;
747 <    for(k=0; k<3; k++){
747 >    for (k = 0; k < 3; k++){
748        ji[i] += rot[i][k] * tempJ[k];
749      }
750    }
# Line 736 | Line 757 | void Integrator::rotate( int axes1, int axes2, double
757    // calculation as:
758    //                transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
759  
760 <  for(i=0; i<3; i++){
761 <    for(j=0; j<3; j++){
762 <      A[3*j + i] = 0.0;
763 <      for(k=0; k<3; k++){
764 <        A[3*j + i] += tempA[i][k] * rot[j][k];
760 >  for (i = 0; i < 3; i++){
761 >    for (j = 0; j < 3; j++){
762 >      A[j][i] = 0.0;
763 >      for (k = 0; k < 3; k++){
764 >        A[j][i] += tempA[i][k] * rot[j][k];
765        }
766      }
767    }
768   }
769 +
770 + template<typename T> void Integrator<T>::calcForce(int calcPot, int calcStress){
771 +  myFF->doForces(calcPot, calcStress);
772 + }
773 +
774 + template<typename T> void Integrator<T>::thermalize(){
775 +  tStats->velocitize();
776 + }
777 +
778 + template<typename T> double Integrator<T>::getConservedQuantity(void){
779 +  return tStats->getTotalE();
780 + }

Diff Legend

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