How can I keep Cu fixed in protein during the MD simulation?

GROMACS version:2020
GROMACS modification: No

Hello all.I performed an MD simulation on a laccase containing 4 Cu ions in water solvent without any substrate.
simulation run without any error, but when i used vmd for visualizing the last trajectory, one of the ions that was complexed by His and cys, has changed its place and is not inside the protein at the end of simulation time. Is this simulation reliable or not? how can i keep cu fixed in protein during the MD simulation?
I used charmm-gui to generate input files for simulation.

thank you.


You can apply harmonic restraint to restrain the distance between the Cu ions and the protein atoms to which they are coordinated

For harmonic restraints you can use [bonds] type 6 in the topology file

See for more details “Distance restraint” in

and topology syntax in

Best regards

Hi @alevilla
I have same problem but with zn ion between 3 cys
I tried add these lines at the end of [bonds]
4954 2989 6
4954 3036 6
4954 3534 6
First atom number is zn and second is cysteines HG

But in the grompp I face the error:
Fatal error:
Invalid pairs type 6

Where is my mistake and how can I solve this?

From the error it looks that the values are located in the [pairs] directive, and not in the [bonds] directive. Can you check that the lines are added at the end of the [bonds] directive?
You also need to specify the values of distance b0(nm) and the force [ kJ(mol)-1(nm)-2 ], you want. As a distance you will use the distance between the Zn and S.


I also have the same problem.
I’ve got b0 and the force but I don’t know in which file I need to put them.
can you please help me with that?

do you mean that you want to use harmonic restraints between 2 atoms?
Then you can add them in the topology file that describe your molecule (note topology files are hierarchical). See links above posts for topology syntax

Under [moleculetype] directive, you find a directive calls [bonds] there you can add lines (i.e. at the end of the directive) corresponding to those harmonic restraints you want to add.

For example you want an harminic restraints between atoms 503 and 900 of a protein (1000 atoms) with a distance of {XX} nm and a force of {YYYY} [ kJ(mol)-1(nm)-2 ]

[ bonds ]
; aj ai func b0 force0

503 900 6 XX YYYY

I hope it helps

Thank you so much for your help.
I had 6 Cu atoms at my pdb file with same number, I changed their number and then pdb2gmx made 6 different ion.itp file for me, and that worked, I hope that was the right thing to do.

Note, that you need to define all the atoms involved in the restraints under the same [moleculetype] directive.
For example if you want the harmonic restraints between a Cu atom and a protein atom, both atoms have to be defined under [atoms] of the same [moleculetype] directive and the restraints under bonds.
To get a topology file that include all atoms under the same [moleculetype] directive you can use the option -merge of gmx pdb2gmx.