| 1 |
+ |
#include <stdlib.h> |
| 2 |
|
#include <math.h> |
| 3 |
+ |
|
| 4 |
|
#include "Atom.hpp" |
| 5 |
|
#include "SRI.hpp" |
| 6 |
|
#include "AbstractClasses.hpp" |
| 26 |
|
// |
| 27 |
|
// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499. |
| 28 |
|
|
| 29 |
< |
template<typename T> NPT<T>::NPT ( SimInfo *theInfo, ForceFields* the_ff): |
| 30 |
< |
T( theInfo, the_ff ) |
| 29 |
> |
NPT::NPT ( SimInfo *theInfo, ForceFields* the_ff): |
| 30 |
> |
Integrator( theInfo, the_ff ) |
| 31 |
|
{ |
| 32 |
|
GenericData* data; |
| 33 |
< |
DoubleData * chiValue; |
| 32 |
< |
DoubleData * integralOfChidtValue; |
| 33 |
< |
|
| 34 |
< |
chiValue = NULL; |
| 35 |
< |
integralOfChidtValue = NULL; |
| 36 |
< |
|
| 33 |
> |
|
| 34 |
|
chi = 0.0; |
| 35 |
|
integralOfChidt = 0.0; |
| 36 |
|
have_tau_thermostat = 0; |
| 43 |
|
|
| 44 |
|
// retrieve chi and integralOfChidt from simInfo |
| 45 |
|
data = info->getProperty(CHIVALUE_ID); |
| 46 |
< |
if(data){ |
| 47 |
< |
chiValue = dynamic_cast<DoubleData*>(data); |
| 46 |
> |
if(data != NULL ){ |
| 47 |
> |
chi = data->getDval(); |
| 48 |
|
} |
| 49 |
|
|
| 50 |
|
data = info->getProperty(INTEGRALOFCHIDT_ID); |
| 51 |
< |
if(data){ |
| 52 |
< |
integralOfChidtValue = dynamic_cast<DoubleData*>(data); |
| 51 |
> |
if(data != NULL ){ |
| 52 |
> |
integralOfChidt = data->getDval(); |
| 53 |
|
} |
| 54 |
|
|
| 58 |
– |
// chi and integralOfChidt should appear by pair |
| 59 |
– |
if(chiValue && integralOfChidtValue){ |
| 60 |
– |
chi = chiValue->getData(); |
| 61 |
– |
integralOfChidt = integralOfChidtValue->getData(); |
| 62 |
– |
} |
| 63 |
– |
|
| 55 |
|
oldPos = new double[3*nAtoms]; |
| 56 |
|
oldVel = new double[3*nAtoms]; |
| 57 |
|
oldJi = new double[3*nAtoms]; |
| 63 |
|
|
| 64 |
|
} |
| 65 |
|
|
| 66 |
< |
template<typename T> NPT<T>::~NPT() { |
| 66 |
> |
NPT::~NPT() { |
| 67 |
|
delete[] oldPos; |
| 68 |
|
delete[] oldVel; |
| 69 |
|
delete[] oldJi; |
| 70 |
|
} |
| 71 |
|
|
| 72 |
< |
template<typename T> void NPT<T>::moveA() { |
| 72 |
> |
void NPT::moveA() { |
| 73 |
|
|
| 74 |
|
//new version of NPT |
| 75 |
|
int i, j, k; |
| 81 |
|
double COM[3]; |
| 82 |
|
|
| 83 |
|
instaTemp = tStats->getTemperature(); |
| 84 |
+ |
|
| 85 |
|
tStats->getPressureTensor( press ); |
| 86 |
|
instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0; |
| 87 |
+ |
|
| 88 |
|
instaVol = tStats->getVolume(); |
| 89 |
|
|
| 90 |
|
tStats->getCOM(COM); |
| 173 |
|
this->scaleSimBox(); |
| 174 |
|
} |
| 175 |
|
|
| 176 |
< |
template<typename T> void NPT<T>::moveB( void ){ |
| 176 |
> |
void NPT::moveB( void ){ |
| 177 |
|
|
| 178 |
|
//new version of NPT |
| 179 |
|
int i, j, k; |
| 180 |
|
DirectionalAtom* dAtom; |
| 181 |
|
double Tb[3], ji[3], sc[3]; |
| 182 |
< |
double vel[3], frc[3]; |
| 182 |
> |
double vel[3], frc[3], qVel[3]; |
| 183 |
|
double mass; |
| 184 |
|
|
| 185 |
|
// Set things up for the iteration: |
| 209 |
|
|
| 210 |
|
for (k=0; k < 4; k++) { |
| 211 |
|
|
| 212 |
+ |
atoms[0]->getVel(vel); |
| 213 |
+ |
|
| 214 |
|
instaTemp = tStats->getTemperature(); |
| 215 |
+ |
|
| 216 |
|
instaPress = tStats->getPressure(); |
| 217 |
|
|
| 218 |
|
// evolve chi another half step using the temperature at t + dt/2 |
| 229 |
|
|
| 230 |
|
getVelScaleB( sc, i ); |
| 231 |
|
|
| 232 |
+ |
for(j=0;j<3;j++) |
| 233 |
+ |
qVel[j] = vel[j]; |
| 234 |
+ |
|
| 235 |
|
// velocity half step |
| 236 |
|
for (j=0; j < 3; j++) |
| 237 |
|
vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]); |
| 267 |
|
|
| 268 |
|
} |
| 269 |
|
|
| 270 |
< |
template<typename T> void NPT<T>::resetIntegrator() { |
| 271 |
< |
chi = 0.0; |
| 272 |
< |
T::resetIntegrator(); |
| 270 |
> |
void NPT::resetIntegrator() { |
| 271 |
> |
chi = 0.0; |
| 272 |
> |
Integrator::resetIntegrator(); |
| 273 |
|
} |
| 274 |
|
|
| 275 |
< |
template<typename T> void NPT<T>::evolveChiA() { |
| 275 |
> |
void NPT::evolveChiA() { |
| 276 |
> |
|
| 277 |
|
chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; |
| 278 |
|
oldChi = chi; |
| 279 |
|
} |
| 280 |
|
|
| 281 |
< |
template<typename T> void NPT<T>::evolveChiB() { |
| 281 |
> |
void NPT::evolveChiB() { |
| 282 |
|
|
| 283 |
|
prevChi = chi; |
| 284 |
|
chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2; |
| 285 |
|
} |
| 286 |
|
|
| 287 |
< |
template<typename T> bool NPT<T>::chiConverged() { |
| 287 |
> |
bool NPT::chiConverged() { |
| 288 |
|
|
| 289 |
|
return ( fabs( prevChi - chi ) <= chiTolerance ); |
| 290 |
|
} |
| 291 |
|
|
| 292 |
< |
template<typename T> int NPT<T>::readyCheck() { |
| 292 |
> |
int NPT::readyCheck() { |
| 293 |
|
|
| 294 |
|
//check parent's readyCheck() first |
| 295 |
< |
if (T::readyCheck() == -1) |
| 295 |
> |
if (Integrator::readyCheck() == -1) |
| 296 |
|
return -1; |
| 297 |
|
|
| 298 |
|
// First check to see if we have a target temperature. |
