It was asked to perform an energy minimization of a PDB file (protein/ligand/water). Unlike online tutorials, I have to maintain the water molecules and just energy to minimize the atoms in the crystal structure. The idea is to obtain a relaxed crystal structure to perform some analyses before MD simulation.
The topologies for the protein, ligand, and water molecules (water molecules were recognized by the forcefield CHARMM36) were created, but I’m not sure about the charge of the protein/ligand. In the tutorials, ions are added to neutralize the system and these ions occupy the position of water molecules added by the software.
However, since I shouldn’t add water molecules (I just have to work with crystal water molecules), there isn’t the possibility of including the ions, so the system is charged, and the Gromacs showed I couldn’t use the PME as coulombtype during energy minimization (it must be cutoff).
Since I have never done something like this, I want to know if it is ok or if there are other options. Could someone help me, please?
Inserting the ions is easy using gmx insert-molecules. I imagine there is lots of vacuum in the system if you only include crystal water molecules, so it should be easy finding places to insert the ions.
But more importantly, do you believe the system to be realistic without bulk water?
I create the topology of protein, ligand, and water as usual;
Then I create the box (using gmx editconf -f complex_final.gro -o complex_box.gro -c -d 1.2 -bt dodecahedron);
Create the PDB file with only 1 atom (for example Na or Cl, depending on the net charges);
Run the gmx insert-molecules -f complex_box.gro -ci atom.pdb -nmol "number of atom repetition" -o output.gro to include the ions randomly;
And finally, run the EM.
Is that right?
One more question about the number of ion repetitions. The protein has a net charge of -1. The ligand is -2 and there’s a NADPH that is -3. Thus, if I understood it, the -nmol must be 6 (using Na.pdb), is it?
Well, it was told me the idea of the project is exactly obtain this answer for the problem they are studying.
Yes, your process and ions count look correct. You might want to create a bounding box around the ion as well (e.g., using gmx editconf -d 0.1 -c, but it should not be necessary, I think.
One more thing you will need to do is to create settings for restraining the water. I would suggest copying the water model itp from the force field folder to where your topol.top file is located, remove the path (if any) pointing to the directory of the water model in the topology file, and add a [ position_restraints ] section at the end of the water model itp file, e.g. (for a 3-site water model):
then run gmx editconf -f output.gro -o final_ion.gro -c -d 0.1 -bt dodecahedron
and run EM
Is it, or do I understand it wrong?
At the end of the gmx pdb2gmx, it add the following lines in the topol.top (I selected tip3p model):
; Include water topology
#include "./charmm36-jul2022.ff/tip3p.itp"
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
Thus, I should delete these lines from topol.top, copy the tip3p.itp to my folder, open it, and at the end include these lines:
I don’t think it really matters, though. What you want to make sure is that the atom coordinates are close to 0, especially if you ever want to specify where you insert it.
Your solution, putting the restraints in the topol.top should work just as well and in that case you don’t need to modify the water model itp file. But do you only want to restrain the O atom (with only one atom restrained), i.e., not the water molecule orientation? Just make sure to remember to keep that restraint section right after including the water model in the topology file.
You don’t have to add the restraints more than once per molecule type.
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
That means that you apply restraints on atom 1 in the tip3p.itp file, or as it says in the comment: “Position restraint for each water oxygen”. If you want to restrain the hydrogens as well (i.e., keep the water molecule orientations) you should:
#ifdef POSRES_WATER
; Position restraint for each water oxygen and hydrogen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
2 1 1000 1000 1000
3 1 1000 1000 1000
#endif
Is there a way to prevent the Na ions from getting closer to the protein/ligand/water? After energy minimization, some Na atoms are closest to the protein, and atoms of the amino acid side chain and water molecules move toward these ions. I’ve tried to include the posre_na, but it didn’t solve it.
Since there is no bulk water to shield the ions they will be attracted to the protein. You can restrain the ions, but if the restrained water molecules move towards the ions, you get an idea of the forces involved. You can of course increase the restraint forces. If the simulation box is large enough and the ions far away from the protein (and each other) it might help. As I indicated above, simulating a protein in vacuum will probably not give good results, it might even unfold. There is no possibility to use implicit solvent in GROMACS anymore. Otherwise, that might have helped.
Ok. I’m trying here to increase the box and the forces…
Let me ask if I increase the box, should I change anything in my mdp file because of this change? Something like my rcoulomb, rvdw, rlist, etc
And one last question. I want to restrain the position of the protein mainchain atoms to perform the EM (only side-chain atoms, ligands and water molecules are “free” to move). Are these steps right?
gmx make_ndx -f complex_ion.gro -o index_mchain.ndx (complex_ion is the gro file with protein/ligand/water/Na ions)
You do not need to change the mdp settings based on the simulation box size. After having increased the box size, you can manually move the ions (set the coordinates in the .gro file) so that they are far away from each other and from the protein and ligand.