255 |
|
template<typename T> double NVT<T>::getConservedQuantity(void){ |
256 |
|
|
257 |
|
double conservedQuantity; |
258 |
< |
double E_NVT; |
259 |
< |
|
260 |
< |
//HNVE |
261 |
< |
conservedQuantity = tStats->getTotalE(); |
262 |
< |
//HNVE |
263 |
< |
|
264 |
< |
E_NVT = (info->getNDF() * kB * targetTemp * |
265 |
< |
(integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0 )) / eConvert; |
258 |
> |
double fkBT; |
259 |
> |
double Energy; |
260 |
> |
double thermostat_kinetic; |
261 |
> |
double thermostat_potential; |
262 |
|
|
263 |
< |
conservedQuantity += E_NVT; |
263 |
> |
fkBT = (double)(info->getNDF() ) * kB * targetTemp; |
264 |
> |
|
265 |
> |
Energy = tStats->getTotalE(); |
266 |
|
|
267 |
< |
//cerr << info->getTime() << "\t" << chi << "\t" << integralOfChidt << "\t" << E_NVT << endl; |
267 |
> |
thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / |
268 |
> |
(2.0 * eConvert); |
269 |
|
|
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 |
|
} |