GROMACS version: 2020
I am trying to run the Energy minmization of the Ribosome (PDB ID 6i0y) with a Protein in the exit channel. I have freezed the Ribosome atoms using freezegrps in mdp file and left the protein free.
But once I run the energy minimization I result in
very high potential energy in +7.65 *e+07 range which is extermely high and
did not reach the condition Fmax<1000.
Using this state I cannot run the Production MD simulation
I tried to removing peptide and running it still only with ribosome the potential energy is very high. May I know how to resolve this issue. Any input would be really helpful.
Some additional questions:
- Does freezing of atoms still influence in the potential Energy?
- Can I eliminate the forces of particular atoms which are in boundary, will it help in solving this issue.
Yes, frozen atoms still contribute to the energy - they are simply not moved during geometry optimization. If you’re freezing atoms that clash with other ones, your potential energy will stay high because the procedure won’t be able to resolve the clashes.
But usually Gromacs gives you extra information on where the highest forces occur, so the reasonable thing is to visualize the system and see what’s happening around there. With many automatic processing tools (and even deep learning models), you can get side or main chains going through aromatic rings or other crazy arrangements, so be sure to go back to the structure and resolve any major issues it might have.
But my focus is to freeze the Ribosome group atoms. Since just freezing still calculates the forces but to avoid the force calculation I was thinking to include energygrp-excl parameter this way it would not calculate the forces and result in high potential energy. Rather only the forces will be calculated from the atoms which are not freezed and not included in energygrp-excl But I am not able to set it in Gromacs 2023.
I tried checking with GROMACS 2019 as few posts suggested but results in issue of installation 403 Forbidden Regression_Test * http://gerrit.gromacs.org/download/regressiontests-2019.tar.gz
May I know is this right way and how to fix this issue of fixing specific atoms, not resulting high potential energy.
Don’t freeze them. Apply a very strong position restraint. I’m not sure that energygrp-excl works with the Verlet scheme. It needed the group scheme, which was removed (unless I am mistaken).
Hi @jalemkul I understand that Position restraints are used to limit the movement of atoms while still allowing them some freedom to vibrate. This is typically implemented by adding a harmonic potential to restrain the atoms around their initial positions.
But my goal is to fix the high potential energy and reduce time both, If I freeze the atoms fully and also eliminate the forces so that they dont result in high potential energy and also help me in not wasting my time calculating those forces.
Freezing does not skip the force calculation and therefore does not save time. It just skips the update of positions. That’s why you have unsolvable high forces. Again I emphasize the only legitimate way to do this in GROMACS is with position restraints.
Okay I get it. I will restrain it but along with that what if I add the if condition in the source code to not calculate the forces of certain atoms based on their atom id. I think this way I will be able to not only fix atoms but also reduce time.
In most systems, the majority of time is spent on calculating the interactions involving solvent anyway. If you want to focus on a small subsystem of the ribosome and save computational time, it might be more reasonable to extract a subset of the full complex and (optionally) restrain the interface with the rest of the ribosome, given it’s sufficiently far away. This has its own challenges and approximations, but for certain purposes the increase in sampling time might outweigh the error resulting from reduced mobility.
Thank you, Yeah I will try cropping the ribosome further to small system which i want to focus.
But, I see that restraining the subset of the atoms does not fully fix the positions there is still movement in the atoms I see during the simulation. Am I doing something wrong let me know:
I need to fix the ribosome atoms fully except for protein in exit channel, I tried using freezegrps but it was leading to very high potential energy, so my goal was to remove the forces on those ribosome atoms for which I am looking into source code still, which is not resolved.
As alternative I found position restraints will help me in fixing the position by applying the counter force so tried including in the topology.top file but did not work as mentioned below
#include "topol_Protein_chain_2.itp"
#ifdef POSRES
#include "posre_Protein_chain_2.itp"
#endif
#include "topol_Protein_chain_7.itp"
#include "topol_RNA_chain_A.itp"
#ifdef POSRES
#include "posre_RNA_chain_A.itp"
#endif
- I run the Energy Minimization on freezegrps led to high potential energy in e+14 and Fmax did not converge.
- so I tried with restraints on ribosome atoms but no restraints on protein chain 7 mentioned above. Potential energy resulted to e+07 which is still high, Fmax did converge.
But the Ribosome atoms displaced from the original position which i wanted to be fixed.
Continuing with the output of Energy Minimization then running nvt.
Lastly, I ran the Production MD. This displaced the ribosome atoms further with increased RMSD in ribosome atoms alone.
Am i doing something wrong?
- So position restraints does not fully fix the atoms position, they still jiggle and also displace from original position?
@jalemkul @milosz.wieczor
Restraints and constraints work very differently on the theoretical level. The former, as you say, just allow for some wiggle a tiny bit around the reference position (with the amount of wiggling controlled by the force constant(s). The latter require special treatment in the algorithms (integration, thermo- and barostats etc), as they remove some degrees of freedom.
Again, whenever you refer to high energy in minimization, I suggest you look at the highest-force atom indicated by mdrun, and check it visually. This will help you understand where the clashes are located, and whether they are acceptable given your strategy to put restraints.
If parts of your system drift away even when they were supposed to be held by position restraints, check the .log to make sure there’s an energy contribution coming from position restraints. If not, that might indicate you failed to turn them on correctly (e.g. because of a misspelled or commented-out directive in .mdp).