I’m running a sequential (staged) multi-mutation alchemical calculation with pmx + GROMACS: I mutate one residue per stage (e.g. T2V), equilibrate both end states, then use the state-B equilibrated structure as the starting point for the next stage’s mutation, and so on through 16 stages.
Problem: starting from stage 2, GROMACS/pmx appears to treat all previously mutated residues (still named T2V in the structure/topology) as active dual-topology (A/B) mutations again, not just the current stage’s target. I confirmed this by checking [ atoms ] in each stage’s topol.top — every prior T2V residue still carries distinct typeB/chargeB entries, not just the new one. As supporting evidence, the peak dH/dλ magnitude for forward window 1 scales roughly linearly with stage number (~2,000 kJ/mol·λ⁻¹ at stage 1 up to ~22,000 at stage 16), consistent with an increasing number of residues being simultaneously coupled to λ rather than just one per stage.
I believe this happens because gmx pdb2gmx -ff charmm36m-mut (and pmx gentop) attach B-state parameters to any residue named T2V, with no way to distinguish “this stage’s new target” from “an already-completed prior mutation.”
Question: What’s the correct/recommended way to “freeze” a previously mutated residue as a plain, single-topology amino acid (e.g. VAL) before starting the next stage’s mutation, so it’s no longer treated as dual-topology? I want to avoid simply renaming T2V→VAL in the coordinate file, since the state-B structure still contains the extra (now non-interacting) dummy atoms unique to the original residue (e.g. THR’s OG1/HG1), which wouldn’t match VAL’s .rtp entry.
Is there a pmx utility/workflow for this, or is the expected approach to manually strip the dummy atoms and regenerate the topology for that residue via the plain (non-mut) force field before proceeding to the next stage? Section 4.3 of Aldeghi et al. (Methods Mol Biol 2019) discusses reusing equilibrium trajectories across a mutation scan, but I haven’t found guidance for this specific chained/sequential case.
you’re right. Every time you run gmx pdb2gmx -ff charmm36m-mut (and pmx gentop) a new hybrid structure + topology is added, and so subsequent stages in your workflow will define more and more morphs.
I am guessing that you are doing the mutations sequentially rather than all at once to give the system more time to equilibrate after each individual mutation and therefore get a more converged result? I can see the reasoning behind that, but please also keep in mind that the overall error propagates over all 16 individual steps, which might become substantial as well.
So if that were the rationale, I would consider as alternative to do all mutations at once. This is a massive perturbation, so you may need long non-equilibrium transitions (of perhaps 10 ns or more) to get sufficient overlap between forward and backward runs, but if it converges you do not need to worry about reverting mutated residues back to normal residues in the structures and topologies, or about error propagation across the 16 stages.
I tried running each window for 15 ns but the error is still coming high
ANALYSIS
========================================================
Number of forward (0->1) trajectories: 26
Number of reverse (1->0) trajectories: 25
Temperature : 300.00 K
Running KS-test…
Forward: gaussian quality = 0.92
—> KS-Test Ok
Reverse: gaussian quality = 0.86
—> KS-Test Ok
Bennett Acceptance Ratio
Running Nelder-Mead Simplex algorithm…
Bootstrap (Std Err): iteration 100/100
Bootstrap (Conv): iteration 100/100
BAR: dG = 939.00 kJ/mol
BAR: Std Err (analytical) = 359.98 kJ/mol
BAR: Std Err (bootstrap) = 2.29 kJ/mol
BAR: Conv = 1.00
BAR: Conv Std Err (bootstrap) = 0.00
you’re introducing a massive perturbation by running a 16-fold mutation. This is probably the best you can expect. You could try with even longer transitions, but given the size of the perturbation, you probably have to live a certain statistical uncertainty in the obtained free energy value.
I suggest to do mutations sequentially: this way you can inspect which mutation is showing lack of convergence and increase transition times for that mutation only.
To perform sequential mutations, you cannot use stateB of the previous mutation for the next unrelated mutant. You need to prepare vanilla (non-hybrid) structures for all your mutants of interest, e.g. 1) WT, 2) Mut1, 3) Mut1_Mut2, 4) Mut1_Mut2_Mut3, … You can introduce mutations with Pymol, Rosetta, FoldX or any other tool built for this purpose. Afterwards, you can do alchemical perturbations of one mutation at a time: i) WT ↔ Mut1; (ii) Mut1 ↔ Mut1_Mut2; (iii) Mut1_Mut2 ↔ Mut1_Mut2_Mut3, …
