ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Symplectic.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Symplectic.cpp (file contents):
Revision 472 by mmeineke, Mon Apr 7 21:16:35 2003 UTC vs.
Revision 497 by chuckv, Mon Apr 14 21:16:37 2003 UTC

# Line 38 | Line 38 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
38    myES = the_es;
39    isFirst = 1;
40  
41 –  std::cerr<< "calling symplectic constructor\n";
42 –
41    molecules = entry_plug->molecules;
42    nMols = entry_plug->n_mol;
43  
# Line 164 | Line 162 | void Symplectic::integrate( void ){
162    double dt2;                       // half the dt
163  
164    double vx, vy, vz;    // the velocities
165 < //  double vx2, vy2, vz2; // the square of the velocities
165 >  double vx2, vy2, vz2; // the square of the velocities
166    double rx, ry, rz;    // the postitions
167    
168    double ji[3];   // the body frame angular momentum
# Line 172 | Line 170 | void Symplectic::integrate( void ){
170    double Tb[3];   // torque in the body frame
171    double angle;   // the angle through which to rotate the rotation matrix
172    double A[3][3]; // the rotation matrix
173 +  double press[9];
174  
175    int time;
176  
# Line 192 | Line 191 | void Symplectic::integrate( void ){
191    StatWriter*  e_out;
192    DumpWriter*  dump_out;
193  
195 –  std::cerr << "about to call new thermo\n";
196 –
194    tStats   = new Thermo( entry_plug );
195    e_out    = new StatWriter( entry_plug );
199 –
200 –  std::cerr << "calling dumpWriter \n";
196    dump_out = new DumpWriter( entry_plug );
202 –  std::cerr << "called dumpWriter \n";
197  
198    Atom** atoms = entry_plug->atoms;
199    DirectionalAtom* dAtom;
# Line 219 | Line 213 | void Symplectic::integrate( void ){
213    
214    calcPot = 0;
215  
216 +  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
217 +    calcStress = 1;
218 +  } else {
219 +    calcStress = 0;
220 +  }
221 +
222    if( n_constrained ){
223  
224      double *Rx = new double[nAtoms];
# Line 235 | Line 235 | void Symplectic::integrate( void ){
235      
236  
237      for( tl=0; tl < n_loops; tl++ ){
238 +
239        
240        for( j=0; j<nAtoms; j++ ){
241  
# Line 272 | Line 273 | void Symplectic::integrate( void ){
273        }
274  
275  
276 < //       for( i=0; i<nAtoms; i++ ){
277 < // //   if( atoms[i]->isDirectional() ){
276 >      for( i=0; i<nAtoms; i++ ){
277 >        if( atoms[i]->isDirectional() ){
278                    
279 < // //     dAtom = (DirectionalAtom *)atoms[i];
279 >          dAtom = (DirectionalAtom *)atoms[i];
280            
281 < // //     // get and convert the torque to body frame
281 >          // get and convert the torque to body frame
282            
283 < // //     Tb[0] = dAtom->getTx();
284 < // //     Tb[1] = dAtom->getTy();
285 < // //     Tb[2] = dAtom->getTz();
283 >          Tb[0] = dAtom->getTx();
284 >          Tb[1] = dAtom->getTy();
285 >          Tb[2] = dAtom->getTz();
286            
287 < // //     dAtom->lab2Body( Tb );
287 >          dAtom->lab2Body( Tb );
288 >          
289 >          // get the angular momentum, and propagate a half step
290            
291 < // //     // get the angular momentum, and propagate a half step
291 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
292 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
293 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
294            
295 < // //     ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
291 < // //     ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
292 < // //     ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
295 >          // get the atom's rotation matrix
296            
297 < // //     // get the atom's rotation matrix
297 >          A[0][0] = dAtom->getAxx();
298 >          A[0][1] = dAtom->getAxy();
299 >          A[0][2] = dAtom->getAxz();
300            
301 < // //     A[0][0] = dAtom->getAxx();
302 < // //     A[0][1] = dAtom->getAxy();
303 < // //     A[0][2] = dAtom->getAxz();
299 <          
300 < // //     A[1][0] = dAtom->getAyx();
301 < // //     A[1][1] = dAtom->getAyy();
302 < // //     A[1][2] = dAtom->getAyz();
303 <          
304 < // //     A[2][0] = dAtom->getAzx();
305 < // //     A[2][1] = dAtom->getAzy();
306 < // //     A[2][2] = dAtom->getAzz();
301 >          A[1][0] = dAtom->getAyx();
302 >          A[1][1] = dAtom->getAyy();
303 >          A[1][2] = dAtom->getAyz();
304            
305 +          A[2][0] = dAtom->getAzx();
306 +          A[2][1] = dAtom->getAzy();
307 +          A[2][2] = dAtom->getAzz();
308            
309 – // //     // use the angular velocities to propagate the rotation matrix a
310 – // //     // full time step
309            
310 +          // use the angular velocities to propagate the rotation matrix a
311 +          // full time step
312            
313 – // //     angle = dt2 * ji[0] / dAtom->getIxx();
314 – // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
313            
314 < // //     angle = dt2 * ji[1] / dAtom->getIyy();
315 < // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
314 >          angle = dt2 * ji[0] / dAtom->getIxx();
315 >          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
316            
317 < // //     angle = dt * ji[2] / dAtom->getIzz();
318 < // //     this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
317 >          angle = dt2 * ji[1] / dAtom->getIyy();
318 >          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
319            
320 < // //     angle = dt2 * ji[1] / dAtom->getIyy();
321 < // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
320 >          angle = dt * ji[2] / dAtom->getIzz();
321 >          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
322            
323 < // //     angle = dt2 * ji[0] / dAtom->getIxx();
324 < // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
323 >          angle = dt2 * ji[1] / dAtom->getIyy();
324 >          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
325            
326 +          angle = dt2 * ji[0] / dAtom->getIxx();
327 +          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
328            
329 < // //     dAtom->setA( A );
330 < // //     dAtom->setJx( ji[0] );
331 < // //     dAtom->setJy( ji[1] );
332 < // //     dAtom->setJz( ji[2] );
333 < // //   }
334 < //       }
329 >          
330 >          dAtom->setA( A );
331 >          dAtom->setJx( ji[0] );
332 >          dAtom->setJy( ji[1] );
333 >          dAtom->setJz( ji[2] );
334 >        }
335 >      }
336        
337 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
338 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
339 +
340        // calculate the forces
341        
342        myFF->doForces(calcPot, calcStress);
# Line 373 | Line 377 | void Symplectic::integrate( void ){
377          atoms[j]->set_vz(Vz[j]);
378        }
379        
380 < //       for( i=0; i< nAtoms; i++ ){
380 >      for( i=0; i< nAtoms; i++ ){
381  
382 < //      if( atoms[i]->isDirectional() ){
382 >        if( atoms[i]->isDirectional() ){
383  
384 < //        dAtom = (DirectionalAtom *)atoms[i];
384 >          dAtom = (DirectionalAtom *)atoms[i];
385            
386 < //        // get and convert the torque to body frame
386 >          // get and convert the torque to body frame
387            
388 < //        Tb[0] = dAtom->getTx();
389 < //        Tb[1] = dAtom->getTy();
390 < //        Tb[2] = dAtom->getTz();
388 >          Tb[0] = dAtom->getTx();
389 >          Tb[1] = dAtom->getTy();
390 >          Tb[2] = dAtom->getTz();
391            
392 < //        dAtom->lab2Body( Tb );
392 >          dAtom->lab2Body( Tb );
393            
394 < //        // get the angular momentum, and complete the angular momentum
395 < //        // half step
394 >          // get the angular momentum, and complete the angular momentum
395 >          // half step
396            
397 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
398 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
399 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
397 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
398 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
399 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
400            
401 < //        dAtom->setJx( ji[0] );
402 < //        dAtom->setJy( ji[1] );
403 < //        dAtom->setJz( ji[2] );
404 < //      }
405 < //       }
401 >          dAtom->setJx( ji[0] );
402 >          dAtom->setJy( ji[1] );
403 >          dAtom->setJz( ji[2] );
404 >        }
405 >      }
406      
407 +
408 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
409 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
410 +
411 +      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
412 +        tStats->getPressureTensor(press);
413 +        myES->NoseHooverAndersonNPT( dt,
414 +                                     tStats->getKinetic(),
415 +                                     press);
416 +      }
417 +
418        time = tl + 1;
419        
420        if( entry_plug->setTemp ){
# Line 408 | Line 423 | void Symplectic::integrate( void ){
423        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
424        if( !((time+1) % status_n) ) {
425          calcPot = 1;
426 <        calcStress = 1;
426 >        calcStress = 1;
427        }
428        if( !(time % status_n) ){
429          e_out->writeStat( time * dt );
430          calcPot = 0;
431 <        calcStress = 0;
431 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
432 >        else calcStress = 0;
433        }
434      }
435    }
# Line 424 | Line 440 | void Symplectic::integrate( void ){
440        kE = 0.0;
441        rot_kE= 0.0;
442        trans_kE = 0.0;
443 <      
443 >
444        for( i=0; i<nAtoms; i++ ){
445          
446          // velocity half step
# Line 450 | Line 466 | void Symplectic::integrate( void ){
466          atoms[i]->set_vy( vy );
467          atoms[i]->set_vz( vz );
468          
469 < //      if( atoms[i]->isDirectional() ){
469 >        if( atoms[i]->isDirectional() ){
470  
471 < //        dAtom = (DirectionalAtom *)atoms[i];
471 >          dAtom = (DirectionalAtom *)atoms[i];
472            
473 < //        // get and convert the torque to body frame
473 >          // get and convert the torque to body frame
474            
475 < //        Tb[0] = dAtom->getTx();
476 < //        Tb[1] = dAtom->getTy();
477 < //        Tb[2] = dAtom->getTz();
475 >          Tb[0] = dAtom->getTx();
476 >          Tb[1] = dAtom->getTy();
477 >          Tb[2] = dAtom->getTz();
478            
479 < //        dAtom->lab2Body( Tb );
479 >          dAtom->lab2Body( Tb );
480            
481 < //        // get the angular momentum, and propagate a half step
481 >          // get the angular momentum, and propagate a half step
482            
483 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
484 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
485 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
483 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
484 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
485 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
486            
487 < //        // get the atom's rotation matrix
487 >          // get the atom's rotation matrix
488            
489 < //        A[0][0] = dAtom->getAxx();
490 < //        A[0][1] = dAtom->getAxy();
491 < //        A[0][2] = dAtom->getAxz();
489 >          A[0][0] = dAtom->getAxx();
490 >          A[0][1] = dAtom->getAxy();
491 >          A[0][2] = dAtom->getAxz();
492            
493 < //        A[1][0] = dAtom->getAyx();
494 < //        A[1][1] = dAtom->getAyy();
495 < //        A[1][2] = dAtom->getAyz();
493 >          A[1][0] = dAtom->getAyx();
494 >          A[1][1] = dAtom->getAyy();
495 >          A[1][2] = dAtom->getAyz();
496            
497 < //        A[2][0] = dAtom->getAzx();
498 < //        A[2][1] = dAtom->getAzy();
499 < //        A[2][2] = dAtom->getAzz();
497 >          A[2][0] = dAtom->getAzx();
498 >          A[2][1] = dAtom->getAzy();
499 >          A[2][2] = dAtom->getAzz();
500            
501            
502 < //        // use the angular velocities to propagate the rotation matrix a
503 < //        // full time step
502 >          // use the angular velocities to propagate the rotation matrix a
503 >          // full time step
504            
505            
506 < //        angle = dt2 * ji[0] / dAtom->getIxx();
507 < //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
506 >          angle = dt2 * ji[0] / dAtom->getIxx();
507 >          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
508            
509 < //        angle = dt2 * ji[1] / dAtom->getIyy();
510 < //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
509 >          angle = dt2 * ji[1] / dAtom->getIyy();
510 >          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
511            
512 < //        angle = dt * ji[2] / dAtom->getIzz();
513 < //        this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
512 >          angle = dt * ji[2] / dAtom->getIzz();
513 >          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
514            
515 < //        angle = dt2 * ji[1] / dAtom->getIyy();
516 < //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
515 >          angle = dt2 * ji[1] / dAtom->getIyy();
516 >          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
517            
518 < //        angle = dt2 * ji[0] / dAtom->getIxx();
519 < //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
518 >          angle = dt2 * ji[0] / dAtom->getIxx();
519 >          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
520            
521            
522 < //        dAtom->setA( A );
523 < //        dAtom->setJx( ji[0] );
524 < //        dAtom->setJy( ji[1] );
525 < //        dAtom->setJz( ji[2] );
526 < //      }
522 >          dAtom->setA( A );
523 >          dAtom->setJx( ji[0] );
524 >          dAtom->setJy( ji[1] );
525 >          dAtom->setJz( ji[2] );
526 >        }
527        }
528 +
529 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
530 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
531        
532 +      
533        // calculate the forces
534        
535        myFF->doForces(calcPot,calcStress);
# Line 529 | Line 549 | void Symplectic::integrate( void ){
549          atoms[i]->set_vy( vy );
550          atoms[i]->set_vz( vz );
551          
552 < //      vx2 = vx * vx;
553 < //      vy2 = vy * vy;
554 < //      vz2 = vz * vz;
552 >        vx2 = vx * vx;
553 >        vy2 = vy * vy;
554 >        vz2 = vz * vz;
555          
556 < //      if( atoms[i]->isDirectional() ){
556 >        if( atoms[i]->isDirectional() ){
557  
558 < //        dAtom = (DirectionalAtom *)atoms[i];
558 >          dAtom = (DirectionalAtom *)atoms[i];
559            
560 < //        // get and convert the torque to body frame
560 >          // get and convert the torque to body frame
561            
562 < //        Tb[0] = dAtom->getTx();
563 < //        Tb[1] = dAtom->getTy();
564 < //        Tb[2] = dAtom->getTz();
562 >          Tb[0] = dAtom->getTx();
563 >          Tb[1] = dAtom->getTy();
564 >          Tb[2] = dAtom->getTz();
565            
566 < //        dAtom->lab2Body( Tb );
566 >          dAtom->lab2Body( Tb );
567            
568 < //        // get the angular momentum, and complete the angular momentum
569 < //        // half step
568 >          // get the angular momentum, and complete the angular momentum
569 >          // half step
570            
571 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
572 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
573 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
571 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
572 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
573 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
574            
575 < //        jx2 = ji[0] * ji[0];
576 < //        jy2 = ji[1] * ji[1];
577 < //        jz2 = ji[2] * ji[2];
575 >          jx2 = ji[0] * ji[0];
576 >          jy2 = ji[1] * ji[1];
577 >          jz2 = ji[2] * ji[2];
578            
579 < //        rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
580 < //          + (jz2 / dAtom->getIzz());
579 >          rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
580 >            + (jz2 / dAtom->getIzz());
581            
582 < //        dAtom->setJx( ji[0] );
583 < //        dAtom->setJy( ji[1] );
584 < //        dAtom->setJz( ji[2] );
585 < //      }
582 >          dAtom->setJx( ji[0] );
583 >          dAtom->setJy( ji[1] );
584 >          dAtom->setJz( ji[2] );
585 >        }
586 >
587        }
588 <      
588 >    
589 >      if (!strcasecmp( entry_plug->ensemble, "NVT"))
590 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
591 >
592 >      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
593 >        tStats->getPressureTensor(press);
594 >        myES->NoseHooverAndersonNPT( dt,
595 >                                     tStats->getKinetic(),
596 >                                     press);
597 >      }
598 >  
599        time = tl + 1;
600        
601        if( entry_plug->setTemp ){
# Line 573 | Line 604 | void Symplectic::integrate( void ){
604        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
605        if( !((time+1) % status_n) ) {
606          calcPot = 1;
607 <        calcStress = 1;
607 >        calcStress = 1;
608        }
609        if( !(time % status_n) ){
610          e_out->writeStat( time * dt );
611          calcPot = 0;
612 <        calcStress = 0;
612 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
613 >        else calcStress = 0;
614        }      
615      }
616    }

Diff Legend

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