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

Diff Legend

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