forcebalance and iamoeba

[Andrish Reddy]

Hi,

I am interested in studying water-ion interactions. I know that iamoeba is an accurate forcefield, but ion parameters are lacking. Is it possible to use forcebalance to generate ion parameters utilizing a direct polarization scheme (taking ion dipole parameters from amoeba as a start)? Where would one begin when taking on such a task?

Thanks,
Andrish

[Lee-Ping Wang]

Hi Andrish,

Thanks for getting in touch. Yes, I believe this is possible.

You would set up your calculation by deciding which force field parameters to optimize, then manually building the "targets" for fitting the parameters. ForceBalance provides some convenient modules / scripts for making targets, but the process is not fully automated like the optimization.

There are a couple of targets I could think of for doing this: (1) Generate QM data for water-ion interaction energies. (2) Gather experimental data on densities / compressibilities of ionic crystals.

Thanks,

- Lee-Ping

[Andrish Reddy]

Hi Lee-Ping,

In terms of editing the iamoeba.xml file to include initial parameters for say an Na+ ion, could I just copy the parameters from the amoeba ff for Na+ and add them under the appropriate directives: <AtomTypes>, <Residues> <Amoeba VdwForce> <Multipole> and <Polarize> ?

I could then select for example the vdw and polarizability as targets to optimize with forcebalance. Is that correct?

Thanks,
Andrish

[Lee-Ping Wang]

Yes, that sounds like the correct approach. When I was optimizing the iAMOEBA water parameters, I allowed the damping parameter to change (the second floating point number in the polarizability line). You might find this to be helpful.

Note that the vdW and polarizability are the parameters that you want to optimize. The "targets" are data sets from experimental measurements or ab initio calculations which the optimization tries to reproduce using the force field.

[Andrish Reddy]

Thanks Lee-Ping,

Yes that is what I meant to say regarding the parameters/targets. I am trying to replicate your results for the density and diffusion constant for iamoeba using openmm 6.2 and the CUDA platform with a stock-clocked gtx 970.

For the diffusion constant I get 2.90 E-5 cm^2/s, where I should get 2.54.
For the density @ 298K & 1bar I get 0.989 where it should be 0.997.

Here is my script:

Code: Select all

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout
Platform.loadPluginsFromDirectory(Platform.getDefaultPluginsDirectory())

# Initial System Parameters
struct = PDBFile('equil_nvt.pdb')
forcefield = ForceField('iamoeba.xml')
system = forcefield.createSystem(struct.topology, nonbondedMethod=PME,nonbondedCutoff=1*nanometer,ewaldErrorTolerance=0.0005,
useDispersionCorrection=True, vdwCutoff=1*nanometer, constraints=HBonds,rigidWater=False,polarization='direct')
system.addForce(MonteCarloBarostat(1*bar,298*kelvin,12))
integrator = LangevinIntegrator(298*kelvin, 1/picoseconds ,1*femtoseconds)
integrator.setConstraintTolerance(0.000001)
platform = Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed','CudaUseCpuPme': 'true', 'CudaDeviceIndex': '0'}
sim = Simulation(struct.topology, system, integrator,platform,properties)
sim.context.setPositions(struct.positions)

# Equilibrate NPT
print('Equilibrating NPT...')
sim.reporters.append(StateDataReporter(stdout, 500, progress=True,remainingTime=True, step=True, temperature=True, volume=True,potentialEnergy=True, density=True, totalSteps=1700000))
sim.step(200000) #200ps equilibration NPT
	

# Production Run
print('Running Production...')
sim.reporters.append(PDBReporter('run.pdb', 5000)) #Trajectory saved every 5ps
sim.reporters.append(StateDataReporter('run.log', 1, step=True,potentialEnergy=True,kineticEnergy=True,totalEnergy=True,temperature=True,volume=True,density=True))
sim.step(1500000) #1.5ns production run

Here are average energy values over the run:

P.Energy = -1.9783e+04
K.Energy = 4325.8
T.Energy = -1.5457e+04
Temp = 297.55
Vol = 15.125
Density = 0.98904

I have tried a smaller integration step (0.5fs) and changing the barostat frequency and constraint tolerance. I still cannot replicate your values. I am using a small box of 500 waters equilibrated via NVT (attached). I would appreciate your advice!

Thanks,
Andrish

[Lee-Ping Wang]

Hi Andresh,

I would set constraints=None and use a 0.5 fs time step. Everything else looks fine to me. Let me know if this brings your result into agreement. I typically use vdwCutoff=0.9*nanometers and nonbondedCutoff=0.7*nanometers (this is the real space PME cutoff), but I don't think that should affect the results.

Note that to reproduce the diffusion coefficient, you need to correct for finite size effects. This means that your raw simulated value will be lower than 2.54 - probably closer to 2.2. One way to apply this correction is to simulate several box sizes and fit a straight line of the simulated diffusion coefficient vs. 1/(box vector length).

Another aspect is that because you are using Langevin dynamics, you might get slightly slower diffusion coefficients due to the friction force acting on the atoms. I would recommend saving positions and velocities from your NPT simulation (e.g. every 10 ps for 1 ns, 100 initial conditions total), and then running short (~30 ps) NVE simulations from these. Now you can calculate the diffusion coefficient for each of the short trajectories and obtain your average value.

Thanks,

- Lee-Ping

PS: Once you reproduce the results, you can use the following multi-timestep integrator to speed things up.

Code: Select all

