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

Diff Legend

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