# Assessing the stability of a mutant protein withGROMACS **URL:** \ **Category:** User discussions\ **Tags:** mdrun, simulation-setup\ **Created:** [April 16, 2024, 1:40pm UTC](https://gromacs.bioexcel.eu/t/assessing-the-stability-of-a-mutant-protein-withgromacs/8899 "2024-04-16T13:40:09Z")\ **Posts on this page:** 7\ **Page:** 1
**Author:** ![tkoh](https://avatars.discourse-cdn.com/v4/letter/t/a9adbd/32.png) [@tkoh](https://gromacs.bioexcel.eu/u/tkoh)\ **Post date:** [April 16, 2024, 1:40pm UTC](https://gromacs.bioexcel.eu/t/assessing-the-stability-of-a-mutant-protein-withgromacs/8899/1 "2024-04-16T13:40:09Z")
GROMACS version: 2021.4-Ubuntu-2021.4-2 GROMACS modification: No Hi, So I am trying to assess the stability of a mutant protein with respect to its native conformation, in essence obtain a ddg value as well as possibly a break down of the components that contribute to the change in stability. I have gone through the free energy of solvation tutorial, the lysozome tutorial, as well as the umbrella sampling tutorial. A superior recommended I use umbrella sampling but I don’t understand how the pulling concept would be applied to my single chain protein. It might be asking for a lot but if someone who perhaps has experience in obtaining ddg values from gromacs or have done similar simulations could assist in directing me to a tutorial or protocol that may provide me with some guidance. I did want to make use of thermodynamic integration for the MD simulations, however I was unable to find a tutorial to guide this process. ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [April 18, 2024, 3:54pm UTC](https://gromacs.bioexcel.eu/t/assessing-the-stability-of-a-mutant-protein-withgromacs/8899/2 "2024-04-18T15:54:58Z")
Indeed, alchemical simulations (TI) are a better option than US for estimating the effects of mutation on a relative free energy between conformations. You’ll have to have a model for the reference (here: unfolded) state, and it can be just a tiny peptide with the same transition, but you need to close the thermodynamic cycle to have a meaningful estimate. Here’s a great paper that explains the key concept and choices: [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8388617/](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8388617/) The PMX webserver will help you create alchemical amino acids: [http://pmx.mpibpc.mpg.de/](http://pmx.mpibpc.mpg.de/) In Fig 8 of that paper you will find the most popular strategies for running alchemical simulations. Depending on your choice, you will find different tools and tutorials; I should soon have a tutorial for case (D), I have been postponing it for a while, but there are excellent ones for a similar case (solvation free energy): [http://www.mdtutorials.com/gmx/free\_energy/index.html](http://www.mdtutorials.com/gmx/free_energy/index.html) [https://tutorials.gromacs.org/free-energy-of-solvation.html](https://tutorials.gromacs.org/free-energy-of-solvation.html) ---
**Author:** ![tkoh](https://avatars.discourse-cdn.com/v4/letter/t/a9adbd/32.png) [@tkoh](https://gromacs.bioexcel.eu/u/tkoh)\ **Post date:** [July 3, 2024, 7:17am UTC](https://gromacs.bioexcel.eu/t/assessing-the-stability-of-a-mutant-protein-withgromacs/8899/3 "2024-07-03T07:17:12Z")
Hi, So I have made use of PMX for the unfolded state and used their tripeptide database for the CHARMM36 ff. However when I attempt the command pdb2gmx on the tripeptide: gmx\_mpi pdb2gmx -f tripep\_V2C.pdb -o conf.pdb -p topol.top -ff charmm36m-mut -water tip3p I get the following error: :-) GROMACS - gmx pdb2gmx, 2023.2 (-: Executable: /usr/local/gromacs/bin/gmx\_mpi Data prefix: /usr/local/gromacs Working dir: /home/tammara/Documents/gmx\_test Command line: gmx\_mpi pdb2gmx -f tripep\_V2C.pdb -o conf.pdb -p topol.top -ff charmm36m-mut -water tip3p Using the Charmm36m-mut force field in directory /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff going to rename /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.r2b Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.r2b Reading tripep\_V2C.pdb… Read ‘’, 42 atoms Analyzing pdb file Splitting chemical chains based on TER records or chain id changing. There are 1 chains and 0 blocks of water and 5 residues with 42 atoms chain #res #atoms 1 ‘X’ 5 42 All occupancies are one All occupancies are one Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/atomtypes.atp Reading residue database… (Charmm36m-mut) Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.rtp Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/mutres.rtp Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/mutres\_dna.rtp Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.hdb Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.n.tdb Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.c.tdb Back Off! I just backed up topol.top to ./#topol.top.8# Processing chain 1 ‘X’ (42 atoms, 5 residues) Identified residue ACE1 as a starting terminus. Identified residue CT35 as a ending terminus. 8 out of 8 lines of specbond.dat converted successfully Start terminus ACE-1: NH3+ End terminus CT3-5: COO- * * * Program: gmx pdb2gmx, version 2023.2 Source file: src/gromacs/gmxpreprocess/pdb2top.cpp (line 1115) Fatal error: atom N not found in buiding block 1ACE while combining tdb and rtp ## For more information and tips for troubleshooting, please check the GROMACS website at [Common Errors — GROMACS webpage https://www.gromacs.org documentation](http://www.gromacs.org/Documentation/Errors) * * * MPI\_ABORT was invoked on rank 0 in communicator MPI\_COMM\_WORLD with errorcode 1. NOTE: invoking MPI\_ABORT causes Open MPI to kill all MPI processes. You may or may not see output from other processes, depending on exactly when Open MPI kills them. To that effect I realised I had to use the -ter option to specify the capped terminal ends, and so I added the option and selected none for both ends, however I get the following: gmx\_mpi pdb2gmx -f tripep\_V2C.pdb -o conf.pdb -p topol.top -ff charmm36m-mut -water tip3p -ter :-) GROMACS - gmx pdb2gmx, 2023.2 (-: Executable: /usr/local/gromacs/bin/gmx\_mpi Data prefix: /usr/local/gromacs Working dir: /home/tammara/Documents/gmx\_test Command line: gmx\_mpi pdb2gmx -f tripep\_V2C.pdb -o conf.pdb -p topol.top -ff charmm36m-mut -water tip3p -ter Using the Charmm36m-mut force field in directory /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff going to rename /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.r2b Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.r2b Reading tripep\_V2C.pdb… Read ‘’, 42 atoms Analyzing pdb file Splitting chemical chains based on TER records or chain id changing. There are 1 chains and 0 blocks of water and 5 residues with 42 atoms chain #res #atoms 1 ‘X’ 5 42 All occupancies are one All occupancies are one Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/atomtypes.atp Reading residue database… (Charmm36m-mut) Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.rtp Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/mutres.rtp Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/mutres\_dna.rtp Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.hdb Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.n.tdb Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.c.tdb Back Off! I just backed up topol.top to ./#topol.top.9# Processing chain 1 ‘X’ (42 atoms, 5 residues) Identified residue ACE1 as a starting terminus. Identified residue CT35 as a ending terminus. 8 out of 8 lines of specbond.dat converted successfully Select start terminus type for ACE-1 0: NH3+ 1: NH2 2: 5TER 3: None 3 Start terminus ACE-1: None Select end terminus type for CT3-5 0: COO- 1: COOH 2: CT2 3: 3TER 4: None 4 End terminus CT3-5: None Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/merged.arn Opening force field file /usr/local/lib/python3.10/dist-packages/pmx/data/mutff/charmm36m-mut.ff/mutres.arn Checking for duplicate atoms… Generating any missing hydrogen atoms and/or adding termini. Now there are 5 residues with 42 atoms Making bonds… atom HV is missing in residue V2C 3 in the pdb file atom HV is missing in residue V2C 3 in the pdb file You might need to add atom HV to the hydrogen database of building block V2C in the file mutres.hdb (see the manual) You might need to add atom HV to the hydrogen database of building block V2C in the file mutres.hdb (see the manual) * * * Program: gmx pdb2gmx, version 2023.2 Source file: src/gromacs/gmxpreprocess/pdb2top.cpp (line 1567) Fatal error: There were 1 missing atoms in molecule Protein\_chain\_X, if you want to use this incomplete topology anyhow, use the option -missing ## For more information and tips for troubleshooting, please check the GROMACS website at [Common Errors — GROMACS webpage https://www.gromacs.org documentation](http://www.gromacs.org/Documentation/Errors) * * * MPI\_ABORT was invoked on rank 0 in communicator MPI\_COMM\_WORLD with errorcode 1. ## NOTE: invoking MPI\_ABORT causes Open MPI to kill all MPI processes. You may or may not see output from other processes, depending on exactly when Open MPI kills them. ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [July 3, 2024, 11:34am UTC](https://gromacs.bioexcel.eu/t/assessing-the-stability-of-a-mutant-protein-withgromacs/8899/4 "2024-07-03T11:34:33Z")
Arguably working with the local version of pmx is a bit more tricky - in principle, the webserver should give you working topologies so that you don’t have to go through pdb2gmx in order to obtain them. What you’re getting are pretty common errors that come from missing hydrogen database entries for alchemical residues, or automatic selection of terminal patches. ---
**Author:** ![tkoh](https://avatars.discourse-cdn.com/v4/letter/t/a9adbd/32.png) [@tkoh](https://gromacs.bioexcel.eu/u/tkoh)\ **Post date:** [September 9, 2024, 1:15pm UTC](https://gromacs.bioexcel.eu/t/assessing-the-stability-of-a-mutant-protein-withgromacs/8899/5 "2024-09-09T13:15:23Z")
Hi, So I have gotten pretty far with these free energy simulations however something strange is happening when I get to the non-equilibrium part of the MD. When I try to extract 50 snapshots from the 5ns equilibrium simulation at every 100ps starting from 100ps, I only generate 47 frames instead of 50. I have gone through my MD parameters but I don’t see anything that particularly stands out. Below is the command and the output, along with the output from gmx check on the .trr file. Attached is the mdp files used for equil and non-equil simulations. [f\_equil.mdp](https://gromacs.bioexcel.eu/uploads/short-url/zndAo7VGSi8ePLO2aBGKoFZdpZx.mdp) (4.0 KB) The command and out: Command line: gmx\_mpi trjconv -f equil.trr -s equil.tpr -sep -b 100 -o frame\_.gro Will write gro: Coordinate file in Gromos-87 format Reading file equil.tpr, VERSION 2023.2 (single precision) Reading file equil.tpr, VERSION 2023.2 (single precision) Group 0 ( System) has 3003 elements Group 1 ( Protein) has 45 elements Group 2 ( Protein-H) has 21 elements Group 3 ( C-alpha) has 3 elements Group 4 ( Backbone) has 11 elements Group 5 ( MainChain) has 15 elements Group 6 ( MainChain+Cb) has 16 elements Group 7 ( MainChain+H) has 19 elements Group 8 ( SideChain) has 26 elements Group 9 ( SideChain-H) has 6 elements Group 10 ( Prot-Masses) has 45 elements Group 11 ( non-Protein) has 2958 elements Group 12 ( Water) has 2958 elements Group 13 ( SOL) has 2958 elements Group 14 ( non-Water) has 45 elements Select a group: trr version: GMX\_trn\_file (single precision) Skipping frame 0 time 0.000 Skipping frame 1 time 40.000 Skipping frame 2 time 80.000 Reading frame 3 time 120.000 Reading frame 4 time 160.000 → frame 0 time 120.000 Reading frame 5 time 200.000 → frame 1 time 160.000 Reading frame 6 time 240.000 → frame 2 time 200.000 Reading frame 7 time 280.000 → frame 3 time 240.000 Reading frame 8 time 320.000 → frame 4 time 280.000 Reading frame 9 time 360.000 → frame 5 time 320.000 Reading frame 10 time 400.000 → frame 6 time 360.000 Reading frame 11 time 440.000 → frame 7 time 400.000 Reading frame 12 time 480.000 → frame 8 time 440.000 Reading frame 13 time 520.000 → frame 9 time 480.000 Reading frame 14 time 560.000 → frame 10 time 520.000 Reading frame 15 time 600.000 → frame 11 time 560.000 Reading frame 16 time 640.000 → frame 12 time 600.000 Reading frame 17 time 680.000 → frame 13 time 640.000 Reading frame 18 time 720.000 → frame 14 time 680.000 Reading frame 19 time 760.000 → frame 15 time 720.000 Reading frame 20 time 800.000 → frame 16 time 760.000 Reading frame 30 time 1200.000 → frame 26 time 1160.000 Reading frame 40 time 1600.000 → frame 36 time 1560.000 Reading frame 50 time 2000.000 → frame 46 time 1960.000 Last frame 50 time 2000.000 → frame 47 time 2000.000 Last written: frame 47 time 2000.000 The trr file check: gmx\_mpi check -f equil.trr ``` :-) GROMACS - gmx check, 2023.2 (-: ``` Executable: /usr/local/gromacs/bin/gmx\_mpi Data prefix: /usr/local/gromacs Working dir: /home/user/Documents/test\_nonequil\_process Command line: gmx\_mpi check -f equil.trr Checking file equil.trr trr version: GMX\_trn\_file (single precision) Reading frame 0 time 0.000 # Atoms 3003 Last frame 50 time 2000.000 Item #frames Timestep (ps) Step 51 40 Time 51 40 Lambda 51 40 Coords 51 40 Velocities 51 40 Forces 51 40 Box 51 40 GROMACS reminds you: “Make the Floor Burn” (2 Unlimited) [f\_nonequil.mdp](https://gromacs.bioexcel.eu/uploads/short-url/odamTO30Y5tUcExMtXIbAegxwbo.mdp) (4.2 KB) ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [September 9, 2024, 1:52pm UTC](https://gromacs.bioexcel.eu/t/assessing-the-stability-of-a-mutant-protein-withgromacs/8899/6 "2024-09-09T13:52:01Z")
You say “When I try to extract 50 snapshots from the 5ns equilibrium simulation at every 100 ps starting from 100 ps, I only generate 47 frames instead of 50.” From the output you posted it seems like the simulation is 2 ns and the trajectory is written every 40ps. As you can also see from the output, starting from 100ps means that you discard the frames at 0, 40, 80 ps. So, what you do get is 48 frames (0 to 47), instead of 51 (0 to 50). If you expect the simulation to write every 100 ps and run for 5 ns, you need to check your mdp file and also that it matches the tpr file. Edit: Checking your f\_equil.mdp it can be seen that the source of the problem is that `dt = 0.0008`. It seems like you assume that it is 0.002. ---
**Author:** ![tkoh](https://avatars.discourse-cdn.com/v4/letter/t/a9adbd/32.png) [@tkoh](https://gromacs.bioexcel.eu/u/tkoh)\ **Post date:** [September 10, 2024, 6:25am UTC](https://gromacs.bioexcel.eu/t/assessing-the-stability-of-a-mutant-protein-withgromacs/8899/7 "2024-09-10T06:25:26Z")
You are absolutely right, thank you for pointing that out.