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 466 by gezelter, Mon Apr 7 14:30:36 2003 UTC vs.
Revision 472 by mmeineke, Mon Apr 7 21:16:35 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 +
43    molecules = entry_plug->molecules;
44    nMols = entry_plug->n_mol;
45  
# Line 184 | Line 186 | void Symplectic::integrate( void ){
186    int status_n = (int)( statusTime / dt );
187    int vel_n    = (int)( thermalTime / dt );
188  
189 <  int calcPot;
189 >  int calcPot, calcStress;
190  
191 <   Thermo *tStats = new Thermo( entry_plug );
191 >  Thermo *tStats;
192 >  StatWriter*  e_out;
193 >  DumpWriter*  dump_out;
194  
195 <  StatWriter*  e_out    = new StatWriter( entry_plug );
192 <  DumpWriter*  dump_out = new DumpWriter( entry_plug );
195 >  std::cerr << "about to call new thermo\n";
196  
197 +  tStats   = new Thermo( entry_plug );
198 +  e_out    = new StatWriter( entry_plug );
199 +
200 +  std::cerr << "calling dumpWriter \n";
201 +  dump_out = new DumpWriter( entry_plug );
202 +  std::cerr << "called dumpWriter \n";
203 +
204    Atom** atoms = entry_plug->atoms;
205    DirectionalAtom* dAtom;
206    dt2 = 0.5 * dt;
207  
208    // initialize the forces the before the first step
209  
210 <  myFF->doForces(1,0);
210 >  myFF->doForces(1,1);
211    
212    if( entry_plug->setTemp ){
213      
# Line 262 | Line 272 | void Symplectic::integrate( void ){
272        }
273  
274  
275 <      for( i=0; i<nAtoms; i++ ){
276 <        if( atoms[i]->isDirectional() ){
275 > //       for( i=0; i<nAtoms; i++ ){
276 > // //   if( atoms[i]->isDirectional() ){
277                    
278 <          dAtom = (DirectionalAtom *)atoms[i];
278 > // //     dAtom = (DirectionalAtom *)atoms[i];
279            
280 <          // get and convert the torque to body frame
280 > // //     // get and convert the torque to body frame
281            
282 <          Tb[0] = dAtom->getTx();
283 <          Tb[1] = dAtom->getTy();
284 <          Tb[2] = dAtom->getTz();
282 > // //     Tb[0] = dAtom->getTx();
283 > // //     Tb[1] = dAtom->getTy();
284 > // //     Tb[2] = dAtom->getTz();
285            
286 <          dAtom->lab2Body( Tb );
286 > // //     dAtom->lab2Body( Tb );
287            
288 <          // get the angular momentum, and propagate a half step
288 > // //     // get the angular momentum, and propagate a half step
289            
290 <          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;
290 > // //     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;
293            
294 <          // get the atom's rotation matrix
294 > // //     // get the atom's rotation matrix
295            
296 <          A[0][0] = dAtom->getAxx();
297 <          A[0][1] = dAtom->getAxy();
298 <          A[0][2] = dAtom->getAxz();
289 <          
290 <          A[1][0] = dAtom->getAyx();
291 <          A[1][1] = dAtom->getAyy();
292 <          A[1][2] = dAtom->getAyz();
296 > // //     A[0][0] = dAtom->getAxx();
297 > // //     A[0][1] = dAtom->getAxy();
298 > // //     A[0][2] = dAtom->getAxz();
299            
300 <          A[2][0] = dAtom->getAzx();
301 <          A[2][1] = dAtom->getAzy();
302 <          A[2][2] = dAtom->getAzz();
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();
307            
299 –          // use the angular velocities to propagate the rotation matrix a
300 –          // full time step
308            
309 + // //     // use the angular velocities to propagate the rotation matrix a
310 + // //     // full time step
311            
303 –          angle = dt2 * ji[0] / dAtom->getIxx();
304 –          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
312            
313 <          angle = dt2 * ji[1] / dAtom->getIyy();
314 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
313 > // //     angle = dt2 * ji[0] / dAtom->getIxx();
314 > // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
315            
316 <          angle = dt * ji[2] / dAtom->getIzz();
317 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
316 > // //     angle = dt2 * ji[1] / dAtom->getIyy();
317 > // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
318            
319 <          angle = dt2 * ji[1] / dAtom->getIyy();
320 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
319 > // //     angle = dt * ji[2] / dAtom->getIzz();
320 > // //     this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
321            
322 <          angle = dt2 * ji[0] / dAtom->getIxx();
323 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
322 > // //     angle = dt2 * ji[1] / dAtom->getIyy();
323 > // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
324            
325 + // //     angle = dt2 * ji[0] / dAtom->getIxx();
326 + // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
327            
328 <          dAtom->setA( A );
329 <          dAtom->setJx( ji[0] );
330 <          dAtom->setJy( ji[1] );
331 <          dAtom->setJz( ji[2] );
332 <        }
333 <      }
328 >          
329 > // //     dAtom->setA( A );
330 > // //     dAtom->setJx( ji[0] );
331 > // //     dAtom->setJy( ji[1] );
332 > // //     dAtom->setJz( ji[2] );
333 > // //   }
334 > //       }
335        
336        // calculate the forces
337        
338 <      myFF->doForces(calcPot, 0);
338 >      myFF->doForces(calcPot, calcStress);
339        
340        // move b
341  
# Line 363 | Line 373 | void Symplectic::integrate( void ){
373          atoms[j]->set_vz(Vz[j]);
374        }
375        
376 <      for( i=0; i< nAtoms; i++ ){
376 > //       for( i=0; i< nAtoms; i++ ){
377  
378 <        if( atoms[i]->isDirectional() ){
378 > //      if( atoms[i]->isDirectional() ){
379  
380 <          dAtom = (DirectionalAtom *)atoms[i];
380 > //        dAtom = (DirectionalAtom *)atoms[i];
381            
382 <          // get and convert the torque to body frame
382 > //        // get and convert the torque to body frame
383            
384 <          Tb[0] = dAtom->getTx();
385 <          Tb[1] = dAtom->getTy();
386 <          Tb[2] = dAtom->getTz();
384 > //        Tb[0] = dAtom->getTx();
385 > //        Tb[1] = dAtom->getTy();
386 > //        Tb[2] = dAtom->getTz();
387            
388 <          dAtom->lab2Body( Tb );
388 > //        dAtom->lab2Body( Tb );
389            
390 <          // get the angular momentum, and complete the angular momentum
391 <          // half step
390 > //        // get the angular momentum, and complete the angular momentum
391 > //        // half step
392            
393 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
394 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
395 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
393 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
394 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
395 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
396            
397 <          dAtom->setJx( ji[0] );
398 <          dAtom->setJy( ji[1] );
399 <          dAtom->setJz( ji[2] );
400 <        }
401 <      }
397 > //        dAtom->setJx( ji[0] );
398 > //        dAtom->setJy( ji[1] );
399 > //        dAtom->setJz( ji[2] );
400 > //      }
401 > //       }
402      
403        time = tl + 1;
404        
# Line 396 | Line 406 | void Symplectic::integrate( void ){
406          if( !(time % vel_n) ) tStats->velocitize();
407        }
408        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
409 <      if( !((time+1) % status_n) ) calcPot = 1;
410 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
409 >      if( !((time+1) % status_n) ) {
410 >        calcPot = 1;
411 >        calcStress = 1;
412 >      }
413 >      if( !(time % status_n) ){
414 >        e_out->writeStat( time * dt );
415 >        calcPot = 0;
416 >        calcStress = 0;
417 >      }
418      }
419    }
420    else{
# Line 433 | Line 450 | void Symplectic::integrate( void ){
450          atoms[i]->set_vy( vy );
451          atoms[i]->set_vz( vz );
452          
453 <        if( atoms[i]->isDirectional() ){
453 > //      if( atoms[i]->isDirectional() ){
454  
455 <          dAtom = (DirectionalAtom *)atoms[i];
455 > //        dAtom = (DirectionalAtom *)atoms[i];
456            
457 <          // get and convert the torque to body frame
457 > //        // get and convert the torque to body frame
458            
459 <          Tb[0] = dAtom->getTx();
460 <          Tb[1] = dAtom->getTy();
461 <          Tb[2] = dAtom->getTz();
459 > //        Tb[0] = dAtom->getTx();
460 > //        Tb[1] = dAtom->getTy();
461 > //        Tb[2] = dAtom->getTz();
462            
463 <          dAtom->lab2Body( Tb );
463 > //        dAtom->lab2Body( Tb );
464            
465 <          // get the angular momentum, and propagate a half step
465 > //        // get the angular momentum, and propagate a half step
466            
467 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
468 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
469 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
467 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
468 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
469 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
470            
471 <          // get the atom's rotation matrix
471 > //        // get the atom's rotation matrix
472            
473 <          A[0][0] = dAtom->getAxx();
474 <          A[0][1] = dAtom->getAxy();
475 <          A[0][2] = dAtom->getAxz();
473 > //        A[0][0] = dAtom->getAxx();
474 > //        A[0][1] = dAtom->getAxy();
475 > //        A[0][2] = dAtom->getAxz();
476            
477 <          A[1][0] = dAtom->getAyx();
478 <          A[1][1] = dAtom->getAyy();
479 <          A[1][2] = dAtom->getAyz();
477 > //        A[1][0] = dAtom->getAyx();
478 > //        A[1][1] = dAtom->getAyy();
479 > //        A[1][2] = dAtom->getAyz();
480            
481 <          A[2][0] = dAtom->getAzx();
482 <          A[2][1] = dAtom->getAzy();
483 <          A[2][2] = dAtom->getAzz();
481 > //        A[2][0] = dAtom->getAzx();
482 > //        A[2][1] = dAtom->getAzy();
483 > //        A[2][2] = dAtom->getAzz();
484            
485            
486 <          // use the angular velocities to propagate the rotation matrix a
487 <          // full time step
486 > //        // use the angular velocities to propagate the rotation matrix a
487 > //        // full time step
488            
489            
490 <          angle = dt2 * ji[0] / dAtom->getIxx();
491 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
490 > //        angle = dt2 * ji[0] / dAtom->getIxx();
491 > //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
492            
493 <          angle = dt2 * ji[1] / dAtom->getIyy();
494 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
493 > //        angle = dt2 * ji[1] / dAtom->getIyy();
494 > //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
495            
496 <          angle = dt * ji[2] / dAtom->getIzz();
497 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
496 > //        angle = dt * ji[2] / dAtom->getIzz();
497 > //        this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
498            
499 <          angle = dt2 * ji[1] / dAtom->getIyy();
500 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
499 > //        angle = dt2 * ji[1] / dAtom->getIyy();
500 > //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
501            
502 <          angle = dt2 * ji[0] / dAtom->getIxx();
503 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
502 > //        angle = dt2 * ji[0] / dAtom->getIxx();
503 > //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
504            
505            
506 <          dAtom->setA( A );
507 <          dAtom->setJx( ji[0] );
508 <          dAtom->setJy( ji[1] );
509 <          dAtom->setJz( ji[2] );
510 <        }
506 > //        dAtom->setA( A );
507 > //        dAtom->setJx( ji[0] );
508 > //        dAtom->setJy( ji[1] );
509 > //        dAtom->setJz( ji[2] );
510 > //      }
511        }
512        
513        // calculate the forces
514        
515 <      myFF->doForces(calcPot,0);
515 >      myFF->doForces(calcPot,calcStress);
516        
517        for( i=0; i< nAtoms; i++ ){
518          
# Line 516 | Line 533 | void Symplectic::integrate( void ){
533   //      vy2 = vy * vy;
534   //      vz2 = vz * vz;
535          
536 <        if( atoms[i]->isDirectional() ){
536 > //      if( atoms[i]->isDirectional() ){
537  
538 <          dAtom = (DirectionalAtom *)atoms[i];
538 > //        dAtom = (DirectionalAtom *)atoms[i];
539            
540 <          // get and convert the torque to body frame
540 > //        // get and convert the torque to body frame
541            
542 <          Tb[0] = dAtom->getTx();
543 <          Tb[1] = dAtom->getTy();
544 <          Tb[2] = dAtom->getTz();
542 > //        Tb[0] = dAtom->getTx();
543 > //        Tb[1] = dAtom->getTy();
544 > //        Tb[2] = dAtom->getTz();
545            
546 <          dAtom->lab2Body( Tb );
546 > //        dAtom->lab2Body( Tb );
547            
548 <          // get the angular momentum, and complete the angular momentum
549 <          // half step
548 > //        // get the angular momentum, and complete the angular momentum
549 > //        // half step
550            
551 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
552 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
553 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
551 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
552 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
553 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
554            
555 <          jx2 = ji[0] * ji[0];
556 <          jy2 = ji[1] * ji[1];
557 <          jz2 = ji[2] * ji[2];
555 > //        jx2 = ji[0] * ji[0];
556 > //        jy2 = ji[1] * ji[1];
557 > //        jz2 = ji[2] * ji[2];
558            
559 <          rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
560 <            + (jz2 / dAtom->getIzz());
559 > //        rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
560 > //          + (jz2 / dAtom->getIzz());
561            
562 <          dAtom->setJx( ji[0] );
563 <          dAtom->setJy( ji[1] );
564 <          dAtom->setJz( ji[2] );
565 <        }
562 > //        dAtom->setJx( ji[0] );
563 > //        dAtom->setJy( ji[1] );
564 > //        dAtom->setJz( ji[2] );
565 > //      }
566        }
567        
568        time = tl + 1;
# Line 554 | Line 571 | void Symplectic::integrate( void ){
571          if( !(time % vel_n) ) tStats->velocitize();
572        }
573        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
574 <      if( !((time+1) % status_n) ) calcPot = 1;
575 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
574 >      if( !((time+1) % status_n) ) {
575 >        calcPot = 1;
576 >        calcStress = 1;
577 >      }
578 >      if( !(time % status_n) ){
579 >        e_out->writeStat( time * dt );
580 >        calcPot = 0;
581 >        calcStress = 0;
582 >      }      
583      }
584    }
585  

Diff Legend

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