108 |
|
dAtom->setA( A ); |
109 |
|
} |
110 |
|
} |
111 |
+ |
|
112 |
+ |
if (nConstrained){ |
113 |
+ |
constrainA(); |
114 |
+ |
} |
115 |
|
|
116 |
|
// Finally, evolve chi a half step (just like a velocity) using |
117 |
|
// temperature at time t, not time t+dt/2 |
194 |
|
} |
195 |
|
} |
196 |
|
|
197 |
+ |
if (nConstrained){ |
198 |
+ |
constrainB(); |
199 |
+ |
} |
200 |
+ |
|
201 |
|
if (fabs(prevChi - chi) <= chiTolerance) break; |
202 |
|
} |
203 |
|
|
255 |
|
template<typename T> double NVT<T>::getConservedQuantity(void){ |
256 |
|
|
257 |
|
double conservedQuantity; |
258 |
< |
double E_NVT; |
258 |
> |
double fkBT; |
259 |
> |
double Energy; |
260 |
> |
double thermostat_kinetic; |
261 |
> |
double thermostat_potential; |
262 |
|
|
263 |
< |
//HNVE |
264 |
< |
conservedQuantity = tStats->getTotalE(); |
265 |
< |
//HNVE |
255 |
< |
|
256 |
< |
E_NVT = (info->getNDF() * kB * targetTemp * |
257 |
< |
(integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0 )) / eConvert; |
263 |
> |
fkBT = (double)(info->getNDF() ) * kB * targetTemp; |
264 |
> |
|
265 |
> |
Energy = tStats->getTotalE(); |
266 |
|
|
267 |
< |
conservedQuantity += E_NVT; |
267 |
> |
thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / |
268 |
> |
(2.0 * eConvert); |
269 |
|
|
270 |
< |
//cerr << info->getTime() << "\t" << chi << "\t" << integralOfChidt << "\t" << E_NVT << endl; |
270 |
> |
thermostat_potential = fkBT * integralOfChidt / eConvert; |
271 |
|
|
272 |
+ |
conservedQuantity = Energy + thermostat_kinetic + thermostat_potential; |
273 |
+ |
|
274 |
+ |
cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic << |
275 |
+ |
"\t" << thermostat_potential << "\t" << conservedQuantity << endl; |
276 |
+ |
|
277 |
|
return conservedQuantity; |
278 |
|
} |