Umbrella Sampling Method For Ligand Water System

1.Software and Method Abstract

In this article, we will mainly introduce how to perform umbrella sampling, free energy calculations, and construct a free energy landscape (FEL) for a complex system such as ligand-water using the GROMACS software. GROMACS is a widely used molecular dynamics simulation software that can simulate various systems, including biomolecular systems, nanomaterials, homogeneous systems, and multi-molecular self-assembly. It also supports algorithms such as adaptive biasing force (ABF) and umbrella sampling. However, umbrella sampling is just one form of enhanced sampling in molecular dynamics simulations. Personally, I feel that CF-SMD (Constant Force Steered Molecular Dynamics) may be more suitable for investigating the dissociation processes between proteins and peptides/ligands. Next, we will discuss the relevant background and specific details of the simulation.

2.Software compile and installation guide

The software we mainly needed is as follows:

1.Gromacs (installation and compile guide):

For more details, please visit the gromacs official websites follows:
https://manual.gromacs.org/2023.2/install-guide/index.html

And i strongly recommand you to compile gromacs with gpu(CUDA) support.

tar xfz gromacs-2023.2.tar.gz
cd gromacs-2023.2
mkdir build
cd build
cmake .. -DGMX_BUILD_OWN_FFTW=ON -DREGRESSIONTEST_DOWNLOAD=ON -DGMX_GPU=CUDA
make -j
make check
sudo make install -j
source /usr/local/gromacs/bin/GMXRC

2.Packmol:

For MacOS Users:

brew install packmol

For CentOS or Ubuntu Users:

sudo dnf install packmol
/
sudo apt install packmol

Because packmol is only a serial software but it is really fast and customizable, we prefer not compile this software from the source code but install it from dnf/apt in linux or brew in the MacOS

3.Background

Organosulfates (OSs) are well-known and ubiq-uitous constituents of atmospheric aerosol particles and have been used as secondary organic aerosol markers in many field studies. Hence, it is imperative to understand the formation of OS species in the atmosphere. Recently, hydroxy acids (HAs) and hydroxy acid sulfates have been extensively detected in the atmospheric environment. However, the reaction mechanism of HAs to form OSs is much less understood. In this work, we have mainly investigated the reaction of typical α-HAs, including glycolic acid (GA) and lactic acid (LA), and SO3 at the liquid aerosol surface using quantum chemistry calculations and Born−Oppenheimer molecular dynamics simulations. The OH group orientation of α-HAs at the air−water interface is found to exert a significant impact on the formation of OSs. The OH group pointing to the gas phase is obviously beneficial to the formation of OSs. Two key factors are discovered important to the reaction of α-HAs adsorbed on the liquid surface with SO3: (a) the exposure position of the active site to the gas phase and (b) the reactivity of the exposed site to the attracted SO3 molecule. Moreover, we found that the air−water interface exerts a significant influence on the physicochemical behaviors of GA and LA, especially on their OH group orientation, and thus leads to their different properties for the SO3 colliding reaction. The presented reaction mechanism provides a new feasible pathway for the production of OSs at the liquid aerosol surface, which may have important impacts on the formation of organic aerosols.

Paper Reference DOI: 10.1021/jacs.2c05807

As follows is the initial structure which can be built with the software packmol:

4.Model Build Ways and Parameters Setting

First, we need to construct a water box with dimensions of [40, 40, 40]. We will add 5 nm vacuum layers above and below the box, and apply periodic boundary conditions. The main purpose of adding vacuum layers to the water box is to prevent the influence of mirror images due to the periodicity of the system during the simulation. We use the software VMD to implement the setup of the vacuum layers. Additionally, we have written a script which can be used in the software VMD to facilitate obtaining the central reference atom index of a specific object when setting parameters for umbrella sampling in later stages. The script that is based on tcl is as follows:

# sqr(x-20)+sqr(y-20)+sqr(z-20) < sqr(2)
proc atmlab {range id} {
set sel [atomselect $id $range]
set k 0
foreach i [$sel list] {
label add Atoms $id/$i
label textformat Atoms $k { %i }
label textoffset Atoms $k { -0.11 -0.0055 }
incr k
}
$sel delete
label textsize 1.2
}

proc lab {} {
atmlab all [molinfo top]
}

