So you are not using lambda-coupling free-energy functionality?
I’m using the CHARMM-GUI .mdp script for minimization with free-energy functionality and soft-core minimization.
The script is as follows:
define = -DFLEXIBLE
integrator = steep
emtol = 700.0
tinit = 0.0
nsteps = 5000
nstlog = 100
nstenergy = 100
nstxout-compressed = 1000
compressed-x-precision = 100
cutoff-scheme = Verlet
nstlist = 20
rlist = 1.2
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
epsilon_r = 2.5
coulombtype = pme
rcoulomb = 1.2
vdw_type = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 0.9
rvdw = 1.2
; GENERATE VELOCITIES FOR STARTUP RUN:
gen_vel = yes
gen_temp = 310
gen_seed = 0154403735
;soft-core-minimization so that single precision GROMACS works here
; Free energy parameters
free-energy = yes
init-lambda = 0.01
sc-alpha = 4
sc-power = 2
sc-coul = yes
nstdhdl = 0
couple-moltype = system
; we are changing both the vdw and the charge. In the initial state, both are on
couple-lambda0 = vdw-q
; in the final state, both are off.
couple-lambda1 = none
couple-intramol = yes
refcoord_scaling = all
I didn’t know about this protocol from charmm-gui.
One explanation might be that an atom get displaced too much in a step. At which step do you get the error?
What is it you want to achieve with that “soft-core-minimization” scheme? You are scaling down the interactions marginally, by 1%.
This error happened after 41 steps.
"
Step 35, time 0.035 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000205, max 0.007230 (between atoms 1951943 and 1951942)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1951943 1951942 30.3 0.1400 0.1410 0.1400
Wrote pdb files with previous and current coordinates
Step= 35, Dmax= 3.6e-01 nm, Epot= inf Fmax= 6.63678e+08, atom= 64485
Step 36, time 0.036 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000389, max 0.052175 (between atoms 331206 and 331204)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
331205 331204 70.7 0.1400 0.1427 0.1400
331206 331204 79.6 0.1400 0.1327 0.1400
Step= 36, Dmax= 1.8e-01 nm, Epot= 7.09476e+08 Fmax= 1.97399e+09, atom= 2185505
Step 37, time 0.037 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.003645, max 0.485409 (between atoms 324357 and 324355)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
324356 324355 90.0 0.1400 0.1680 0.1400
324357 324355 90.0 0.1402 0.2080 0.1400
Step= 37, Dmax= 2.1e-01 nm, Epot= 6.31489e+08 Fmax= 1.62423e+08, atom= 38075
Step= 38, Dmax= 2.6e-01 nm, Epot= 5.41662e+08 Fmax= 1.16308e+07, atom= 134329
Step= 39, Dmax= 3.1e-01 nm, Epot= 2.25710e+08 Fmax= 6.04571e+06, atom= 45487
Step= 40, Dmax= 3.7e-01 nm, Epot= 1.90998e+08 Fmax= 1.06678e+09, atom= 45486
Step= 41, Dmax= 4.4e-01 nm, Epot= 1.71392e+08 Fmax= 8.63684e+08, atom= 47503
Program: gmx mdrun, version 2021.5-EasyBuild-4.7.1
Source file: src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp (line 854)
MPI rank: 25 (out of 64)
Fatal error:
There are 2 perturbed non-bonded pair interactions beyond the pair-list cutoff
of 1.26 nm, 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.
"
Without that soft-core minimization, I got infinite forces error and was suggested to use soft-core potentials with the free energy code to avoid infinite forces.
"
Step= 14, Dmax= 1.2e-06 nm, Epot= 3.40388e+20 Fmax= inf, atom= 1342438
Energy minimization has stopped because the force on at least one atom is not
finite. This usually means atoms are overlapping. Modify the input
coordinates to remove atom overlap or use soft-core potentials with the free
energy code to avoid infinite forces.
You could also be lucky that switching to double precision is sufficient to
obtain finite forces.
writing lowest energy coordinates.
Steepest Descents converged to machine precision in 15 steps,
but did not reach the requested Fmax < 700.
Potential Energy = 3.4038730e+20
Maximum force = inf on atom 1342438
Norm of force = inf
"
So, it’s quite apparent that you have atoms overlapping. I’ve never tried using the free-energy code with soft-core potentials to solve that during energy-minimization. I have used similar methods (according to the second suggestion below) in order to grow in small molecules in a system during equilibration. If you want to try it further I would suggest trying:
init-lambda = 0.9
So that you start in a nearly decoupled state. However, even if you finish that energy minimization you should probably run another, without free-energy perturbations, as follow-up. But the risk is that that will also crash. What you could try is:
init-lambda = 0.9
delta-lambda = -0.00018
I am not sure if that works properly together with energy minimization. But if it does, you would turn on the interactions over time (or rather over the minimization steps). Keep in mind that growing in large molecules like this can results in unphysical configurations and I would not be surprised if the simulations would still crash, even if you might manage to finish the energy minimization step.
Thank you very much for your suggestion. I tried as you suggested but it still crashed and I got the same error:
“Fatal error:
There are 1 perturbed non-bonded pair interactions beyond the pair-list cutoff
of 2.2 nm, 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.”
So you mean that it crashed during the equilibration even after energy minimization with my last suggested setting? I.e., when growing in the ligand using delta-lambda during the energy minimization? Then you could try using a similar free-energy “growing-in” phase during your equilibration.
But I think it would be best for you to go back and check how you generate your initial configuration. Apparently you’ve got atoms overlapping in your system. Fix that and everything will hopefully work better.
Hi @MagnusL, I couldn’t get past the minimization step due to the same error even after trying your last suggested editing.
I generated my initial configuration with CHARMM-GUI. I’ve tried to start over from the beginning, generating another initial configuration, but still got the same error.
I checked the coordinates of atom with infinite force (atom 1342438) with other atoms around it. Their coordinates are close but not exactly overlapped.
ATOM ***** W PW 37815 27.507 309.791 18.107 1.00 0.00 WAT
ATOM ***** WP PW 37815 27.207 309.391 19.407 1.00 0.00 WAT
ATOM ***** WM PW 37815 27.107 310.391 16.907 1.00 0.00 WAT
(***** is because if the atom number > 99999 it is shown like that in the .pdb file of the initial structure)
I tried to change their coordinates a bit then run minimization again without applying the free energy functionality, the infinite force happened on another atom. I repeat the process for that atom and the infinite force happened again on another atom… I checked their coordinates and they’re also not exactly overlapped.
Could you please let me know a better way to check and fix the atom overlapping?
Thank you.
Since moving the atoms a little bit helped it seems like there was some overlap, but I agree that the ones you listed should not cause any problem, but it’s difficult to know which one is atom 1342438. How large is the periodic box of the system? Is the overlap across the periodic boundary?
Hi @MagnusL,
Atom 1342438 is this one:
ATOM ***** W PW 37815 27.507 309.791 18.107 1.00 0.00 WAT
The simulation box size is 750.477 x 750.477 x 189.191 Å, and the overlap is not across the periodic boundary.
That’'s a very large simulation box. Even if you meant dimensions in A it’s large, but that’s beside the point.
Everything indicates a problematic input configuration, but if the atoms with infinite (or extremely high) force are not involved in overlaps and do not have extremely extended bonds I don’t have any more suggestions at the moment.
I’m sorry it’s in Å. I’ve edited my reply above with the correct unit now.
Thank you very much for your suggestions so far. If you think of any other reasons or solutions, please let me know. I appreciate your help and support.
Hello!
You mentioned that additional, solvent-free simulation need to be made when couple-intramol=yes, is this referring to decoupling simulation in a vacuum?
I am calculating the ligand binding free energy of the protein-ligand system. Do I also need to conduct decoupling simulation of the ligand of the complex system in the vacuum?
Thank you!
If you use couple-intramol=yes, then yes.