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

Diff Legend

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