FES from biased trajectories, too!

Vladas Oleinikovas         uccavol@ucl.ac.uk
Guillermo Pérez-Hernández  guille.perez@fu-berlin.de

We do that by passing correct weights for each conformation in our trajectory.

We exemplify this with metadynamics trajectory of the ever popular test system - alanine dipeptide. The 2ns simulation used phi and psi angles as CVs, and the corresponding weight was derived using Tiwary and Parrinello reweighting algorithm (JPCB 2014, doi:10.1021/jp504920s).

Try commenting out the weights to notice the effect on the estimated depth of the basins.

import molpx
import numpy as np
%matplotlib ipympl

# Topology
top = molpx._molpxdir(join='notebooks/data/ala2.pdb')

# Generated during the metadynamics run (extra reversed copy to have a list of files)
MD_trajfiles = [molpx._molpxdir(join='notebooks/data/ala2.meta.xtc'),
               molpx._molpxdir(join='notebooks/data/ala2.meta.reversed.xtc')]

# COLVAR files containing CV projections and logarithm of the corresponding weight
colvar_files = [molpx._molpxdir(join='notebooks/data/ala2.meta.CV.txt'),
                molpx._molpxdir(join='notebooks/data/ala2.meta.CV.reversed.txt')]

# Load logweights and turn them into proper weights
weights = [np.exp(np.loadtxt(ifile)[:,7]) for ifile in colvar_files]

mpx_wdg_box = molpx.visualize.FES(MD_trajfiles,
                                  top,
                                  colvar_files,
                                  proj_idxs=[1,2],
                                  nbins=50,
                                  n_sample=200,
                                  proj_labels='CV',
                                  weights=weights, # this is important for biased trajectories!
                                 )
mpx_wdg_box

We can also visualize the weights together with the CV-trajectories

mpx_wdg_box = molpx.visualize.traj(MD_trajfiles,
                                  top,
                                  colvar_files,
                                  max_frames=10000,
                                  proj_idxs=[1, 2, 7],
                                  panel_height=1,
                                  proj_labels=['$CV_1$','$CV_2$','log(weights)'],
                                  )
mpx_wdg_box