Question about Custom Thermostat and Barostat
- Kailong Mao
- Posts: 16
- Joined: Fri Oct 02, 2015 8:56 am
Question about Custom Thermostat and Barostat
Hi Friends:
I've been trying to implement a custom BUSSI thermostat in OpenMM as a Custom Force. Basically, what I'd like to do is to scale the velocities of all atoms in the system based on the current temperature. Is there a way to access/calculate the temperature (or the kinetic energy) of the system directly from a custom force class? Also, I found a CUDA implementation of the Andersen thermostat in andersenThermostat.cu. I wonder if I can just modify this file to implement my own thermostat. How would you access/calculate the kinetic energy from within this file?
Thank you very much!
I've been trying to implement a custom BUSSI thermostat in OpenMM as a Custom Force. Basically, what I'd like to do is to scale the velocities of all atoms in the system based on the current temperature. Is there a way to access/calculate the temperature (or the kinetic energy) of the system directly from a custom force class? Also, I found a CUDA implementation of the Andersen thermostat in andersenThermostat.cu. I wonder if I can just modify this file to implement my own thermostat. How would you access/calculate the kinetic energy from within this file?
Thank you very much!
- Peter Eastman
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Question about Custom Thermostat and Barostat
Hi,
Doing this with any of the custom force classes would be hard. None of them is really designed for this sort of thing.
I suspect you could do it pretty easily with a CustomIntegrator. Use it to implement whatever integration method you want, but add a section where it rescales the velocities. A CustomIntegrator can easily compute the kinetic energy:
The other approach would be to write a plugin. You could then write a Force subclass that implements the calculation similar to the way AndersenThermostat does. If you think you might want to do that, you should start by reading the developer guide: http://docs.openmm.org/6.3.0/developerguide/index.html
Peter
Doing this with any of the custom force classes would be hard. None of them is really designed for this sort of thing.
I suspect you could do it pretty easily with a CustomIntegrator. Use it to implement whatever integration method you want, but add a section where it rescales the velocities. A CustomIntegrator can easily compute the kinetic energy:
Code: Select all
integrator.addComputeSum("ke", "m*v*v/2)
Peter
- Kailong Mao
- Posts: 16
- Joined: Fri Oct 02, 2015 8:56 am
Re: Question about Custom Thermostat and Barostat
Thank you very much Peter! I also wonder if creating and using my own CustomIntegrator will make my code a lot slower than using the original integrator with the Andersen thermostat implemented in CUDA. We hope to use our GPU's as much as possible.
- Peter Eastman
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Question about Custom Thermostat and Barostat
No, it should have almost no effect at all on speed. A CustomIntegrator still does all its expensive calculations on the GPU. And even the "expensive" parts of integration are usually very cheap compared to computing forces.
Peter
Peter
- Kailong Mao
- Posts: 16
- Joined: Fri Oct 02, 2015 8:56 am
Re: Question about Custom Thermostat and Barostat
Thank you very much Peter! Also, to implement a custom barostat, would it be possible to change the box size using a Custom Integrator?
- Peter Eastman
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Question about Custom Thermostat and Barostat
No, CustomIntegrator doesn't have any way to change the box size. If you want to do that, you'll need a plugin.
Peter
Peter
- Kailong Mao
- Posts: 16
- Joined: Fri Oct 02, 2015 8:56 am
Re: Question about Custom Thermostat and Barostat
Thank you! One last question: Based on what I read from the documentation, I can access individual forces f0, f1, f2,... etc from within CustomIntegrator. But I can't find what these forces are (harmonic bond forces? Coulomb interaction? van der Waals interaction?). Is there documentation somewhere that gives me this information? How do I know how many fx's there are in total? I assume that f = f0 + f1 + f2 + ...fn. Or can I only access some forces, but not all of them that together make up f.
- Peter Eastman
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Question about Custom Thermostat and Barostat
The numbers indicate force groups. By default everything is in group 0, but you can call setForceGroup() on any Force object to assign it to a different group. NonbondedForce also has a setReciprocalSpaceForceGroup() method, if you want to put the direct space and reciprocal space calculations into different groups. The most common reason for doing this with a CustomIntegrator is to implement multiple time step algorithms, where some forces are computed more often than others.
Peter
Peter
- Kailong Mao
- Posts: 16
- Joined: Fri Oct 02, 2015 8:56 am
Re: Question about Custom Thermostat and Barostat
Thanks a lot! I wrote the code for a CustomIntegrator and found out it works correctly on the "Reference" platform but not on "CUDA". The problem is that after I update one variable, for example, the kinetic energy, and want to update the temperature using the kinetic energy I just calculated, the code would still use the default kinetic energy I set when I added the global variable.
Here is how I add the variables to CustomIntegrator:
Here is how I add the computations (I removed all the irrelevant steps):
Every 10 steps, I make the integrator print out all the global variables with the following code:
So, I would assume that "ke" is calculated before "temp".
Using the reference platform, here is what I got after 20 steps. (I deleted all the irrelevant info)
Using CUDA, this is what I got:
So, the default value of "ke" is set to 0.1. It seems that it never gets updated (using CUDA). I wonder if there is a special update or sync step that I need to do after every computation if I want to use this code on CUDA?
By the way, here is how I set the platform:
Thanks a lot!
Here is how I add the variables to CustomIntegrator:
Code: Select all
self.integrator.addGlobalVariable("ke", 0.1)
self.integrator.addGlobalVariable("temp", 220.0)
Code: Select all
self.integrator.addComputeSum("ke", "m*v*v/2.0")
self.integrator.addComputeGlobal("temp", "2.0*ke/(df *gasconst)")
Code: Select all
totalVariables = self.integrator.getNumGlobalVariables();
for i in range(totalVariables):
name = self.integrator.getGlobalVariableName(i)
value = self.integrator.getGlobalVariable(i)
print(str(name) + " -> " + str(value))
Using the reference platform, here is what I got after 20 steps. (I deleted all the irrelevant info)
Code: Select all
ke -> 20705.7288735
temp -> 187.277879616
========================
ke -> 22046.7747921
temp -> 199.407287745
Code: Select all
ke -> 17099.9335938
temp -> 0.000904473736524
========================
ke -> 16971.2304688
temp -> 0.000904473736524
========================
By the way, here is how I set the platform:
Code: Select all
simulation = Simulation(pdb.topology, system, integrator, Platform.getPlatformByName("CUDA"))
- Peter Eastman
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Question about Custom Thermostat and Barostat
Could you post the full code for this? Or if it's too large you can email it to me. I'll try to reproduce it and figure out what's going on.
Peter
Peter