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 614 by mmeineke, Tue Jul 15 17:57:04 2003 UTC vs.
Revision 778 by mmeineke, Fri Sep 19 20:00:27 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 +  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 184 | 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 <    
191 <    tStats->velocitize();
181 >  if (info->setTemp){
182 >    thermalize();
183    }
184 <  
194 <  dumpOut->writeDump( 0.0 );
195 <  statOut->writeStat( 0.0 );
196 <  
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,
209 <          "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 <
215 <    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 );
221 <      
222 <    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,
245 <            "successfully took a time step." );
239 >    strcpy(checkPointMsg, "successfully took a time step.");
240      MPIcheckPoint();
241   #endif // is_mpi
248 –
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 <
259 <
260 <      
250 > template<typename T> void Integrator<T>::integrateStep(int calcPot,
251 >                                                       int calcStress){
252    // Position full step, and velocity half step
262 –
253    preMove();
254 +
255    moveA();
265 –  if( nConstrained ) constrainA();
256  
257 <  
257 >
258 >
259 >
260   #ifdef IS_MPI
261 <  strcpy( checkPointMsg, "Succesful moveA\n" );
261 >  strcpy(checkPointMsg, "Succesful moveA\n");
262    MPIcheckPoint();
263   #endif // is_mpi
272 –  
264  
265 +
266    // calc forces
267  
268 <  myFF->doForces(calcPot,calcStress);
268 >  calcForce(calcPot, calcStress);
269  
270   #ifdef IS_MPI
271 <  strcpy( checkPointMsg, "Succesful doForces\n" );
271 >  strcpy(checkPointMsg, "Succesful doForces\n");
272    MPIcheckPoint();
273   #endif // is_mpi
282 –  
274  
275 +
276    // finish the velocity  half step
277 <  
277 >
278    moveB();
279 <  if( nConstrained ) constrainB();
280 <  
279 >
280 >
281 >
282   #ifdef IS_MPI
283 <  strcpy( checkPointMsg, "Succesful moveB\n" );
283 >  strcpy(checkPointMsg, "Succesful moveB\n");
284    MPIcheckPoint();
285   #endif // is_mpi
293 –  
294 –
286   }
287  
288  
289 < void Integrator::moveA( void ){
299 <  
289 > template<typename T> void Integrator<T>::moveA(void){
290    int i, j;
291    DirectionalAtom* dAtom;
292    double Tb[3], ji[3];
303 –  double A[3][3], I[3][3];
304 –  double angle;
293    double vel[3], pos[3], frc[3];
294    double mass;
295  
296 <  for( i=0; i<nAtoms; i++ ){
296 >  for (i = 0; i < nAtoms; i++){
297 >    atoms[i]->getVel(vel);
298 >    atoms[i]->getPos(pos);
299 >    atoms[i]->getFrc(frc);
300  
310 –    atoms[i]->getVel( vel );
311 –    atoms[i]->getPos( pos );
312 –    atoms[i]->getFrc( frc );
313 –
301      mass = atoms[i]->getMass();
302  
303 <    for (j=0; j < 3; j++) {
303 >    for (j = 0; j < 3; j++){
304        // velocity half step
305 <      vel[j] += ( dt2 * frc[j] / mass ) * eConvert;
305 >      vel[j] += (dt2 * frc[j] / mass) * eConvert;
306        // position whole step
307        pos[j] += dt * vel[j];
308      }
309  
310 <    atoms[i]->setVel( vel );
311 <    atoms[i]->setPos( pos );
310 >    atoms[i]->setVel(vel);
311 >    atoms[i]->setPos(pos);
312  
313 <    if( atoms[i]->isDirectional() ){
313 >    if (atoms[i]->isDirectional()){
314 >      dAtom = (DirectionalAtom *) atoms[i];
315  
328 –      dAtom = (DirectionalAtom *)atoms[i];
329 –          
316        // get and convert the torque to body frame
331 –      
332 –      dAtom->getTrq( Tb );
333 –      dAtom->lab2Body( Tb );
317  
318 +      dAtom->getTrq(Tb);
319 +      dAtom->lab2Body(Tb);
320 +
321        // get the angular momentum, and propagate a half step
322  
323 <      dAtom->getJ( ji );
323 >      dAtom->getJ(ji);
324  
325 <      for (j=0; j < 3; j++)
325 >      for (j = 0; j < 3; j++)
326          ji[j] += (dt2 * Tb[j]) * eConvert;
341 –      
342 –      // use the angular velocities to propagate the rotation matrix a
343 –      // full time step
327  
328 <      dAtom->getA(A);
346 <      dAtom->getI(I);
347 <    
348 <      // rotate about the x-axis      
349 <      angle = dt2 * ji[0] / I[0][0];
350 <      this->rotate( 1, 2, angle, ji, A );
328 >      this->rotationPropagation( dAtom, ji );
329  
330 <      // rotate about the y-axis
331 <      angle = dt2 * ji[1] / I[1][1];
332 <      this->rotate( 2, 0, angle, ji, A );
355 <      
356 <      // rotate about the z-axis
357 <      angle = dt * ji[2] / I[2][2];
358 <      this->rotate( 0, 1, angle, ji, A);
359 <      
360 <      // rotate about the y-axis
361 <      angle = dt2 * ji[1] / I[1][1];
362 <      this->rotate( 2, 0, angle, ji, A );
363 <      
364 <       // rotate about the x-axis
365 <      angle = dt2 * ji[0] / I[0][0];
366 <      this->rotate( 1, 2, angle, ji, A );
367 <      
330 >      dAtom->setJ(ji);
331 >    }
332 >  }
333  
334 <      dAtom->setJ( ji );
335 <      dAtom->setA( A  );
371 <          
372 <    }    
334 >  if (nConstrained){
335 >    constrainA();
336    }
337   }
338  
339  
340 < void Integrator::moveB( void ){
340 > template<typename T> void Integrator<T>::moveB(void){
341    int i, j;
342    DirectionalAtom* dAtom;
343    double Tb[3], ji[3];
344    double vel[3], frc[3];
345    double mass;
346  
347 <  for( i=0; i<nAtoms; i++ ){
348 <
349 <    atoms[i]->getVel( vel );
387 <    atoms[i]->getFrc( frc );
347 >  for (i = 0; i < nAtoms; i++){
348 >    atoms[i]->getVel(vel);
349 >    atoms[i]->getFrc(frc);
350  
351      mass = atoms[i]->getMass();
352  
353      // velocity half step
354 <    for (j=0; j < 3; j++)
355 <      vel[j] += ( dt2 * frc[j] / mass ) * eConvert;
394 <    
395 <    atoms[i]->setVel( vel );
396 <
397 <    if( atoms[i]->isDirectional() ){
354 >    for (j = 0; j < 3; j++)
355 >      vel[j] += (dt2 * frc[j] / mass) * eConvert;
356  
357 <      dAtom = (DirectionalAtom *)atoms[i];
357 >    atoms[i]->setVel(vel);
358  
359 +    if (atoms[i]->isDirectional()){
360 +      dAtom = (DirectionalAtom *) atoms[i];
361 +
362        // get and convert the torque to body frame      
363  
364 <      dAtom->getTrq( Tb );
365 <      dAtom->lab2Body( Tb );
364 >      dAtom->getTrq(Tb);
365 >      dAtom->lab2Body(Tb);
366  
367        // get the angular momentum, and propagate a half step
368  
369 <      dAtom->getJ( ji );
369 >      dAtom->getJ(ji);
370  
371 <      for (j=0; j < 3; j++)
371 >      for (j = 0; j < 3; j++)
372          ji[j] += (dt2 * Tb[j]) * eConvert;
412 –      
373  
374 <      dAtom->setJ( ji );
374 >
375 >      dAtom->setJ(ji);
376      }
377    }
378 +
379 +  if (nConstrained){
380 +    constrainB();
381 +  }
382   }
383  
384 < void Integrator::preMove( void ){
384 > template<typename T> void Integrator<T>::preMove(void){
385    int i, j;
386    double pos[3];
387  
388 <  if( nConstrained ){
388 >  if (nConstrained){
389 >    for (i = 0; i < nAtoms; i++){
390 >      atoms[i]->getPos(pos);
391  
392 <    for(i=0; i < nAtoms; i++) {
393 <
427 <      atoms[i]->getPos( pos );
428 <
429 <      for (j = 0; j < 3; j++) {        
430 <        oldPos[3*i + j] = pos[j];
392 >      for (j = 0; j < 3; j++){
393 >        oldPos[3 * i + j] = pos[j];
394        }
432 –
395      }
396 <  }  
396 >  }
397   }
398  
399 < void Integrator::constrainA(){
400 <
439 <  int i,j,k;
399 > template<typename T> void Integrator<T>::constrainA(){
400 >  int i, j, k;
401    int done;
402    double posA[3], posB[3];
403    double velA[3], velB[3];
# Line 451 | Line 412 | void Integrator::constrainA(){
412    double gab;
413    int iteration;
414  
415 <  for( i=0; i<nAtoms; i++){    
415 >  for (i = 0; i < nAtoms; i++){
416      moving[i] = 0;
417 <    moved[i]  = 1;
417 >    moved[i] = 1;
418    }
419  
420    iteration = 0;
421    done = 0;
422 <  while( !done && (iteration < maxIteration )){
462 <
422 >  while (!done && (iteration < maxIteration)){
423      done = 1;
424 <    for(i=0; i<nConstrained; i++){
465 <
424 >    for (i = 0; i < nConstrained; i++){
425        a = constrainedA[i];
426        b = constrainedB[i];
468 –      
469 –      ax = (a*3) + 0;
470 –      ay = (a*3) + 1;
471 –      az = (a*3) + 2;
427  
428 <      bx = (b*3) + 0;
429 <      by = (b*3) + 1;
430 <      bz = (b*3) + 2;
428 >      ax = (a * 3) + 0;
429 >      ay = (a * 3) + 1;
430 >      az = (a * 3) + 2;
431  
432 <      if( moved[a] || moved[b] ){
433 <        
434 <        atoms[a]->getPos( posA );
435 <        atoms[b]->getPos( posB );
436 <        
437 <        for (j = 0; j < 3; j++ )
432 >      bx = (b * 3) + 0;
433 >      by = (b * 3) + 1;
434 >      bz = (b * 3) + 2;
435 >
436 >      if (moved[a] || moved[b]){
437 >        atoms[a]->getPos(posA);
438 >        atoms[b]->getPos(posB);
439 >
440 >        for (j = 0; j < 3; j++)
441            pab[j] = posA[j] - posB[j];
484 –        
485 –        //periodic boundary condition
442  
443 <        info->wrapVector( pab );
443 >        //periodic boundary condition
444  
445 <        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
445 >        info->wrapVector(pab);
446  
447 <        rabsq = constrainedDsqr[i];
492 <        diffsq = rabsq - pabsq;
447 >        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
448  
449 <        // the original rattle code from alan tidesley
450 <        if (fabs(diffsq) > (tol*rabsq*2)) {
496 <          rab[0] = oldPos[ax] - oldPos[bx];
497 <          rab[1] = oldPos[ay] - oldPos[by];
498 <          rab[2] = oldPos[az] - oldPos[bz];
449 >        rabsq = constrainedDsqr[i];
450 >        diffsq = rabsq - pabsq;
451  
452 <          info->wrapVector( rab );
452 >        // the original rattle code from alan tidesley
453 >        if (fabs(diffsq) > (tol * rabsq * 2)){
454 >          rab[0] = oldPos[ax] - oldPos[bx];
455 >          rab[1] = oldPos[ay] - oldPos[by];
456 >          rab[2] = oldPos[az] - oldPos[bz];
457  
458 <          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
458 >          info->wrapVector(rab);
459  
460 <          rpabsq = rpab * rpab;
460 >          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
461  
462 +          rpabsq = rpab * rpab;
463  
507 –          if (rpabsq < (rabsq * -diffsq)){
464  
465 +          if (rpabsq < (rabsq * -diffsq)){
466   #ifdef IS_MPI
467 <            a = atoms[a]->getGlobalIndex();
468 <            b = atoms[b]->getGlobalIndex();
467 >            a = atoms[a]->getGlobalIndex();
468 >            b = atoms[b]->getGlobalIndex();
469   #endif //is_mpi
470 <            sprintf( painCave.errMsg,
471 <                     "Constraint failure in constrainA at atom %d and %d.\n",
472 <                     a, b );
473 <            painCave.isFatal = 1;
474 <            simError();
475 <          }
470 >            sprintf(painCave.errMsg,
471 >                    "Constraint failure in constrainA at atom %d and %d.\n", a,
472 >                    b);
473 >            painCave.isFatal = 1;
474 >            simError();
475 >          }
476  
477 <          rma = 1.0 / atoms[a]->getMass();
478 <          rmb = 1.0 / atoms[b]->getMass();
477 >          rma = 1.0 / atoms[a]->getMass();
478 >          rmb = 1.0 / atoms[b]->getMass();
479  
480 <          gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
480 >          gab = diffsq / (2.0 * (rma + rmb) * rpab);
481  
482            dx = rab[0] * gab;
483            dy = rab[1] * gab;
484            dz = rab[2] * gab;
485  
486 <          posA[0] += rma * dx;
487 <          posA[1] += rma * dy;
488 <          posA[2] += rma * dz;
486 >          posA[0] += rma * dx;
487 >          posA[1] += rma * dy;
488 >          posA[2] += rma * dz;
489  
490 <          atoms[a]->setPos( posA );
490 >          atoms[a]->setPos(posA);
491  
492 <          posB[0] -= rmb * dx;
493 <          posB[1] -= rmb * dy;
494 <          posB[2] -= rmb * dz;
492 >          posB[0] -= rmb * dx;
493 >          posB[1] -= rmb * dy;
494 >          posB[2] -= rmb * dz;
495  
496 <          atoms[b]->setPos( posB );
496 >          atoms[b]->setPos(posB);
497  
498            dx = dx / dt;
499            dy = dy / dt;
500            dz = dz / dt;
501  
502 <          atoms[a]->getVel( velA );
502 >          atoms[a]->getVel(velA);
503  
504 <          velA[0] += rma * dx;
505 <          velA[1] += rma * dy;
506 <          velA[2] += rma * dz;
504 >          velA[0] += rma * dx;
505 >          velA[1] += rma * dy;
506 >          velA[2] += rma * dz;
507  
508 <          atoms[a]->setVel( velA );
508 >          atoms[a]->setVel(velA);
509  
510 <          atoms[b]->getVel( velB );
510 >          atoms[b]->getVel(velB);
511  
512 <          velB[0] -= rmb * dx;
513 <          velB[1] -= rmb * dy;
514 <          velB[2] -= rmb * dz;
512 >          velB[0] -= rmb * dx;
513 >          velB[1] -= rmb * dy;
514 >          velB[2] -= rmb * dz;
515  
516 <          atoms[b]->setVel( velB );
516 >          atoms[b]->setVel(velB);
517  
518 <          moving[a] = 1;
519 <          moving[b] = 1;
520 <          done = 0;
521 <        }
518 >          moving[a] = 1;
519 >          moving[b] = 1;
520 >          done = 0;
521 >        }
522        }
523      }
524 <    
525 <    for(i=0; i<nAtoms; i++){
569 <      
524 >
525 >    for (i = 0; i < nAtoms; i++){
526        moved[i] = moving[i];
527        moving[i] = 0;
528      }
# Line 574 | Line 530 | void Integrator::constrainA(){
530      iteration++;
531    }
532  
533 <  if( !done ){
534 <
535 <    sprintf( painCave.errMsg,
536 <             "Constraint failure in constrainA, too many iterations: %d\n",
581 <             iteration );
533 >  if (!done){
534 >    sprintf(painCave.errMsg,
535 >            "Constraint failure in constrainA, too many iterations: %d\n",
536 >            iteration);
537      painCave.isFatal = 1;
538      simError();
539    }
540  
541   }
542  
543 < void Integrator::constrainB( void ){
544 <  
590 <  int i,j,k;
543 > template<typename T> void Integrator<T>::constrainB(void){
544 >  int i, j, k;
545    int done;
546    double posA[3], posB[3];
547    double velA[3], velB[3];
# Line 601 | Line 555 | void Integrator::constrainB( void ){
555    double gab;
556    int iteration;
557  
558 <  for(i=0; i<nAtoms; i++){
558 >  for (i = 0; i < nAtoms; i++){
559      moving[i] = 0;
560      moved[i] = 1;
561    }
562  
563    done = 0;
564    iteration = 0;
565 <  while( !done && (iteration < maxIteration ) ){
612 <
565 >  while (!done && (iteration < maxIteration)){
566      done = 1;
567  
568 <    for(i=0; i<nConstrained; i++){
616 <      
568 >    for (i = 0; i < nConstrained; i++){
569        a = constrainedA[i];
570        b = constrainedB[i];
571  
572 <      ax = (a*3) + 0;
573 <      ay = (a*3) + 1;
574 <      az = (a*3) + 2;
572 >      ax = (a * 3) + 0;
573 >      ay = (a * 3) + 1;
574 >      az = (a * 3) + 2;
575  
576 <      bx = (b*3) + 0;
577 <      by = (b*3) + 1;
578 <      bz = (b*3) + 2;
576 >      bx = (b * 3) + 0;
577 >      by = (b * 3) + 1;
578 >      bz = (b * 3) + 2;
579  
580 <      if( moved[a] || moved[b] ){
580 >      if (moved[a] || moved[b]){
581 >        atoms[a]->getVel(velA);
582 >        atoms[b]->getVel(velB);
583  
584 <        atoms[a]->getVel( velA );
585 <        atoms[b]->getVel( velB );
586 <          
633 <        vxab = velA[0] - velB[0];
634 <        vyab = velA[1] - velB[1];
635 <        vzab = velA[2] - velB[2];
584 >        vxab = velA[0] - velB[0];
585 >        vyab = velA[1] - velB[1];
586 >        vzab = velA[2] - velB[2];
587  
588 <        atoms[a]->getPos( posA );
589 <        atoms[b]->getPos( posB );
588 >        atoms[a]->getPos(posA);
589 >        atoms[b]->getPos(posB);
590  
591 <        for (j = 0; j < 3; j++)
591 >        for (j = 0; j < 3; j++)
592            rab[j] = posA[j] - posB[j];
642 –          
643 –        info->wrapVector( rab );
644 –        
645 –        rma = 1.0 / atoms[a]->getMass();
646 –        rmb = 1.0 / atoms[b]->getMass();
593  
594 <        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
649 <          
650 <        gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
594 >        info->wrapVector(rab);
595  
596 <        if (fabs(gab) > tol) {
597 <          
654 <          dx = rab[0] * gab;
655 <          dy = rab[1] * gab;
656 <          dz = rab[2] * gab;
657 <        
658 <          velA[0] += rma * dx;
659 <          velA[1] += rma * dy;
660 <          velA[2] += rma * dz;
596 >        rma = 1.0 / atoms[a]->getMass();
597 >        rmb = 1.0 / atoms[b]->getMass();
598  
599 <          atoms[a]->setVel( velA );
599 >        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
600  
601 <          velB[0] -= rmb * dx;
665 <          velB[1] -= rmb * dy;
666 <          velB[2] -= rmb * dz;
601 >        gab = -rvab / ((rma + rmb) * constrainedDsqr[i]);
602  
603 <          atoms[b]->setVel( velB );
604 <          
605 <          moving[a] = 1;
606 <          moving[b] = 1;
607 <          done = 0;
608 <        }
603 >        if (fabs(gab) > tol){
604 >          dx = rab[0] * gab;
605 >          dy = rab[1] * gab;
606 >          dz = rab[2] * gab;
607 >
608 >          velA[0] += rma * dx;
609 >          velA[1] += rma * dy;
610 >          velA[2] += rma * dz;
611 >
612 >          atoms[a]->setVel(velA);
613 >
614 >          velB[0] -= rmb * dx;
615 >          velB[1] -= rmb * dy;
616 >          velB[2] -= rmb * dz;
617 >
618 >          atoms[b]->setVel(velB);
619 >
620 >          moving[a] = 1;
621 >          moving[b] = 1;
622 >          done = 0;
623 >        }
624        }
625      }
626  
627 <    for(i=0; i<nAtoms; i++){
627 >    for (i = 0; i < nAtoms; i++){
628        moved[i] = moving[i];
629        moving[i] = 0;
630      }
631 <    
631 >
632      iteration++;
633    }
684 –  
685 –  if( !done ){
634  
635 <  
636 <    sprintf( painCave.errMsg,
637 <             "Constraint failure in constrainB, too many iterations: %d\n",
638 <             iteration );
635 >  if (!done){
636 >    sprintf(painCave.errMsg,
637 >            "Constraint failure in constrainB, too many iterations: %d\n",
638 >            iteration);
639      painCave.isFatal = 1;
640      simError();
641 <  }
694 <
641 >  }
642   }
643  
644 < void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
645 <                         double A[3][3] ){
644 > template<typename T> void Integrator<T>::rotationPropagation
645 > ( DirectionalAtom* dAtom, double ji[3] ){
646  
647 <  int i,j,k;
647 >  double angle;
648 >  double A[3][3], I[3][3];
649 >
650 >  // use the angular velocities to propagate the rotation matrix a
651 >  // full time step
652 >
653 >  dAtom->getA(A);
654 >  dAtom->getI(I);
655 >  
656 >  // rotate about the x-axis      
657 >  angle = dt2 * ji[0] / I[0][0];
658 >  this->rotate( 1, 2, angle, ji, A );
659 >  
660 >  // rotate about the y-axis
661 >  angle = dt2 * ji[1] / I[1][1];
662 >  this->rotate( 2, 0, angle, ji, A );
663 >  
664 >  // rotate about the z-axis
665 >  angle = dt * ji[2] / I[2][2];
666 >  this->rotate( 0, 1, angle, ji, A);
667 >  
668 >  // rotate about the y-axis
669 >  angle = dt2 * ji[1] / I[1][1];
670 >  this->rotate( 2, 0, angle, ji, A );
671 >  
672 >  // rotate about the x-axis
673 >  angle = dt2 * ji[0] / I[0][0];
674 >  this->rotate( 1, 2, angle, ji, A );
675 >  
676 >  dAtom->setA( A  );    
677 > }
678 >
679 > template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
680 >                                                double angle, double ji[3],
681 >                                                double A[3][3]){
682 >  int i, j, k;
683    double sinAngle;
684    double cosAngle;
685    double angleSqr;
# Line 709 | Line 691 | void Integrator::rotate( int axes1, int axes2, double
691  
692    // initialize the tempA
693  
694 <  for(i=0; i<3; i++){
695 <    for(j=0; j<3; j++){
694 >  for (i = 0; i < 3; i++){
695 >    for (j = 0; j < 3; j++){
696        tempA[j][i] = A[i][j];
697      }
698    }
699  
700    // initialize the tempJ
701  
702 <  for( i=0; i<3; i++) tempJ[i] = ji[i];
703 <  
702 >  for (i = 0; i < 3; i++)
703 >    tempJ[i] = ji[i];
704 >
705    // initalize rot as a unit matrix
706  
707    rot[0][0] = 1.0;
# Line 728 | Line 711 | void Integrator::rotate( int axes1, int axes2, double
711    rot[1][0] = 0.0;
712    rot[1][1] = 1.0;
713    rot[1][2] = 0.0;
714 <  
714 >
715    rot[2][0] = 0.0;
716    rot[2][1] = 0.0;
717    rot[2][2] = 1.0;
718 <  
718 >
719    // use a small angle aproximation for sin and cosine
720  
721 <  angleSqr  = angle * angle;
721 >  angleSqr = angle * angle;
722    angleSqrOver4 = angleSqr / 4.0;
723    top = 1.0 - angleSqrOver4;
724    bottom = 1.0 + angleSqrOver4;
# Line 748 | Line 731 | void Integrator::rotate( int axes1, int axes2, double
731  
732    rot[axes1][axes2] = sinAngle;
733    rot[axes2][axes1] = -sinAngle;
734 <  
734 >
735    // rotate the momentum acoording to: ji[] = rot[][] * ji[]
736 <  
737 <  for(i=0; i<3; i++){
736 >
737 >  for (i = 0; i < 3; i++){
738      ji[i] = 0.0;
739 <    for(k=0; k<3; k++){
739 >    for (k = 0; k < 3; k++){
740        ji[i] += rot[i][k] * tempJ[k];
741      }
742    }
# Line 766 | Line 749 | void Integrator::rotate( int axes1, int axes2, double
749    // calculation as:
750    //                transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
751  
752 <  for(i=0; i<3; i++){
753 <    for(j=0; j<3; j++){
752 >  for (i = 0; i < 3; i++){
753 >    for (j = 0; j < 3; j++){
754        A[j][i] = 0.0;
755 <      for(k=0; k<3; k++){
756 <        A[j][i] += tempA[i][k] * rot[j][k];
755 >      for (k = 0; k < 3; k++){
756 >        A[j][i] += tempA[i][k] * rot[j][k];
757        }
758      }
759    }
760   }
761 +
762 + template<typename T> void Integrator<T>::calcForce(int calcPot, int calcStress){
763 +  myFF->doForces(calcPot, calcStress);
764 + }
765 +
766 + template<typename T> void Integrator<T>::thermalize(){
767 +  tStats->velocitize();
768 + }
769 +
770 + template<typename T> double Integrator<T>::getConservedQuantity(void){
771 +  return tStats->getTotalE();
772 + }

Diff Legend

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