1. Preparations: .mdp & .gro & .top

By Zijian Download gromacs ``` wget ftp://ftp.gromacs.org/pub/gromacs/gromacs-2020-beta1.tar.gz ``` and install it Clone Bussi's plumed2 repo ``` git clone https://github.com/GiovanniBussi/plumed2 -b v2.1-hrex ``` and install it ###Compare partial_tempering with scale=1.0, scale=0.5 to non-scaled force field Produce un-scaled `tpr` file ``` gmx grompp -f grompp.mdp -c conf.gro -p topol.top -o topol-unscaled.tpr -maxwarn 10 ``` Produce `processed.top` ``` gmx grompp -f grompp.mdp -c conf.gro -p topol.top -pp -maxwarn 10 ``` Choose the "hot" atoms appending "_" in `processed.top`, e.g. change ``` [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls_138 1 ALAB CB 1 0.000 12.011 2 opls_140 1 ALAB HB1 2 0.000 1.008 3 opls_140 1 ALAB HB2 3 0.000 1.008 4 opls_140 1 ALAB HB3 4 0.000 1.008 5 opls_140 1 ALAB HB4 5 0.000 1.008 ``` to ``` [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls_138_ 1 ALAB CB 1 0.000 12.011 2 opls_140_ 1 ALAB HB1 2 0.000 1.008 3 opls_140_ 1 ALAB HB2 3 0.000 1.008 4 opls_140_ 1 ALAB HB3 4 0.000 1.008 5 opls_140_ 1 ALAB HB4 5 0.000 1.008 ``` Process `top` file with plumed ``` plumed partial_tempering 1.0 < processed.top > topol-1.0.top ``` ``` plumed partial_tempering 0.5 < processed.top > topol-0.5.top ``` Produce scaled `tpr` files ``` gmx grompp -f grompp.mdp -c conf.gro -p topol-1.0.top -o topol-1.0.tpr -maxwarn 10 ``` ``` gmx grompp -f grompp.mdp -c conf.gro -p topol-0.5.top -o topol-0.5.tpr -maxwarn 10 ``` Run MD simulations ``` gmx mdrun -s topol-unscaled.tpr -rerun rerun.trr -g md-unscaled.log ``` ``` gmx mdrun -s topol-1.0.tpr -rerun rerun.trr -g md-1.0.log ``` ``` gmx mdrun -s topol-0.5.tpr -rerun rerun.trr -g md-0.5.log ``` Compare `md-unscaled.log`, `md-1.0.log` and `md-0.5.log`. By Yuzhi: # 1. Preparations: .mdp & .gro & .top In a typical free-energy calculation, the files .gro & .top are often the same. (Of course need time to generate) The .mdp files differ in only one line: ``` init-lambda-state = X ``` The .tpr files can be generated by `grompp` For halmitonian replica exchange (HREX) methods, it's necesary to perform "scaling" on force field parameters. The `partial_tempering` in `plumed` can help. # 2. MD (Plumed + Gromacs) plumed2.1-hrex (Bussi) + groamac 4.6.7 are needed. See http://bbs.keinsci.com/thread-760-1-1.html. and http://sobereva.com/247 for installation guides. Here I take direct ethanol solvation free energy as an example. This follows directly from http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy (well designed website and tutorial) The MD process can run in such ways: ``` mdrun -nt 1 -deffnm ethanol.X -dhdl ethanol.X.dhdl.xvg ``` The `*.dhdl.xvg` files are all we need for further analysis. HREX can be run by ``` mpiexec -np $NP mdrun_mpi_d -defnm -multi 9 -replex 1000 -nex 100000 -deffnm ethanol. -dhdl ethanol.dhdl. ``` If GROMACS has been patched with PLUMED it should accept the -hrex option in mdrun. By typing `gmx_mpi mdrun -h`, one can verify whether he/she can perform HREX. Notice that HREX requires for a `mpi` verison of gromacs. The `*.dhdl.xvg` are also all we need (containing the MD outputs) ## 3. Analysis (MBar) pymbar 0.3.0 git clone https://github.com/choderalab/pymbar.git alchemistry: git clone https://github.com/choderalab/pymbar-examples.git (The latest version has changes into https://github.com/MobleyLab/alchemical-analysis (alchemlyb) compatible with Python3) In this part, I directly follow the example provided by pymbar, to compute the hydration free energy of 3-methylindole using a beta version of Gromacs 4.6. Actually it's nearly the same with the above ethanol example. Given the files `*.dhdl.xvg` obtained from the MD step and stored in `data/3-methylindole-11steps`, the script `alchemical_analysis.py` can used to analysis: ``` python alchemical_analysis.py -d data/3-methylindole-11steps -q xvg -p dhdl -u kJ ``` The free energy change can be refered to https://github.com/choderalab/pymbar-examples/blob/master/alchemical-free-energy/gromacs/output_11steps/dF_state.pdf The total screen outputs are https://github.com/choderalab/pymbar-examples/blob/master/alchemical-free-energy/gromacs/output_11steps/screen_printout.txt It supports for these six methods: TI, TI-CUBIC, DEXP, IEXP, BAR, MBAR. DEXP and IEXP are exponential averaging in both directions, insertion and deletion. To do-list: 1. this is an example for Relative Binding Affinity http://www.alchemistry.org/wiki/Example:_Relative_Binding_Affinity 2. Benchmarks and more complicated test cases. 3. The verison of gromacs? (4.6.7 is too old) 4. Further reading and discussions.