I am new to molecular simulation, so please excuse what is likely a very silly question. I am using the python package openmmtools to run a simple molecular simulation problem using the metropolis algorithm. Though my question may not be programming related at all.

To know if the algorithm is doing anything, I first load a default system of particles (alanine-dipeptide in vacuum), add some uniform random disturbance to every position, then let the algorithm run to see if it brings me back to the starting system position.

It looks like the metropolis moves decrease energy, as expected, but do not decrease the euclidean distance of each particle back to its starting position. The following graph is what I see when I measure these quantities. The big initial jump is when I add a disturbance to the positions of every particle in the system. I have even subtracted the mean of the particle positions from each particle to make sure that it isn't just that the entire system is shifting away from the starting position.

distance and energy during metropolis algorithm for simulation of 22 molecule system

Why do the Euclidean distances of each particle not decrease with the energy. In fact, they seem to increase as the simulation continues, which means the particles are growing farther and farther apart, which doesn't seem right.

I expect is something very silly on my end. I hope the confusion is clear.

Here's the code that gave me this graph, just in case.

from openmmtools import mcmc, testsystems, states, cache
from simtk import openmm
from simtk import unit
import numpy as np
import matplotlib.pyplot as plt

# user parameters
timesteps = 100
draw_init_final_positions = False
plot_distance_energy = True

class AddGaussianVector(mcmc.MetropolizedMove):
    def __init__(self, **kwargs):
        super(AddGaussianVector, self).__init__(**kwargs)

    def _propose_positions(self, initial_positions):
        # displacement = unit.Quantity(np.array([1, 1, 1]), initial_positions.unit)
        mean, var = 0, .1
        x_prop, y_prop, z_prop = np.random.normal(mean, var), np.random.normal(mean, var), np.random.normal(mean, var)
        displacement = unit.Quantity(np.array([x_prop, y_prop, z_prop]), initial_positions.unit)
        return initial_positions + displacement

def distance(x, correct):
    for k in range(3):
        x[:, k] -= np.mean(x[:, k])
        correct[:, k] -= np.mean(correct[:, k])
    noise = np.subtract(x, correct)
    return np.sqrt(np.sum(np.square(noise)))

platform = openmm.Platform.getPlatformByName('CUDA')
cache.global_context_cache.platform = platform

# Create the initial state (thermodynamic and microscopic) for an alanine dipeptide system in vacuum.
alanine = testsystems.AlanineDipeptideVacuum()
sampler_state = states.SamplerState(alanine.positions)
thermodynamic_state = states.ThermodynamicState(alanine.system, 300 * unit.kelvin)
e, d = [], []
context_cache = cache.global_context_cache
context, unused_integrator = context_cache.get_context(thermodynamic_state)
sampler_state.apply_to_context(context, ignore_velocities=True)
d.append(distance(sampler_state.positions.copy(), sampler_state.positions.copy()))

# draw initial sampler_state
if draw_init_final_positions:
    ax = plt.axes(projection='3d')
    xdata = sampler_state.positions[:, 0]
    ydata = sampler_state.positions[:, 1]
    zdata = sampler_state.positions[:, 2]
    ax.scatter3D(xdata, ydata, zdata, cmap='Greens')
    plt.title('Initial particle conformation')

noise = np.zeros((sampler_state.n_particles, 3))
correct_state = sampler_state.positions.copy()
for i, molecule in enumerate(sampler_state.positions):

    # add noise to correct state
    noise_low, noise_high = -1, 1
    sampler_state.positions[i][0] += unit.quantity.Quantity(value= np.random.uniform(noise_low, noise_high), unit=unit.nanometer)
    sampler_state.positions[i][1] += unit.quantity.Quantity(value= np.random.uniform(noise_low, noise_high), unit=unit.nanometer)
    sampler_state.positions[i][2] += unit.quantity.Quantity(value= np.random.uniform(noise_low, noise_high), unit=unit.nanometer)

print('Correct_state: ', correct_state)
print('Noisey_state', sampler_state.positions.copy())
print('INITIAL DIST AND ENERGY: {} nm {} energy'.format(d[0], e[0]))

sampler_state.apply_to_context(context, ignore_velocities=True)

# Create an update MCMC move that brings us back to the initial configuration.
current_state = sampler_state.positions.copy()
total_accepted, total_proposed = 0, 0

print('TOTAL PARTICLES: ', sampler_state.n_particles)
for ii in range(timesteps):
    for jj in range(sampler_state.n_particles):
        move = AddGaussianVector(atom_subset=[jj])
        move.apply(thermodynamic_state, sampler_state)
        if move.n_accepted == 1:
            total_accepted += 1
        total_proposed += 1

    current_state = sampler_state.positions.copy()
    # measure euclidean distance
    d.append(distance(current_state, correct_state))

    # measure energy
    sampler_state.apply_to_context(context, ignore_velocities=True)

fig, ((ax1),(ax2)) = plt.subplots(2, 1, sharex=True)
ax1.plot(d, label='Euclidean dist')
ax1.set_title('Euclidean Distance to ground state')
ax1.set_ylabel('Euclidean Distance (nm)')
ax2.plot(e, label='Energy')
ax2.set_title('Energy of conformation')
ax2.text(2, 10, r'Original Energy {}'.format(e[0]), fontsize=15)
  • 1
    $\begingroup$ Is it rotating? $\endgroup$ Mar 7 at 4:16
  • 1
    $\begingroup$ This is more suited for Matter Modelling SE imho: mattermodeling.stackexchange.com $\endgroup$
    – S R Maiti
    Mar 7 at 10:57
  • $\begingroup$ @karstenTheis If it were rotating, I would think the euclidean distance of all particles to 0,0,0 would remain constant right? I see, I'll try posting it to Matter Modelling! $\endgroup$
    – Squishy
    Mar 7 at 22:56
  • $\begingroup$ This was answered in Matter Modelling SE mattermodeling.stackexchange.com/questions/4458/… $\endgroup$
    – Squishy
    Mar 11 at 19:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.