I cannot follow the martini tutorial Parametrization of a new small molecule

User discussions
1

GROMACS version: 2024.3
GROMACS modification: No
Dear all, I am trying to parametrize FMN in martini3 forcefield to do a coarse grained simualtion. For that I am following the Parametrization of a new small molecule tutorial from martini oficial site. But I run into this error:

Segmentation fault (core dumped) gmx mdrun -v -deffnm 2-eq >> mdrun.log 2.

I tried both with my ligand and with the example files that are included thanks to the tutorial. I always get the same error. I have searched for the root of the problem and apparently at the step 0 of equilibration the system has achieved somehow a temperature of 700k and a negative pressure of - 60 bar approx, which raises a warning “Warning: pressure scaling more than 1%, mu: 0.489707 0.489707 0.489707”.

I don’t know if this warning is the one causing the segmentation fault, but it certaily seemed suspicious both the temperature and the pressure conditions. The minimization also seems supicious, since only 15 steps are done, and it exits with this comment:

" Step Time
14 14.00000

Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 100 (which may not be possible for your system).
It stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters."

I am a little bit at lost here, specifically since following a tutorial should make things easier, but still I fail to do the very firs step of this tutorial. So any help or comment would be greatly appreciated, Thank you all for your attention,

Best,
Alex

Should I maybe try with an older GROMACS version?

2

Are you running the energy minimization on GPU? If so, see if it runs longer on CPU. How large are the forces that are listed in the energy minimization output if you run with the -v flag?

3

Hi Magnus, thanks a lot for your reply, it is great knowing that there is someone out there hepling.

I have run gromax disabeling the GPU and it still does the same 14 steps, and in both cases the forces listed in the energy minimization output are :

with GPU enabled
Potential Energy = -6.7915320e+04
Maximum force = 1.5379228e+05 on atom 4328
Norm of force = 3.7540321e+03

without GPU
Potential Energy = -6.8119086e+04
Maximum force = 1.5379244e+05 on atom 4328
Norm of force = 3.7540325e+03

After seeing these results I have tried to solvate the system differently to the steps described by the tutorial, first by generating a buffer region with the command :gmx editconf -f ENAP_LigParGen.pdb -o boxed.gro -c -d 1.0 -bt cubic, then solvating the box, gmx solvate -cp boxed.gro -cs spc216.gro -o solvated.gro -p topol.top and finally neutralizing the system. When I run the simulation now, it runs for 154 steps, with the steep minimization. The resulting minimization gives :

Potential Energy = -4.3677281e+04
Maximum force = 9.6763250e+04 on atom 31
Norm of force = 2.0883371e+03

I changed the solvation method because I thought that maybe two or more solvent atoms were very close to the ligand and that was messing up the minimization.

Again thanks a lot for your comments, and help,
Best,
Alex

4

Hi,

If you are using only their example files, I’m a bit surprised it does not work as expected. It might be worth trying an older version of GROMACS (2023.5 or even an unsupported version, such as 2022.6). If it works in an older version of GROMACS it would be very good if you could open an issue (Issues · GROMACS / GROMACS · GitLab). If it does not work in an older version of GROMACS then I think you would have to contact the Martini developers.

5

Hi Magnus,
I uninstalled the version I had of GROMACS, and I installed 2 different versions at different locations the 2023.6 and the 2024.3, paying attention to all the flags in the cmake command, then I tried running the files as in the tutorial, and both the 2024 and 2023 version fail to minimize the system further than the 14th step. I don’t know what I am doing wrong, but will keep trying things, see if I find what is failing, I will post anything related to the solution. Thank you so much for your help !

6

I’d recommend contacting the Martini developers then, i.e., if it’s using only the input from their tutorial.

If you are using your own input it’s a completely different matter of course.

7

Dear all, I finally found the error that was preventing me to continue in the tutorial, I was having trouble with the simulation box size now when I change it from 3.8 nm (tutorial size) to 5.0 nm in the tutorial it works. However when I try with my molecule, it will work for some time and then suddenly interrupt with this error:

step 77000, remaining wall clock time: 102 s [LAPTOP-ITIGCFT1:256174] *** Process received signal ***
[LAPTOP-ITIGCFT1:256174] Signal: Segmentation fault (11)
[LAPTOP-ITIGCFT1:256174] Signal code: Address not mapped (1)
[LAPTOP-ITIGCFT1:256174] Failing at address: 0x56352a9ba508

I am now confused on what to do next, I guess it is probably a problem with the itp or .pdb file (then transformed to gromacs, according to the tutorial) since this problem is only happening with my files. These files were generated with LigParGene, for FMN, so the should be correct … I will add here the comands I am runing just to make sure that I am not mistaken there.

#!/bin/bash

shopt -s expand_aliases
alias gmx=“gmx_mpi”

#------------------------------------------------------------------------------------------------#

USE:

./prepare_1mol_AA_system.sh solute.gro box_solvent.gro NAME

DESCRIPTION:

Generates a box with a molecule solvated in a solvent of choice. The script needs 4 inputs:

1. SOLUTE.gro,

2. BOX_SOLVENT.gro

3. solvent name

4. solvent number of atoms

BEFORE:

Make sure your topology file is correctly set up.

COMMENTS:

-

#------------------------------------------------------------------------------------------------#

Load GROMACS 2016.5 or later

source /usr/local/gromacs-2024/bin/GMXRC

Check if the files have been passed to the script

