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 594 by mmeineke, Fri Jul 11 22:34:48 2003 UTC vs.
Revision 787 by mmeineke, Thu Sep 25 19:27:15 2003 UTC

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

Diff Legend

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