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

Diff Legend

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