simtk.org

Hi,

I am running the example file: ExampleAdenylateMobilitiesVTK.cpp.
I tried to change the thermostat to the Nose-Hoover by adding the following line in the main:
Force::Thermostat (forces, matter,SimTK_BOLTZMANN_CONSTANT_MD,293.15,1.0);
but it failed.
How is the correct way of implementing it?
And, what is the correct form of using "clacCurrentTemperature(state)"?
I appreciate your help.
Regards

Hi, Mohammad.

Force::Thermostat is one of the forces supported by a GeneralForceSubsystem; it is not part of the DuMMForceFieldSubsystem which is what that example calls "forces". So you need to add a GeneralForceSubsystem and then put the thermostat in there:
// At top of file.
GeneralForceSubsystem genForces(system);
//...
Force::Thermostat thermo(genForces, matter,SimTK_BOLTZMANN_CONSTANT_MD,293.15,1.0);

(Be sure to remove the VelocityRescalingThermostat event handler.)

Then above I named the Thermostat force element "thermo" so I can ask it for the current temperature via thermo.getCurrentTemperature(state), e.g.:
std::cout << "temp=" << thermo.getCurrentTemperature(state) << std::endl;

(You would have to put the above line in an event reporter that you have the TimeStepper call periodically -- please repost if you aren't sure how to do that.)

Be sure to check the doxygen documentation for Simbody's Force::Temperature element -- it doesn't have a method calcCurrentTemperature().

Best regards,
Sherm

Also, Molmodel has a class NoseHooverThermostat which is a slightly nicer interface to Simbody's Force::Thermostat (it derives from Force::Thermostat). It knows that Molmodel uses MD units so you don't have to give it Boltzmann's constant:

NoseHooverThermostat thermo(genForces,matter,293.14,1.0).

Sherm

Dear Sherm,

I implemented
std::cout << "temp=" << thermo.getCurrentTemperature(state) << std::endl;
in the eventReporter. But the error
thermo undeclared indentifier
appears. I just guess that I need to use some methods to get the thermostat then add line.
I appreciate your time and help.
Regards

Yes, typically you just pass a reference to the objects you need in a reporter to the constructor for that reporter, then keep the reference as a data member of the reporter class so you can access it when you need it.

Sherm

Hi Sherm,

Great, Got it to work.

Thanks for your help.