Software package collection for dynamical Girsanov path reweighting
Dynamical reweighting methods are becoming more and more important to study complex mechanisms of chemical processes, whose direct simulations of the full state-space at the correct statistical weight are impossible. In order to gain a better understanding of the applicability of Girsanov path reweighting to chemical processes in complex environments, a broad network between already high-quality methods must be created. The article presents the basis for how Girsanov path reweighting can be integrated into any simulation and analysis software. We also provide a ready-to-use simulation program that enables on-the-fly estimation of reweighting factors and a unified procedure for estimating a reweighted dynamical model from this data.
- simulation programs: openmmtools
- reporter function: ReweightingReporter
- analysis package: deeptime
-
download existing software with reweighting modification
a) go to the respective GitHub repository
b) clone the respective package by
- copying the "HTTPS clone URL" link that appears after selecting the code icon and
- run
git clone [url]
in the directory you want to download it to
-
create new environment (e.g. RWGHTSoftware) with
conda create -n RWGHTSoftware python=3.11
-
install the package and all corresponding dependencies, e.g.
conda install -c anaconda matplotlib
,conda install -c anaconda spyder
a) openmmtools: install the openmmtools extension in cloned openmmtools folder via
pip insatll .
, further possible dependencies areconda install -c conda-forge openmm
conda install -c conda-forge mdtraj
conda install -c conda-forge netcdf4
conda install -c conda-forge mpiplus
conda install -c omnia pymbar
conda install -c numba numba
b) deeptime: install the deeptime extension in cloned deeptime folder via
pip insatll .
-
create scripts for simulation and analysis with Girsanov reweighting or use the setup and templates in reweightingtools
Example studies are discussed in Implementation of Girsanov reweighting in OpenMM and Deeptime.
a) The propagation of a system, cf. Fig. 3 and Fig. 5 b-d, according to a Langevin symmetric splittig scheme at a unbiased (openMM_simulation()
) or biased potential (openMM_biased_simulation()
) can be realized with the help of some top-level functions of the reweightingtools simulation-api. A detailed description of the setup can be found in the respective function documentation, e.g. here for the biased simulation:
# Import
from reweightingtools.simulation import openMM_biased_simulation, restraints_harmonic_force
from openmm.unit import nanometer
# MD input
forcefield='amber99.ff'
equisteps= 5000
nsteps= 5000000000
integrator_scheme= 'LangevinWithGirsanov'
integrator_splitting='R V O V R'
nstxout=100
solvent='HOH'
restraints_kwargs=dict()
restraints_kwargs['multiple_restraints']=True
restraints_kwargs['atom_name']=[('CA'),('CL')]
restraints_kwargs['force_constants']=[200000.0,500000.0]
restraints_kwargs['restraint_x']=[True, True]
restraints_kwargs['restraint_y']=[True, True]
restraints_kwargs['restraint_z']=[False, True]
# Run openMM simulation
openMM_biased_simulation(forcefield,
equisteps,
nsteps,
integrator_scheme,
integrator_splitting,
nstxout,
gro_input='ClCaCl',
plumed_file='plumedCOLVAR.dat',
externalForce_file='bias_05z.txt',
restraints_func= [restraints_harmonic_force,restraints_harmonic_force],
restraints_kwargs=restraints_kwargs,
fix_system=True,
#restart=True,
dcd_output=False,
pdp_output=False,
data_output=False
)
b) Dynamic properties of a system, such as implicit timescales (deeptime_implied_timescales()
) or dominant eigenvectors (deeptime_eigenvectors()
) of a Markov state models (MSM) (cf. Fig. 4 or Fig. 5 b-d) can be analyzed with the help of some top-level functions of the reweightingtools analysis-api. For more customization of the analysis see the deeptime documentation. Here an example for evaluating implied timescales:
# Import
from reweightingtools.analysis.discretisation import *
from reweightingtools.analysis.api import *
# Set input
cwd='/home/project/'
analysisMD=cwd+'analysis_MD'
gridsize=100
reversible=True
nstxout=100
lagtimes=[1,3,5,7,9,10]
max2D, min2D = np.array([ 1.5 , 3.5]), np.array([-3.5, -1.5])
traj = np.load(MD_directory+'/traj.npy')
dtraj = trajectory2D_MBP(traj,
gridsize,
min2D,
max2D)
np.save(analysisMD+'/discretized_'+str(gridsize),dtraj)
# Prepare directory, eg. /discretized, containing discretized
# trajectory, g and M factor output from openMM simulation using
# the LangevinSplittingGirsanov integrator
deeptime_implied_timescales(lagtimes=lagtimes,
gridsize=gridsize,
reversible=reversible,
reweighting=True,
its_file='/its_rwght',
analysisMD=analysisMD,
dtraj_file='/discretized',
gF_file='/g_factors',
MF_file='/M_factors',
save_its=True)
# Visualize reweighted implied timescales
its_rwght = np.load(analysisMD+'/its_rwght_'+str(gridsize)+'.npy')
i= 0 # slowest process
plt.plot(np.array(lagtimes),its_rwght[:,i],label='rwght')
To get the first two dominant eigenvectors of a from reweightingtools.analysis.api import *
. Prepare a directory analysisMD:str='analysis_MD'
, with the discretized trajectory dtraj_file:str='/discretized'
, as well with both gF_file:str='/g_factors'
, MF_file:str='/M_factors'
) from openMM simulation using LangevinSplittingGirsanov
integrator.
deeptime_eigenvectors(lagtime=150,
gridsize=36,
number_eigenvectors=2,
reversible=True,
reweighting=True,
analysisMD:str='analysis_MD',
eigenvector_file:str='/evecs',
dtraj_file:str='/discretized',
gF_file:str='/g_factors',
MF_file:str='/M_factors'
)