Issue with Copper Ion Restraints: Two Ions Moving Despite Position Restraints

User discussions
1

GROMACS version:2022
Hello GROMACS Community,

I’m encountering an issue with my simulation involving a protein containing four copper ions. I’ve applied position restraints to ensure the ions remain fixed, but during the simulation, two of them are exhibiting unexpected movement.
I’ve checked the simulation, and despite the high force constants, two of the copper ions are not staying in place. The equilibration time is set to 1 ns, and I’ve visually inspected the trajectory.

If anyone has insights into why this might be happening or suggestions for troubleshooting, I would greatly appreciate your input.

Thank you!

; Position restraints for copper ions in GROMACS

[ position_restraints ]
; i funct fcx fcy fcz
3585 1 100000 100000 100000
3586 1 100000 100000 100000
3587 1 100000 100000 100000
3588 1 100000 100000 100000

; LINES STARTING WITH ‘;’ ARE COMMENTS
title = Minimization ; Title of run
define = -DPOSRES-strong ; position restrain the protein

; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 10.0 kJ/mol
emstep = 0.0001 ; Energy step size
nsteps = 500000 ; Maximum number of (minimization) steps to perform
constraints = none

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet
ns_type = grid ; Method to determine neighbor list (simple, grid)
rlist = 1.0 ; Cut-off for making neighbor list (short range forces)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; long range electrostatic cut-off
rvdw = 1.0 ; long range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions

2

There are two quick checks to see if Gromacs correctly reads position restraints, one is to run gmx grompp without the -r option (it should complain about the lack of reference coordinates), the other is to look for the respective energy term in the .log file produced by the run. If either check fails, that means you somehow didn’t manage to set up the restraints correctly.

You only share an .mdp for minimization, do you set the same defines in your equilibration .mdp too?

3

Thanks for your reply, I haven’t used -r. Also, I don’t understand what exactly do you mean by fail in the log files? Is it mentioned in the file?