.. code:: ipython3 # Execute this first # # * trigger notebook styling # * check if notebook had been modified since its distribution # # Note: executing any cells before this modifies the notebook. # %run src/init_notebooks.py hide_toggle() check_notebook() A Beginner’s Guide to Perform Molecular Dynamics Simulation of a Membrane Protein using GROMACS ----------------------------------------------------------------------------------------------- **Authors :** Farzaneh Jalalypour, Maryam Majdolhosseini, Alessandra Villa | **Goal :** Learn step-by-step how to run a molecular dynamics simulation of a simple membrane protein using GROMACS | **Reading time :** 40 minutes | **Software :** GROMACS 2023.2, python modules (numpy, matplotlib, nglviewer, MDTraj, and pandas) | **Optional Software :** Visualization software `VMD `__, Xmgrace plotting tool | **Tutorial Source :** tutorials.gromacs.org | **Version :** Final **Small talk with the audience** This tutorial is designed for beginner GROMACS users and assumes you understand the basics of MD simulations and are familiar with GROMACS files and how to generate and execute them. Don’t worry if you are not yet acquainted. We have already developed an `“Introduction to Molecular Dynamics” tutorial `__ that covers all you need to know to get started with this one. Please note that this tutorial is for educational purposes and and won’t go through every step in details. For your system, you must change the settings accordingly (see `GROMACS documentation `__). Preparations to run this notebook --------------------------------- .. code:: ipython3 # Change to the data directory # Note that executing this command twice will result in an error you can ignore %cd data Model system for this tutorial ------------------------------ For this tutorial, we will utilize maltoporin channel, a trimeric protein located at the outer membrane of Gram-negative bacteria that facilitates the translocation of the polysaccharide maltodextrin. You can download the 3D structure of maltoporin channel (PDB code 1MAL) from RCSB website https://www.rcsb.org/structure/1MAL or use the previously downloaded file, “1mal.pdb” in the “input” directory. **Visualize the structure** .. code:: ipython3 import nglview as ng view = ng.show_structure_file("input/1mal.pdb") view #click and drag to rotate, zoom with your mouse wheel #for more information on this viewer have a look at https://github.com/nglviewer/nglview System preparation for this tutorial ------------------------------------ **Building the protein-membrane system in CHARMM-GUI** We are now ready to embed the protein structure in the membrane in the proper location and orientation and construct the membrane composition we desire. To do this, we utilized the CHARMM-GUI input Generator, a handy web-based tool to generate GROMACS inputs for the protein-membrane system. In order to access CHARMM-GUI (http://www.charmm-gui.org), a web browser such as Chrome, Firefox, etc. is required. In this Tutorial, rather than providing an in-depth discussion of how to create input via CHARMM-GUI, we will focus on how to perform a correct simulation using GROMACS. However, you can find useful information on `Membrane Builder Intro `__ and in this `Membrane Proteins Tutorial `__ . CHARMM-GUI generates all required input files as a zip file. You can unzip the file using the command below, or directly use the previously unzipped and modified folder named “charmm-gui-1MAL”. This folder contains a diverse set of files (not available on reference directory). .. code:: ipython3 ## command to unzip the zip file. The command is executable when you remove the hashtags #import shutil #shutil.unpack_archive('input/charmm-gui-1MAL.zip') #print("Archive file unpacked successfully.") What we require, though, is in a “gromacs” folder. Let’s take a closer look: .. code:: ipython3 %ls input/charmm-gui-1MAL/gromacs/ In a prior tutorial titled `“Introduction to Molecular Dynamics” `__, you may get an explanation for each file including structure and topology files. Some of these mdp files may need to be changed to fit the system as further discussed in the next section. Now lets create a new folder, called “run”. Here we will perform the minimization and equilibration steps. Ensure that you are always in the correct working directory, you can use the ``pwd`` command, which stands for “print working directory”. The command can be used to determine what directory you are currently working in. .. code:: ipython3 # Creat working directory named run %mkdir run # Change to the directory "run" %cd run .. code:: ipython3 # Make sure you are in the right folder. This will show you the path :) !pwd Modified MDP files ================== Before starting the simulations we need the topology, parameter, and structure files. We will copy them from the ``CHARMM-GUI-1MAL`` folder to the ``run`` directory. Please note that mdp files are customized to fit the requirements of the system and they are located in the ``input/mdp`` folder. .. code:: ipython3 ##copy required files to a folder !cp ../input/charmm-gui-1MAL/gromacs/step5_input.gro . !cp ../input/charmm-gui-1MAL/gromacs/step5_input.pdb . !cp ../input/charmm-gui-1MAL/gromacs/topol.top . !cp ../input/charmm-gui-1MAL/gromacs/index.ndx . !cp -r ../input/charmm-gui-1MAL/gromacs/toppar . #updated mdp files !cp ../input/mdp/*.mdp . Because GROMACS always assumes a default value for each parameter, simulations may be conducted with an entirely empty mdp file. Hence, you must double-check each option to ensure that the default settings are acceptable. If not, you can modify the mdp file. GROMACS is user-friendly and prepared an intriguing list of possibilities that is clearly explained in `mdp options webpage `__. To make your life easier, in all mdp files of this tutorial, a small explanation is included next to each parameter. Now we have all required files to start minimization and equilibration parts: .. code:: ipython3 %ls +-----------------------+-----------------------+-----------------------+ | Step 6.0 | EM | Energy minimisation | +-----------------------+-----------------------+-----------------------+ | Step 6.1 | eq | constant number (N), | | | uilibration_NVT_step1 | constant-volume (V), | | | | and | | | | constant-temperature | | | | (T), | +-----------------------+-----------------------+-----------------------+ | Step 6.2 | eq | | | | uilibration_NVT_step2 | | +-----------------------+-----------------------+-----------------------+ | Step 6.3 | eq | constant number (N), | | | uilibration-NPT-step1 | constant-pressure | | | | (P), and | | | | constant-temperature | | | | (T), | +-----------------------+-----------------------+-----------------------+ | Step 6.4 | eq | | | | uilibration-NPT-step2 | | +-----------------------+-----------------------+-----------------------+ | Step 6.5 | eq | | | | uilibration-NPT-step3 | | +-----------------------+-----------------------+-----------------------+ | Step 6.6 | eq | | | | uilibration-NPT-step4 | | +-----------------------+-----------------------+-----------------------+ | Step 7 | production | | +-----------------------+-----------------------+-----------------------+ Let’s start! Although it will take time, you can run the commands below one by one to generate files in every step. Fortunately, we have prepared output of runs ready to use if you prefer the quick route (see ``input/reference`` folder). Energy minimisation =================== Energy minimisation (EM) is used to relax the protein-membrane system by removing steric conflicts and improper geometry, e.g. to reduce excessive forces caused by too-close particle interaction in the system. We will first generate a binary (.tpr) file which assembles the structure(.gro), topology(.top), and simulation parameters(.mdp). Then we run energy minimization using ``gmx mdrun``. But first we go through the parameters in the minimization.mdp file: .. code:: ipython3 #revised minimization.mdp !cat ../input/mdp/step6.0_minimization.mdp .. code:: ipython3 ## EM-step1 > to generate binary file !gmx grompp -f step6.0_minimization.mdp -o minimization.tpr -c step5_input.gro -r step5_input.gro -p topol.top **Note:** From now on, if you intend to execute commands on your system, refer to cells labeled with an odd number ``##{number]##``. Remove one hashtag from each line, and avoid running cells marked with even numbers. Alternatively, if you wish to utilize pre-prepared outputs, only run cells with even numbers. .. code:: ipython3 ##1## EM-step2 > to invoke mdrun #!gmx mdrun -v -deffnm minimization ##The command is executable when you remove one of the hashtags Five files (tpr,gro,trr,edr,log) are generated after you run the minimization command. You can always find the previously prepared files for each step in the reference folder: .. code:: ipython3 ##2## EM-step2 > To transfer minimization output from the reference directory to working directory !cp ../reference/minimization* . ##remove one hashtag to transfer pre-prepared minimization files to working directory .. code:: ipython3 #folder of pre-prepared outputs %ls minimization* Let’s check if the minimization criteria are reached by looking at the last lines of the log file: .. code:: ipython3 !tail -19 ../reference/minimization.log There are two very important factors to evaluate to determine if energy minimization was successful. The first is the potential energy (printed at the end of the EM process, even without -v). The second important feature is the maximum force, Fmax, the target for which was set in minim.mdp - “emtol = 1000.0” - indicating a target Fmax of no greater than 1000 kJ/(mol nm). It is possible to arrive at a reasonable Epot with Fmax > emtol. If this happens, your system may not be stable enough for simulation. Evaluate why it may be happening, and perhaps change your minimization parameters (integrator, emstep, etc). .. code:: ipython3 ## Minimized structure import nglview as ng view = ng.show_structure_file("../reference/minimization.gro") ## Lipids view.add_representation("ball+stick",selection="POPC") ## water view.add_representation("licorice",selection="TIP3P", opacity=0.9) view.center() view #Click and drag to rotate, zoom with your mouse wheel #For more information on this viewer have a look at https://github.com/nglviewer/nglview In alternative, remove the comment character (#) to use VMD: .. code:: ipython3 #!vmd ../reference/minimization.gro Let’s do a bit of analysis. The minimization.edr file contains all of the energy terms that GROMACS collects during EM. To extract energy information you can use the tool ``gmx energy`` .. code:: ipython3 !printf "Potential\n0\n" | gmx energy -f minimization.edr -o potential.xvg -xvg none .. code:: ipython3 import pandas as pd df = pd.read_csv('potential.xvg', sep='\s+', header=None, names=['step','energy(kJ Mol-1)']) df.plot('step') Note: to analyse or visualize simulation data in python or jupyter notebooks, we can output a simplified xvg format from gmx-analysis tools with the option -xvg none Equilibration run - NVT ======================= **position restraints** For free proteins in solution, e.g. cytosolic proteins, position restraints are not as crucial as membrane proteins, where different parts of a system must be restrained in different levels and released gradually in several steps of minimization and equilibration. Otherwise you will end up in serious problems such as undesired change in system size and conformation. To set position restraints we need to use the ``define`` option in the simulation parameter file, .mdp, which switch on by -DPOSRES flag. Have a look at the .mdp files for the following runs: .. code:: ipython3 !head -4 ../input/mdp/step6.1_equilibration_NVT_step1.mdp !head -4 ../input/mdp/step6.2_equilibration_NVT_step2.mdp !head -4 ../input/mdp/step6.3_equilibration_NPT_step1.mdp !head -4 ../input/mdp/step6.4_equilibration_NPT_step2.mdp !head -4 ../input/mdp/step6.5_equilibration_NPT_step3.mdp !head -4 ../input/mdp/step6.6_equilibration_NPT_step4.mdp !head -4 ../input/mdp/step7_production_revised.mdp Let’s use chain A (segment PROA) as an example to demonstrate how it works. The position restraint is defined in the last part of the itp file. ``DPOSRES_FC_BB`` and ``POSRES_FC_SC`` specify the type of restraints that apply to the different part of the protein. BB stands for the protein ``BackBone`` and SC for the protein ``Side Chains``! Please execute the command below to see the position restraints set for chain A: .. code:: ipython3 !sed -n '7,33p' ../input/charmm-gui-1MAL/gromacs/topol.top .. code:: ipython3 !sed -n '59914,59932p' ../input/charmm-gui-1MAL/gromacs/toppar/PROA.itp Hint: The number of atoms may be misleading. Please follow the PDB file to ensure which atoms are selected. Here you can find the atom numbers for the first two residues: .. code:: ipython3 !head -n 30 ../input/charmm-gui-1MAL/gromacs/step5_input.pdb As you can see, the 1st and 5th atoms are N and CA protein backbone atoms of the first residue, respectively, and the value is defined by ``POSRES_FC_BB`` in PROA.itp file. The 7th and 9th are atom indexes of CB and CG which belong to the protein side chain and specify with ``POSRES_FC_SC``. Take some moment to find the rest of the atoms and understand the concept by comparing itp and pdb files. The P atoms of the lipid headgroups are constrained by ``DPOSRES_FC_LIPID``. As the value is only determined for the Z direction, POPC lipid molecules can only move in the X-Y plane of the bilayer, maintaining the bilayer’s thickness. ``i`` is the atom index of the P atom of the first POPC. .. code:: ipython3 !sed -n '1271,1276p' ../input/charmm-gui-1MAL/gromacs/toppar/POPC.itp In POPC lipid molecules, two dihedral angles are typically restrained to preserve the cis-isomer. Here, the angles formed by C1-C3-C2-O21 and C28-C29-C210-C211 atoms are described by the atom indexes of the first POPC 25-36-28-30 and 60-63-65-67, respectively, and the value is defined by ``DIHRES_FC``. .. code:: ipython3 !sed -n '1278,1285p' ../input/charmm-gui-1MAL/gromacs/toppar/POPC.itp Equilibration run - NVT ======================= To perform NVT equilibration we need two steps. In each step we will remove the restraints gradually to avoid undesired change in system size and conformation. As discussed before, by removing the first hashtag (#) from the beginning of each line in odd-numbered cells, you may “uncomment” and execute them on your system, but be warned that it will take a while. As an alternative, you can transfer the pre-prepared outputs by executing even-numbered cells to be able to continue. .. code:: ipython3 ##1## equilibration-NVT6.1-step1 > To run #!gmx grompp -f step6.1_equilibration_NVT_step1.mdp -o step6.1_equilibration_NVT_step1.tpr -c minimization.gro -r step5_input.gro -p topol.top -n index.ndx #!gmx mdrun -v -deffnm step6.1_equilibration_NVT_step1 As an alternative, you can transfer the pre-prepared outputs by executing even-numbered cells to be able to continue. .. code:: ipython3 ##2## equilibration-NVT6.1-step1 > To transfer !cp ../reference/step6.1_equilibration_NVT_step1* . Now we run step 2 using results from step 1 as input .. code:: ipython3 ##3## equilibration-NVT6.2-step2 > To run #!gmx grompp -f step6.2_equilibration_NVT_step2.mdp -o step6.2_equilibration_NVT_step2.tpr -c step6.1_equilibration_NVT_step1.gro -r step5_input.gro -p topol.top -n index.ndx #!gmx mdrun -v -deffnm step6.2_equilibration_NVT_step2 As an alternative, you can transfer the pre-prepared outputs by executing even-numbered cells to be able to continue. .. code:: ipython3 ##4## equilibration-NVT6.2-step2 > To transfer !cp ../reference/step6.2_equilibration_NVT_step2* . .. code:: ipython3 #Equilibration - NVT outputs %ls ../reference/*NVT_step* It is crucial to evaluate every step to ensure that the simulation is going right. For this purpose, we should check the temperature using **gmx energy** tool. To plot the data, we used python data analysis library (pandas), but there are several alternatives out there such as Xmgrace tool. **Evaluate temperature** .. code:: ipython3 ## check the temperature > To run !echo 17|gmx energy -f step6.2_equilibration_NVT_step2.edr -o NVT_S2-temp.xvg -xvg none .. code:: ipython3 import pandas as pd df = pd.read_csv('NVT_S2-temp.xvg', sep='\s+', header=None, names=['time (ps)','temperature (K)']) df.plot('time (ps)') .. code:: ipython3 ##Alternative##Plotting data by Xmgrace tool #!xmgrace ../run/NVT_S2-temp.xvg **Evaluate the final structure** It is also strongly recommended to visualize the final structure of each step: .. code:: ipython3 import nglview as ng view = ng.show_structure_file("../reference/step6.1_equilibration_NVT_step1.gro") ## Lipids view.add_representation("ball+stick",selection="POPC") ## water view.add_representation("licorice",selection="TIP3P", opacity=0.9) view.center() view #click and drag to rotate, zoom with your mouse wheel #for more information on this viewer have a look at https://github.com/nglviewer/nglview .. code:: ipython3 import nglview as ng view = ng.show_structure_file("../reference/step6.2_equilibration_NVT_step2.gro") ## Lipids view.add_representation("ball+stick",selection="POPC") ## water view.add_representation("licorice",selection="TIP3P", opacity=0.9) view.center() view **Inside the mdp files** As you may realize, the mdp files for EM and NVT are not the same. The ``integrator`` has been updated to md: .. code:: ipython3 !sed -n '8,8p' ../input/mdp/step6.1_equilibration_NVT_step1.mdp !sed -n '7,7p' ../input/mdp/step6.0_minimization.mdp And a temperature coupling section has been added: .. code:: ipython3 !sed -n '42,55p' ../input/mdp/step6.1_equilibration_NVT_step1.mdp .. code:: ipython3 ## For more details and compare at the entire files remove one hashtag #!cat ../input/mdp/step6.0_minimization.mdp #!cat ../input/mdp/step6.1_equilibration_NVT_step1.mdp We used the v-rescale thermostat, which is the modified version of Berendsen having a stochastic term that assures the generation of a correct canonical ensemble. Sometimes you need to define groups, for instance to couple them separately. Here, SOLU_MEMB stands for protein and membrane, whereas SOLV stands for water and ions. If you’re wondering about how we suggest these names, go through the index files. You can utilize pre-defined groups in the index file or even alter the index file to establish more groups if necessary. Typically, CHARMM-GUI provides five categories in the index file, as shown below: ============= ============ ================ Group Index Molecules ============= ============ ================ [ SOLU ] 1-19260 protein [ MEMB ] 19261-56378 POPC [ SOLV ] 56379-138195 ions and water [ SOLU_MEMB ] 1-56378 protein and POPC [ SYSTEM ] 1-138195 whole ============= ============ ================ The interpretation of this information is as follows: .. code:: ipython3 !printf "q\n" | gmx make_ndx -f step5_input.pdb -n index.ndx -o You may also notice a variation in values if you look at the first lines of NVT mdp files. The sole difference between the prior stage (step6.1_equilibration_NVT_step1) and the final stage of NVT (step6.2_equilibration_NVT_step2) is that position restraints are gradually decreased and relaxed: .. code:: ipython3 !head -6 ../input/mdp/step6.1_equilibration_NVT_step1.mdp !head -4 ../input/mdp/step6.2_equilibration_NVT_step2.mdp Equilibration run - NPT ======================= Pressure equilibration takes place in an NPT ensemble. To keep the pressure stable, we turn the barostat on. We chose the semisotropic pressure coupling, which is isotropic in the x and y directions but not in the z, and is suitable for membrane simulations. Because Berendsen Barostat is outdated and Parrinello-Rahman is better to utilize ‘only’ in production runs when the structure is ‘properly equilibrated’, the c-rescale algorithm is appropriate for this part: **Pressure coupling is on** pcoupl = c-rescale pcoupltype = semiisotropic tau_p = 5.0 compressibility = 4.5e-5 4.5e-5 ref_p = 1.0 1.0 refcoord_scaling = com In order to minimize unintended changes in system size and conformation, NPT equilibration is carried out in four phases, with each step including a progressive removal of restraint. .. code:: ipython3 ##9## equilibration-NPT6.3-step1 > To run #!gmx grompp -f step6.3_equilibration_NPT_step1.mdp -o step6.3_equilibration_NPT_step1.tpr -c step6.2_equilibration_NVT_step2.gro -r step5_input.gro -p topol.top -n index.ndx #!gmx mdrun -v -deffnm step6.3_equilibration_NPT_step1 In alternative: .. code:: ipython3 ##10## equilibration-NPT6.3-step1 > To transfer !cp ../reference/step6.3_equilibration_NPT_step1* . In the next phases, we gradually reduce the constraints: .. code:: ipython3 ##11## equilibration-NPT6.4-step2 > To run #!gmx grompp -f step6.4_equilibration_NPT_step2.mdp -o step6.4_equilibration_NPT_step2.tpr -c step6.3_equilibration_NPT_step1.gro -r step5_input.gro -p topol.top -n index.ndx #!gmx mdrun -v -deffnm step6.4_equilibration_NPT_step2 In alternative: .. code:: ipython3 ##12## equilibration-NPT6.4-step2 > To transfer !cp ../reference/step6.4_equilibration_NPT_step2* . .. code:: ipython3 ##13## equilibration-NPT6.5-step3 > To run #!gmx grompp -f step6.5_equilibration_NPT_step3.mdp -o step6.5_equilibration_NPT_step3.tpr -c step6.4_equilibration_NPT_step2.gro -r step5_input.gro -p topol.top -n index.ndx #!gmx mdrun -v -deffnm step6.5_equilibration_NPT_step3 In alternative: .. code:: ipython3 ##14## equilibration-NPT6.5-step3 > To transfer !cp ../reference/step6.5_equilibration_NPT_step3* . .. code:: ipython3 ##15## equilibration-NPT6.6-step4 > To run #!gmx grompp -f step6.6_equilibration_NPT_step4.mdp -o step6.6_equilibration_NPT_step4.tpr -c step6.5_equilibration_NPT_step3.gro -r step5_input.gro -p topol.top -n index.ndx #!gmx mdrun -v -deffnm step6.6_equilibration_NPT_step4 In alternative: .. code:: ipython3 ##16## equilibration-NPT6.6-step4 > To transfer !cp ../reference/step6.6_equilibration_NPT_step4* . Here you can find the output of NPT steps: .. code:: ipython3 ##Equilibration - NPT outputs %ls ../reference/*NPT_step* **Evaluate the Pressure** ``Attention``! I hope you have not forgotten to evaluate the temperature in each step so far. In NPT equilibration steps, we need to check at extra parameter, ``pressure``: .. code:: ipython3 ## check the pressure > To run !echo 18|gmx energy -f step6.6_equilibration_NPT_step4.edr -o NPT_S4_Press.xvg -xvg none .. code:: ipython3 import pandas as pd import statistics df = pd.read_csv('NPT_S4_Press.xvg', sep='\s+', header=None, names=['time(ps)','pressure(bar)']) df.plot('time(ps)') .. code:: ipython3 ##Alternative##Plotting data by Xmgrace tool #!xmgrace NPT_S4_Press.xvg As you can see, pressure fluctuates more than temperature which is totally ok :) The system’s temperature should rise to the desired level (~303 K) and stay steady for the rest of the equilibration. While, pressure is the parameter which fluctuates a lot during the MD simulation, and its average should be statistically comparable with the pre-defined pressure (1 bar). Production ========== **Invoke mdrun for production (one single simulation)** The system is now properly equilibrated and reached the appropriate temperature and pressure. The position restrictions may now be lifted, to start the production MD and begin data gathering.For production run, we will first generate the ``binary tpr`` filr using the following mdp parameters: .. code:: ipython3 !cat step7_production_revised.mdp .. code:: ipython3 ##19.1## step7-production > To run ##generate tpr file #!gmx grompp -f step7_production_revised.mdp -o step7_production.tpr -c step6.6_equilibration_NPT_step4.gro -t step6.6_equilibration_NPT_step4.cpt -p topol.top -n index.ndx The completion of a production run (cell 19.2) might take hours or days. Dont execute this command here and try it on the terminal of the system you choose to perform the simulation. **!Warning:** This run will generate roughly 9028 Mb of data .. code:: ipython3 ##19.2## step7-production(one single simulation) > To run #!gmx mdrun -s step7_production -cpi A 100ns simulation was performed. Every 10000 of the 50000000 steps were recorded as the trajectory frames. The trr file is compressed into an xtc file. ======================================================================================================================================================== For convenience, we saved every 500 frames of the final trajectory and put the outputs in the reference folder. .. code:: ipython3 ##20## step7-production > To transfer !cp ../reference/step7_production* . .. code:: ipython3 ##step7-production outputs %ls *step7_production* **Trajectory Visualization** .. code:: ipython3 import nglview as ng import mdtraj as md trajectory = md.load("step7_production_traj_comp_skip500.xtc", top="step5_input.pdb") view = ng.show_mdtraj(trajectory) view.clear() view.add_cartoon('protein',color_scheme='residueindex') view Now, place the mouse cursor over the protein and press the play button. You will see the protein wiggling and jiggling! .. code:: ipython3 view2 = ng.show_mdtraj(trajectory) view2.add_cartoon('protein',color_scheme='residueindex') view2.add_spacefill('not protein', opacity=0.2) view2.center() view2.camera = 'orthographic' view2 In this Part, lipids (gray and red spheres), and ions (green and purpule spheres) are illustrated. You can find the name of the atoms by simply placing the cursor on it. Try to find chloride ions. -------------- **HINT** **How to generate replicas?** You may need more than one replica. You can generate several tpr files using the random seed option for the velocity and temperature in the ‘mdp’ file. **How to run several replicas?** You don’t need to perform simulations separately, because GROMACS has a simple solution for that: **step 1:** make folders (here we set four replicas) **#mkdir Rep1 Rep2 Rep3 Rep4** **step 2:** copy each tpr file in a different folder and execute mdrun by command below: **#mpirun -np 4 gmx_mpi mdrun -v -s step7_production.tpr -multidir Rep1 Rep2 Rep3 Rep4 -cpi** -------------- **Quiz** Why do we need to run several replicas? Do you know how they differ from one another? Exactly the same? Analysis ======== **RMSD analysis** .. code:: ipython3 !echo 4 4 |gmx rms -s step7_production.tpr -f step7_production_traj_comp_skip500.xtc -o production_rmsd.xvg -bin production.rmsd.dat -tu ns -fit rot+trans -xvg none .. code:: ipython3 import pandas as pd df = pd.read_csv('production_rmsd.xvg', sep='\s+', header=None, names=['time(ns)','RMSD(nm)']) df.plot('time(ns)') The root mean square deviation (RMSD) is a practical parameter to compare the backbones of a protein from its initial to final state, which illustrates the dynamics of structure during the simulation. Here, the protein structure is pretty stable and fluctuates within expected margins(1-2 Å). .. code:: ipython3 ##Alternative##Plotting data by Xmgrace tool #!xmgrace production_rmsd.xvg **Energy** .. code:: ipython3 ##check the energy !echo 13|gmx energy -f step7_production.edr -o production_Etot.xvg -xvg none .. code:: ipython3 import pandas as pd df = pd.read_csv('production_Etot.xvg', sep='\s+', header=None, names=['time(ps)','Energy(kJ mol-1)']) df.plot('time(ps)') .. code:: ipython3 ##Alternative##Plotting data by Xmgrace tool #!xmgrace production_Etot.xvg **Temperature** .. code:: ipython3 #check the temperature !echo 15|gmx energy -f step7_production.edr -o production_temp.xvg -xvg none .. code:: ipython3 import pandas as pd df = pd.read_csv('production_temp.xvg', sep='\s+', header=None, names=['time(ps)','Temperature*(K)']) df.plot('time(ps)') .. code:: ipython3 ##Alternative##Plotting data by Xmgrace tool #!xmgrace production_temp.xvg **Pressure** .. code:: ipython3 #check the pressure !echo 16|gmx energy -f step7_production.edr -o production_press.xvg -xvg none .. code:: ipython3 import pandas as pd df = pd.read_csv('production_press.xvg', sep='\s+', header=None, names=['time(ps)','Pressure(bar)']) df.plot('time(ps)') .. code:: ipython3 ##Alternative##Plotting data by Xmgrace tool #!xmgrace production_press.xvg Easy ha? Now try the tutorial on your protein and let us know if it is useful. if you have any question write in the forum and we will answer as soon as possible: http://forums.gromacs.org *References* 1. Yvonnesdotter, L., Rovšnik, U., Blau, C., Lycksell, M., Howard, R.J. and Lindahl, E., 2023. Automated simulation-based membrane-protein refinement into cryo-EM data. Biophysical Journal.