MdEngFlow

Example demonstrating a free energy estimate for comparing the mutants derived from LigandMPNN.

This is a simplified version of the workflow used in the project, but should give you a good starting point. Including the MD simulation and analysis steps batchly or parallelly.

installation

requirements

• conda or microconda
• @TODO: add env.yaml file

usage

file prep:

see files in wdir dir.

• provide the PDB file that pre-equilibrates
  • where the ligand pose has been determined,like confirming from the quantum chemistry
  • the water molecules should be relaxed fully
• ff dir
  • prepare the force field files for the tleap script
  • rewrite the tleap script inclued in the /MdEngFlow/script/generate_md_system.py

     tmp_leap = """
     source oldff/leaprc.ff14SB
     source leaprc.gaff2
     source leaprc.water.tip3p
     addAtomTypes {{                                 
             {{ "M1"  "Fe" "sp3"}}
             {{ "Y1"  "S" "sp3" }}
             {{ "Y2"  "N" "sp3" }}
             {{ "Y3"  "N" "sp3" }}
             {{ "Y4"  "N" "sp3" }}
             {{ "Y5"  "N" "sp3" }}
             {{ "Y6"  "O" "sp3" }}
     }}
     CM1 = loadmol2 {ffdir}/CM1.mol2
     HM1 = loadmol2 {ffdir}/HM1.mol2
     FE1 = loadmol2 {ffdir}/FE1.mol2
     OE1 = loadmol2 {ffdir}/OE1.mol2
     DOK = loadmol2 {ffdir}/DOK.mol2
     loadamberparams {ffdir}/HEM.frcmod
     loadamberparams frcmod.ions1lm_126_tip3p
     loadamberparams frcmod.ions234lm_126_tip3p
     loadamberparams {ffdir}/CpdI_OFE_mcpbpy.frcmod
     loadamberparams {ffdir}/DOK.frcmod
    
     mol = loadpdb {pdbfile}
    
     bond mol.434.SG mol.498.FE
     bond mol.433.C mol.434.N
     bond mol.434.C mol.435.N
     bond mol.497.NA mol.498.FE
     bond mol.497.NB mol.498.FE
     bond mol.497.NC mol.498.FE
     bond mol.497.ND mol.498.FE
     bond mol.498.FE mol.499.O1
    
     # solvatebox mol TIP3PBOX 10.0 0.80
     # addions mol Na+ 0
     # addions mol Cl- 0
     setbox mol vdw
     savepdb mol {tmpprefix}.pdb
     saveamberparm mol {tmpprefix}.prmtop {tmpprefix}.inpcrd
     quit
     """

• config.yaml This file contains the parameters for the MD simulation. You can adjust these as needed.
  • mutations
  • list format:

      list: 
        [ 
        ['T1A_P2A_P3A'],    # residue 1 Thr mutation to Ala, residue 2 Pro mutation to Ala, residue 3 Pro mutation to Ala      
        ['T1A'],            # residue 1 Thr mutation to Ala
        ['P2A'],            # residue 2 Pro mutation to Ala
        ]

  • @TODO: add more explanation for the other parameters in the config.yaml file.

commands

• MD + MMPBSA

cd widr
CUDA_VISIBLE_DEVICES=0 python ../engine config.yaml > run.log 2>&1

It will generate xxx.log files and output dir.

• only MD

cd widr
CUDA_VISIBLE_DEVICES=0 python ../engine config.yaml --skip_postprocess_mm > run.log 2>&1

• Only the energy calculation, in the outputdirectory which has the md trajectory, if u want to run the binding affinity calculation (MMPBSA/MMGBSA)

cd widr
CUDA_VISIBLE_DEVICES=0 python ../engine config.yaml --skip-simulation > run.log 2>&1

results

.
├── config.yaml
├── demo.pdb
├── ff
│   ├── CM1.mol2
│   ├── CpdI_OFE_mcpbpy.frcmod
│   ├── DOK.frcmod
│   ├── DOK.mol2
│   ├── FE1.mol2
│   ├── HEM.frcmod
│   ├── HM1.mol2
│   └── OE1.mol2
├── leap.log
├── output
│   └── T1A_P2A_P3A
│       ├── md
│       │   ├── colvars.colvar
│       │   ├── finished_step_list
│       │   ├── index.ndx
│       │   ├── input.gro
│       │   ├── input.pdb
│       │   ├── input.top
│       │   ├── mdout.mdp
│       │   ├── mdrun.sh
│       │   ├── min.colvars_prev.state
│       │   ├── min.colvars.state
│       │   ├── min.colvars.traj
│       │   ├── min.edr
│       │   ├── min.gro
│       │   ├── min.log
│       │   ├── min.mdp
│       │   ├── min.tpr
│       │   ├── pdt1.edr
│       │   ├── pdt1.log
│       │   ├── pdt1.mdp
│       │   ├── pdt1.tpr
│       │   ├── pdt1.xtc
│       │   └── posre.itp
│       └── prep
│           ├── T1A_P2A_P3A.amb2gmx
│           │   ├── acpype.log
│           │   ├── em.mdp
│           │   ├── md.mdp
│           │   ├── posre_T1A_P2A_P3A.itp
│           │   ├── rungmx.sh
│           │   ├── T1A_P2A_P3A_GMX.gro
│           │   └── T1A_P2A_P3A_GMX.top
│           ├── T1A_P2A_P3A.gro
│           ├── T1A_P2A_P3A.in
│           ├── T1A_P2A_P3A.inpcrd
│           ├── T1A_P2A_P3A.in.result
│           ├── T1A_P2A_P3A.pdb
│           ├── T1A_P2A_P3A.prmtop
│           └── T1A_P2A_P3A.top
└── run.log

tidbits