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

Diff Legend

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