Freeze groups error in grompp/mdrun

User discussions
1

GROMACS version: 2020.1 for Ubuntu
GROMACS modification: No

Hi,

I am having an issue with the freeze groups functionality. I am simulating two proteins in water. I want to freeze the position of the larger of the two proteins. When I use grompp to set up the production run:

gmx grompp -f md-equilibration/md.mdp -c pr-equilibration/pep_actin_pr.gro -r pr-equilibration/pep_actin_pr.gro -p common/pep.top -o md-equilibration/pep_actin_md.tpr -n system.ndx -maxwarn 2

It produces a tpr file but this is the output of the command:

Back Off! I just backed up md-equilibration/pep_actin_md.tpr to md-equilibration/#pep_actin_md.tpr.3#
free(): invalid next size (fast)
Aborted (core dumped)

That last part is concerning. When I try to use mdrun with this tpr file I end up with a segmentation fault. Here is my mdp file (cannot attach because new user…):

integrator = md
dt = 0.002
tinit = 0 ; ps
nsteps = 150000000 ; 300 ns
nstxout = 10000
nstvout = 10000
nstlog = 1000
nstenergy = 1000
nstxtcout = 1000
nstlist = 10
ns_type = grid
vdwtype = Cut-off
optimize_fft = yes
fourierspacing = 0.125
pme_order = 6
ewald_rtol = 1e-5
tcoupl = v-rescale
tc-grps = System ;Protein_chain_A Protein_chain_B NA SOL
tau_t = 0.1 ;0.1 0.1 0.1
ref_t = 300 ;300 300 300
Pcoupl = berendsen
pcoupltype = isotropic
tau_p = 2
compressibility = 0.45e-4
ref_p = 1.0
gen_vel = yes
gen_temp = 300
gen_seed = 173529
constraints = h-bonds
cutoff-scheme = Verlet
vdw-modifier = force-switch
rlist = 1.2
rvdw = 1.2
rvdw-switch = 1.0
coulombtype = PME
rcoulomb = 1.2
DispCorr = no

; FREEZING ACTIN GROUPS COORDINATES
freezegrps = Protein_chain_A Protein_chain_B NA SOL
freezedim = Y Y Y N N N N N N N N N

I have also tried listing the tc-grps explicitly without using the group System and I get a similar error (or a seg fault). It is also worth mentioning that my energy minimization and pressure equilibration steps without freezing any groups works without error. And I can start a production simulation without error if I remove the freeze groups sections from the above mdp file. That leads me to believe this feature is explicitly the reason for the errors I am receiving. I am unsure why this is creating an issue.

Thank you for any help,

Brian

2

Hi,
I understood that you want to freeze one group ( that is X, Y, and/or Z position will not be updated), thus
freezegrps = the group that is to be frozen
(see Molecular dynamics parameters (.mdp options) — GROMACS 2020.5 documentation)
Best regards
Alessandra

3

Hi Alessandra,

That option is included in my mdp file that I provided in the original post. My issue is that the inclusion of freezegrps causes errors. If I take that option out of the mdp file, everything works. The error messages I am getting do not provide a clear path for debugging. Any help is appreciated.

Thanks,

Brian

4

What is the purpose of -maxwarn 2? This is a very dangerous practice and you may be overriding important warnings.

Note that freezing is extremely unphysical, and without energygrp-excl, you can get spurious contributions to the virial that lead to instability, as the manual notes. I think such exclusions may be incompatible with the Verlet scheme, though. I guess the real question is: why do you (think you) need freezing? Would a strong restraint accomplish the same thing?

5

The maxwarn doesn’t have to be there. My grompp command only produces 2 notes which are inconsequential (one is how much data the production run will produce). So I don’t think this negatively affecting anything. Especially since I can run the simulation without the freezegrps.

I am mostly interested in this because my system is quite large now ~10 nm^3. I expect the smaller of my two proteins to bind to the larger. By freezing the larger protein, I am hoping to reduce both the real and simulation time required to observe how the unfrozen protein behaves in the presence of this larger protein.

I tried adding energy groups and exclusions as below:

; FREEZING ACTIN GROUPS COORDINATES
freezegrps = Protein_chain_A Protein_chain_B NA SOL
freezedim = Y Y Y N N N N N N N N N
energygrps = Protein_chain_B NA SOL
energygrp-excl = Protein_chain_A Protein_chain_A

But this lead to a segmentation fault. Could you explain the strong restraint you mention? *Edited to add: Also tried putting Protein_chain_A in energygrps and that still produced a segmentation fault.

Thanks for the help.

Brian

6

You just need to use define to implement the restraints, using the keyword in the topology that pertains to chain A. Restraints won’t speed up the calculation, but they’re not as crude and unstable as freezing (which is really in place for modeling surfaces and rigid materials).

7

Thank you for the help! This solution is definitely sufficient for what I want to achieve.