Page 1 of 2
Question about Custom Thermostat and Barostat
Posted: Thu Jan 28, 2016 2:16 pm
by zebra_123456
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!
Re: Question about Custom Thermostat and Barostat
Posted: Thu Jan 28, 2016 2:33 pm
by peastman
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:
Code: Select all
integrator.addComputeSum("ke", "m*v*v/2)
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
Re: Question about Custom Thermostat and Barostat
Posted: Thu Jan 28, 2016 2:44 pm
by zebra_123456
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.
Re: Question about Custom Thermostat and Barostat
Posted: Thu Jan 28, 2016 3:29 pm
by peastman
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
Re: Question about Custom Thermostat and Barostat
Posted: Mon Feb 01, 2016 1:35 pm
by zebra_123456
Thank you very much Peter! Also, to implement a custom barostat, would it be possible to change the box size using a Custom Integrator?
Re: Question about Custom Thermostat and Barostat
Posted: Mon Feb 01, 2016 1:52 pm
by peastman
No, CustomIntegrator doesn't have any way to change the box size. If you want to do that, you'll need a plugin.
Peter
Re: Question about Custom Thermostat and Barostat
Posted: Tue Feb 02, 2016 11:23 am
by zebra_123456
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.
Re: Question about Custom Thermostat and Barostat
Posted: Tue Feb 02, 2016 12:23 pm
by peastman
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
Re: Question about Custom Thermostat and Barostat
Posted: Wed Feb 03, 2016 12:33 pm
by zebra_123456
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:
Code: Select all
self.integrator.addGlobalVariable("ke", 0.1)
self.integrator.addGlobalVariable("temp", 220.0)
Here is how I add the computations (I removed all the irrelevant steps):
Code: Select all
self.integrator.addComputeSum("ke", "m*v*v/2.0")
self.integrator.addComputeGlobal("temp", "2.0*ke/(df *gasconst)")
Every 10 steps, I make the integrator print out all the global variables with the following code:
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))
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)
Code: Select all
ke -> 20705.7288735
temp -> 187.277879616
========================
ke -> 22046.7747921
temp -> 199.407287745
Using CUDA, this is what I got:
Code: Select all
ke -> 17099.9335938
temp -> 0.000904473736524
========================
ke -> 16971.2304688
temp -> 0.000904473736524
========================
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:
Code: Select all
simulation = Simulation(pdb.topology, system, integrator, Platform.getPlatformByName("CUDA"))
Thanks a lot!
Re: Question about Custom Thermostat and Barostat
Posted: Wed Feb 03, 2016 3:29 pm
by peastman
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