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

Diff Legend

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