Water is not settled when em run

User discussions
1

GROMACS version: 2019.4
GROMACS modification: Yes

Hello Gromacs users,
I have tried to pulling + Umbrella sampling with DOPC bilayer system.
Before I pull a molecule, I made the dopc+water system that em and equilibration step is done.
And I just put the molecule in the water region, followed by the energy minimization.
But I encountered the error below:

step 9: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates
Step=    9, Dmax= 8.6e-02 nm, Epot=          inf Fmax= 6.22838e+06, atom= 7
Step=   10, Dmax= 4.3e-02 nm, Epot= -1.37278e+05 Fmax= 2.93889e+05, atom= 40
Step=   11, Dmax= 5.2e-02 nm, Epot= -2.20640e+05 Fmax= 3.01888e+05, atom= 2
Step=   12, Dmax= 6.2e-02 nm, Epot= -2.30162e+05 Fmax= 6.42174e+05, atom= 5
Step=   13, Dmax= 7.4e-02 nm, Epot= -2.49762e+05 Fmax= 9.33565e+04, atom= 37
Step=   14, Dmax= 8.9e-02 nm, Epot= -2.68927e+05 Fmax= 5.14408e+04, atom= 17551
step 15: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates
Step=   15, Dmax= 1.1e-01 nm, Epot=          inf Fmax= 9.33565e+04, atom= 37^MStep=   16, Dmax= 5.3e-02 nm, Epot= -2.76287e+05 Fmax= 3.41605e+04, atom= 17551
step 17: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates
Step=   17, Dmax= 6.4e-02 nm, Epot=          inf Fmax= 5.14408e+04, atom= 17551^MStep=   18, Dmax= 3.2e-02 nm, Epot= -2.80820e+05 Fmax= 3.08559e+04, atom= 17482
Step=   19, Dmax= 3.9e-02 nm, Epot= -2.83108e+05 Fmax= 6.51723e+04, atom= 17482
Step=   20, Dmax= 4.6e-02 nm, Epot= -2.86219e+05 Fmax= 2.73405e+04, atom= 17482
Step=   21, Dmax= 5.5e-02 nm, Epot= -2.87274e+05 Fmax= 7.82276e+04, atom= 17482
step 22: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates

Below is the result of em step(em.gro file) that the bonds in the molecule is dissociated.
(The other water and dopc molecules are omitted in figure.)
Screenshot from 2020-07-24 09-52-35

So I tried to find the source for the problem that raise error.

  1. Make only ligand system
  2. Solvate ligand only system
  3. Energy minimization system

Below is the em-processed ligand only system

Screenshot from 2020-07-24 13-16-23
Screenshot from 2020-07-24 13-16-23884Ă—750 77.1 KB

As you can see, the ligand only system does not have problem, meaning that the ligand topology file does not have problem.

Could you know what I’m missing?
Thank you for your reply in advance :)

2

If you use water model from gromacs’s share directory, it has flexible force field for minimization/equilibration purpose. You can use those with define=-DFLEXIBLE option in mdp file. First, run minimization with define option, and then run without it ( rigid water).

1 Like
3

Thank you :)

4

Flexible water won’t solve a totally broken ligand. Did you happen to get any warnings from grompp about non-matching names? The bonded connectivity of the ligand in the “bad” image looks totally different from that of the correct system, implying that something is out of order in the topology with respect to the coordinates.

1 Like
5

Because I overwrote the result file I uploaded to the other file, I try it again without flexible water model and tell whether there are the grompp non-matching names lines.

But there is no error when I give the flexible water option

6

I tried to energy minimization without flexible water and there is no warning about the non-matching.

7

Check the starting coordinates. Is the ligand in contact with/embedded within the DOPC? Depending on how you prepared the coordinates, it may be possible that there’s a lipid alkyl chain sticking through the ligand, leading to massive distortion.

8

As I prepared the DOPC bilayer and then put the ligand by merely concatenating the dopc .gro file and ligand .gro file, It is correct ,as you said, that there is the ligand in contact with the DOPC molecule( which already suffered from energy minimization step.)
Just like below:

FILENAME="COMPLEX/complex_0.gro"
lig_num=$(tail -n +2 LIGAND/lig_newbox_0.gro | head -1)
dopc_num=$(tail -n +2 DOPC/dopc72.gro | head -1)
total_num=$((lig_num+dopc_num))

echo 'LIGAND-DOPC COMPLEX' > $FILENAME
echo ' '$total_num >> $FILENAME
tail -n +3 LIGAND/lig_newbox_0.gro | head -n -1 >> $FILENAME
tail -n +3 DOPC/dopc72.gro >> $FILENAME

So as I wanted to prevent clashes between ligand and DOPC molecule, I would run the energy minimization.
But during the pulling step, there is massive distortion as you mentioned when I pulled the ligand along the z reaction coordinates from middle inner bilayer to upper direction.
At first, I think it is because of the charges from ligand and DOPC molecules as the polar ester groups in DOPC are interacting with the charged ligand

Screenshot from 2020-07-28 15-02-47
Screenshot from 2020-07-28 15-02-47470Ă—619 114 KB

(The ligand is represented as bond, the DOPC molecule is represented as vdw, and the yellow/blue ball is the N4, P8 atom in a DOPC molecule, respectively.)

What do you think about the reason for the distortion??? Is it from the the ligand stick to the DOPC molecule or the strong charge interaction?

Thank you for your reply :)

9

You need to prepare the coordinates more carefully, with e.g. Packmol.

10

Thank you :)