| 4 |
|
#include "Integrator.hpp" |
| 5 |
|
#include "Thermo.hpp" |
| 6 |
|
#include "ReadWrite.hpp" |
| 7 |
+ |
#include "ForceField.hpp" |
| 8 |
+ |
#include "simError.h" |
| 9 |
|
|
| 8 |
– |
|
| 10 |
|
extern "C"{ |
| 11 |
|
|
| 12 |
|
void v_constrain_a_( double &dt, int &n_atoms, double* mass, |
| 35 |
|
isFirst = 1; |
| 36 |
|
|
| 37 |
|
srInteractions = entry_plug->sr_interactions; |
| 37 |
– |
longRange = entry_plug->longRange; |
| 38 |
|
nSRI = entry_plug->n_SRI; |
| 39 |
|
|
| 40 |
|
// give a little love back to the SimInfo object |
| 120 |
|
double dt2; // half the dt |
| 121 |
|
|
| 122 |
|
double vx, vy, vz; // the velocities |
| 123 |
< |
double vx2, vy2, vz2; // the square of the velocities |
| 123 |
> |
// double vx2, vy2, vz2; // the square of the velocities |
| 124 |
|
double rx, ry, rz; // the postitions |
| 125 |
|
|
| 126 |
|
double ji[3]; // the body frame angular momentum |
| 129 |
|
double angle; // the angle through which to rotate the rotation matrix |
| 130 |
|
double A[3][3]; // the rotation matrix |
| 131 |
|
|
| 132 |
+ |
int time; |
| 133 |
|
|
| 134 |
|
double dt = entry_plug->dt; |
| 135 |
|
double runTime = entry_plug->run_time; |
| 142 |
|
int status_n = (int)( statusTime / dt ); |
| 143 |
|
int vel_n = (int)( thermalTime / dt ); |
| 144 |
|
|
| 145 |
+ |
ForceFields *ff = entry_plug-> |
| 146 |
+ |
|
| 147 |
|
Thermo *tStats = new Thermo( entry_plug ); |
| 148 |
|
|
| 149 |
|
StatWriter* e_out = new StatWriter( entry_plug ); |
| 150 |
|
DumpWriter* dump_out = new DumpWriter( entry_plug ); |
| 151 |
|
|
| 149 |
– |
|
| 152 |
|
Atom** atoms = entry_plug->atoms; |
| 153 |
|
DirectionalAtom* dAtom; |
| 154 |
|
dt2 = 0.5 * dt; |
| 156 |
|
// initialize the forces the before the first step |
| 157 |
|
|
| 158 |
|
|
| 159 |
+ |
|
| 160 |
|
for(i = 0; i < nAtoms; i++){ |
| 161 |
|
atoms[i]->zeroForces(); |
| 162 |
|
} |
| 173 |
|
tStats->velocitize(); |
| 174 |
|
} |
| 175 |
|
|
| 176 |
+ |
dump_out->writeDump( 0.0 ); |
| 177 |
+ |
e_out->writeStat( 0.0 ); |
| 178 |
+ |
|
| 179 |
|
if( n_constrained ){ |
| 180 |
|
|
| 181 |
|
double *Rx = new double[nAtoms]; |
| 364 |
|
dAtom->setJz( ji[2] ); |
| 365 |
|
} |
| 366 |
|
} |
| 367 |
< |
|
| 367 |
> |
|
| 368 |
> |
time = tl + 1; |
| 369 |
> |
|
| 370 |
|
if( entry_plug->setTemp ){ |
| 371 |
< |
if( !(tl % vel_n) ) tStats->velocitize(); |
| 371 |
> |
if( !(time % vel_n) ) tStats->velocitize(); |
| 372 |
|
} |
| 373 |
< |
if( !(tl % sample_n) ) dump_out->writeDump( tl * dt ); |
| 374 |
< |
if( !(tl % status_n) ) e_out->writeStat( tl * dt ); |
| 373 |
> |
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 374 |
> |
if( !(time % status_n) ) e_out->writeStat( time * dt ); |
| 375 |
|
} |
| 376 |
|
} |
| 377 |
|
else{ |
| 494 |
|
atoms[i]->set_vy( vy ); |
| 495 |
|
atoms[i]->set_vz( vz ); |
| 496 |
|
|
| 497 |
< |
vx2 = vx * vx; |
| 498 |
< |
vy2 = vy * vy; |
| 499 |
< |
vz2 = vz * vz; |
| 497 |
> |
// vx2 = vx * vx; |
| 498 |
> |
// vy2 = vy * vy; |
| 499 |
> |
// vz2 = vz * vz; |
| 500 |
|
|
| 501 |
|
if( atoms[i]->isDirectional() ){ |
| 502 |
|
|
| 529 |
|
dAtom->setJz( ji[2] ); |
| 530 |
|
} |
| 531 |
|
} |
| 532 |
< |
|
| 532 |
> |
|
| 533 |
> |
time = tl + 1; |
| 534 |
> |
|
| 535 |
|
if( entry_plug->setTemp ){ |
| 536 |
< |
if( !(tl % vel_n) ) tStats->velocitize(); |
| 536 |
> |
if( !(time % vel_n) ) tStats->velocitize(); |
| 537 |
|
} |
| 538 |
< |
if( !(tl % sample_n) ) dump_out->writeDump( tl * dt ); |
| 539 |
< |
if( !(tl % status_n) ) e_out->writeStat( tl * dt ); |
| 538 |
> |
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 539 |
> |
if( !(time % status_n) ) e_out->writeStat( time * dt ); |
| 540 |
|
} |
| 541 |
|
} |
| 542 |
|
|