solute=“$1”
solvent_box=“$2”
solvent_name=“$3”
solvent_atoms="4" size_1={#solute}
size_2=${#solvent_box}
if [ $size_1 = 0 ] || [$size_2 = 0 ] ; then
echo “”
echo “Missing file of solute or solvent box. Check that. e.g.:”
echo “”
echo “(WATER) ./prepare_1mol_AA_system.sh benz.gro spc216.gro SOL 3”
exit
fi

Insert molecule into a box

gmx insert-molecules -ci $solute -nmol 1 -box 5.0 5.0 5.0 -rot xyz -seed 0 -o inserted.gro

Center the molecule in the box explicitly

gmx editconf -f inserted.gro -o centered.gro -c

Solvate the system with the solvent box

gmx solvate -cp centered.gro -cs $solvent_box -o initial.gro

cp system_EMPTY.top system.top
solvent_lines=$(grep solvent_name initial.gro | wc -l) solvent_molecules=(expr $solvent_lines / $solvent_atoms )
echo “$solvent_name $solvent_molecules” >> system.top

gmx grompp -p system.top -c initial.gro -f em.mdp -o 1-min.tpr -po 1-min.mdp
mpirun -np 2 gmx_mpi mdrun -v -deffnm 1-min -ntomp 2 -gpu_id 0 -pin on >> mdrun.log 2>&1

gmx grompp -p system.top -c 1-min.gro -f eq1.mdp -o 2-eq1.tpr -po 2-eq1.mdp
gmx mdrun -v -deffnm 2-eq1 >> mdrun.log 2>&1
gmx grompp -p system.top -c 2-eq1.gro -f eq2.mdp -o 2-eq2.tpr -po 2-eq2.mdp
gmx mdrun -v -deffnm 2-eq2 >> mdrun.log 2>&1

gmx grompp -p system.top -c 2-eq.gro -f run.mdp -o 3-run.tpr -po 3-run.mdp
gmx mdrun -v -deffnm 3-run >> mdrun.log 2>&1

then from the terminal : bash prepare_1mol_AA_system.sh ENAP_LigParGen.pdb spc216.gro SOL 3

As always, thanks for your attention

8

Which step of gmx mdrun is crashing? Is it one of the equilibration runs 2-eq1 or 2-eq2)? Does it work if you try a lower timestep dt in the mdp file?

The equilibration stages in the tutorial might be tailored to the example system. It is not certain that the same settings (except the force field specific parameters) can be used for all systems. You might need to lower the integration time step (as I suggested above) and/or equilibrate for a lot longer (depending on in which simulation stage the crash is occuring). You might also have to change the thermostat and barostat (v-rescale and c-rescale with tau-t = 1.0 and tau-p in the range of 5.0-10.0 are usually decent choices).

9

The simulation crashes during the first equilibration, I have now tried with the decreased time step and increased tau, but it still crashes:

"step 63000, remaining wall clock time: 231 s
step 63100, remaining wall clock time: 231 s
step 63200, remaining wall clock time: 231 s [LAPTOP-ITIGCFT1:264920] *** Process received signal ***
[LAPTOP-ITIGCFT1:264920] Signal: Segmentation fault (11)
[LAPTOP-ITIGCFT1:264920] Signal code: Address not mapped (1)
[LAPTOP-ITIGCFT1:264920] Failing at address: 0x5589ac4a3ef8 "

I will add my eq1.mdp file so it is more clear what I am doing:

;
; **** GMX general run parameters for an atomistic simulation ****
;
; Run parameters are according to [R. Alessandri, P.C.T. Souza, et al. JCTC, 2019, 15, 5448].
;
; You may want to double-check whether:
;
; coulombtype
; Tcoupl
; Pcoupl
;
; suit your needs.
; This is the equilibration mdp.
;
; RA
;

integrator = md ;
dt = 0.0005 ; ps
nsteps = 500000 ; total time: 0.250 ns
tinit = 0 ; initial time, ps
nstcomm = 100 ; freq. for cm-motion removal
ld_seed = -1

; Bond constraints
constraints = h-bonds ; constrain hidrogen bonds
constraint_algorithm = lincs ; default

; X/V/F/E outputs
nstxout = 500000 ; pos out — 1000 ps
nstvout = 500000 ; vel out — 1000 ps
nstfout = 0 ; force out — no
nstlog = 10000 ; energies to log (20 ps)
nstenergy = 10000 ; energies to energy file
nstxtcout = 10000 ; xtc, 10 ps
xtc_precision = 1000

; Neighbour list
ns_type = grid ; neighlist type
nstlist = 5 ; Freq. to update neighbour list
rlist = 0.8 ; nm (cutoff for short-range NL)

; Coulomb interactions
;coulombtype = Reaction-field ;
;epsilon_rf = 80 ; water
;epsilon_rf = 4.8 ; CHCl3 (CRC Handbook)
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
rcoulomb = 1.4 ; nm (direct space sum cut-off)
optimize_fft = yes ; optimal FFT plan for the grid

; van der Waals interactions
vdwtype = Cut-off ; Van der Waals interactions
rvdw = 1.4 ; nm (LJ cut-off)
DispCorr = EnerPres ; use dispersion correction

; Temperature coupling
;Tcoupl = Berendsen ; This is ignored with sd integrator
;tc-grps = System ;
;tau_t = 2 ; ps, recommended value for sd
;ref_t = 298.15 ; K

tcoupl = V-rescale ; Stochastic velocity rescaling thermostat
tc-grps = System ; Apply temperature coupling to the entire system
tau_t = 1 ; Time constant for temperature coupling
ref_t = 298.15 ; Reference temperature (K)

; Energy monitoring
energygrps = System

; Pressure coupling
Pcoupl = no

; Generate velocites in the beginning
continuation = no ; continue from npt equilibration
gen_vel = yes ; continue from npt equilibration
gen_temp = 298.0
gen_seed = -1 ; -1 = the seed is calculated from the process ID number

10

Thank you alexativic. I was having the same problem.