GROMACS version:
GROMACS modification: Yes/No
Here post your question
Hello,
I have run my protein with a metal and ligand for 30 ns. I need to run the md simulation 30 ns. I ran into the segmentation fault error when performing nvt or npt. I need npt to run md. I adjusted the energy minimization steps, but that doesn’t work. I looked online for help and found the pressure coupling in npt.mdp can be adjusted.
pcoupl = Parrinello-Rahman
in the npt.mdp file with
pcoupl = C-rescale
This doesn’t fix the problem. I created my protein with the calcium ion, disulfide bond in Amber. There was no problem “building” in Amber. I tried using Amber forcefield and Charm forcefield in the topology, and I tried using either dry protein or solvate protein generated by acpype to run nvt, npt, and md, but I still get a segmentation fault.
Hi Ryan,
I theory you should not get a segmentation fault at all, but sometimes it happens if the system is exploding. Since it happens already in your first equilibration stage (nvt) I would presume that the system is not energy minimized well enough.
Do you get any messages before the segmentation fault? Often those can give a hint which atoms are moving fast, or which bonds are extended too far. It’s also good to check the output from the energy minimization run to see which atoms have the highest force. You say that you adjusted the energy minimization steps, but that it did not work. Does the energy minimization converge? Sometimes it’s good to force the energy minimization to run on CPU (if it’s running on GPU with the default settings), by using the -nb cpu option.
I had the energy minimization step up to 5000 at the start. I decreased the number steps to the energy minimization knowing that at some point it will not converge. Then, increased the steps to find the balance between where energy minimization converges and doesn’t converge, and use that energy minimization to run the md nvt run, which still gave me a segmentation fault. I can get the energy minimization to converge at above 250 steps. I haven’t tried forcing the energy minimization to run on CPU.
I re-ran the energy minimization again. The system converges with no errors. It converged 298 steps. I used steepest descent minimization. I ran the nvt step to generate nvt.tpr. The only warning I got was:
WARNING 1 [file mdp file]:
Malformed define option -POSRES
I have my position restraint stated in the mdp and top files as POSRES to restraint the system at 1000. I think gmx doesn’t recognize my position restraint in the mdp and top file?
I changed the position restriants to DPOSRES. I tried running the mdrun em at 5000 steps with Verlet cut-off and cpu on my laptop.
gmx mdrun -nb cpu -deffnm em, the command line
I get a segmentation fault running this command, but without -nb cpu it converges at 298 steps.
Yes, I ran the gmx mdrun -nb cpu -deffnm em -v with 5000 steps. It converged at 294 steps with a maximum force of:
Steepest Descents converged to machine precision in 294 steps,
but did not reach the requested Fmax < 1000.
Potential Energy = -4.2708806e+05
Maximum force = 1.1458638e+05 on atom 3896
Norm of force = 1.2184215e+03
I expect problems with the input structure. I would start with checking atom 3896 and its surroundings. I would not be surprised if it’s overlapping with (or very close to) another atom or if there is, e.g., a bond nearby passing through a ring structure, or if a chemical bond is far too long.
Okay, thank you. I looked at atom 3896. It is a water molecule used in the metal optimization in Amber and Games. From what I found there are two water molecules that optimize the metal when running GAMES and AMBER. There was no problem running GAMES and Amber, and no problems converting the amber format to GROMACS in gro and top files. The distances in the GROMACS format from the metal to the water molecules is 2.402 and 3.077 angstroms.
Do you mean the connectivities/bonds? I’m afraid not. It’s not stored in the gro file.
If the water molecules has been optimized in GAMESS it is possible that you need to make water molecules flexible during energy minimization. See if adding define = -DFLEXIBLE to the mdp file helps.
Otherwise, how large is the system? Would you be able to post the gro file here (or on a file sharing service)?
I was thinking of manipulating the coordinates in the gro file to adjust the distance, but that won’t change the bond parameters generated by GAMES. I tried running the simulation with “define = -DFLEXIBLE” in the mdp file and it outputs a segmentation fault. I can post it as a top file here. solv_GMX_gro.top (1.3 MB)
I had a look at the system. I don’t see any apparent problems. But the distance between the water hydrogen 3897 to the metal is only 2.30 Å. What metal is this? It says CA, but I guess it’s not a Ca2+ ion, right? Especially since the QM optimized water is orienting its hydrogen atom towards it. What is its VdW radius? Is it correctly parameterized in the force field?
I’m still quite sure it’s the input structure that’s the problem. Without the topology and mdp parameters I’m afraid I cannot help you much more at the moment. Using define = -DFLEXIBLE was mainly to see if applying constraints from an unfavorable structure caused any problems. Since it did not help there’s no reason to retain it.
My enzyme has no metal when you download it. According to scientific articles, they performed matchmake in Chimera to copy the coordinates of the metal ion in another enzyme and “paste it” relative to my enzyme. The metal, CA, is a calcium +2 ion, vdw radius 231. It was parameterized by GAMES, then AMBER MCPB.py builds the GAMES output into the enzyme. When I check the Amber generated topology, .prmtop, it shows the correct atomic number for calcium, 20. It is parameterized in Amber forcefield. If I add:
; Include forcefield parameters #include “amber99sb.ff/forcefield.itp”
to the topology file I get warnings saying it already recognizes the bonds and I can suppress the warning with maxwarn, or there is no warnings if I remove the forcefield parameters from the topology. Both will generate the em.tpr, but with the forcefield parameters leads to a segmentation fault appears when running gmx mdrun -deffnm -nb cpu em -v. I didn’t run GAMES on my computer. The only problem with the GAMES output was naming of the two water molecules. I also checked where I got the metal ion from, and it said Calcium ion 2+. 4eb0_solv_GMX.top (3.0 MB) em.mdp (1.2 KB) nvt.mdp (1017 Bytes) npt.mdp (1.2 KB) md.mdp (1.2 KB)
Having a closer look at this I’m afraid I can’t help you. The input coordinates energy minimize to a state that will crash. Moving the two QM optimized water molecules a bit still did not help.
I’m not an expert at parameterizing these kinds of systems. Unless you are 100% sure about what you are doing, I would start over and use common MM force field parameters and skip the QM step. Currently the partial charge on your Ca2+ ion and the neighboring waters, as well as some of the nearby protein side chains, will have their partial charges set by QM. But these will move. I think you should either stick to MD simulations (with conventional parameters) or run QM/MM for this part of the system. The latter option may be better for coordinated ions.
If you are completely sure that your parameters and input coordinates are correct then I’m afraid you’ll have to find another MD software to simulate this system.
I just noticed that there are covalent bonds between the Ca2+ ions and the coordinated waters. This will at least prevent them from leaving. But apparently something needs to be done about the input coordinates and/or force field parameters.