How to generate a MOL2 ligand file with the correct docking pose, full hydrogens, and physiological

User discussions
1

GROMACS version: 2026.2
GROMACS modification: No
Hello everyone,

I am currently preparing a molecular dynamics simulation for a protein–ligand complex. The ligand pose was obtained from molecular docking using AutoDock Vina, and the coordinates of the ligand in the binding pocket seem to be correct.

However, I have encountered a problem during ligand preparation for MD. The ligand file extracted from the docking pose is in PDB/PDBQT format, but the atom names are very generic, for example C, C, C, O, O, instead of unique names such as C1, C2, C3, O1, O2, O3. In addition, the ligand does not contain all hydrogens because the non-polar hydrogens were merged during docking preparation.

I would like to run the MD simulation under physiological conditions, around pH 7.4. Therefore, I need to generate a MOL2 file that satisfies the following conditions:

  1. It preserves the correct docking pose/conformation obtained from AutoDock Vina.
  2. It contains all hydrogens, including non-polar hydrogens.
  3. It has the correct protonation state at physiological pH.
  4. It has atom names that are consistent with the topology files used in GROMACS.

I also have the original MOL2 file of the ligand with the correct chemical structure. I used this MOL2 file to generate the ligand parameter files, including the .str and .itp files. However, when I try to merge the ligand with the receptor to build the complex in GROMACS, I get an error indicating that the number of hydrogen atoms does not match. On the other hand, if I use the ini.pdb file generated during the topology preparation, the ligand no longer has the correct docking pose.

Could anyone suggest a proper workflow to generate a ligand MOL2/PDB file that keeps the docking pose but also has the correct full-hydrogen structure and matches the generated topology?

Any advice would be greatly appreciated.

Thank you very much.

2

Hi, I did something similar recently, although maybe this is not exactly what you want.

I extracted the pdb with vina using

vina_split --input ${lig_name}_${input_file}_docking_results.pdbqt --ligand ${lig_name}_${input_file}_docking_results_output > ${lig_name}_${input_file}_docking_results_output.log
for file in ${lig_name}_${input_file}_docking_results_output*.pdbqt; do obabel "$file" -O "${file%.pdbqt}.pdb" -p 7.5 -xhn; done

So I copied out the top x docking poses, and converted them to .pdb using obabel
Then I uploaded the charmm-gui to get .gro formatted itps with unique ligands (i may of converted them again to .sdf to make cgen work better)

Then I used gromolgist to combine these with my protein of interest.
eg.
python ../../gromologist_ligand_receptor_prep.py -f {lig1} -c {protein}

Here is the python script, should be sefl explanatory but lmk if you need any help.

Gromology Automation Script (gromology_script.py)

## Gromology Automation Script (gromology_script.py)
import argparse
import gromologist as gml

def main():
    # Set up argument parser
    parser = argparse.ArgumentParser(description='Automate molecular topology and structure preparation')
    parser.add_argument('-f', '--molecule', required=True, help='Molecule name (e.g., AMX)')
    parser.add_argument('-c', '--input_file', required=True, help='PDB ID (e.g., 6i1f)')
    
    args = parser.parse_args()
    lig_name = args.molecule
    input_file = args.input_file

    # Generate filenames using f-strings
    top_file = f"{lig_name}_{input_file}_docked.top"
    output_pdb = f"{lig_name}_{input_file}_docked.pdb"
    docking_results = f"{lig_name}_{input_file}_docking_results_output-charmm.pdb"
    prot_pdb = f"{input_file}_prod_prot_only.pdb"
    charrmm_itp = f"{lig_name}_{input_file}_charmm36.itp"
    lig_itp = f"{lig_name}_{input_file}.itp"

    try:
        # Process topology file
        print(f"Processing topology file: {top_file}")
        t = gml.Top(top_file)
        t.add_parameters_from_file(charrmm_itp, prefix_type='Y')
        t.add_molecule_from_file(lig_itp, molcount=1, prefix_type='Y')
        t.save_top(top_file)

        # Process PDB file
        print(f"Processing structure file: {prot_pdb}")
        p = gml.Pdb(prot_pdb)
        p.insert_from(docking_results)
        p.save_pdb(output_pdb)

        print("\nProcess completed successfully!")
        print(f"Generated files: {top_file}, {output_pdb}")

    except FileNotFoundError as e:
        print(f"Error: Required file not found - {e.filename}")
    except Exception as e:
        print(f"An error occurred: {str(e)}")

if __name__ == "__main__":
    main()

This should keep all the hydrogens in your ligand intact.