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 637 by gezelter, Thu Jul 17 21:50:01 2003 UTC vs.
Revision 733 by tim, Wed Aug 27 19:23:29 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){
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 <        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 <      }
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 );
164 >  tStats = new Thermo(info);
165 >  statOut = new StatWriter(info);
166 >  dumpOut = new DumpWriter(info);
167  
168    atoms = info->atoms;
169    DirectionalAtom* dAtom;
# Line 184 | Line 173 | void Integrator::integrate( void ){
173  
174    // initialize the forces before the first step
175  
176 <  myFF->doForces(1,1);
176 >  calcForce(1, 1);
177    
178 <  if( info->setTemp ){
179 <    
191 <    tStats->velocitize();
178 >  if (info->setTemp){
179 >    thermalize();
180    }
181 +
182 +  calcPot = 0;
183 +  calcStress = 0;
184 +  currSample = sampleTime;
185 +  currThermal = thermalTime;
186 +  currStatus = statusTime;
187    
194 –  dumpOut->writeDump( 0.0 );
195 –  statOut->writeStat( 0.0 );
196 –  
188    calcPot     = 0;
189    calcStress  = 0;
190 <  currSample  = sampleTime;
191 <  currThermal = thermalTime;
192 <  currStatus  = statusTime;
202 <  currTime    = 0.0;;
190 >  currSample  = sampleTime + info->getTime();
191 >  currThermal = thermalTime+ info->getTime();
192 >  currStatus  = statusTime + info->getTime();
193  
194 +  dumpOut->writeDump(info->getTime());
195 +  statOut->writeStat(info->getTime());
196  
197    readyCheck();
198  
199   #ifdef IS_MPI
200 <  strcpy( checkPointMsg,
209 <          "The integrator is ready to go." );
200 >  strcpy(checkPointMsg, "The integrator is ready to go.");
201    MPIcheckPoint();
202   #endif // is_mpi
203  
204 <  while( currTime < runTime ){
205 <
215 <    if( (currTime+dt) >= currStatus ){
204 >  while (info->getTime() < runTime){
205 >    if ((info->getTime() + dt) >= currStatus){
206        calcPot = 1;
207        calcStress = 1;
208      }
209  
210 <    integrateStep( calcPot, calcStress );
221 <      
222 <    currTime += dt;
223 <    info->setTime(currTime);
210 >    integrateStep(calcPot, calcStress);
211  
212 <    if( info->setTemp ){
213 <      if( currTime >= currThermal ){
214 <        tStats->velocitize();
215 <        currThermal += thermalTime;
212 >    info->incrTime(dt);
213 >
214 >    if (info->setTemp){
215 >      if (info->getTime() >= currThermal){
216 >        thermalize();
217 >        currThermal += thermalTime;
218        }
219      }
220  
221 <    if( currTime >= currSample ){
222 <      dumpOut->writeDump( currTime );
221 >    if (info->getTime() >= currSample){
222 >      dumpOut->writeDump(info->getTime());
223        currSample += sampleTime;
224      }
225  
226 <    if( currTime >= currStatus ){
227 <      statOut->writeStat( currTime );
226 >    if (info->getTime() >= currStatus){
227 >      statOut->writeStat(info->getTime());
228        calcPot = 0;
229        calcStress = 0;
230        currStatus += statusTime;
231      }
232  
233   #ifdef IS_MPI
234 <    strcpy( checkPointMsg,
246 <            "successfully took a time step." );
234 >    strcpy(checkPointMsg, "successfully took a time step.");
235      MPIcheckPoint();
236   #endif // is_mpi
249 –
237    }
238  
239 <  dumpOut->writeFinal(currTime);
239 >  dumpOut->writeFinal(info->getTime());
240  
241    delete dumpOut;
242    delete statOut;
243   }
244  
245 < void Integrator::integrateStep( int calcPot, int calcStress ){
246 <
260 <
261 <      
245 > template<typename T> void Integrator<T>::integrateStep(int calcPot,
246 >                                                       int calcStress){
247    // Position full step, and velocity half step
263 –
248    preMove();
249 +
250    moveA();
266 –  if( nConstrained ) constrainA();
251  
252 <  
252 >  if (nConstrained){
253 >    constrainA();
254 >  }
255 >
256 >
257   #ifdef IS_MPI
258 <  strcpy( checkPointMsg, "Succesful moveA\n" );
258 >  strcpy(checkPointMsg, "Succesful moveA\n");
259    MPIcheckPoint();
260   #endif // is_mpi
273 –  
261  
262 +
263    // calc forces
264  
265 <  myFF->doForces(calcPot,calcStress);
265 >  calcForce(calcPot, calcStress);
266  
267   #ifdef IS_MPI
268 <  strcpy( checkPointMsg, "Succesful doForces\n" );
268 >  strcpy(checkPointMsg, "Succesful doForces\n");
269    MPIcheckPoint();
270   #endif // is_mpi
283 –  
271  
272 +
273    // finish the velocity  half step
274 <  
274 >
275    moveB();
276 <  if( nConstrained ) constrainB();
277 <  
276 >
277 >  if (nConstrained){
278 >    constrainB();
279 >  }
280 >
281   #ifdef IS_MPI
282 <  strcpy( checkPointMsg, "Succesful moveB\n" );
282 >  strcpy(checkPointMsg, "Succesful moveB\n");
283    MPIcheckPoint();
284   #endif // is_mpi
294 –  
295 –
285   }
286  
287  
288 < void Integrator::moveA( void ){
300 <  
288 > template<typename T> void Integrator<T>::moveA(void){
289    int i, j;
290    DirectionalAtom* dAtom;
291    double Tb[3], ji[3];
# Line 306 | Line 294 | void Integrator::moveA( void ){
294    double vel[3], pos[3], frc[3];
295    double mass;
296  
297 <  for( i=0; i<nAtoms; i++ ){
297 >  for (i = 0; i < nAtoms; i++){
298 >    atoms[i]->getVel(vel);
299 >    atoms[i]->getPos(pos);
300 >    atoms[i]->getFrc(frc);
301  
311 –    atoms[i]->getVel( vel );
312 –    atoms[i]->getPos( pos );
313 –    atoms[i]->getFrc( frc );
314 –
302      mass = atoms[i]->getMass();
303  
304 <    for (j=0; j < 3; j++) {
304 >    for (j = 0; j < 3; j++){
305        // velocity half step
306 <      vel[j] += ( dt2 * frc[j] / mass ) * eConvert;
306 >      vel[j] += (dt2 * frc[j] / mass) * eConvert;
307        // position whole step
308        pos[j] += dt * vel[j];
309      }
310  
311 <    atoms[i]->setVel( vel );
312 <    atoms[i]->setPos( pos );
311 >    atoms[i]->setVel(vel);
312 >    atoms[i]->setPos(pos);
313  
314 <    if( atoms[i]->isDirectional() ){
314 >    if (atoms[i]->isDirectional()){
315 >      dAtom = (DirectionalAtom *) atoms[i];
316  
329 –      dAtom = (DirectionalAtom *)atoms[i];
330 –          
317        // get and convert the torque to body frame
332 –      
333 –      dAtom->getTrq( Tb );
334 –      dAtom->lab2Body( Tb );
318  
319 +      dAtom->getTrq(Tb);
320 +      dAtom->lab2Body(Tb);
321 +
322        // get the angular momentum, and propagate a half step
323  
324 <      dAtom->getJ( ji );
324 >      dAtom->getJ(ji);
325  
326 <      for (j=0; j < 3; j++)
326 >      for (j = 0; j < 3; j++)
327          ji[j] += (dt2 * Tb[j]) * eConvert;
328 <      
328 >
329        // use the angular velocities to propagate the rotation matrix a
330        // full time step
331  
332        dAtom->getA(A);
333        dAtom->getI(I);
334 <    
334 >
335        // rotate about the x-axis      
336        angle = dt2 * ji[0] / I[0][0];
337 <      this->rotate( 1, 2, angle, ji, A );
337 >      this->rotate(1, 2, angle, ji, A);
338  
339        // rotate about the y-axis
340        angle = dt2 * ji[1] / I[1][1];
341 <      this->rotate( 2, 0, angle, ji, A );
342 <      
341 >      this->rotate(2, 0, angle, ji, A);
342 >
343        // rotate about the z-axis
344        angle = dt * ji[2] / I[2][2];
345 <      this->rotate( 0, 1, angle, ji, A);
346 <      
345 >      this->rotate(0, 1, angle, ji, A);
346 >
347        // rotate about the y-axis
348        angle = dt2 * ji[1] / I[1][1];
349 <      this->rotate( 2, 0, angle, ji, A );
350 <      
351 <       // rotate about the x-axis
349 >      this->rotate(2, 0, angle, ji, A);
350 >
351 >      // rotate about the x-axis
352        angle = dt2 * ji[0] / I[0][0];
353 <      this->rotate( 1, 2, angle, ji, A );
368 <      
353 >      this->rotate(1, 2, angle, ji, A);
354  
355 <      dAtom->setJ( ji );
356 <      dAtom->setA( A  );
357 <          
358 <    }    
355 >
356 >      dAtom->setJ(ji);
357 >      dAtom->setA(A);
358 >    }
359    }
360   }
361  
362  
363 < void Integrator::moveB( void ){
363 > template<typename T> void Integrator<T>::moveB(void){
364    int i, j;
365    DirectionalAtom* dAtom;
366    double Tb[3], ji[3];
367    double vel[3], frc[3];
368    double mass;
369  
370 <  for( i=0; i<nAtoms; i++ ){
371 <
372 <    atoms[i]->getVel( vel );
388 <    atoms[i]->getFrc( frc );
370 >  for (i = 0; i < nAtoms; i++){
371 >    atoms[i]->getVel(vel);
372 >    atoms[i]->getFrc(frc);
373  
374      mass = atoms[i]->getMass();
375  
376      // velocity half step
377 <    for (j=0; j < 3; j++)
378 <      vel[j] += ( dt2 * frc[j] / mass ) * eConvert;
395 <    
396 <    atoms[i]->setVel( vel );
397 <
398 <    if( atoms[i]->isDirectional() ){
377 >    for (j = 0; j < 3; j++)
378 >      vel[j] += (dt2 * frc[j] / mass) * eConvert;
379  
380 <      dAtom = (DirectionalAtom *)atoms[i];
380 >    atoms[i]->setVel(vel);
381  
382 +    if (atoms[i]->isDirectional()){
383 +      dAtom = (DirectionalAtom *) atoms[i];
384 +
385        // get and convert the torque to body frame      
386  
387 <      dAtom->getTrq( Tb );
388 <      dAtom->lab2Body( Tb );
387 >      dAtom->getTrq(Tb);
388 >      dAtom->lab2Body(Tb);
389  
390        // get the angular momentum, and propagate a half step
391  
392 <      dAtom->getJ( ji );
392 >      dAtom->getJ(ji);
393  
394 <      for (j=0; j < 3; j++)
394 >      for (j = 0; j < 3; j++)
395          ji[j] += (dt2 * Tb[j]) * eConvert;
413 –      
396  
397 <      dAtom->setJ( ji );
397 >
398 >      dAtom->setJ(ji);
399      }
400    }
401   }
402  
403 < void Integrator::preMove( void ){
403 > template<typename T> void Integrator<T>::preMove(void){
404    int i, j;
405    double pos[3];
406  
407 <  if( nConstrained ){
407 >  if (nConstrained){
408 >    for (i = 0; i < nAtoms; i++){
409 >      atoms[i]->getPos(pos);
410  
411 <    for(i=0; i < nAtoms; i++) {
412 <
428 <      atoms[i]->getPos( pos );
429 <
430 <      for (j = 0; j < 3; j++) {        
431 <        oldPos[3*i + j] = pos[j];
411 >      for (j = 0; j < 3; j++){
412 >        oldPos[3 * i + j] = pos[j];
413        }
433 –
414      }
415 <  }  
415 >  }
416   }
417  
418 < void Integrator::constrainA(){
419 <
440 <  int i,j,k;
418 > template<typename T> void Integrator<T>::constrainA(){
419 >  int i, j, k;
420    int done;
421    double posA[3], posB[3];
422    double velA[3], velB[3];
# Line 452 | Line 431 | void Integrator::constrainA(){
431    double gab;
432    int iteration;
433  
434 <  for( i=0; i<nAtoms; i++){    
434 >  for (i = 0; i < nAtoms; i++){
435      moving[i] = 0;
436 <    moved[i]  = 1;
436 >    moved[i] = 1;
437    }
438  
439    iteration = 0;
440    done = 0;
441 <  while( !done && (iteration < maxIteration )){
463 <
441 >  while (!done && (iteration < maxIteration)){
442      done = 1;
443 <    for(i=0; i<nConstrained; i++){
466 <
443 >    for (i = 0; i < nConstrained; i++){
444        a = constrainedA[i];
445        b = constrainedB[i];
469 –      
470 –      ax = (a*3) + 0;
471 –      ay = (a*3) + 1;
472 –      az = (a*3) + 2;
446  
447 <      bx = (b*3) + 0;
448 <      by = (b*3) + 1;
449 <      bz = (b*3) + 2;
447 >      ax = (a * 3) + 0;
448 >      ay = (a * 3) + 1;
449 >      az = (a * 3) + 2;
450  
451 <      if( moved[a] || moved[b] ){
452 <        
453 <        atoms[a]->getPos( posA );
454 <        atoms[b]->getPos( posB );
455 <        
456 <        for (j = 0; j < 3; j++ )
451 >      bx = (b * 3) + 0;
452 >      by = (b * 3) + 1;
453 >      bz = (b * 3) + 2;
454 >
455 >      if (moved[a] || moved[b]){
456 >        atoms[a]->getPos(posA);
457 >        atoms[b]->getPos(posB);
458 >
459 >        for (j = 0; j < 3; j++)
460            pab[j] = posA[j] - posB[j];
485 –        
486 –        //periodic boundary condition
461  
462 <        info->wrapVector( pab );
462 >        //periodic boundary condition
463  
464 <        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
464 >        info->wrapVector(pab);
465  
466 <        rabsq = constrainedDsqr[i];
493 <        diffsq = rabsq - pabsq;
466 >        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
467  
468 <        // the original rattle code from alan tidesley
469 <        if (fabs(diffsq) > (tol*rabsq*2)) {
497 <          rab[0] = oldPos[ax] - oldPos[bx];
498 <          rab[1] = oldPos[ay] - oldPos[by];
499 <          rab[2] = oldPos[az] - oldPos[bz];
468 >        rabsq = constrainedDsqr[i];
469 >        diffsq = rabsq - pabsq;
470  
471 <          info->wrapVector( rab );
471 >        // the original rattle code from alan tidesley
472 >        if (fabs(diffsq) > (tol * rabsq * 2)){
473 >          rab[0] = oldPos[ax] - oldPos[bx];
474 >          rab[1] = oldPos[ay] - oldPos[by];
475 >          rab[2] = oldPos[az] - oldPos[bz];
476  
477 <          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
477 >          info->wrapVector(rab);
478  
479 <          rpabsq = rpab * rpab;
479 >          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
480  
481 +          rpabsq = rpab * rpab;
482  
508 –          if (rpabsq < (rabsq * -diffsq)){
483  
484 +          if (rpabsq < (rabsq * -diffsq)){
485   #ifdef IS_MPI
486 <            a = atoms[a]->getGlobalIndex();
487 <            b = atoms[b]->getGlobalIndex();
486 >            a = atoms[a]->getGlobalIndex();
487 >            b = atoms[b]->getGlobalIndex();
488   #endif //is_mpi
489 <            sprintf( painCave.errMsg,
490 <                     "Constraint failure in constrainA at atom %d and %d.\n",
491 <                     a, b );
492 <            painCave.isFatal = 1;
493 <            simError();
494 <          }
489 >            sprintf(painCave.errMsg,
490 >                    "Constraint failure in constrainA at atom %d and %d.\n", a,
491 >                    b);
492 >            painCave.isFatal = 1;
493 >            simError();
494 >          }
495  
496 <          rma = 1.0 / atoms[a]->getMass();
497 <          rmb = 1.0 / atoms[b]->getMass();
496 >          rma = 1.0 / atoms[a]->getMass();
497 >          rmb = 1.0 / atoms[b]->getMass();
498  
499 <          gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
499 >          gab = diffsq / (2.0 * (rma + rmb) * rpab);
500  
501            dx = rab[0] * gab;
502            dy = rab[1] * gab;
503            dz = rab[2] * gab;
504  
505 <          posA[0] += rma * dx;
506 <          posA[1] += rma * dy;
507 <          posA[2] += rma * dz;
505 >          posA[0] += rma * dx;
506 >          posA[1] += rma * dy;
507 >          posA[2] += rma * dz;
508  
509 <          atoms[a]->setPos( posA );
509 >          atoms[a]->setPos(posA);
510  
511 <          posB[0] -= rmb * dx;
512 <          posB[1] -= rmb * dy;
513 <          posB[2] -= rmb * dz;
511 >          posB[0] -= rmb * dx;
512 >          posB[1] -= rmb * dy;
513 >          posB[2] -= rmb * dz;
514  
515 <          atoms[b]->setPos( posB );
515 >          atoms[b]->setPos(posB);
516  
517            dx = dx / dt;
518            dy = dy / dt;
519            dz = dz / dt;
520  
521 <          atoms[a]->getVel( velA );
521 >          atoms[a]->getVel(velA);
522  
523 <          velA[0] += rma * dx;
524 <          velA[1] += rma * dy;
525 <          velA[2] += rma * dz;
523 >          velA[0] += rma * dx;
524 >          velA[1] += rma * dy;
525 >          velA[2] += rma * dz;
526  
527 <          atoms[a]->setVel( velA );
527 >          atoms[a]->setVel(velA);
528  
529 <          atoms[b]->getVel( velB );
529 >          atoms[b]->getVel(velB);
530  
531 <          velB[0] -= rmb * dx;
532 <          velB[1] -= rmb * dy;
533 <          velB[2] -= rmb * dz;
531 >          velB[0] -= rmb * dx;
532 >          velB[1] -= rmb * dy;
533 >          velB[2] -= rmb * dz;
534  
535 <          atoms[b]->setVel( velB );
535 >          atoms[b]->setVel(velB);
536  
537 <          moving[a] = 1;
538 <          moving[b] = 1;
539 <          done = 0;
540 <        }
537 >          moving[a] = 1;
538 >          moving[b] = 1;
539 >          done = 0;
540 >        }
541        }
542      }
543 <    
544 <    for(i=0; i<nAtoms; i++){
570 <      
543 >
544 >    for (i = 0; i < nAtoms; i++){
545        moved[i] = moving[i];
546        moving[i] = 0;
547      }
# Line 575 | Line 549 | void Integrator::constrainA(){
549      iteration++;
550    }
551  
552 <  if( !done ){
553 <
554 <    sprintf( painCave.errMsg,
555 <             "Constraint failure in constrainA, too many iterations: %d\n",
582 <             iteration );
552 >  if (!done){
553 >    sprintf(painCave.errMsg,
554 >            "Constraint failure in constrainA, too many iterations: %d\n",
555 >            iteration);
556      painCave.isFatal = 1;
557      simError();
558    }
586 –
559   }
560  
561 < void Integrator::constrainB( void ){
562 <  
591 <  int i,j,k;
561 > template<typename T> void Integrator<T>::constrainB(void){
562 >  int i, j, k;
563    int done;
564    double posA[3], posB[3];
565    double velA[3], velB[3];
# Line 602 | Line 573 | void Integrator::constrainB( void ){
573    double gab;
574    int iteration;
575  
576 <  for(i=0; i<nAtoms; i++){
576 >  for (i = 0; i < nAtoms; i++){
577      moving[i] = 0;
578      moved[i] = 1;
579    }
580  
581    done = 0;
582    iteration = 0;
583 <  while( !done && (iteration < maxIteration ) ){
613 <
583 >  while (!done && (iteration < maxIteration)){
584      done = 1;
585  
586 <    for(i=0; i<nConstrained; i++){
617 <      
586 >    for (i = 0; i < nConstrained; i++){
587        a = constrainedA[i];
588        b = constrainedB[i];
589  
590 <      ax = (a*3) + 0;
591 <      ay = (a*3) + 1;
592 <      az = (a*3) + 2;
590 >      ax = (a * 3) + 0;
591 >      ay = (a * 3) + 1;
592 >      az = (a * 3) + 2;
593  
594 <      bx = (b*3) + 0;
595 <      by = (b*3) + 1;
596 <      bz = (b*3) + 2;
628 <
629 <      if( moved[a] || moved[b] ){
594 >      bx = (b * 3) + 0;
595 >      by = (b * 3) + 1;
596 >      bz = (b * 3) + 2;
597  
598 <        atoms[a]->getVel( velA );
599 <        atoms[b]->getVel( velB );
600 <          
634 <        vxab = velA[0] - velB[0];
635 <        vyab = velA[1] - velB[1];
636 <        vzab = velA[2] - velB[2];
598 >      if (moved[a] || moved[b]){
599 >        atoms[a]->getVel(velA);
600 >        atoms[b]->getVel(velB);
601  
602 <        atoms[a]->getPos( posA );
603 <        atoms[b]->getPos( posB );
602 >        vxab = velA[0] - velB[0];
603 >        vyab = velA[1] - velB[1];
604 >        vzab = velA[2] - velB[2];
605  
606 <        for (j = 0; j < 3; j++)
606 >        atoms[a]->getPos(posA);
607 >        atoms[b]->getPos(posB);
608 >
609 >        for (j = 0; j < 3; j++)
610            rab[j] = posA[j] - posB[j];
643 –          
644 –        info->wrapVector( rab );
645 –        
646 –        rma = 1.0 / atoms[a]->getMass();
647 –        rmb = 1.0 / atoms[b]->getMass();
611  
612 <        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
650 <          
651 <        gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
612 >        info->wrapVector(rab);
613  
614 <        if (fabs(gab) > tol) {
615 <          
655 <          dx = rab[0] * gab;
656 <          dy = rab[1] * gab;
657 <          dz = rab[2] * gab;
658 <        
659 <          velA[0] += rma * dx;
660 <          velA[1] += rma * dy;
661 <          velA[2] += rma * dz;
614 >        rma = 1.0 / atoms[a]->getMass();
615 >        rmb = 1.0 / atoms[b]->getMass();
616  
617 <          atoms[a]->setVel( velA );
617 >        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
618  
619 <          velB[0] -= rmb * dx;
666 <          velB[1] -= rmb * dy;
667 <          velB[2] -= rmb * dz;
619 >        gab = -rvab / ((rma + rmb) * constrainedDsqr[i]);
620  
621 <          atoms[b]->setVel( velB );
622 <          
623 <          moving[a] = 1;
624 <          moving[b] = 1;
625 <          done = 0;
626 <        }
621 >        if (fabs(gab) > tol){
622 >          dx = rab[0] * gab;
623 >          dy = rab[1] * gab;
624 >          dz = rab[2] * gab;
625 >
626 >          velA[0] += rma * dx;
627 >          velA[1] += rma * dy;
628 >          velA[2] += rma * dz;
629 >
630 >          atoms[a]->setVel(velA);
631 >
632 >          velB[0] -= rmb * dx;
633 >          velB[1] -= rmb * dy;
634 >          velB[2] -= rmb * dz;
635 >
636 >          atoms[b]->setVel(velB);
637 >
638 >          moving[a] = 1;
639 >          moving[b] = 1;
640 >          done = 0;
641 >        }
642        }
643      }
644  
645 <    for(i=0; i<nAtoms; i++){
645 >    for (i = 0; i < nAtoms; i++){
646        moved[i] = moving[i];
647        moving[i] = 0;
648      }
649 <    
649 >
650      iteration++;
651    }
685 –  
686 –  if( !done ){
652  
653 <  
654 <    sprintf( painCave.errMsg,
655 <             "Constraint failure in constrainB, too many iterations: %d\n",
656 <             iteration );
653 >  if (!done){
654 >    sprintf(painCave.errMsg,
655 >            "Constraint failure in constrainB, too many iterations: %d\n",
656 >            iteration);
657      painCave.isFatal = 1;
658      simError();
659 <  }
695 <
659 >  }
660   }
661  
662 < void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
663 <                         double A[3][3] ){
664 <
665 <  int i,j,k;
662 > template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
663 >                                                double angle, double ji[3],
664 >                                                double A[3][3]){
665 >  int i, j, k;
666    double sinAngle;
667    double cosAngle;
668    double angleSqr;
# Line 710 | Line 674 | void Integrator::rotate( int axes1, int axes2, double
674  
675    // initialize the tempA
676  
677 <  for(i=0; i<3; i++){
678 <    for(j=0; j<3; j++){
677 >  for (i = 0; i < 3; i++){
678 >    for (j = 0; j < 3; j++){
679        tempA[j][i] = A[i][j];
680      }
681    }
682  
683    // initialize the tempJ
684  
685 <  for( i=0; i<3; i++) tempJ[i] = ji[i];
686 <  
685 >  for (i = 0; i < 3; i++)
686 >    tempJ[i] = ji[i];
687 >
688    // initalize rot as a unit matrix
689  
690    rot[0][0] = 1.0;
# Line 729 | Line 694 | void Integrator::rotate( int axes1, int axes2, double
694    rot[1][0] = 0.0;
695    rot[1][1] = 1.0;
696    rot[1][2] = 0.0;
697 <  
697 >
698    rot[2][0] = 0.0;
699    rot[2][1] = 0.0;
700    rot[2][2] = 1.0;
701 <  
701 >
702    // use a small angle aproximation for sin and cosine
703  
704 <  angleSqr  = angle * angle;
704 >  angleSqr = angle * angle;
705    angleSqrOver4 = angleSqr / 4.0;
706    top = 1.0 - angleSqrOver4;
707    bottom = 1.0 + angleSqrOver4;
# Line 749 | Line 714 | void Integrator::rotate( int axes1, int axes2, double
714  
715    rot[axes1][axes2] = sinAngle;
716    rot[axes2][axes1] = -sinAngle;
717 <  
717 >
718    // rotate the momentum acoording to: ji[] = rot[][] * ji[]
719 <  
720 <  for(i=0; i<3; i++){
719 >
720 >  for (i = 0; i < 3; i++){
721      ji[i] = 0.0;
722 <    for(k=0; k<3; k++){
722 >    for (k = 0; k < 3; k++){
723        ji[i] += rot[i][k] * tempJ[k];
724      }
725    }
# Line 767 | Line 732 | void Integrator::rotate( int axes1, int axes2, double
732    // calculation as:
733    //                transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
734  
735 <  for(i=0; i<3; i++){
736 <    for(j=0; j<3; j++){
735 >  for (i = 0; i < 3; i++){
736 >    for (j = 0; j < 3; j++){
737        A[j][i] = 0.0;
738 <      for(k=0; k<3; k++){
739 <        A[j][i] += tempA[i][k] * rot[j][k];
738 >      for (k = 0; k < 3; k++){
739 >        A[j][i] += tempA[i][k] * rot[j][k];
740        }
741      }
742    }
743   }
744 +
745 + template<typename T> void Integrator<T>::calcForce(int calcPot, int calcStress){
746 +  myFF->doForces(calcPot, calcStress);
747 + }
748 +
749 + template<typename T> void Integrator<T>::thermalize(){
750 +  tStats->velocitize();
751 + }

Diff Legend

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