Drude Force Field

[moradza]

Is there any tutorial available for charm36 with water and ions using openmm. Can someone share the link? Thanks,

[peastman]

Are you asking about CHARMM36 (which is a point charge force field) or the Drude force field? Assuming you mean the latter, I'm currently working on converting the latest version of the Drude force field to OpenMM's format. In the mean time, the easiest way of doing it is to build your system with CHARMM and then load the psf file into OpenMM.

[moradza]

I mean the Drude force field. Should I change anything particular in the integrator and other details for the Drude force field?

[peastman]

Yes, you should use one of the integrators that's specifically designed for Drude particles. Usually that means DrudeLangevinIntegrator.

Once I have the new version of the force field converted, I'll put together some proper documentation on using it.

[moradza]

Do you have an example of the SWM4 with ions? I am new to open-mm and would appreciate if I get pdb, psf, py, and xml files for it?

[moradza]

I was wondering if the sigma and epsilon are normal for water based on charm_polar_2013.xml:

forces = { force.__class__.__name__ : force for force in system.getForces() }
nbforce = forces['NonbondedForce']
for cnt, at in enumerate(modeller.topology.atoms()):
    print(at.residue, cnt)
    if at.name[0]=='D':
        system.setParticleMass(cnt, mass=0.4)
    [charge, sigma, epsilon] = nbforce.getParticleParameters(cnt)
    print(cnt, sigma, epsilon, charge)

<Residue 0 (HOH) of chain 0> 0 0 1.0 nm 0.0 kJ/mol 1.71639 e <Residue 0 (HOH) of chain 0> 1 1 1.0 nm 0.0 kJ/mol 0.55733 e <Residue 0 (HOH) of chain 0> 2 2 1.0 nm 0.0 kJ/mol 0.55733 e <Residue 0 (HOH) of chain 0> 3 3 1.0 nm 0.0 kJ/mol -1.11466 e <Residue 0 (HOH) of chain 0> 4 4 1.0 nm 0.0 kJ/mol -1.71639 e

[peastman]

That force field doesn't rely on standard combining rules. It therefore uses a CustomNonbondedForce to implement Lennard-Jones, with lookup tables for sigma and epsilon. That's why the NonbondedForce has all the epsilons set to 0. Those interactions are computed by a different force.

[moradza]

Okay, can I access those parameters? More importantly, since it is running fine, can I be sure it is correct?

[peastman]

The CustomNonbondedForce defines a per-particle parameter called type and two tabulated functions called sigma and epsilon. For any pair of atoms, use the atoms' types as indices into the tables.

>>> force = [f for f in system.getForces() if isinstance(f, CustomNonbondedForce)][0]
>>> for i in range(force.getNumTabulatedFunctions()):
...   print(i, force.getTabulatedFunctionName(i))
... 
0 epsilon
1 sigma
>>> epsilon = force.getTabulatedFunction(0)
>>> sigma = force.getTabulatedFunction(1)
>>> print(force.getParticleParameters(0))
(87.0,)
>>> print(force.getParticleParameters(1))
(67.0,)
>>> xsize, ysize, values = sigma.getFunctionParameters()
>>> values[87+xsize*67]
0.204016

Particle 0 has type 87. Particle 1 has type 67. For their interaction, sigma is 0.204106. For more details see http://docs.openmm.org/latest/api-python/generated/simtk.openmm.openmm.CustomNonbondedForce.html#simtk.openmm.openmm.CustomNonbondedForce and http://docs.openmm.org/latest/api-c++/generated/OpenMM.Discrete2DFunction.html#_CPPv2N6OpenMM18Discrete2DFunctionE.

[moradza]

Thanks @peastman for the help. I have one more question, I want to simulate fix carbon nanotube but I am not sure which type or residue should I select. Can you help me with this as well or point to some references.

[peastman]

I don't think the standard force field has parameters for carbon nanotubes. Someone may have developed some, but it's not a subject I know about.

[moradza]

Because what I am looking to do is to freeze CNT, I don't need bonds, angle, and dihedral. What I need is adding a new particle type with only sigma and epsilon (no charge or polarizability) and using it with combination rule for other interactions. I appreciate if you help with it?

[peastman]

No idea. Parametrizing a new class of molecules is a research project. Chances are good that someone has done it, but I'm not familiar with that work. You'll need to search the literature to see if you can find suitable parameters.

[moradza]

Thanks, I found sigma and epsilon. Now, I am trying to figure out how to add them to the charm_polar_2013.xml