How to get the pressure components using "-rerun" command with a different itp file?

User discussions
1

I generated a trajectory of 300 ns of a system using bonded and non-bonded force filed parameters, and I got the pressure components (Pxy, Pyz,…) from the existing edr file. For example, I want to generate the pressure components obtained only from the bonds or only from Coulomb interactions by rerunning the existing xtc file with a different itp file. How can I get the edr file consisting of the pressure components using the new force filed parameters?

I have tried “-rerun” command and generated the edr file. But I am not able to get the new pressure components. Can anyone help me in this topic?

2

I think that in most cases it is not meaningful to look at a decomposition of pressure contributions, as several terms in force fields are balanced against each other.

Rerun can not compute anything to do with velocities. So we can not compute kinetic energies and the constraint virial. If you are interested in only the force component of the virial, without the constraint contribution, the fix below force the 2023 release will get you that. Apply with git am

commit 56a2ed53b882798c9aae1abb667d1d469e4268bc (HEAD -> bh-rerun-add-force-virial)
Author: Berk Hess <hess@kth.se>
Date:   Fri Nov 3 09:58:34 2023 +0100

    Add force virial output to rerun

diff --git a/src/gromacs/mdlib/energyoutput.cpp b/src/gromacs/mdlib/energyoutput.cpp
index 6fbd4238b0..965532ecf9 100644
--- a/src/gromacs/mdlib/energyoutput.cpp
+++ b/src/gromacs/mdlib/energyoutput.cpp
@@ -320,10 +320,10 @@ EnergyOutput::EnergyOutput(ener_file*                fp_ene,
             ienthalpy_ = get_ebin_space(ebin_, 1, &enthalpyEnergyFieldName, unit_energy);
         }
     }
+    ivir_ = get_ebin_space(
+            ebin_, virialEnergyFieldNames.size(), virialEnergyFieldNames.data(), unit_energy);
     if (bPres_)
     {
-        ivir_ = get_ebin_space(
-                ebin_, virialEnergyFieldNames.size(), virialEnergyFieldNames.data(), unit_energy);
         ipres_  = get_ebin_space(ebin_, asize(pres_nm), pres_nm, unit_pres_bar);
         isurft_ = get_ebin_space(ebin_, asize(surft_nm), surft_nm, unit_surft_bar);
     }
@@ -911,9 +911,9 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
             add_ebin(ebin_, ienthalpy_, 1, &enthalpy, bSum);
         }
     }
+    add_ebin(ebin_, ivir_, 9, vir[0], bSum);
     if (bPres_)
     {
-        add_ebin(ebin_, ivir_, 9, vir[0], bSum);
         add_ebin(ebin_, ipres_, 9, pres[0], bSum);
         tmp = (pres[ZZ][ZZ] - (pres[XX][XX] + pres[YY][YY]) * 0.5) * box[ZZ][ZZ];
         add_ebin(ebin_, isurft_, 1, &tmp, bSum);
diff --git a/src/gromacs/mdlib/stat.cpp b/src/gromacs/mdlib/stat.cpp
index fc67a089ab..154288a4f8 100644
--- a/src/gromacs/mdlib/stat.cpp
+++ b/src/gromacs/mdlib/stat.cpp
@@ -161,6 +161,7 @@ void global_stat(const gmx_global_stat&   gs,
     bool bTemp       = ((flags & CGLO_TEMPERATURE) != 0);
     bool bEner       = ((flags & CGLO_ENERGY) != 0);
     bool bPres       = ((flags & CGLO_PRESSURE) != 0);
+    bool bForceVir   = bEner;
     bool bConstrVir  = ((flags & CGLO_CONSTRAINT) != 0);
     bool bEkinAveVel = (inputrec.eI == IntegrationAlgorithm::VV
                         || (inputrec.eI == IntegrationAlgorithm::VVAK && bPres));
@@ -230,7 +231,7 @@ void global_stat(const gmx_global_stat&   gs,
         }
     }
 
-    if (bPres)
+    if (bForceVir)
     {
         ifv = add_binr(rb, DIM * DIM, fvir[0]);
     }
@@ -327,7 +328,7 @@ void global_stat(const gmx_global_stat&   gs,
             }
         }
     }
-    if (bPres)
+    if (bForceVir)
     {
         extract_binr(rb, ifv, DIM * DIM, fvir[0]);
     }
diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp
index 9f56f1f5e5..c35d679047 100644
--- a/src/gromacs/mdrun/rerun.cpp
+++ b/src/gromacs/mdrun/rerun.cpp
@@ -778,7 +778,7 @@ void gmx::LegacySimulator::do_rerun()
                                                                 state->nhpres_xi,
                                                                 state->nhpres_vxi }),
                                              state->fep_state,
-                                             total_vir,
+                                             force_vir, // We have only the force component of the virial
                                              pres,
                                              ekind,
                                              mu_tot,
1 Like