def MTSVVVRIntegrator(temperature, collision_rate, timestep, system, ninnersteps=4):
    """
    Create a multiple timestep velocity verlet with velocity randomization (VVVR) integrator.
    
    ARGUMENTS

    temperature (Quantity compatible with kelvin) - the temperature
    collision_rate (Quantity compatible with 1/picoseconds) - the collision rate
    timestep (Quantity compatible with femtoseconds) - the integration timestep
    system (simtk.openmm.System) - system whose forces will be partitioned
    ninnersteps (int) - number of inner timesteps (default: 4)

    RETURNS

    integrator (openmm.CustomIntegrator) - a VVVR integrator

    NOTES
    
    This integrator is equivalent to a Langevin integrator in the velocity Verlet discretization with a
    timestep correction to ensure that the field-free diffusion constant is timestep invariant.  The inner
    velocity Verlet discretization is transformed into a multiple timestep algorithm.

    REFERENCES

    VVVR Langevin integrator: 
    * http://arxiv.org/abs/1301.3800
    * http://arxiv.org/abs/1107.2967 (to appear in PRX 2013)    
    
    TODO

    Move initialization of 'sigma' to setting the per-particle variables.
    
    """
    # Multiple timestep Langevin integrator.
    for i in system.getForces():
        if i.__class__.__name__ in ["NonbondedForce", "CustomNonbondedForce", "AmoebaVdwForce", "AmoebaMultipoleForce"]:
            # Slow force.
            logger.info(i.__class__.__name__ + " is a Slow Force\n")
            i.setForceGroup(1)
        else:
            logger.info(i.__class__.__name__ + " is a Fast Force\n")
            # Fast force.
            i.setForceGroup(0)

    kB = BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA
    kT = kB * temperature
    
    integrator = CustomIntegrator(timestep)
    
    integrator.addGlobalVariable("dt_fast", timestep/float(ninnersteps)) # fast inner timestep
    integrator.addGlobalVariable("kT", kT) # thermal energy
    integrator.addGlobalVariable("a", np.exp(-collision_rate*timestep)) # velocity mixing parameter
    integrator.addGlobalVariable("b", np.sqrt((2/(collision_rate*timestep)) * np.tanh(collision_rate*timestep/2))) # timestep correction parameter
    integrator.addPerDofVariable("sigma", 0) 
    integrator.addPerDofVariable("x1", 0) # position before application of constraints

    #
    # Pre-computation.
    # This only needs to be done once, but it needs to be done for each degree of freedom.
    # Could move this to initialization?
    #
    integrator.addComputePerDof("sigma", "sqrt(kT/m)")

    # 
    # Velocity perturbation.
    #
    integrator.addComputePerDof("v", "sqrt(a)*v + sqrt(1-a)*sigma*gaussian")
    integrator.addConstrainVelocities();
    
    #
    # Symplectic inner multiple timestep.
    #
    integrator.addUpdateContextState(); 
    integrator.addComputePerDof("v", "v + 0.5*b*dt*f1/m")
    for innerstep in range(ninnersteps):
        # Fast inner symplectic timestep.
        integrator.addComputePerDof("v", "v + 0.5*b*dt_fast*f0/m")
        integrator.addComputePerDof("x", "x + v*b*dt_fast")
        integrator.addComputePerDof("x1", "x")
        integrator.addConstrainPositions();        
        integrator.addComputePerDof("v", "v + 0.5*b*dt_fast*f0/m + (x-x1)/dt_fast")
    integrator.addComputePerDof("v", "v + 0.5*b*dt*f1/m") # TODO: Additional velocity constraint correction?
    integrator.addConstrainVelocities();

    #
    # Velocity randomization
    #
    integrator.addComputePerDof("v", "sqrt(a)*v + sqrt(1-a)*sigma*gaussian")
    integrator.addConstrainVelocities();

    return integrator

[Andrish Reddy]

Thanks for the quick response

I will remove restraints and let you know how it went.

Andrish

[Andrish Reddy]

Hi Lee-Ping,

I was able to get a better density approximation by removing the constraints. I am trying to reproduce the diffusion constant by the method you suggested. I wanted to know whether it's fine to pull the velocities from getVelocities in my NPT system and calling setVelocities without correcting them for NVE? I'm using NPT with the LangevinIntegrator and will be doing the NVE with the VerletIntegrator. I assume that since they are both leap-frog integrators I won't have to adjust the velocities? If I use the multi-step integrator that you suggested, then I would have to adjust the velocities, correct?

Thanks,
Andrish

[Lee-Ping Wang]

Hi Andrish,

You are correct that the velocities do not correspond to the current time. In the LangevinIntegrator and VerletIntegrator, the current velocities are half a time step behind the coordinates (tx - 0.5*timestep). I believe that in the MTS-VVVR integrator, the velocities are current. Thus, if you use MTS-VVVR to obtain initial conditions for a VerletIntegrator simulation, you need to propagate the velocities backward by half a time step.

I have some code here for writing restart files that you could look at. It basically makes sure all velocities are current before writing the restart file, and for leapfrog integrators it propagates the velocities backward by half a time step before proceeding.

https://github.com/leeping/OpenMM-MD/bl ... MD.py#L782

Thanks,

- Lee-Ping

[Andrish Reddy]

Hi Lee-Ping,

Thanks for your help thus far. I have one more question as to the information required by the qdata file. If I am optimizing the ion parameters: vdw-epsilon, vdw-sigma, polarizability, and thole, then how many snapshots would I need from the qm run? Also what would you suggest in terms of water-ion cluster sizes for a range of ions (NA+, CL-, K+, etc)?

Regards,
Andrish