GROMACS version: 2020.2
GROMACS modification: No
I am trying to execute the free energy solvation calculation as described here “Free Energy Calculations”. In the job.sh file “http://www.mdtutorials.com/gmx/free_energy/Files/job.sh”, there is a loop from 0 to 20. the “init_lambda_state” in the mdp changes from 0 - 20. If I choose to execute the files simultaneously (not consecutively), will I that be okay?
From my understanding of how free energy simulation work, yeah. Each simulation is independent from one another. The reason the init_lambda_state changes is because each simulation is simulating a different lambda point. init_lambda_state=0 simulates the first lambda point init_lambda_state=1 simulates the next and so on. job.sh makes it to where you don’t have to run each simulation by hand. You could run them simultaneously, as to whether your computer can handle all that is a different discussion.
Thank you @Noob. Regarding the mdp settings, I have modified the mdp provided in the tutorial + livecomms paper by @jalemkul. Please note that the molecule is a non-ionic surfactant. Are my settings okay as it is?
; Free energy control stuff
free_energy = yes ; instructs grompp to read the relevant free energy settings
init_lambda_state = 0 ; determines the value of λ, which is actually a vector
delta_lambda = 0 ; allows λ to change as a function of time, in what are called "slow-growth" free energy calculations
calc_lambda_neighbors = 1 ; only immediate neighboring windows
couple-moltype = ML4 ; name of moleculetype to decouple
couple-lambda0 = vdw-q ; vdw-q (meaning both vdW and Coulombic interactions are fully on)
couple-lambda1 = vdw-q ; set none if only vdW
couple-intramol = yes ; determines whether intramol nonbonded interactns be transformed as a function of λ, set to no for small and yes for larger molecules
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; We are not transforming any bonded or restrained interactions
bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated tempering here
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none in this case)
sc-power = 1
sc-sigma = 0.3
nstdhdl = 10 ; interval (in number of time steps) for writing the values of ∂H(λ)/∂λ to the output file
From the vdw-lambdas and coul_lambdas it seems like your intention is to decouple VdW interactions, but keep Coulomb interactions. This is generally not a good idea. The other way around, i.e., to decouple Coulomb interactions while VdW interactions are unchanged is usually better. On the other hand, with the same settings for couple-lambda0 and couple-lambda1 both end-states will be the same, so you are not actually doing any perturbations.
I would recommend either:
couple-lambda0 = none
couple-lambda1 = vdw
and do Coulomb interactions separately afterwards. Or:
couple-lambda0 = vdw
couple-lambda1 = vdwq
and modify your lambda vectors accordingly.
Regarding running all simulations in parallel, yes, that works. They are all independent. However, it might be more efficient to first run a simulation with lambda going from one state to the other (using delta-lambda) to find suitable snapshots to start your 21 simulations from. That can shorten your equilibration times (the time you need to discard in the beginning of each simulation when you analyze the data) for each lambda state and/or reduce the risk of crashes due to large instantaneous state changes.
Thank you! I will settle for option 2.
couple-lambda0 = vdw
couple-lambda1 = vdwq
Could you please advise on how to set the lambda vectors? How about the following?
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
In that case you are only changing the Coulomb interactions. So, the vdw_lambdas line is not important. The exact lambda point distribution can be very system dependent, but a rough guide could be.
couple-lambda0 = vdw
couple-lambda1 = vdwq
coul-lambdas = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
sc-alpha = 0
Then you follow this by a set of simulations with:
couple-lambda0 = none
couple-lambda1 = vdw
vdw-lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
sc-alpha = 0.5
Or you can do both transformations in one simulation, e.g.:
couple-lambda0 = none
couple-lambda1 = vdwq
fep-lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
sc-alpha = 0.5
sc-coul = yes
Which alternative is most efficient may vary from system to system.
Thanks for answers! I will start with the option that apply both transformations in one simulation and I have updated the mdp file for minimization task as follows. Kindly advise on whether it’s good.
integrator = steep
nsteps = 5000
; EM criteria and other stuff
emtol = 100
emstep = 0.01
niter = 20
nbfgscorr = 10
; Output control
nstlog = 1
nstenergy = 1
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 1
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
; Temperature and pressure coupling are off during EM
tcoupl = no
pcoupl = no
; Free energy control stuff
free_energy = yes ; instructs grompp to read the relevant free energy settings
init_lambda_state = 0 ; determines the value of λ, which is actually a vector
delta_lambda = 0 ; allows λ to change as a function of time, in what are called "slow-growth" free energy calculations
calc_lambda_neighbors = 1 ; only immediate neighboring windows
couple-moltype = ML4 ; name of moleculetype to decouple
couple-lambda0 = none ;
couple-lambda1 = vdw-q ; vdw-q (meaning both vdW and Coulombic interactions are fully on)
couple-intramol = yes ; determines whether intramol nonbonded interactns be transformed as a function of λ, set to no for small and yes for larger molecules
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
fep-lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated tempering here
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = yes ;
sc-power = 1
sc-sigma = 0.3
nstdhdl = 10 ; interval (in number of time steps) for writing the values of ∂H(λ)/∂λ to the output file
; No velocities during EM
gen_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Do not constrain the starting configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
For energy-minimization you can set free_energy = no. You are not interested in sampling the free energy differences at that point. You would probably also want to equilibrate the system before you start the free energy calculations.
Okay. Before proceeding to the energy calc, I have pre-equilibrated the system as follows: minimization —> 200 ps NVT —> 200ps NPT —> 1ns NPT prod. free_energy = no .
For the solvation energy calc, using the gro file from 1ns NPT prod, I plan to do 100 ps NVT —> 100 ps NPT —> 5 ns NPT production. this will have the free_energy = yes
That sounds reasonable, yes. You probably don’t need to do 100 ps NVT after the 1ns NPT prod.
Thank you!
I tried starting an NVT simulation for a system with 2636 atoms (1 surfactant and 856 water) in box size 2.98160 nm * 2.98160 nm * 2.9816 nm and got the following error. Is there a way to speed up the simulation using an HPC with the following computing resources 1 node 48 ntasks without experiencing the same problem?
Initializing Domain Decomposition on 48 ranks
Dynamic load balancing: auto
Minimum cell size due to atom displacement: 0.700 nm
Initial maximum distances in bonded interactions:
two-body bonded interactions: 0.419 nm, LJC-14 q, atoms 29 12
multi-body bonded interactions: 0.419 nm, Ryckaert-Bell., atoms 12 29
Minimum cell size due to bonded interactions: 0.461 nm
Maximum distance for 13 constraints, at 120 deg. angles, all-trans: 0.218 nm
Estimated maximum distance required for P-LINCS: 0.218 nm
Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25
Guess for relative PME load: 0.55
Using 0 separate PME ranks, as guessed by mdrun
Optimizing the DD grid for 48 cells with a minimum initial size of 0.875 nm
The maximum allowed number of cells is: X 3 Y 3 Z 3
-------------------------------------------------------
Program: mdrun_mpi, version 2020.2
Source file: src/gromacs/domdec/domdec.cpp (line 2277)
MPI rank: 0 (out of 48)
Fatal error:
There is no domain decomposition for 48 ranks that is compatible with the
given box and a minimum cell size of 0.875 nm
Change the number of ranks or mdrun option -rdd or -dds
Look in the log file for details on the domain decomposition
With such a small system you need fewer tasks. I.e. divide the system into fewer domains. You could run 24 tasks, but I would suggest 12 to be on the safe side. E.g.,
mpirun -np 12 gmx_mpi mdrun
Check the output/log file to make sure that you get 4 OMP threads per rank. Otherwise, you must specify that explicitly:
mpirun -np 12 gmx_mpi mdrun -ntomp 4
But you will probably notice that such a small system does not scale very well to 24 physical cores, 48 threads with hyper threading.
The simulation is somehow slow using 12 MPI processes (4 OpenMP threads per MPI process). If I increase box size, do you think I can solve the problem? If yes, what size should I be looking at? Any other suggestions to possibly speed it ups?
It will not go faster with a larger box (more atoms), but it will scale better with more hardware resources. I don’t see why your system would be slow, unless you are calculating energies (nstcalcenergy) or writing outputs very often (nstlog, nstenergy, nstxout, nstxout-compressed etc).
For a 100ps simulation with 2fs and My current setting is:
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 0
Why do you want to write that often? That will certainly slow things down. Unless you have special reasons to keep lots of data (from equilibration) I’d suggest:
; Output control
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 10000
nstenergy = 10000
nstxout-compressed = 10000
No specifics. I was just using the default from the equilibration steps.
The simulation is working perfectly now. The slow speed was due to my HPC not optimally using the available processors. The GROMACS was compiled mvapich2 not IntelMPI. I have since set the MV2_ENABLE_AFFINITY=0 to ensure that the MPI tasks don’t get locked to single core.
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 2500000 ; 2 ns
nstcomm = 100
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 0
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = PME
rcoulomb = 1.2
; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps = system
tau_t = 1.0
ref_t = 298
; Pressure coupling is on for NPT
Pcoupl = Parrinello-Rahman
tau_p = 1.0
compressibility = 4.5e-05
ref_p = 1.0
; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
couple-moltype = ML4 ; name of moleculetype to decouple
couple-lambda0 = none ; turn off everything
couple-lambda1 = vdw-q ; vdw-q (meaning both vdW and Coulombic interactions are fully on)
couple-intramol = yes
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
fep-lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
; We are not transforming any restrained interactions
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated tempering here
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = yes ; transformation to the Columbic interaction of a molecule
sc-power = 1
sc-sigma = 0.3
nstdhdl = 10
; Do not generate velocities
gen_vel = no
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation = yes
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
Using the production settings shown above, I ran a simulation and got the following result:
Final results in kJ/mol:
point 0 - 1, DG -46.08 +/- 0.03
point 1 - 2, DG -42.00 +/- 0.04
point 2 - 3, DG -36.28 +/- 0.01
point 3 - 4, DG -28.74 +/- 0.02
point 4 - 5, DG -19.19 +/- 0.02
point 5 - 6, DG -7.45 +/- 0.02
point 6 - 7, DG 5.82 +/- 0.05
point 7 - 8, DG 18.67 +/- 0.05
point 8 - 9, DG 29.43 +/- 0.09
point 9 - 10, DG 28.86 +/- 0.33
point 10 - 11, DG 13.63 +/- 0.14
point 11 - 12, DG 9.68 +/- 0.06
point 12 - 13, DG 9.85 +/- 0.06
point 13 - 14, DG 11.07 +/- 0.10
point 14 - 15, DG 12.80 +/- 0.08
point 15 - 16, DG 14.36 +/- 0.06
point 16 - 17, DG 16.09 +/- 0.15
point 17 - 18, DG 16.61 +/- 0.17
point 18 - 19, DG 12.40 +/- 0.12
point 19 - 20, DG -15.55 +/- 0.24total 0 - 20, DG 3.95 +/- 0.64
System was has big -ve values at lower points but not at higher lambda. Any advice how to resolve this? I was expecting to get a negative average value
PS: Before proceeding to the energy calc, I have pre-equilibrated the system as follows: minimization —> 200 ps NVT —> 200ps NPT —> 10ns NPT prod. free_energy = no . For the solvation energy calc, using the gro file from 10ns NPT prod, I plan to do 100 ps NVT —> 100 ps NPT —> 5 ns NPT production (MDP setting as above).
the perl script from Justin’s tutorial was used to create the mdp for the different lambdas.