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 643 by mmeineke, Mon Jul 21 21:27:40 2003 UTC vs.
Revision 829 by gezelter, Tue Oct 28 16:03:37 2003 UTC

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

Diff Legend

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