Question about Custom Thermostat and Barostat

[Kailong Mao]

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!

[Peter Eastman]

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

[Kailong Mao]

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]

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

[Kailong Mao]

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]

No, CustomIntegrator doesn't have any way to change the box size. If you want to do that, you'll need a plugin.

Peter

[Kailong Mao]

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]

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

[Kailong Mao]

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!

[Peter Eastman]

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