I am trying to use pymol and open babel in a pipeline to eventually perform docking of protein and ligands in a complex. Currently, I have two methods: (1)Open Babel only and (2) PyMOL+Open Babel.

In (1), the pipeline only uses open babel to eventually output the PDBQT files (as is needed for AutoDock Vina).

In (2), the pipeline uses PyMOL to generate the separated PDB files (of the protein and ligand) and then subsequently Open Babel to generate the PDBQT files.

Initially, I was using (1) and works without any issues until I came across PDB files with missing elements in the "element" column, this caused issues for Open Babel especially when it came to sanitising the molecules which is why I incorporated PyMOL to interpret the structure files in an attempt to "fix" the PDB file prior to generating the protein and ligand PDBs (which would be put through the Vina pre-processing steps without issues).

I noticed the outputs for (1) and (2) generated different ligand protonated states when I inspected the final PDBQT files. I have included screenshots of the ligand outputs for both methods.

Ligand output from (1) Ligand output from (2)

The code and processing steps are completely identical for both (1) and (2) where hydrogens are added for pH 7.4 (using Open Babel) and gasteiger charges are added (using Open Babel). The only difference was using PyMOL in the beginning of (2) to produce the separated files of the Protein.pdb and Ligand.pdb from the Complex.pdb.

I am running PyMol through an API on my linux cluster via: I am not sure if it is an Open Babel issue or a PyMol issue.

Versions being used: Open Babel: 3.1.0 PyMOL: 3.1.3.1 Incentive Product (698d66e4c6), 2025-02-05

    ##########################################################################
# Function: Extract protein and ligand from a complex using PyMOL.
#   - The protein (all atoms not matching resn LIGAND) is saved as protein.pdb.
#   - The ligand (atoms with resn LIGAND) is saved as ${LIGAND}.mol2.
##########################################################################
extract_ligand_as_mol2() {
  local dir="$1"
  echo "Extracting protein and ligand from complex in $dir using PyMOL..."
  
  # Find the first complex PDB file in the directory.
  local pdb_file
  pdb_file=$(find "$dir" -maxdepth 1 -type f -name "*.pdb" | head -n 1)
  if [[ -z "$pdb_file" ]]; then
    echo "No pdb file found in $dir"
    return 1
  fi
  
  # Define output files.
  local protein_file="$dir/protein.pdb"
  local ligand_mol2_file="$dir/${LIGAND}.mol2"
  
  # Create a temporary PyMOL script.
  local tmp_pymol_script
  tmp_pymol_script=$(mktemp /tmp/pymol_script.XXXXXX.py)
  
  cat <<EOF > "$tmp_pymol_script"
from pymol import cmd
cmd.reinitialize()
cmd.load("$pdb_file", "complex")
# Select the protein: everything NOT with the residue name matching LIGAND.
cmd.select("protein", "not resn $LIGAND")
cmd.save("$protein_file", "protein")
# Select the ligand: only atoms with residue name matching LIGAND.
cmd.select("ligand", "resn $LIGAND")
cmd.save("$ligand_mol2_file", "ligand", format="mol2")
cmd.quit()
EOF