Error when running umbralla sampling

GROMACS version:2021.2
GROMACS modification: Yes
Here post your question
Dear GROMACS users,
I’m try to run umbrella sampling to study the the permeation of nitric oxide(NO) in DPPC and cholesterol bilayer system. Prior to insertion of the two NO molecules, the system was equilibrated for about 50 ns. Then, the two solute molecules were inserted at different depths after which a crude minimization was run to remove bad contacts. However, upon proceeding the system to umbrella sampling run, I keep getting this error below with other details:

"Step 267, time 0.534 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000075, max 0.003388 (between atoms 8620 and 8622)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
8620 8621 57.2 0.1111 0.1111 0.1111
12962 12963 70.4 0.1111 0.1111 0.1111

Step 268, time 0.536 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000041, max 0.001670 (between atoms 12962 and 12964)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
8620 8621 47.5 0.1111 0.1111 0.1111
12962 12963 72.6 0.1111 0.1111 0.1111

Step 269, time 0.538 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000074, max 0.002686 (between atoms 12962 and 12964)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1063 1064 42.5 0.1111 0.1111 0.1111
8620 8621 35.9 0.1111 0.1111 0.1111
12962 12963 55.0 0.1111 0.1111 0.1111

Step 270, time 0.54 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.144456, max 6.443329 (between atoms 12959 and 12961)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1063 1064 32.9 0.1111 0.1111 0.1111
8620 8621 39.9 0.1111 0.1112 0.1111
12962 12963 90.0 0.1111 0.3651 0.1111
12959 12961 90.0 0.1111 0.8270 0.1111
Wrote pdb files with previous and current coordinates

Step 271, time 0.542 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.767628, max 24.837708 (between atoms 1063 and 1066)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1063 1064 65.5 0.1111 0.1048 0.1111
1063 1066 90.0 0.1111 2.8706 0.1111
1086 1087 90.0 0.1111 2.7116 0.1111
8620 8621 90.0 0.1112 0.1585 0.1111
8623 8624 90.0 0.1111 0.2467 0.1111
12962 12963 90.0 0.3651 1.4864 0.1111
12962 12964 45.9 0.1090 0.1164 0.1111
12959 12960 82.5 0.1132 0.1180 0.1111
12959 12961 90.0 0.8270 0.1613 0.1111
Wrote pdb files with previous and current coordinates

Step 272, time 0.544 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 15.565684, max 554.185791 (between atoms 8629 and 8631)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1063 1064 90.0 0.1050 0.3954 0.1111
1063 1065 90.0 0.1103 0.3808 0.1111
1063 1066 90.0 2.8745 0.1174 0.1111
1086 1087 90.0 2.7134 0.1122 0.1111
8614 8615 90.0 0.1113 0.5110 0.1111
8614 8616 90.0 0.1113 0.5223 0.1111
8620 8621 90.0 0.1586 5.6136 0.1111
8620 8622 90.0 0.1126 0.3078 0.1111
12953 12955 54.2 0.1113 0.1111 0.1111
8623 8624 90.0 0.2468 53.7077 0.1111
8623 8625 90.0 0.1115 2.1322 0.1111
12962 12963 90.0 1.4886 0.3751 0.1111
12962 12964 90.0 0.1165 0.1786 0.1111
12959 12960 90.0 0.1181 0.1192 0.1111
12959 12961 90.0 0.1616 0.8942 0.1111
8629 8631 90.0 0.1113 61.6811 0.1111
Wrote pdb files with previous and current coordinates


Program: gmx mdrun, version 2021.2
Source file: src/gromacs/ewald/pme_redistribute.cpp (line 305)
MPI rank: 1 (out of 4)

Fatal error:
1 particles communicated to PME rank 1 are more than 2/3 times the cut-off out
of the domain decomposition cell of their charge group in dimension x.
This usually means that your system is not well equilibrated."

I’ve tried to inspect my system for possible source of the error but it seems my system look fine. I would be glad if anyone could help me in troubleshooting this error. The snap short of the system of inserting the solutes is attached. I have also included the parameters used for minimization and umbrella sampling:
umbrella sampling parameters:
#!/bin/bash
#SBATCH --nodes=1 # number of nodes
#SBATCH --gres=gpu:2 # request 2 GPUs per node (Graham)
#SBATCH --ntasks-per-node=4 # request 4 MPI tasks per node
#SBATCH --cpus-per-task=8 # 8 OpenMP threads per MPI process
#SBATCH --mem=4G # Request all available memory in the node
#SBATCH --time=48:00:00 # time limit (D-HH:MM:ss)

mpirun or srun also work

module purge
module load StdEnv/2020 gcc/9.3.0 cuda/11.0 openmpi/4.0.3 gromacs/2021.2
export OMP_NUM_THREADS="${SLURM_CPUS_PER_TASK:-1}"

echo ";title = Umbrella sampling
integrator = md
dt = 0.002
nsteps = 15000000
nstxout = 50000
nstvout = 50000
nstfout = 50000
nstcalcenergy = 100
nstenergy = 1000
nstlog = 1000
;
cutoff-scheme = Verlet
nstlist = 20
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
coulombtype = pme
rcoulomb = 1.2
;
tcoupl = Nose-Hoover
tc_grps = MEMB SOLV
tau_t = 1.0 1.0
ref_t = 323 323
;
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
tau_p = 5.0
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0
;
constraints = h-bonds
constraint_algorithm = LINCS
continuation = no
;
nstcomm = 100
comm_mode = linear
comm_grps = MEMB SOLV
;
refcoord_scaling = com
; Pull code
pull = yes
pull_ncoords = 2 ;
pull_ngroups = 3 ;
pull_group1_name = MEMB
pull_group2_name = r_3784
pull_group3_name = r_3785
pull_coord1_type = umbrella ;
pull_coord1_geometry = direction ;
pull_coord2_geometry = direction ;
pull_coord1_groups = 1 2
pull-coord2-groups = 1 3
pull-group1-pbcatom = 7594
pull-pbc-ref-prev-step-com = Yes
pull_coord1_dim = N N Y
pull_coord2_dim = N N Y
pull-coord1-vec = 0.0 0.0 3.76
pull-coord1-rate = 0.0 ;
pull-coord1-k = 500 ; kJ mol^-1 nm^-2
pull-coord2-vec = 0.0 0.0 1.46
pull-coord2-rate = 0.0 ;
pull-coord2-k = 500 ; kJ mol^-1 nm^-2
pull_coord1_start = yes ;
pull_coord2_start = yes ;

echo "Starting run at: date"
gmx_mpi mdrun -v -deffnm $SCRATCH/20%/Run1/1RL

minimization parameter:
integrator = steep
emtol = 100
nsteps = 5000
nstlist = 10
cutoff-scheme = Verlet
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
coulombtype = pme
rcoulomb = 1.2
;
constraints = h-bonds
constraint_algorithm = LINCS"

Thanks
gro

2021-06-23T04:00:00Z

Hi,
Just a suggestion to avoid thw issue, in place to insert the NO molecule directly in the bilayer you can pull in the molecule in the membrane using the pull code.
kind regards
Alessandra

Thanks Avilla for your response. I would consider that. Could you please share with me a sample of such pull code?