ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Integrator.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Integrator.cpp (file contents):
Revision 567 by mmeineke, Wed Jun 25 21:12:14 2003 UTC vs.
Revision 1127 by tim, Tue Apr 20 16:56:40 2004 UTC

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

Diff Legend

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