Packmol build progress:

First of all, we use the following configuration files to build the basic water box:

Packmol tutorial reference: https://m3g.github.io/packmol/

tolerance 2.0
add_box_sides 1.2
output mix.pdb
structure water.pdb 
  number 2165
  inside box 0. 0. 0. 40. 40. 40.
end structure

The final water box is shown in the following figure:

Of course, after building the water box, we should allocate and build the vacuum layer through the two commands pbc set {x y z} and [atomselect top all] moveby {0 0 n} in vmd, and then write the corresponding topological file and perform energy minimization and nvt simulation pre-balancing operations and then the same operation for the ligand-water complex reaction system and the following topology file is a example for everyone who want to preform this kind simulation:

Water_slab topology example:

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/tip3p.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

[ system ]
Water_Slab_Nvt

[ molecules ]
; Compound        nmols
SOL  2165

Ligand_Water complex topology example:

The ligand topology is generated by acpype with its atom charge which is calculated by ORCA based on density function with def2-TZVP basis-set.

#include "oplsaa.ff/forcefield.itp"

[ atomtypes ]
;name   bond_type     mass     charge   ptype   sigma         epsilon       Amb
 c3       c3          0.00000  0.00000   A     3.39967e-01   4.57730e-01 ; 1.91  0.1094
 hc       hc          0.00000  0.00000   A     2.64953e-01   6.56888e-02 ; 1.49  0.0157
 h1       h1          0.00000  0.00000   A     2.47135e-01   6.56888e-02 ; 1.39  0.0157
 oh       oh          0.00000  0.00000   A     3.06647e-01   8.80314e-01 ; 1.72  0.2104
 ho       ho          0.00000  0.00000   A     0.00000e+00   0.00000e+00 ; 0.00  0.0000
 c        c           0.00000  0.00000   A     3.39967e-01   3.59824e-01 ; 1.91  0.0860
 o        o           0.00000  0.00000   A     2.95992e-01   8.78640e-01 ; 1.66  0.2100

#include "la_GMX.itp"
#include "oplsaa.ff/tip3p.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

[ system ]
equil_COM

[ molecules ]
; Compound        nmols
SOL  2165
la   1 

Next, I will provide an mdp file of gromacs for umbrella sampling for your reference. It is worth noting that umbrella sampling does not require a long time. Different times are set according to different systems, and the simulation time of about 800ps is set to our current system, and the step length is 2fs.

define = -DPOSRES_WATER 
integrator               = md
dt                       = 0.002     ; 2 fs
nsteps                   = 400000   ; 400 ps

nstenergy                = 500
nstlog                   = 500
nstxout-compressed       = 200

continuation             = yes
constraint-algorithm     = lincs
constraints              = h-bonds

cutoff-scheme            = Verlet

coulombtype              = PME
rcoulomb                 = 1.2

vdwtype                  = Cut-off
rvdw                     = 1.2
DispCorr                 = EnerPres

tcoupl                   = v-rescale
tc-grps                  = System
tau-t                    = 0.1
ref-t                    = 298.15

pcoupl                   = Parrinello-Rahman
pcoupltype               = isotropic
tau_p                    = 2.0
compressibility          = 4.46e-5
ref_p                    = 1.0
refcoord_scaling = com

;
freezegrps  = UNK 
freezedim   = Y Y N
;pull code

;Pull code
pull                 = yes 
pull_ngroups         = 2 
pull_ncoords         = 1 
pull_group1_name     = UNK
pull_group2_name     = SOL 
pull_coord1_type     = umbrella 
pull_coord1_geometry = distance
pull_coord1_groups   = 1 2 
pull_coord1_dim      = N N Y
;pull-coord1-vec      = 0 0 1
pull_coord1_rate     = -0.005
pull_coord1_k        = 1000
pull_coord1_start    = yes

Meanwhile, when your reference group is SOL and WATER or any other reference system, you can set the reference atoms in the reference system by configuring the following two parameters: “pull-pbc-ref-prev-step-com=yes“ and “*pull-group1-pbcatom=Number(A)*“. After performing umbrella sampling, you can extract all trajectory frames within 800ps and filter out 32 structures as windows for another round of umbrella sampling simulations. The resulting structures are shown below.

Keith's Blog