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 1108 by tim, Wed Apr 14 15:37:41 2004 UTC

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

Diff Legend

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