5 |
|
#include "Thermo.hpp" |
6 |
|
#include "ReadWrite.hpp" |
7 |
|
#include "ForceFields.hpp" |
8 |
+ |
#include "ExtendedSystem.hpp" |
9 |
|
#include "simError.h" |
10 |
|
|
11 |
|
extern "C"{ |
31 |
|
|
32 |
|
|
33 |
|
|
34 |
< |
Symplectic::Symplectic( SimInfo* the_entry_plug, ForceFields* the_ff ){ |
34 |
> |
Symplectic::Symplectic( SimInfo* the_entry_plug, ForceFields* the_ff, |
35 |
> |
ExtendedSystem* the_es ){ |
36 |
|
entry_plug = the_entry_plug; |
37 |
|
myFF = the_ff; |
38 |
+ |
myES = the_es; |
39 |
|
isFirst = 1; |
40 |
|
|
41 |
|
molecules = entry_plug->molecules; |
51 |
|
mass = new double[entry_plug->n_atoms]; |
52 |
|
for(int i = 0; i < entry_plug->n_atoms; i++){ |
53 |
|
mass[i] = entry_plug->atoms[i]->getMass(); |
54 |
< |
} |
52 |
< |
|
53 |
< |
|
54 |
> |
} |
55 |
|
|
56 |
|
// check for constraints |
57 |
|
|
184 |
|
int status_n = (int)( statusTime / dt ); |
185 |
|
int vel_n = (int)( thermalTime / dt ); |
186 |
|
|
187 |
< |
int calcPot; |
187 |
> |
int calcPot, calcStress; |
188 |
|
|
189 |
|
Thermo *tStats = new Thermo( entry_plug ); |
190 |
|
|
197 |
|
|
198 |
|
// initialize the forces the before the first step |
199 |
|
|
200 |
< |
myFF->doForces(1,0); |
200 |
> |
myFF->doForces(1,1); |
201 |
|
|
202 |
|
if( entry_plug->setTemp ){ |
203 |
|
|
325 |
|
|
326 |
|
// calculate the forces |
327 |
|
|
328 |
< |
myFF->doForces(calcPot, 0); |
328 |
> |
myFF->doForces(calcPot, calcStress); |
329 |
|
|
330 |
|
// move b |
331 |
|
|
396 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
397 |
|
} |
398 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
399 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
400 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
399 |
> |
if( !((time+1) % status_n) ) { |
400 |
> |
calcPot = 1; |
401 |
> |
calcStress = 1; |
402 |
> |
} |
403 |
> |
if( !(time % status_n) ){ |
404 |
> |
e_out->writeStat( time * dt ); |
405 |
> |
calcPot = 0; |
406 |
> |
calcStress = 0; |
407 |
> |
} |
408 |
|
} |
409 |
|
} |
410 |
|
else{ |
502 |
|
|
503 |
|
// calculate the forces |
504 |
|
|
505 |
< |
myFF->doForces(calcPot,0); |
505 |
> |
myFF->doForces(calcPot,calcStress); |
506 |
|
|
507 |
|
for( i=0; i< nAtoms; i++ ){ |
508 |
|
|
561 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
562 |
|
} |
563 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
564 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
565 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
564 |
> |
if( !((time+1) % status_n) ) { |
565 |
> |
calcPot = 1; |
566 |
> |
calcStress = 1; |
567 |
> |
} |
568 |
> |
if( !(time % status_n) ){ |
569 |
> |
e_out->writeStat( time * dt ); |
570 |
> |
calcPot = 0; |
571 |
> |
calcStress = 0; |
572 |
> |
} |
573 |
|
} |
574 |
|
} |
575 |
|
|