test particle insertion method

[Alexej Jerschow]

Can somebody suggest a good way to implement a test-particle insertion method (e.g. Widom's method) to calculate chemical potentials. In particular, I am looking to evaluate the chemical potential for small molecules in different solvents. What would be the right way to go about this, so that it is also reasonably efficient in order to perform many test insertions.
Thank you
Alexej

[Peter Eastman]

I'm not familiar with that algorithm. Can you describe how it works?

[Alexej Jerschow]

One does for example an NVT simulation, and then at random snapshots inserts a particle (molecule) at a random position and calculates the free energy difference between the N and the N+1 configuration. The particle is then removed again and this insertion process is repeated many times to get good statistics (e.g. up to 100,000 times). The chemical potential \mu is then given by $\beta\mu = \ln \langle \exp ( -\beta \Delta U ) \rangle$, where the average is done over all the test insertions. The method is described for example in Section 7.2.1 in Dan Frankel's book.

I believe GROMACS implements it by taking an existing N-particle simulation as input and then rerunning over the trajectory with these random insertions of a molecule with random orientations and averaging over the energies due to the random insertions ('tpi' parameter https://manual.gromacs.org/documentatio ... tions.html).

[Peter Eastman]

It sounds like the question is how to insert a particle, calculate the energy, and then remove it again? Since it's an isolated particle, I assume it isn't involved in any bonded interactions, only nonbonded ones?

The simplest way of doing that is to include the extra particle in the system from the beginning, but set its parameters so it doesn't interact with anything else. Something like this:

Code: Select all

extra_particle = system.addParticle(1.0)
nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]
nonbonded.addParticle(0.0, 1.0, 0.0)

To compute the energy difference from adding the particle, you set its parameters to make it interact and call updateParametersInContext() to copy the updated particle to the existing context. Something like this:

Code: Select all

# Compute initial energy and positions

state = context.getState(getEnergy=True, getPositions=True)
e1 = state.getPotentialEnergy()
pos = state.getPositions()

# Add the particle

nonbonded.setParticleParameters(extra_particle, charge, sigma, epsilon)
nonbonded.updateParametersInContext(context)
pos[extra_particle] = test_position
context.setPositions(pos)

# Compute final energy

e2 = context.getState(getEnergy=True).getPotentialEnergy()

# Remove the particle again

nonbonded.setParticleParameters(extra_particle, 0.0, 1.0, 0.0)
nonbonded.updateParametersInContext(context)

[Alexej Jerschow]

Super! Thank you so much! I will try this.

[Alexej Jerschow]

Thank you. In terms of efficiency, it seems that it will be best to first run a full simulation without the particle (or with the nonbonding interactions turned off), and then loop through the trajectory, and recompute the change in energy for a subset of time steps and random positions of the test particle.

I did not find examples about how one could pick out particular time steps from a finished trajectory. What would be the best way to do that?

Thank you
Alexej