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

Diff Legend

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