Firstly, just to check that you are comparing the right values. With a negative average value you mean a negative DG going from point 0-20, right? I.e. going from a state with inter- and intramolecular interactions turned off (an eliminated state) to full interactions with itself and its surroundings.
This does not correspond to the solvation free energy of the molecule, which would not involve eliminating the molecule completely (decoupling its intramolecular interactions). If you want to get the solvation free energy you will need to do a similar decoupling of the molecule in vacuum. I usually just put the molecule of interest in a box of the same dimensions and turn off pressure scaling. You will need some NVT equilibration, of course. The simulations will be quick, but will not scale well. So, you can run as many in parallel as you have cores (or virtual threads) on your computer.
By the way, also here I’d recommend increasing the output writing interval significantly. I don’t think you will need 25,000 uncompressed frames with coordinates and velocities. I’d recommend similar settings as I wrote above. Even 250 compressed coordinate frames per lambda point is a lot.
Yes. My initial understanding was that the computed DG from the previous simulation translates to free energy of solvation and a -ve DG would indicate a thermodynamically stable system.
The goal is to set up the simulation to enable the calculation of the solvation free energy.
If I understand correctly, I should re-run same simulation but remove the water molecules inside the box. I plan to do the following:
step 1: Using the md.gro (containing H20) from equilibration step, remove the water molecules. do 1ns NVT free_energy = no .
Use same setting as from solvation energy calc free_energy = yes containing H20,
step 2: Extract the gro file from step 1 (1ns NVT),
step 3: do 100 ps NVT —> 100 ps NPT —> 5 ns NPT production
step 4: DG_solv = DG (with solvent) - DG (without solvent)`
For reference, If your molecule is small enough (and without significant intramolecular interactions) you can run the original simulations with couple-intramol = no. Then you don’t have to run the vacuum simulations - the output you get from the decoupling in solvent will correspond to the solvation free energy.
But in this case, as you are eliminating the molecule itself as well, you need to calculate the free energy of eliminating it in vacuum.
Is this version correct? I changed the NPT to NVT as suggested
step 1: Using the md.gro (containing H20) from equilibration step, remove the water molecules. do 1ns NVT free_energy = no .
Use same setting as from solvation energy calc free_energy = yes containing H20,
step 2: Extract the gro file from step 1 (1ns NVT),
step 3: do 100 ps NVT —> 5 ns NVT production
step 4: DG_solv = DG (with solvent) - DG (without solvent)`
My molecule (non-ionic) has 68 atoms with 2 oxygens and the rest being C & H. Do you think that’s small enough? Are you saying the DG obtained from couple-intramol = no is basically the DG_solv ?
The last line 'But in this case…" is not clear. Kindly help clarify
I’m sorry I was unclear. What I meant with “But in this case …” was when running with couple-intramol = yes, i.e. eliminating the molecule. With couple-intramol = no the obtained DG should correspond to DG_solv, indeed, with the usual caveats for limited sampling.
Thank you very much. I will implement what we discussed and drop questions if any. Based on @hess comment that Large, flexible molecules tend to fold onto them selves and get stuck in one conformation with couple-intramol=no. It’s better to avoid that. I will run both conditions (with couple-intramol = yes and couple-intramol = no) and compare the results since I have some other molecules that I intend to simulate.
For couple-intramol = no, I intend to use same settings as yes except that I set intramol = no as shown below:
; 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, in this case only vdW
couple-lambda1 = vdw-q ; vdw-q (meaning both vdW and Coulombic interactions are fully on)
couple-intramol = no
; 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 ; transformation to the Columbic interaction of a molecule
sc-power = 1
sc-sigma = 0.3
nstdhdl = 10
Fatal error:
There are 2 perturbed, excluded non-bonded pair interactions beyond the
pair-list cut-off, which is not supported. This can happen because the system
is unstable or because intra-molecular interactions at long distances are
excluded. If the latter is the case, you can try to increase nstlist or rlist
to avoid this.The error is likely triggered by the use of couple-intramol=no
and the maximal distance in the decoupled molecule exceeding rlist.
What you can do to circumvent that is to set
verlet-buffer-tolerance = -1
rlist = 1.6
The actual rlist value you need might be different (hopefully not larger). Higher rlist will make the simulation slower and the maximum depends on your box size.
the above settings results in the following error:
ERROR 1 [file topol.top, line 701]:
ERROR: The cut-off length is longer than half the shortest box vector or
longer than the smallest box diagonal element. Increase the box size or
decrease rlist.
Setting gen_seed to -86654465
Velocities were taken from a Maxwell distribution at 300 K
There were 2 notes
Program: gmx grompp, version 2020.2
Source file: src/gromacs/gmxpreprocess/grompp.cpp (line 1928)
Yes, as I hinted at. You might need to find a suitable value that is high enough to work for your decoupled molecule and small enough to work with your simulation box. Sometimes that is not trivial, not possible or not worth it, especially since large molecules are expected to be more efficiently sampled with intramolecular decoupling.
But in your case you could see if rlist = 1.45 would work, but I suspect it will be too small for your molecule.
Yes, you are right. How did you come about the list of 1.45? Is it 2.98160 / 2?
mdrun_mpi -deffnm nvt0
Reading file nvt0.tpr, VERSION 2020.2 (single precision)
Can not increase nstlist because verlet-buffer-tolerance is not set or used
Using 12 MPI processes
Using 4 OpenMP threads per MPI process
starting mdrun ‘ML4 in water’
50000 steps, 100.0 ps.
Not all bonded interactions have been properly assigned to the domain decomposition cells
A list of missing interactions:
LJC Pairs NB of 1887 missing 1
Molecule type ‘ML4’
the first 10 missing interactions, except for exclusions:
LJC Pairs NB atoms 25 66 global 25 66
Program: mdrun_mpi, version 2020.2
Source file: src/gromacs/domdec/domdec_topology.cpp (line 421)
MPI rank: 0 (out of 12)
Fatal error:
1 of the 2429 bonded interactions could not be calculated because some atoms
involved moved further apart than the multi-body cut-off distance (0.59632 nm)
or the two-body cut-off distance (1.45 nm), see option -rdd, for pairs and
tabulated bonds also see option -ddcheck
Yes, I just based it on 2.98160 / 2 and with some margin for box size changes due to pressure coupling. When the molecule is fully decoupled the box is expected to shrink.
I’m a bit surprised that you get the error from the bonded interactions. I wonder if that has to do with your domain decomposition. Would you get the (similar, but not identical) error I reported above if you use 1 or 2 MPI ranks instead?
It is difficult to know if the result is OK. One thing you can do is change the beginning (-b) and end (-e) times when you do the bar analysis. A higher -b means that you discard more as equilibration. See how much difference you get when you discard the first 10%, 25% vs 50% of the simulation. Also, see how much your results change when you finish (-e) the analysis after 50% and 75% of the total run time. You could also analyse time windows of, e.g. 1 ns, but with increasing starting times. If the differences have a clear trend with increasing amount of data and/or longer equilibration times, you might want to sample longer and/or discard more as equilibration.
