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 1234 by tim, Fri Jun 4 03:15:31 2004 UTC

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

Diff Legend

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