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
test particle insertion method
- Peter Eastman
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: test particle insertion method
I'm not familiar with that algorithm. Can you describe how it works?
- Alexej Jerschow
- Posts: 22
- Joined: Tue Mar 10, 2020 7:31 am
Re: test particle insertion method
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).
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
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: test particle insertion method
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:
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:
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)
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
- Posts: 22
- Joined: Tue Mar 10, 2020 7:31 am
Re: test particle insertion method
Super! Thank you so much! I will try this.
- Alexej Jerschow
- Posts: 22
- Joined: Tue Mar 10, 2020 7:31 am
Re: test particle insertion method
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
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