From 76f3e75093c890e2f3760316c2cf819929530717 Mon Sep 17 00:00:00 2001
From: Amy He <ayhe1768@gmail.com>
Date: Wed, 6 Nov 2024 19:30:29 -0800
Subject: [PATCH 01/18] remove repeated meeko installation from env setup in
 tutorials; remove TOC from tutorials for heading formatting

---
 docs/source/tutorial1.rst | 11 -----------
 docs/source/tutorial2.rst | 11 -----------
 docs/source/tutorial3.rst |  4 ----
 3 files changed, 26 deletions(-)

diff --git a/docs/source/tutorial1.rst b/docs/source/tutorial1.rst
index 50acc810..019dd963 100644
--- a/docs/source/tutorial1.rst
+++ b/docs/source/tutorial1.rst
@@ -5,10 +5,6 @@ Basic Docking
 
 This tutorial provides practice examples and a step-by-step guide for the two basic procedures, **Ligand Preparation** and **Receptor Preparation**, with Meeko for molecular docking and virtual screening with `AutoDock Vina <https://github.com/ccsb-scripps/AutoDock-Vina>`_ and `AutoDock-GPU <https://github.com/ccsb-scripps/AutoDock-GPU>`_. It is based on, but not a full version of the tutorial materials in `Forlilab tutorials <https://github.com/forlilab/tutorials>`_. 
 
-.. contents::
-   :local:
-   :depth: 2
-
 Prerequisites and Environment Setup
 ===================================
 
@@ -32,13 +28,6 @@ Install the required Python packages through ``conda-forge``
 Install the additional packages and data from GitHub repositories
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-- (Python package) Meeko 
-
-.. code-block:: bash
-
-   git clone --single-branch --branch develop https://github.com/forlilab/Meeko.git
-   cd Meeko; pip install --use-pep517 -e .; cd ..
-
 - (Python package) scrubber 
 
 .. code-block:: bash
diff --git a/docs/source/tutorial2.rst b/docs/source/tutorial2.rst
index b65e6e7a..27739bb5 100644
--- a/docs/source/tutorial2.rst
+++ b/docs/source/tutorial2.rst
@@ -8,10 +8,6 @@ This is a reactive docking example that uses the AutoDock-GPU executable to gene
 
 Follow the instructions to set up the environment and run this command-line example on your own device (Linux, MacOS or WSL). To run this example in a Colab notebook, see :ref:`colab_examples`. 
 
-.. contents::
-   :local:
-   :depth: 2
-
 Introduction
 ============
 
@@ -44,13 +40,6 @@ Install the required Python packages through ``conda-forge``
 Install the additional packages and data from GitHub repositories
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-- (Python package) Meeko 
-
-.. code-block:: bash
-
-   git clone --single-branch --branch develop https://github.com/forlilab/Meeko.git
-   cd Meeko; pip install --use-pep517 -e .; cd ..
-
 - (Python package) scrubber 
 
 .. code-block:: bash
diff --git a/docs/source/tutorial3.rst b/docs/source/tutorial3.rst
index ce4c602e..e2c0a1e5 100644
--- a/docs/source/tutorial3.rst
+++ b/docs/source/tutorial3.rst
@@ -8,10 +8,6 @@ This is a tethered (two-point attached covalent) docking example that uses the A
 
 Follow the instructions to set up the environment and run this command-line example on your own device (Linux, MacOS or WSL). To run this example in a Colab notebook, see :ref:`colab_examples`. 
 
-.. contents::
-   :local:
-   :depth: 2
-
 Introduction
 ============
 

From cf2129864fe01059822910e4cb2d5c9dccd76719 Mon Sep 17 00:00:00 2001
From: Amy He <ayhe1768@gmail.com>
Date: Wed, 6 Nov 2024 20:56:58 -0800
Subject: [PATCH 02/18] add usage to mk_export.py

---
 docs/source/cli_export_result.rst | 63 +++++++++++++++++++++++++++----
 docs/source/tutorial2.rst         |  6 +--
 docs/source/tutorial3.rst         |  2 +-
 3 files changed, 60 insertions(+), 11 deletions(-)

diff --git a/docs/source/cli_export_result.rst b/docs/source/cli_export_result.rst
index 641fbd21..eb0847b2 100644
--- a/docs/source/cli_export_result.rst
+++ b/docs/source/cli_export_result.rst
@@ -1,8 +1,11 @@
 mk_export.py
 ============
 
+About
+-----
+
 Convert docking results to SDF
-------------------------------
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 AutoDock-GPU and Vina write docking results in the PDBQT format. The DLG output
 from AutoDock-GPU contains docked poses in PDBQT blocks, plus additional information.
@@ -13,12 +16,12 @@ from RDKit molecules.
 
 .. code-block:: bash
 
-    mk_export.py molecule.pdbqt -o molecule.sdf
-    mk_export.py vina_results.pdbqt -o vina_results.sdf
-    mk_export.py autodock-gpu_results.dlg -o autodock-gpu_results.sdf
+    mk_export.py molecule.pdbqt -s molecule.sdf
+    mk_export.py vina_results.pdbqt -s vina_results.sdf
+    mk_export.py autodock-gpu_results.dlg -s autodock-gpu_results.sdf
 
 Why this matters
-----------------
+~~~~~~~~~~~~~~~~
 
 Making RDKit molecules from SMILES is safer than guessing bond orders
 from the coordinates, specially because the PDBQT lacks hydrogens bonded
@@ -31,9 +34,8 @@ but because this is a nearly impossible task.
     obabel -:"C1C=CCO1" -o pdbqt --gen3d | obabel -i pdbqt -o smi
     [C]1=[C][C]=[C]O1
 
-
 Caveats
--------
+~~~~~~~
 
 If docking does not use explicit Hs, which it often does not, the
 exported positions of hydrogens are calculated from RDKit. This can
@@ -41,3 +43,50 @@ be annoying if a careful forcefield minimization is employed before
 docking, as probably rigorous Hs positions will be replaced by the
 RDKit geometry rules, which are empirical and much simpler than most
 force fields.
+
+Usage
+-----
+
+.. code-block:: bash
+
+    mk_export.py [OPTIONS] docking_results_filename(s)
+
+Positional Argument
+~~~~~~~~~~~~~~~~~~~
+
+.. option:: docking_results_filename
+
+   One or more docking output files in either PDBQT format (from Vina) or DLG format (from AD-GPU).
+
+Options
+~~~~~~~
+
+.. option:: -s, --write_sdf <output_SDF_filename>
+
+   Specify the output SDF filename. Defaults to the input filename with a suffix defined by ``--suffix``.
+
+.. option:: -p, --write_pdb <output_PDB_filename>
+
+   Specify the output PDB filename. Defaults to the input filename with a suffix defined by ``--suffix``.
+
+.. option:: --suffix <suffix>
+
+   Set a suffix for output filenames that are not explicitly specified. The default suffix is ``_docked``.
+
+.. option:: -j, --read_json <input_JSON_filename>
+
+   Provide a receptor file generated by ``mk_prepare_receptor`` with the ``-j/--write_json`` option.
+
+.. option:: --all_dlg_poses
+
+   (Flag) Write all poses from AutoDock-GPU DLG output files, instead of only the lead of each cluster.
+
+.. option:: -k, --keep_flexres_sdf
+
+   (Flag) Include flexible residues, if any, in the SDF output.
+
+.. option:: -, --redirect_stdout
+
+   (Flag) Instead of writing an SDF file, print it directly to the standard output (STDOUT).
+
+
diff --git a/docs/source/tutorial2.rst b/docs/source/tutorial2.rst
index 27739bb5..59be7092 100644
--- a/docs/source/tutorial2.rst
+++ b/docs/source/tutorial2.rst
@@ -282,9 +282,9 @@ And to run the docking calculation, the ligand PDBQT file (``AMP.pdbqt``), the f
 
 If you're running these calculations on Google T4 backends, here are the pre-compiled executables of autogrid4 and adgpu specifically for T4: 
 
-- autodock-gpu v1.5.3
-`autodock_gpu_128wi <https://github.com/rwxayheee/Colabs/blob/acd2972f4afbf8c5299ebf0686534f466bf6f81b/Compiled_for_Colab/AutoDock-GPU_v1.5.3/autodock_gpu_128wi>`_
-`adgpu_analysis <https://github.com/rwxayheee/Colabs/blob/acd2972f4afbf8c5299ebf0686534f466bf6f81b/Compiled_for_Colab/AutoDock-GPU_v1.5.3/adgpu_analysis>`_
+- autodock-gpu v1.5.3-2e658c3
+`autodock_gpu_128wi <https://github.com/rwxayheee/Colabs/blob/b18a4c61f68647838ac987aacec8b66b0ab3426c/Compiled_for_Colab/AutoDock-GPU_v1.5.3_develop/autodock_gpu_128wi>`_
+`adgpu_analysis <https://github.com/rwxayheee/Colabs/blob/b18a4c61f68647838ac987aacec8b66b0ab3426c/Compiled_for_Colab/AutoDock-GPU_v1.5.3_develop/adgpu_analysis>`_
 
 - autogrid v4.2.6
 `autogrid4 <https://github.com/rwxayheee/Colabs/blob/acd2972f4afbf8c5299ebf0686534f466bf6f81b/Compiled_for_Colab/AutoGird_v4.2.6/autogrid4>`_
diff --git a/docs/source/tutorial3.rst b/docs/source/tutorial3.rst
index e2c0a1e5..c830324c 100644
--- a/docs/source/tutorial3.rst
+++ b/docs/source/tutorial3.rst
@@ -193,4 +193,4 @@ It is also possible to export the docking poses to a multi-model PDB file with u
    --default_altloc A -f $flexres \
    --box_enveloping "LIG.pdb" --padding 8.0 
 
-   mk_export.py HIE_AMP.dlg -s 3kgd_HIE_AMP_adgpu_out.sdf -k 3kgd_receptorH.json -p 3kgd_HIE_AMP_adgpu_out.pdb
\ No newline at end of file
+   mk_export.py HIE_AMP.dlg -s 3kgd_HIE_AMP_adgpu_out.sdf -j 3kgd_receptorH.json -p 3kgd_HIE_AMP_adgpu_out.pdb
\ No newline at end of file

From a731b350eef899a3a84877dd49728c277e1d2380 Mon Sep 17 00:00:00 2001
From: Amy He <ayhe1768@gmail.com>
Date: Wed, 6 Nov 2024 21:49:02 -0800
Subject: [PATCH 03/18] list options in cli docs

---
 docs/source/cli_export_result.rst |  18 +--
 docs/source/cli_lig_prep.rst      |  85 ++++++++++++-
 docs/source/cli_rec_prep.rst      | 193 ++++++++++++++++++++++++++++--
 3 files changed, 275 insertions(+), 21 deletions(-)

diff --git a/docs/source/cli_export_result.rst b/docs/source/cli_export_result.rst
index eb0847b2..1c35736f 100644
--- a/docs/source/cli_export_result.rst
+++ b/docs/source/cli_export_result.rst
@@ -61,25 +61,25 @@ Positional Argument
 Options
 ~~~~~~~
 
-.. option:: -s, --write_sdf <output_SDF_filename>
+.. option:: --suffix <suffix>
 
-   Specify the output SDF filename. Defaults to the input filename with a suffix defined by ``--suffix``.
+   Set a suffix for output filenames that are not explicitly specified. The default suffix is ``_docked``. 
 
-.. option:: -p, --write_pdb <output_PDB_filename>
+.. option:: -s, --write_sdf <output_SDF_filename>
 
-   Specify the output PDB filename. Defaults to the input filename with a suffix defined by ``--suffix``.
+   Specify the output SDF filename. Defaults to the input filename with a suffix defined by ``--suffix``. 
 
-.. option:: --suffix <suffix>
+.. option:: -j, --read_json <input_JSON_filename>
 
-   Set a suffix for output filenames that are not explicitly specified. The default suffix is ``_docked``.
+   Provide a receptor file generated by ``mk_prepare_receptor`` with the ``-j/--write_json`` option. Currently only effective when used with ``-p, --write_pdb``. 
 
-.. option:: -j, --read_json <input_JSON_filename>
+.. option:: -p, --write_pdb <output_PDB_filename>
 
-   Provide a receptor file generated by ``mk_prepare_receptor`` with the ``-j/--write_json`` option.
+   Specify the output PDB filename. Defaults to the input filename with a suffix defined by ``--suffix``. Must be used together with ``-j, --read_json``. 
 
 .. option:: --all_dlg_poses
 
-   (Flag) Write all poses from AutoDock-GPU DLG output files, instead of only the lead of each cluster.
+   (Flag) Write all poses from AutoDock-GPU DLG output files, instead of only the lead of each cluster. Currently only effective for ``-s, --write_sdf``. 
 
 .. option:: -k, --keep_flexres_sdf
 
diff --git a/docs/source/cli_lig_prep.rst b/docs/source/cli_lig_prep.rst
index 91fe3224..b246c043 100644
--- a/docs/source/cli_lig_prep.rst
+++ b/docs/source/cli_lig_prep.rst
@@ -3,8 +3,11 @@ mk_prepare_ligand.py
 
 Command line tool to prepare small organic molecules.
 
+About
+-----
+
 Write PDBQT files
------------------
+~~~~~~~~~~~~~~~~~
 
 AutoDock-GPU and Vina read molecules in the PDBQT format. These can be prepared
 by Meeko from SD files, or from Mol2 files, but SDF is strongly preferred.
@@ -14,3 +17,83 @@ by Meeko from SD files, or from Mol2 files, but SDF is strongly preferred.
     mk_prepare_ligand.py -i molecule.sdf -o molecule.pdbqt
     mk_prepare_ligand.py -i multi_mol.sdf --multimol_outdir folder_for_pdbqt_files
 
+Usage
+-----
+
+.. code-block:: bash
+
+   python mk_prepare_ligand.py [OPTIONS]
+
+Positional Argument
+~~~~~~~~~~~~~~~~~~~
+
+.. option:: -i, --mol <input_molecule_filename>
+
+   The input molecule file, in formats such as MOL2, SDF, etc.
+
+Options
+~~~~~~~
+
+.. option:: -c, --config_file <config_file>
+
+   Configure `MoleculePreparation` from a JSON file. Command-line arguments will override settings in the file.
+
+   **Example:**
+
+   .. code-block:: bash
+
+      python mk_prepare_ligand.py -c config.json -i ligand.sdf
+
+.. option:: -v, --verbose
+
+   (Flag) Print detailed information about the molecule setup process.
+
+.. option:: --name_from_prop <property_name>
+
+   Set the molecule name using a specified RDKit or SDF property.
+
+   **Example:**
+
+   .. code-block:: bash
+
+      python mk_prepare_ligand.py -i ligand.sdf --name_from_prop compound_name
+
+.. option:: -o, --out <output_pdbqt_filename>
+
+   Specify the output PDBQT filename. Only compatible with single-molecule input.
+
+.. option:: --multimol_outdir <output_directory>
+
+   Specify the directory to write PDBQT output files for multi-molecule inputs. Incompatible with `-o/--out` and `-`/`--`.
+
+.. option:: --multimol_prefix <prefix>
+
+   Replace the internal molecule name in multi-molecule input with the specified prefix. Incompatible with `-o/--out` and `-`/`--`.
+
+.. option:: -z, --multimol_targz
+
+   (Flag) Compress output files into a `.tar.gz` archive.
+
+.. option:: --multimol_targz_size <size>
+
+   Define the number of PDBQT files per `.tar.gz` archive. Default is 10000.
+
+.. option:: -, --
+
+   (Flag) Redirect output to standard output (STDOUT) instead of writing a file. Ignored if `-o/--out` is specified. Only compatible with single-molecule input.
+
+Molecule Preparation Options
+----------------------------
+
+.. option:: --rigid_macrocycles
+
+   (Flag) Keep macrocycles rigid in their input conformation.
+
+.. option:: --macrocycle_allow_A
+
+   (Flag) Allow bond break with atom type A, which will be retyped as carbon (C).
+
+.. option:: --keep_chorded_rings
+
+   (Flag) Retain all rings from exhaustive ring perception.
+
diff --git a/docs/source/cli_rec_prep.rst b/docs/source/cli_rec_prep.rst
index 2fadd173..aa0ecae2 100644
--- a/docs/source/cli_rec_prep.rst
+++ b/docs/source/cli_rec_prep.rst
@@ -1,3 +1,9 @@
+mk_prepare_receptor
+===================
+
+About
+-----
+
 The input structure is matched against templates to
 guarantee chemical correctness and identify problems with the input structures.
 This allows the user to identify and fix problems, resulting in a molecular
@@ -27,28 +33,22 @@ Residue name is primary key unless user overrides.
 
 Currently not supported: capped residues from charmm-gui.
 
-mk_prepare_receptor
-===================
-
 Basic usage
------------
+~~~~~~~~~~~
 
 .. code-block:: bash
 
     mk_prepare_receptor -i examples/system.pdb --write_pdbqt prepared.pdbqt
 
-
-
-
 Protonation states
-------------------
+~~~~~~~~~~~~~~~~~~
 
 
 Adding templates
-----------------
+~~~~~~~~~~~~~~~~
 
 Write flags
------------
+~~~~~~~~~~~
 
 The option flags starting with ``--write`` in  ``mk_prepare_receptor`` can
 be used both with an argument to specify the outpuf filename: 
@@ -77,7 +77,7 @@ in which case the specified filenames have priority over the default basename.
 .. _templates:
 
 Templates
----------
+~~~~~~~~~
 
 The templates contain SMILES strings that are used to create the RDKit
 molecules that constitute every residue in the processed model. In this way,
@@ -102,3 +102,174 @@ is a single carbon atom. Our template SMILES would be "[H]C[H]". The RDKit
 molecule will have three atoms and the carbon will have two implicit hydrogens.
 The implicit hydrogens correspond to bonds to adjacent residues in the
 processed polymer.
+
+Usage
+-----
+
+.. code-block:: bash
+
+   python mk_prepare_receptor.py [OPTIONS]
+
+Options
+~~~~~~~
+
+Input/Output Options
+~~~~~~~~~~~~~~~~~~~~
+
+.. option:: --read_pdb <PDB_FILENAME>
+
+   Read a PDB file directly (not in PDBQT format) without using ProDy.
+
+.. option:: -i, --read_with_prody <MACROMOL_FILENAME>
+
+   Read a PDB or mmCIF file using ProDy (if installed). ProDy can be installed from PyPI or conda-forge.
+
+.. option:: -o, --output_basename <basename>
+
+   Specify a default basename for output files created by `--write` options when no filename is specified.
+
+.. option:: -p, --write_pdbqt <PDBQT_FILENAME> [*]
+
+   Output PDBQT files with `_rigid` or `_flex` suffixes for flexible residues. Defaults to `--output_basename` if no filename is provided.
+
+.. option:: -j, --write_json <JSON_FILENAME> [*]
+
+   Save the receptor's parameterized configuration to JSON format. Defaults to `--output_basename` if unspecified.
+
+.. option:: --write_pdb <PDB_FILENAME> [*]
+
+   Save the prepared receptor in PDB format. Must specify the filename.
+
+.. option:: -g, --write_gpf <GPF_FILENAME> [*]
+
+   Output an AutoGrid input file (GPF). Defaults to `--output_basename` if not specified.
+
+.. option:: -v, --write_vina_box <VINA_BOX_FILENAME> [*]
+
+   Generate a configuration file for Vina with grid box dimensions. Defaults to `--output_basename` if not specified.
+
+Receptor Perception Options
+---------------------------
+
+.. option:: -n, --set_template <template>
+
+   Assign templates to residues, e.g., `A:5,7=CYX,B:17=HID`.
+
+.. option:: -d, --delete_residues <residues>
+
+   Specify residues to delete, e.g., `A:350,B:15,16,17`.
+
+.. option:: -b, --blunt_ends <positions>
+
+   Blunt end definitions, e.g., `A:123,200=2,A:1=0`.
+
+.. option:: --add_templates <JSON_FILENAME> [*]
+
+   Load additional templates from one or more JSON files.
+
+.. option:: --mk_config <JSON_FILENAME>
+
+   Specify a JSON configuration file for receptor preparation.
+
+.. option:: -a, --allow_bad_res
+
+   (Flag) Ignore residues with missing atoms instead of raising an error.
+
+.. option:: --default_altloc <location>
+
+   Define a default alternate location for residues, overridden by `--wanted_altloc`.
+
+.. option:: --wanted_altloc <location>
+
+   Specify alternate locations for particular residues, e.g., `:5=B,B:17=A`.
+
+.. option:: -f, --flexres <residues>
+
+   Define flexible residues by chain ID and residue number, e.g., `-f ":42,B:23"`.
+
+Grid Box Options
+----------------
+
+.. option:: --box_size <X Y Z>
+
+   Set the size of the grid box in Angstroms (x, y, z).
+
+.. option:: --box_center <X Y Z>
+
+   Define the center of the grid box in Angstroms (x, y, z).
+
+.. option:: --box_center_off_reactive_res
+
+   (Flag) Shift the grid box center 5 Å along the CA-CB bond from CB. Applicable only when there is one reactive flexible residue.
+
+.. option:: --box_enveloping <FILENAME>
+
+   Adjust the grid box to enclose atoms in the specified file. Supported formats: `.sdf`, `.mol`, `.mol2`, `.pdb`, `.pdbqt`.
+
+.. option:: --padding <value>
+
+   Set padding around atoms specified in `--box_enveloping` (in Å).
+
+Reactive Options
+----------------
+
+.. option:: -r, --reactive_flexres <residues>
+
+   Define reactive flexible residues by chain ID and residue number, e.g., `-r ":42,B:23"`. Maximum of 8 reactive residues.
+
+.. option:: --reactive_name <residue:atom>
+
+   Specify the reactive atom name for a residue type, e.g., `--reactive_name "TRP:NE1"`. Can be repeated for multiple assignments.
+
+.. option:: -s, --reactive_name_specific <residue:atom>
+
+   Specify the reactive atom for an individual residue by residue ID, e.g., `-s "A:42=NE2"`. The residue becomes reactive.
+
+.. option:: --r_eq_12 <value>
+
+   Set the equilibrium distance (r_eq) for reactive atoms in 1-2 interactions. Default is 1.8 Å.
+
+.. option:: --eps_12 <value>
+
+   Set the epsilon value for reactive atoms in 1-2 interactions. Default is 2.5.
+
+.. option:: --r_eq_13_scaling <factor>
+
+   Scale r_eq for 1-3 interactions across reactive atoms. Default is 0.5.
+
+.. option:: --r_eq_14_scaling <factor>
+
+   Scale r_eq for 1-4 interactions across reactive atoms. Default is 0.5.
+
+Examples
+--------
+
+Basic usage to read a PDB file and specify output basename:
+
+.. code-block:: bash
+
+   python mk_prepare_receptor.py --read_pdb receptor.pdb -o output_basename
+
+Using ProDy to read a PDB file and output PDBQT files:
+
+.. code-block:: bash
+
+   python mk_prepare_receptor.py -i receptor.pdb --write_pdbqt receptor_output.pdbqt
+
+Saving receptor configuration as JSON:
+
+.. code-block:: bash
+
+   python mk_prepare_receptor.py -j receptor.json -i receptor.pdb
+
+Defining grid box size and center:
+
+.. code-block:: bash
+
+   python mk_prepare_receptor.py --box_size 20 20 20 --box_center 10 10 10
+
+Defining reactive flexible residues and adjusting grid box:
+
+.. code-block:: bash
+
+   python mk_prepare_receptor.py -r ":42,B:23" --box_center_off_reactive_res

From 1b88ecf65da8c48ddef51bf0dbc6e9561a3b496a Mon Sep 17 00:00:00 2001
From: Amy He <ayhe1768@gmail.com>
Date: Wed, 6 Nov 2024 22:41:31 -0800
Subject: [PATCH 04/18] explain selection and assignment language

---
 docs/source/cli_rec_prep.rst | 120 ++++++++++++++++++++++-------------
 1 file changed, 76 insertions(+), 44 deletions(-)

diff --git a/docs/source/cli_rec_prep.rst b/docs/source/cli_rec_prep.rst
index aa0ecae2..8007bc99 100644
--- a/docs/source/cli_rec_prep.rst
+++ b/docs/source/cli_rec_prep.rst
@@ -1,5 +1,22 @@
-mk_prepare_receptor
-===================
+mk_prepare_receptor.py
+======================
+
+A command-line script for receptor preparation from a PDB/CIF file, setting up various properties (flexible/reactive residues, box specifications, etc.), and writing useful input files (PDBQT, GPF, box configuration file for Vina) for docking. 
+
+Basic usage
+-----------
+
+.. code-block:: bash
+
+    mk_prepare_receptor.py -i examples/system.pdb --write_pdbqt prepared.pdbqt
+
+This is equivalent to: 
+
+.. code-block:: bash
+
+    mk_prepare_receptor.py -i examples/system.pdb -o prepared -p
+
+Read more about the syntax in :ref:`Write flags` and :ref:`Options`. 
 
 About
 -----
@@ -33,19 +50,34 @@ Residue name is primary key unless user overrides.
 
 Currently not supported: capped residues from charmm-gui.
 
-Basic usage
-~~~~~~~~~~~
-
-.. code-block:: bash
+.. _templates:
 
-    mk_prepare_receptor -i examples/system.pdb --write_pdbqt prepared.pdbqt
+Templates
+~~~~~~~~~
 
-Protonation states
-~~~~~~~~~~~~~~~~~~
+The templates contain SMILES strings that are used to create the RDKit
+molecules that constitute every residue in the processed model. In this way,
+the chemistry of the processed model is fully defined by the templates,
+and the only thing that is preserved from the input are the atom positions
+and the connectivity between residues.
 
+The SMILES strings contain all atoms that exist in the final model,
+and none that do not exist. This also applies to hydrogens,
+meaning that the SMILES are expected to have real hydrogens. Note that
+real hydrogens are different from explicit hydrogens. Real hydrogens will be
+represented as an actual atom in an RDKit molecule, while explicit hydrogens
+are a just property of heavy atoms. In the SMILES, real hydrogens are defined
+with square brackets "[H]" and explicit hydrogens without, e.g. "[nH]" to set
+the number of explicit hydrogens on an aromatic nitrogen to one.
 
-Adding templates
-~~~~~~~~~~~~~~~~
+Residues that are part of a polymer, which is often all of them, will have
+bonds to adjacent residues. The heavy atoms involved in the bonds will miss
+a real hydrogen and have an implicit (or explicit) one instead. As an
+example, consider modeling an alkyl chain as a polymer, in which the monomer
+is a single carbon atom. Our template SMILES would be "[H]C[H]". The RDKit
+molecule will have three atoms and the carbon will have two implicit hydrogens.
+The implicit hydrogens correspond to bonds to adjacent residues in the
+processed polymer. 
 
 Write flags
 ~~~~~~~~~~~
@@ -72,53 +104,53 @@ It is also possible to combine the two types of usage:
     --write_json
     --write_vina_box box_for_myenzyme.txt
 
-in which case the specified filenames have priority over the default basename.
+in which case the specified filenames have priority over the default basename. 
 
-.. _templates:
+Residue selection and assignment language
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-Templates
-~~~~~~~~~
+Meeko uses the **chain ID** and **residue number** to identify a residue. The arguments involving selection of residues: 
 
-The templates contain SMILES strings that are used to create the RDKit
-molecules that constitute every residue in the processed model. In this way,
-the chemistry of the processed model is fully defined by the templates,
-and the only thing that is preserved from the input are the atom positions
-and the connectivity between residues.
+.. option:: -d, --delete_residues <residues>
 
-The SMILES strings contain all atoms that exist in the final model,
-and none that do not exist. This also applies to hydrogens,
-meaning that the SMILES are expected to have real hydrogens. Note that
-real hydrogens are different from explicit hydrogens. Real hydrogens will be
-represented as an actual atom in an RDKit molecule, while explicit hydrogens
-are a just property of heavy atoms. In the SMILES, real hydrogens are defined
-with square brackets "[H]" and explicit hydrogens without, e.g. "[nH]" to set
-the number of explicit hydrogens on an aromatic nitrogen to one.
+.. option:: -f, --flexres <residues>
 
-Residues that are part of a polymer, which is often all of them, will have
-bonds to adjacent residues. The heavy atoms involved in the bonds will miss
-a real hydrogen and have an implicit (or explicit) one instead. As an
-example, consider modeling an alkyl chain as a polymer, in which the monomer
-is a single carbon atom. Our template SMILES would be "[H]C[H]". The RDKit
-molecule will have three atoms and the carbon will have two implicit hydrogens.
-The implicit hydrogens correspond to bonds to adjacent residues in the
-processed polymer.
+.. option:: -r, --reactive_flexres <residues>
+
+use the compact selection language that specify residues efficiently. The chain ID and the residue number(s) are separated by a colon (``:``) delimiter. Each residue number is combined with the most recent chain ID that precedes it, resulting in an expanded list of chain-residue pairs. 
+
+For an input like ``A:5,7,BB:12C``, this selection language represents: ``residues (number) 5 and 7 in Chain A`` and ``residue (number) 12C in Chain BB``. 
+
+The arguments involving assignment of residues to properties: 
+
+.. option:: -n, --set_template <template>
+
+.. option:: -b, --blunt_ends <positions>
+
+.. option:: --wanted_altloc <location>
+
+.. option:: -s, --reactive_name_specific <residue:atom>
+
+use the residue selection lanaguge described above, followed by an equal sign (``=``) as the delimiter and the assigned value, which could be the name of a residue template, the atom index for the blunt end, the wanted altloc ID, or the atom name of the reactive atom. Each residue selection is comibned with the most recent assignment that precedes it, resulting in an expanded list of residue-assignment pairs. 
+
+For an input like ``"A:5,7=CYX,A:19A,B:17=HID``, this assignment language represents: ``residues (number) 5 in Chain A are set to (template name) CYX`` and ``residue (number) 19 A in Chain A, and residue (number) 17 in Chain B are set to (template name) HID``. 
 
 Usage
 -----
 
 .. code-block:: bash
 
-   python mk_prepare_receptor.py [OPTIONS]
+   mk_prepare_receptor.py [OPTIONS]
 
 Options
 ~~~~~~~
 
 Input/Output Options
-~~~~~~~~~~~~~~~~~~~~
+^^^^^^^^^^^^^^^^^^^^
 
 .. option:: --read_pdb <PDB_FILENAME>
 
-   Read a PDB file directly (not in PDBQT format) without using ProDy.
+   Read a PDB file using the PDB parser in RDKit.
 
 .. option:: -i, --read_with_prody <MACROMOL_FILENAME>
 
@@ -136,10 +168,6 @@ Input/Output Options
 
    Save the receptor's parameterized configuration to JSON format. Defaults to `--output_basename` if unspecified.
 
-.. option:: --write_pdb <PDB_FILENAME> [*]
-
-   Save the prepared receptor in PDB format. Must specify the filename.
-
 .. option:: -g, --write_gpf <GPF_FILENAME> [*]
 
    Output an AutoGrid input file (GPF). Defaults to `--output_basename` if not specified.
@@ -148,8 +176,12 @@ Input/Output Options
 
    Generate a configuration file for Vina with grid box dimensions. Defaults to `--output_basename` if not specified.
 
+.. option:: --write_pdb <PDB_FILENAME> [*]
+
+   Save the prepared receptor in PDB format. Must specify the filename.
+
 Receptor Perception Options
----------------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 .. option:: -n, --set_template <template>
 

From ef3ff7bd63fd0ab93bbad3f0468fed42852de6fa Mon Sep 17 00:00:00 2001
From: Amy He <ayhe1768@gmail.com>
Date: Wed, 6 Nov 2024 22:48:56 -0800
Subject: [PATCH 05/18] finish draft cli_rec_prep

---
 docs/source/cli_rec_prep.rst | 50 ++++++++++++++++--------------------
 1 file changed, 22 insertions(+), 28 deletions(-)

diff --git a/docs/source/cli_rec_prep.rst b/docs/source/cli_rec_prep.rst
index 8007bc99..f6d0729e 100644
--- a/docs/source/cli_rec_prep.rst
+++ b/docs/source/cli_rec_prep.rst
@@ -99,10 +99,7 @@ It is also possible to combine the two types of usage:
 
 .. code-block:: bash
 
-    --output_basename myenzyme
-    --write_pdbqt
-    --write_json
-    --write_vina_box box_for_myenzyme.txt
+    --output_basename myenzyme --write_pdbqt --write_json --write_vina_box box_for_myenzyme.txt
 
 in which case the specified filenames have priority over the default basename. 
 
@@ -111,11 +108,13 @@ Residue selection and assignment language
 
 Meeko uses the **chain ID** and **residue number** to identify a residue. The arguments involving selection of residues: 
 
-.. option:: -d, --delete_residues <residues>
+.. code-block:: bash
 
-.. option:: -f, --flexres <residues>
+    -d, --delete_residues <residues>
 
-.. option:: -r, --reactive_flexres <residues>
+    -f, --flexres <residues>
+
+    -r, --reactive_flexres <residues>
 
 use the compact selection language that specify residues efficiently. The chain ID and the residue number(s) are separated by a colon (``:``) delimiter. Each residue number is combined with the most recent chain ID that precedes it, resulting in an expanded list of chain-residue pairs. 
 
@@ -123,27 +122,22 @@ For an input like ``A:5,7,BB:12C``, this selection language represents: ``residu
 
 The arguments involving assignment of residues to properties: 
 
-.. option:: -n, --set_template <template>
-
-.. option:: -b, --blunt_ends <positions>
+.. code-block:: bash
 
-.. option:: --wanted_altloc <location>
+    -n, --set_template <template>
 
-.. option:: -s, --reactive_name_specific <residue:atom>
+    -b, --blunt_ends <positions>
 
-use the residue selection lanaguge described above, followed by an equal sign (``=``) as the delimiter and the assigned value, which could be the name of a residue template, the atom index for the blunt end, the wanted altloc ID, or the atom name of the reactive atom. Each residue selection is comibned with the most recent assignment that precedes it, resulting in an expanded list of residue-assignment pairs. 
+    --wanted_altloc <location>
 
-For an input like ``"A:5,7=CYX,A:19A,B:17=HID``, this assignment language represents: ``residues (number) 5 in Chain A are set to (template name) CYX`` and ``residue (number) 19 A in Chain A, and residue (number) 17 in Chain B are set to (template name) HID``. 
+    -s, --reactive_name_specific <residue:atom>
 
-Usage
------
+use the residue selection lanaguge described above, followed by an equal sign (``=``) as the delimiter and the assigned value, which could be the name of a residue template, the atom index for the blunt end, the wanted altloc ID, or the atom name of the reactive atom. Each residue selection is comibned with the most recent assignment that precedes it, resulting in a further expanded list of residue-assignment pairs. 
 
-.. code-block:: bash
-
-   mk_prepare_receptor.py [OPTIONS]
+For an input like ``"A:5,7=CYX,A:19A,B:17=HID``, this assignment language represents: ``residues (number) 5 in Chain A are set to (template name) CYX`` and ``residue (number) 19 A in Chain A, and residue (number) 17 in Chain B are set to (template name) HID``. 
 
 Options
-~~~~~~~
+-------
 
 Input/Output Options
 ^^^^^^^^^^^^^^^^^^^^
@@ -199,10 +193,6 @@ Receptor Perception Options
 
    Load additional templates from one or more JSON files.
 
-.. option:: --mk_config <JSON_FILENAME>
-
-   Specify a JSON configuration file for receptor preparation.
-
 .. option:: -a, --allow_bad_res
 
    (Flag) Ignore residues with missing atoms instead of raising an error.
@@ -215,9 +205,9 @@ Receptor Perception Options
 
    Specify alternate locations for particular residues, e.g., `:5=B,B:17=A`.
 
-.. option:: -f, --flexres <residues>
+.. option:: --mk_config <JSON_FILENAME>
 
-   Define flexible residues by chain ID and residue number, e.g., `-f ":42,B:23"`.
+   Specify a JSON configuration file for receptor preparation.
 
 Grid Box Options
 ----------------
@@ -242,8 +232,12 @@ Grid Box Options
 
    Set padding around atoms specified in `--box_enveloping` (in Å).
 
-Reactive Options
-----------------
+Flexible and/or Reactive Options
+--------------------------------
+
+.. option:: -f, --flexres <residues>
+
+   Define flexible residues by chain ID and residue number, e.g., `-f ":42,B:23"`. 
 
 .. option:: -r, --reactive_flexres <residues>
 

From 9c70a4299267c0ee1213733d9b39fcd950ba5923 Mon Sep 17 00:00:00 2001
From: Amy He <ayhe1768@gmail.com>
Date: Wed, 6 Nov 2024 22:50:22 -0800
Subject: [PATCH 06/18] fix heading level

---
 docs/source/cli_rec_prep.rst | 37 ++----------------------------------
 1 file changed, 2 insertions(+), 35 deletions(-)

diff --git a/docs/source/cli_rec_prep.rst b/docs/source/cli_rec_prep.rst
index f6d0729e..588b9b17 100644
--- a/docs/source/cli_rec_prep.rst
+++ b/docs/source/cli_rec_prep.rst
@@ -210,7 +210,7 @@ Receptor Perception Options
    Specify a JSON configuration file for receptor preparation.
 
 Grid Box Options
-----------------
+^^^^^^^^^^^^^^^^
 
 .. option:: --box_size <X Y Z>
 
@@ -233,7 +233,7 @@ Grid Box Options
    Set padding around atoms specified in `--box_enveloping` (in Å).
 
 Flexible and/or Reactive Options
---------------------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 .. option:: -f, --flexres <residues>
 
@@ -266,36 +266,3 @@ Flexible and/or Reactive Options
 .. option:: --r_eq_14_scaling <factor>
 
    Scale r_eq for 1-4 interactions across reactive atoms. Default is 0.5.
-
-Examples
---------
-
-Basic usage to read a PDB file and specify output basename:
-
-.. code-block:: bash
-
-   python mk_prepare_receptor.py --read_pdb receptor.pdb -o output_basename
-
-Using ProDy to read a PDB file and output PDBQT files:
-
-.. code-block:: bash
-
-   python mk_prepare_receptor.py -i receptor.pdb --write_pdbqt receptor_output.pdbqt
-
-Saving receptor configuration as JSON:
-
-.. code-block:: bash
-
-   python mk_prepare_receptor.py -j receptor.json -i receptor.pdb
-
-Defining grid box size and center:
-
-.. code-block:: bash
-
-   python mk_prepare_receptor.py --box_size 20 20 20 --box_center 10 10 10
-
-Defining reactive flexible residues and adjusting grid box:
-
-.. code-block:: bash
-
-   python mk_prepare_receptor.py -r ":42,B:23" --box_center_off_reactive_res

From a4b01bafda5b95daeb66a28f9f7c3f09e6015e55 Mon Sep 17 00:00:00 2001
From: Amy He <ayhe1768@gmail.com>
Date: Wed, 6 Nov 2024 22:52:04 -0800
Subject: [PATCH 07/18] fix heading level 2

---
 docs/source/cli_rec_prep.rst | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/source/cli_rec_prep.rst b/docs/source/cli_rec_prep.rst
index 588b9b17..a58eda9a 100644
--- a/docs/source/cli_rec_prep.rst
+++ b/docs/source/cli_rec_prep.rst
@@ -140,7 +140,7 @@ Options
 -------
 
 Input/Output Options
-^^^^^^^^^^^^^^^^^^^^
+~~~~~~~~~~~~~~~~~~~~
 
 .. option:: --read_pdb <PDB_FILENAME>
 
@@ -175,7 +175,7 @@ Input/Output Options
    Save the prepared receptor in PDB format. Must specify the filename.
 
 Receptor Perception Options
-^^^^^^^^^^^^^^^^^^^^^^^^^^^
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 .. option:: -n, --set_template <template>
 
@@ -210,7 +210,7 @@ Receptor Perception Options
    Specify a JSON configuration file for receptor preparation.
 
 Grid Box Options
-^^^^^^^^^^^^^^^^
+~~~~~~~~~~~~~~~~
 
 .. option:: --box_size <X Y Z>
 
@@ -233,7 +233,7 @@ Grid Box Options
    Set padding around atoms specified in `--box_enveloping` (in Å).
 
 Flexible and/or Reactive Options
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 .. option:: -f, --flexres <residues>
 

From fc427efaaef471dd945b7a7641bad08480ce3e98 Mon Sep 17 00:00:00 2001
From: Amy He <ayhe1768@gmail.com>
Date: Wed, 6 Nov 2024 23:41:08 -0800
Subject: [PATCH 08/18] draft cli_lig_prep

---
 docs/source/cli_lig_prep.rst   | 93 +++++++++++++++++++---------------
 docs/source/colab_examples.rst |  8 +++
 2 files changed, 59 insertions(+), 42 deletions(-)

diff --git a/docs/source/cli_lig_prep.rst b/docs/source/cli_lig_prep.rst
index b246c043..503fd4af 100644
--- a/docs/source/cli_lig_prep.rst
+++ b/docs/source/cli_lig_prep.rst
@@ -1,63 +1,33 @@
 mk_prepare_ligand.py
 ====================
 
-Command line tool to prepare small organic molecules.
+A command-line script for ligand preparation to generate the ligand PDBQT file. Currently support SD files (.sdf), Mol2 files (.mol2) and Mol files (.mol), but SDF is strongly preferred. 
 
-About
------
-
-Write PDBQT files
-~~~~~~~~~~~~~~~~~
-
-AutoDock-GPU and Vina read molecules in the PDBQT format. These can be prepared
-by Meeko from SD files, or from Mol2 files, but SDF is strongly preferred.
+Basic usage
+-----------
 
 .. code-block:: bash
 
+    # write a single ligand PDBQT file from a single-molecule input file
     mk_prepare_ligand.py -i molecule.sdf -o molecule.pdbqt
-    mk_prepare_ligand.py -i multi_mol.sdf --multimol_outdir folder_for_pdbqt_files
-
-Usage
------
-
-.. code-block:: bash
 
-   python mk_prepare_ligand.py [OPTIONS]
-
-Positional Argument
-~~~~~~~~~~~~~~~~~~~
-
-.. option:: -i, --mol <input_molecule_filename>
-
-   The input molecule file, in formats such as MOL2, SDF, etc.
+    # prepare ligand PDBQT files in batch from a multiple-molecule input file
+    mk_prepare_ligand.py -i multi_mol.sdf --multimol_outdir folder_for_pdbqt_files
 
 Options
-~~~~~~~
+-------
 
-.. option:: -c, --config_file <config_file>
+Input/Output Options
+~~~~~~~~~~~~~~~~~~~~
 
-   Configure `MoleculePreparation` from a JSON file. Command-line arguments will override settings in the file.
-
-   **Example:**
-
-   .. code-block:: bash
-
-      python mk_prepare_ligand.py -c config.json -i ligand.sdf
-
-.. option:: -v, --verbose
+.. option:: -i, --mol <input_molecule_filename>
 
-   (Flag) Print detailed information about the molecule setup process.
+   The input molecule file, in formats such as MOL2, SDF, etc. This option is required.
 
 .. option:: --name_from_prop <property_name>
 
    Set the molecule name using a specified RDKit or SDF property.
 
-   **Example:**
-
-   .. code-block:: bash
-
-      python mk_prepare_ligand.py -i ligand.sdf --name_from_prop compound_name
-
 .. option:: -o, --out <output_pdbqt_filename>
 
    Specify the output PDBQT filename. Only compatible with single-molecule input.
@@ -83,7 +53,11 @@ Options
    (Flag) Redirect output to standard output (STDOUT) instead of writing a file. Ignored if `-o/--out` is specified. Only compatible with single-molecule input.
 
 Molecule Preparation Options
-----------------------------
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. option:: -c, --config_file <config_file>
+
+   Configure `MoleculePreparation` from a JSON file. Command-line arguments will override settings in the file.
 
 .. option:: --rigid_macrocycles
 
@@ -97,3 +71,38 @@ Molecule Preparation Options
 
    (Flag) Retain all rings from exhaustive ring perception.
 
+.. option:: --keep_equivalent_rings
+
+   (Flag) Retain rings with equivalent sizes and neighboring atoms.
+
+.. option:: --min_ring_size <size>
+
+   Define the minimum number of atoms required in a ring for it to be considered for opening.
+
+.. option:: -w, --hydrate
+
+   (Flag) Add water molecules to the structure for hydrated docking.
+
+.. option:: --merge_these_atom_types <types> [*]
+
+   Specify a list of atom types to merge. The default is `"H"`.
+
+.. option:: -r, --rigidify_bonds_smarts <SMARTS>
+
+   Provide SMARTS patterns to rigidify specific bonds in the molecule.
+
+.. option:: -b, --rigidify_bonds_indices <i j>
+
+   Specify the indices of two atoms that define a bond in the SMARTS pattern (starting from 1).
+
+.. option:: -a, --flexible_amides
+
+   (Flag) Allow amide bonds to rotate, making them non-planar (not recommended).
+
+.. option:: -p, --atom_type_smarts <JSON_FILENAME>
+
+   Specify SMARTS-based atom typing in JSON format.
+
+.. option:: -v, --verbose
+
+   (Flag) Print detailed information about the molecule setup process.
diff --git a/docs/source/colab_examples.rst b/docs/source/colab_examples.rst
index b92554da..eaa2c2d0 100644
--- a/docs/source/colab_examples.rst
+++ b/docs/source/colab_examples.rst
@@ -72,3 +72,11 @@ The basic docking example is developed to showcase the usage of **import additio
 The reactive docking example is based on reactive docking method that has been developed for high-throughput virtual screenings of reactive species. In this example, a small molecule substrate (Adenosine monophosphate, PDB token AMP) is targeting at the catalytic histidine residue of a hollow protein structure of bacteria RNA 3' cyclase (PDB token 3KGD) to generate the near-attack conformation for the formation of the phosphoamide bond. A docked pose that closely resembles the original position of the ligand is expected among the top-ranked poses. 
 
 
+[AutoDock-GPU] Tethered Docking
+---------------
+
+`Run on Colab! <https://colab.research.google.com/drive/1tf9xOgn6u8eDTeFJtc8GCEGRX-8aR9Bo?usp=sharing>`_
+
+The covalent docking example is based on the two-point attractor and flexible side chain method. In this example, a small molecule substrate (Adenosine monophosphate, PDB token AMP) is attached onto the catalytic histidine residue of a hollow protein structure of bacteria RNA 3' cyclase (PDB token 3KGD) to reproduce the covalent intermediate complex structure. A docked pose that closely resembles the original position of the ligand is expected among the top-ranked poses. 
+
+

From 5744bd4089fb1d5f7a72fd9137d9955335b256ac Mon Sep 17 00:00:00 2001
From: Amy He <ayhe1768@gmail.com>
Date: Wed, 6 Nov 2024 23:44:40 -0800
Subject: [PATCH 09/18] include all options in cli_lig_prep

---
 docs/source/cli_lig_prep.rst | 56 ++++++++++++++++++++++++++++++++++--
 1 file changed, 53 insertions(+), 3 deletions(-)

diff --git a/docs/source/cli_lig_prep.rst b/docs/source/cli_lig_prep.rst
index 503fd4af..688095d8 100644
--- a/docs/source/cli_lig_prep.rst
+++ b/docs/source/cli_lig_prep.rst
@@ -1,7 +1,7 @@
 mk_prepare_ligand.py
 ====================
 
-A command-line script for ligand preparation to generate the ligand PDBQT file. Currently support SD files (.sdf), Mol2 files (.mol2) and Mol files (.mol), but SDF is strongly preferred. 
+A command-line script for ligand preparation that generates the ligand PDBQT file(s). Currently supports SD files (.sdf), Mol2 files (.mol2) and Mol files (.mol), but SDF is strongly preferred as input files. 
 
 Basic usage
 -----------
@@ -103,6 +103,56 @@ Molecule Preparation Options
 
    Specify SMARTS-based atom typing in JSON format.
 
-.. option:: -v, --verbose
+.. option:: -aa, --add_atom_types <JSON>
 
-   (Flag) Print detailed information about the molecule setup process.
+   Specify additional atom types to assign in JSON format, with SMARTS patterns and atom type names.
+
+.. option:: --double_bond_penalty <penalty>
+
+   Set a penalty value; values greater than 100 prevent breaking double bonds.
+
+.. option:: --charge_model <model>
+
+   Choose the charge model: `gasteiger`, `espaloma`, or `zero`. Default is `gasteiger`; `zero` sets all charges to zero.
+
+.. option:: --bad_charge_ok
+
+   (Flag) Allow NaN and Inf charges in the PDBQT output.
+
+.. option:: --add_index_map
+
+   (Flag) Include a map of atom indices from the input to the PDBQT file.
+
+.. option:: --remove_smiles
+
+   (Flag) Exclude SMILES from being written as a remark in the PDBQT output.
+
+Reactive Docking Options
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. option:: --reactive_smarts <SMARTS>
+
+   Provide a SMARTS pattern for defining the reactive group.
+
+.. option:: --reactive_smarts_idx <index>
+
+   Specify the 1-based index of the reactive atom within the SMARTS pattern provided by `--reactive_smarts`.
+
+Covalent Docking (Tethered) Options
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. option:: --receptor <filename>
+
+   Specify the receptor file. Supported formats depend on ProDy availability, such as `.pdb` and `.mmcif`.
+
+.. option:: --rec_residue <residue>
+
+   Specify the residue in the receptor for attachment, e.g., `A:LYS:204`.
+
+.. option:: --tether_smarts <SMARTS>
+
+   Provide a SMARTS pattern defining the ligand atoms used for attachment to the receptor.
+
+.. option:: --tether_smarts_indices <IDX IDX>
+
+   Specify the 1-based indices of the two atoms in the SMARTS pattern that will be attached (default: `1 2`).
\ No newline at end of file

From 89f667d8905e646e05214e0ca0b6f7d462d58205 Mon Sep 17 00:00:00 2001
From: Amy He <ayhe1768@gmail.com>
Date: Wed, 6 Nov 2024 23:51:52 -0800
Subject: [PATCH 10/18] adjust section order in cli_export_result

---
 docs/source/cli_export_result.rst | 34 +++++++++++++++----------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/docs/source/cli_export_result.rst b/docs/source/cli_export_result.rst
index 1c35736f..3a35d263 100644
--- a/docs/source/cli_export_result.rst
+++ b/docs/source/cli_export_result.rst
@@ -1,6 +1,16 @@
 mk_export.py
 ============
 
+A command-line script to export docking poses of ligand (and flexible residues) to an SD file (.sdf), and to export the full receptor structures with updated conformations of flexible residues to a PDB file. Currently supports PDBQT files from AutoDock-Vina and DLG files from AutoDock-GPU that have the REMARK lines containing Smiles and index mapping information. 
+
+Basic usage
+-----------
+
+.. code-block:: bash
+
+    mk_export.py vina_results.pdbqt -s vina_results.sdf
+    mk_export.py autodock-gpu_results.dlg -s autodock-gpu_results.sdf
+
 About
 -----
 
@@ -12,13 +22,7 @@ from AutoDock-GPU contains docked poses in PDBQT blocks, plus additional informa
 Meeko generates RDKit molecules from PDBQT using the SMILES
 string in the REMARK lines. The REMARK lines also have the mapping of atom indices
 between SMILES and PDBQT. SD files with docked coordinates are written
-from RDKit molecules.
-
-.. code-block:: bash
-
-    mk_export.py molecule.pdbqt -s molecule.sdf
-    mk_export.py vina_results.pdbqt -s vina_results.sdf
-    mk_export.py autodock-gpu_results.dlg -s autodock-gpu_results.sdf
+from RDKit molecules. 
 
 Why this matters
 ~~~~~~~~~~~~~~~~
@@ -44,22 +48,18 @@ docking, as probably rigorous Hs positions will be replaced by the
 RDKit geometry rules, which are empirical and much simpler than most
 force fields.
 
-Usage
------
-
-.. code-block:: bash
-
-    mk_export.py [OPTIONS] docking_results_filename(s)
+Options
+-------
 
-Positional Argument
-~~~~~~~~~~~~~~~~~~~
+Positional Argument (Input)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 .. option:: docking_results_filename
 
    One or more docking output files in either PDBQT format (from Vina) or DLG format (from AD-GPU).
 
-Options
-~~~~~~~
+Output Options
+~~~~~~~~~~~~~~
 
 .. option:: --suffix <suffix>
 

From d244f5679d4eee7cbcc9716c076fe36d808ba30c Mon Sep 17 00:00:00 2001
From: diogom <diogom@scripps.edu>
Date: Thu, 7 Nov 2024 16:46:00 -0800
Subject: [PATCH 11/18] wip lig docs

---
 docs/source/cli_lig_prep.rst      | 16 ------
 docs/source/index.rst             | 19 +++++--
 docs/source/lig_prep_advanced.rst | 89 +++++++++++++++++++++++++++++++
 docs/source/lig_prep_basic.rst    | 67 +++++++++++++++++++++++
 4 files changed, 170 insertions(+), 21 deletions(-)
 delete mode 100644 docs/source/cli_lig_prep.rst
 create mode 100644 docs/source/lig_prep_advanced.rst
 create mode 100644 docs/source/lig_prep_basic.rst

diff --git a/docs/source/cli_lig_prep.rst b/docs/source/cli_lig_prep.rst
deleted file mode 100644
index 91fe3224..00000000
--- a/docs/source/cli_lig_prep.rst
+++ /dev/null
@@ -1,16 +0,0 @@
-mk_prepare_ligand.py
-====================
-
-Command line tool to prepare small organic molecules.
-
-Write PDBQT files
------------------
-
-AutoDock-GPU and Vina read molecules in the PDBQT format. These can be prepared
-by Meeko from SD files, or from Mol2 files, but SDF is strongly preferred.
-
-.. code-block:: bash
-
-    mk_prepare_ligand.py -i molecule.sdf -o molecule.pdbqt
-    mk_prepare_ligand.py -i multi_mol.sdf --multimol_outdir folder_for_pdbqt_files
-
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 2027b60f..c7aa894b 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -21,10 +21,16 @@ Python API
 
 Meeko is written in Python and exposes functions and classes that operate on
 RDKit molecules for the ligands, leveraging RDKit's popularity to facilitate
-integration with external software. Command line scripts are also available.
+integration with external software that also interfaces with RDKit.
 
-AutoDock ecosystem
-------------------
+Command line scripts
+--------------------
+
+There are three scripts: ``mk_prepare_ligand.py``, ``mk_prepare_receptor.py``,
+and ``mk_export.py``.
+
+Running a docking
+-----------------
 
 To run a docking, more packages are required besides Meeko:
 
@@ -32,6 +38,9 @@ To run a docking, more packages are required besides Meeko:
  * AutoDock-GPU
  * Ringtail
 
+Check the tutorials page to learn about using meeko with these other packages
+to run molecular docking and virtual screening.
+
 
 
 .. toctree::
@@ -49,8 +58,8 @@ To run a docking, more packages are required besides Meeko:
    :caption: Ligand preparation
 
    Overview <lig_overview>
-   cli_lig_prep
-   In Python <py_lig_prep>
+   Basic usage <lig_prep_basic>
+   Advanced usage <lig_prep_advanced>
 
 .. toctree::
    :maxdepth: 2
diff --git a/docs/source/lig_prep_advanced.rst b/docs/source/lig_prep_advanced.rst
new file mode 100644
index 00000000..aec48097
--- /dev/null
+++ b/docs/source/lig_prep_advanced.rst
@@ -0,0 +1,89 @@
+Advanced ligand preparation
+===========================
+
+Parameterization of RDKit molecules is controlled by several parameters
+that can be used when creating an instance of ``MoleculePreparation``.
+Many of such parameters are exposed to the command line script, 
+``mk_prepare_ligand.py``. Here, we cover the most relevant options, and
+how to use them both from command line and Python scripts.
+
+Configuring MoleculePreparation in Python and command line
+----------------------------------------------------------
+
+In Python, an instance of ``MoleculePreparation`` is configured by passing
+any number of optional parameters during initialization:
+
+.. code-block:: python
+
+   from meeko import MoleculePreparation
+
+   mk_prep = MoleculePreparation(
+       merge_these_atom_types=("H"),
+       charge_model="gasteiger",
+   )
+
+Here we are passing two parameters, one that will merge atoms typed "H" into
+their parent atom, and another to use Gasteiger charges. Both are the default
+parameter, so we would have gotten the same configuration by passing zero
+parameters: ``mk_prep = MoleculePreparation()``.
+
+Parameters can be in a dictionary and use the ``from_config`` constructor:
+
+.. code-block:: python
+   config_dict = {"merge_these_atom_types": (), "charge_model": "gasteiger"}
+   mk_prep = MoleculePreapration.form_config(config_dict)
+
+The configuration dictionary can be written to a JSON file:
+
+.. code-block:: python
+   import json
+   with open("my_config.json", "w") as f:
+       json.dump(config_dict, f)
+
+And the command line script ``mk_prepare_ligand.py`` can read the JSON file
+using the ``-c`` or ``--config_file`` option:
+
+.. code-block:: bash
+   mk_prepare_ligand.py -i mol.sdf -c my_config.json
+
+Here, ``mol.sdf`` is some file in the current working directory.
+Alternatively, each parameter can be passed directly:
+
+.. code-block:: bash
+   mk_prepare_ligand.py -i mol.sdf --charge_model gasteiger --merge_these_atom_types H
+
+
+
+Merging hydrogens (or not)
+--------------------------
+
+By default, atoms of type ``H`` are merged. Merging consists of adding
+their partial charge to the parent atom, and setting the ``is_ignore``
+attribute to ``True`` in the instance of ``MoleculeSetup``, which in turn
+prevents the atom from being written to PDBQT. Since the atoms are merged
+based on their atom type, changing the atom typing can change the set of
+atoms that is merged.
+
+Option name is ``merge_these_atom_types``. To prevent merging of any atom
+types set the variable to an empty tuple in Python:
+
+.. code-block:: python
+   mk_prep = MoleculePreparation(merge_these_atom_types=())
+
+or pass no parameters in command line
+.. code-block:: bash
+   mk_prepare_ligand.py -i mol.sdf --merge_these_atom_types
+
+
+
+Modifying atom types
+--------------------
+
+Atom typing relies on
+
+
+No rigid macrocycles
+--------------------
+
+Hydrated docking
+----------------
diff --git a/docs/source/lig_prep_basic.rst b/docs/source/lig_prep_basic.rst
new file mode 100644
index 00000000..19be4f4f
--- /dev/null
+++ b/docs/source/lig_prep_basic.rst
@@ -0,0 +1,67 @@
+Basic ligand preparation
+========================
+
+This covers the fundamentals both from command line and Python API.
+
+
+Command line script
+-------------------
+
+Upon installing Meeko, ``mk_prepare_ligand.py`` is available in the terminal
+and running it should print a help message. This script writes one or more
+PDBQT files for AutoDock-GPU or AutoDock-Vina. It reads SD files, often referred
+to as SDF, and these can contain multiple molecules. The input molecules must
+have all hydrogens as real atoms, and 3D positions. Mol2 files are supported,
+but the SDF is strongly preferred.
+
+To print the help message pass the ``-h`` option:
+
+.. code-block:: bash
+
+   mk_prepare_ligand.py -h
+
+
+Writing a single PDBQT file:
+
+.. code-block:: bash
+
+    mk_prepare_ligand.py -i molecule.sdf -o molecule.pdbqt
+
+Writing multiple PDBQT files from an SD file with multiple molecules:
+
+.. code-block:: bash
+
+    mk_prepare_ligand.py -i multi_mol.sdf --multimol_outdir folder_for_pdbqt_files
+
+
+Python API
+----------
+
+The following script reads a file named ``molecules.sdf`` using RDKit,
+configures an instance of ``MoleculePreparation`` with default options,
+and creates PDBQT strings for all molecules in the input file:
+
+.. code-block:: python
+
+    from meeko import MoleculePreparation
+    from meeko import PDBQTWriterLegacy
+    from rdkit import Chem
+    
+    input_filename = "molecules.sdf"
+    
+    # iterate over molecules in SD file
+    for mol in Chem.SDMolSupplier(input_filename, removeHs=False):
+        preparator = MoleculePreparation()
+        molsetup_list = preparator.prepare(mol)
+
+        # as of v0.6.0, the preparator always returns one molecule setup,
+        # unless it is configured for reactive docking, which it is not
+        # in this case. So we just grab the first (and only) element
+        molsetup = molsetup_list[0]
+
+        pdbqt_string = PDBQTWriterLegacy.write_string(molsetup)
+
+
+The ``pdbqt_string`` can be written to a file for docking with AutoDock-GPU or
+AutoDock-Vina, or passed directly to Vina within Python using Vina's Python API,
+and avoiding writing PDBQT files to the filesystem.

From b173c42e041f75f5bff4403c00d668df77eeb701 Mon Sep 17 00:00:00 2001
From: diogom <diogom@scripps.edu>
Date: Fri, 8 Nov 2024 00:10:55 -0800
Subject: [PATCH 12/18] wip more lig docs

---
 docs/source/cli_rec_prep.rst      |  6 +--
 docs/source/lig_overview.rst      | 28 ++++++++++++-
 docs/source/lig_prep_advanced.rst | 68 +++++++++++++++++++++++++++----
 docs/source/lig_prep_basic.rst    | 15 ++++---
 4 files changed, 97 insertions(+), 20 deletions(-)

diff --git a/docs/source/cli_rec_prep.rst b/docs/source/cli_rec_prep.rst
index 2fadd173..94c68603 100644
--- a/docs/source/cli_rec_prep.rst
+++ b/docs/source/cli_rec_prep.rst
@@ -10,9 +10,9 @@ files from various sources.
 
 Templates are matched on a per residue basis. Each residue is represented
 as an instance of a PolymerResidue object, which contains:
- - an RDKit molecule that represents the actual state
- - a padded RDKit molecule containing a few atoms from the adjacent residues
- - parameters such as partial charges
+ * an RDKit molecule that represents the actual state
+ * a padded RDKit molecule containing a few atoms from the adjacent residues
+ * parameters such as partial charges
 
 The positions are set by the input, and the connectivity and formal charges
 are defined by the templates. Heavy atoms must match exactly. If heavy atoms
diff --git a/docs/source/lig_overview.rst b/docs/source/lig_overview.rst
index 1fbd47ad..346cd05a 100644
--- a/docs/source/lig_overview.rst
+++ b/docs/source/lig_overview.rst
@@ -1,17 +1,41 @@
 Overview of ligand preparation
 ==============================
 
+Meeko provides both command line scripts and a Python API. Here, we briefly
+describe how data flows under the hood, and some of the key methods and data
+structures. While this information is primarily for those using Meeko from
+Python, it provides context for users of the command line script
+``mk_prepare_ligand.py``. When using this script, the input are SD files
+(.sdf extension), and the output are PDBQT files.
+
+Use of RDKit
+------------
+
 Meeko takes as input an RDKit molecule that has 3D positions and
 all hydrogens as real atoms, and creates an object called MoleculeSetup
 that stores all the parameters, such as atom types, partial charges, or
 rotatable bonds. Adding hydrogens and 3D positions is not performed by Meeko.
 
+Parameterization using SMARTS patterns
+--------------------------------------
+
+Many of the parameters are assigned using SMARTS strings, which are a compact
+and versatile query language for identigying chemical substructures. This makes
+it easier to define custom atom types and other parameters. The SMARTS
+that define the AutoDock parameters are included in Meeko and set as the
+default, This goal of this feature is to facilitate the implementation of
+custom docking workflows.
+
 The MoleculePreparation class stores all the configuration required to
-parameterize a ligand, and containes the methods that take an RDKit molecule
+parameterize a ligand, and contains the methods that take an RDKit molecule
 as input, and return a parameterized MoleculeSetup. The MoleculePreparation
 offers several option to control how molecules are parameterized, and many of
 which are exposed as command line options in ``mk_prepare_ligand.py``.
 
-PDBQT files, or strings in Python, are produced using a parameterized instance
+Preparing the input for docking
+-------------------------------
+
+AutoDock-GPU and Vina use the PDBQT format for input molecules.
+PDBQT files, or strings in Python, are produced from a parameterized instance
 of MoleculeSetup. The methods for that live in the PDBQTWriterLegacy class.
 PDBQT strings can be passed directly to Vina using its Python API.
diff --git a/docs/source/lig_prep_advanced.rst b/docs/source/lig_prep_advanced.rst
index aec48097..fceaf1fa 100644
--- a/docs/source/lig_prep_advanced.rst
+++ b/docs/source/lig_prep_advanced.rst
@@ -5,12 +5,16 @@ Parameterization of RDKit molecules is controlled by several parameters
 that can be used when creating an instance of ``MoleculePreparation``.
 Many of such parameters are exposed to the command line script, 
 ``mk_prepare_ligand.py``. Here, we cover the most relevant options, and
-how to use them both from command line and Python scripts.
+how to use them both from command line and Python scripts. Familiarity with
+Python and with the command line is required.
 
-Configuring MoleculePreparation in Python and command line
-----------------------------------------------------------
+Configuring MoleculePreparation
+-------------------------------
 
-In Python, an instance of ``MoleculePreparation`` is configured by passing
+From Python
+^^^^^^^^^^^
+
+An instance of ``MoleculePreparation`` is configured by passing
 any number of optional parameters during initialization:
 
 .. code-block:: python
@@ -30,26 +34,41 @@ parameters: ``mk_prep = MoleculePreparation()``.
 Parameters can be in a dictionary and use the ``from_config`` constructor:
 
 .. code-block:: python
+
    config_dict = {"merge_these_atom_types": (), "charge_model": "gasteiger"}
    mk_prep = MoleculePreapration.form_config(config_dict)
 
 The configuration dictionary can be written to a JSON file:
 
 .. code-block:: python
+
    import json
+
    with open("my_config.json", "w") as f:
        json.dump(config_dict, f)
 
-And the command line script ``mk_prepare_ligand.py`` can read the JSON file
-using the ``-c`` or ``--config_file`` option:
+
+From command line
+^^^^^^^^^^^^^^^^^
+
+The command line script ``mk_prepare_ligand.py`` can read a JSON file
+using the ``-c`` or ``--config_file`` option. The JSON file corresponding
+to the Python example above is:
+
+.. code-block:: json
+
+   {"merge_these_atom_types": ["H"], "charge_model": "gasteiger"}
+
 
 .. code-block:: bash
+
    mk_prepare_ligand.py -i mol.sdf -c my_config.json
 
 Here, ``mol.sdf`` is some file in the current working directory.
 Alternatively, each parameter can be passed directly:
 
 .. code-block:: bash
+
    mk_prepare_ligand.py -i mol.sdf --charge_model gasteiger --merge_these_atom_types H
 
 
@@ -68,10 +87,13 @@ Option name is ``merge_these_atom_types``. To prevent merging of any atom
 types set the variable to an empty tuple in Python:
 
 .. code-block:: python
+
    mk_prep = MoleculePreparation(merge_these_atom_types=())
 
 or pass no parameters in command line
+
 .. code-block:: bash
+
    mk_prepare_ligand.py -i mol.sdf --merge_these_atom_types
 
 
@@ -79,7 +101,39 @@ or pass no parameters in command line
 Modifying atom types
 --------------------
 
-Atom typing relies on
+Atom typing relies on SMARTS patterns to identify chemical substructures.
+AutoDock4 atom types are set by default. The easiest way to modify typing
+is to add new SMARTS that will superseed the existing ones. For example, let's
+assume we want to type hydrogens bound to aromatic carbons as ``HX``. By default,
+hydrogens bound to carbon are typed ``H``. A SMARTS pattern to matches
+hydrogen bound to carbon is ``"[H][c]"``. From command line:
+
+.. code-block:: bash
+
+   mk_prepare_ligand.py --add_atom_types '[{"smarts": "[H]c", "atype": "HX"}]' -i mol.sdf
+
+We pass a JSON string to ``--add_atom_types`` that is a list of dictionaries. Each
+dictionary has a ``"smarts"`` and ``"atype"`` key, and an optional ``"IDX"`` key
+that can be used to specify a list of atom indices (0-based) of the atoms in the SMARTS
+string that will be typed. By default ``IDX = [0]``.
+
+The equivalent from Python is:
+
+.. code-block:: python
+
+   mk_prep = MoleculePreparation(
+       add_atom_types=[{"smarts": "c[H]", "atype": "HX", "IDX": [1]}],
+   )
+
+Note that we swapped the order of the atoms in the SMARTS, and are now
+explicitly defining the ``"IDX"`` key to type the second atom in the SMARTS.
+
+The full set of atom types can also be specified. This can only be done from
+Python or by passing the equivalent configuration JSON file to ``mk_prepare_ligand.py``.
+The easiest way to do so, is to put all SMARTS in a JSON file. See the default
+file for an example, it is located at ``meeko/data/params/ad4_types.json``.
+
+
 
 
 No rigid macrocycles
diff --git a/docs/source/lig_prep_basic.rst b/docs/source/lig_prep_basic.rst
index 19be4f4f..98de381c 100644
--- a/docs/source/lig_prep_basic.rst
+++ b/docs/source/lig_prep_basic.rst
@@ -1,9 +1,6 @@
 Basic ligand preparation
 ========================
 
-This covers the fundamentals both from command line and Python API.
-
-
 Command line script
 -------------------
 
@@ -51,12 +48,9 @@ and creates PDBQT strings for all molecules in the input file:
     
     # iterate over molecules in SD file
     for mol in Chem.SDMolSupplier(input_filename, removeHs=False):
-        preparator = MoleculePreparation()
-        molsetup_list = preparator.prepare(mol)
+        mk_prep = MoleculePreparation()
+        molsetup_list = mk_prep(mol)
 
-        # as of v0.6.0, the preparator always returns one molecule setup,
-        # unless it is configured for reactive docking, which it is not
-        # in this case. So we just grab the first (and only) element
         molsetup = molsetup_list[0]
 
         pdbqt_string = PDBQTWriterLegacy.write_string(molsetup)
@@ -65,3 +59,8 @@ and creates PDBQT strings for all molecules in the input file:
 The ``pdbqt_string`` can be written to a file for docking with AutoDock-GPU or
 AutoDock-Vina, or passed directly to Vina within Python using Vina's Python API,
 and avoiding writing PDBQT files to the filesystem.
+
+Note that calling ``mk_prep`` returns a list of molecule setups.
+As of v0.6.0, this list constains only one element  unless ``mk_prep`` is
+configured for reactive docking, which is not the case in this example. This is
+why we are considering the first (and only) molecule setup in the list.

From 689512cc438ea57c35917ce5debc9b90ea013708 Mon Sep 17 00:00:00 2001
From: diogom <diogom@scripps.edu>
Date: Fri, 8 Nov 2024 15:30:40 -0800
Subject: [PATCH 13/18] address some sphinx docs warnigns

---
 docs/source/cli_rec_prep.rst | 7 ++++---
 docs/source/tutorial3.rst    | 6 +++---
 2 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/docs/source/cli_rec_prep.rst b/docs/source/cli_rec_prep.rst
index a58eda9a..10796405 100644
--- a/docs/source/cli_rec_prep.rst
+++ b/docs/source/cli_rec_prep.rst
@@ -33,9 +33,10 @@ files from various sources.
 
 Templates are matched on a per residue basis. Each residue is represented
 as an instance of a PolymerResidue object, which contains:
- - an RDKit molecule that represents the actual state
- - a padded RDKit molecule containing a few atoms from the adjacent residues
- - parameters such as partial charges
+
+ * an RDKit molecule that represents the actual state
+ * a padded RDKit molecule containing a few atoms from the adjacent residues
+ * parameters such as partial charges
 
 The positions are set by the input, and the connectivity and formal charges
 are defined by the templates. Heavy atoms must match exactly. If heavy atoms
diff --git a/docs/source/tutorial3.rst b/docs/source/tutorial3.rst
index c830324c..49ab6601 100644
--- a/docs/source/tutorial3.rst
+++ b/docs/source/tutorial3.rst
@@ -46,7 +46,7 @@ The ligand of this example will be the covalent conjugate ``HIE_AMP``, where AMP
 To prepare HIE-AMP (1-) as an covalent flexible residue, we will hold this SDF file for further mapping with the specific catalytic residue in receptor structure. In fact, the SDF file can be re-used for different Histidine residues in different receptor structures. 
 
 Receptor Preparation
-===================
+====================
 
 The preparation of a rigid receptor consists of two steps. The receptor structure is first sourced from a PDB file and sent to ``reduce2.py`` for hydrogen addition and optimization, and then, the conversion to a tangible receptor PDBQT file is done by ``mk_prepare_receptor.py``.
 
@@ -139,7 +139,7 @@ For output control: We are expecting at least two types of files, the receptor P
         3kgd_receptorH.box.pdb <-- PDB file to visualize the grid box
 
 Covalent Ligand Preparation
-========================
+===========================
 
 In this step, we will use mk_prepare_ligand.py to generate the PDBQT file for the covalent ligand. Along with the previously generated 3D conformer of the covalent ligand (``HIE_AMP.sdf``), which may be at an arbitrary position, here a reference protein PDB file (``3kgd_receptor.pdb``) will be used to source the positions of the attractor atoms, Cα and Cβ, to keep them unchanged in docking. The reference PDB file does not have to be the full receptor, but it must contain the target residue that matches exactly with ``rec_residue``. Additionally, a SMARTS pattern ``tether_smarts`` is required. Together with the 1-based ``tether_smarts_indices``, they are used to locate the attractor atoms that correspond to Cα and Cβ of a hisitidine (His) residue:
 
@@ -193,4 +193,4 @@ It is also possible to export the docking poses to a multi-model PDB file with u
    --default_altloc A -f $flexres \
    --box_enveloping "LIG.pdb" --padding 8.0 
 
-   mk_export.py HIE_AMP.dlg -s 3kgd_HIE_AMP_adgpu_out.sdf -j 3kgd_receptorH.json -p 3kgd_HIE_AMP_adgpu_out.pdb
\ No newline at end of file
+   mk_export.py HIE_AMP.dlg -s 3kgd_HIE_AMP_adgpu_out.sdf -j 3kgd_receptorH.json -p 3kgd_HIE_AMP_adgpu_out.pdb

From c0cde3f64dd35b52225c784939698f1a96ccd648 Mon Sep 17 00:00:00 2001
From: diogom <diogom@scripps.edu>
Date: Sun, 10 Nov 2024 10:40:34 -0800
Subject: [PATCH 14/18] done with lig docs, changed some idx to 0-based

---
 docs/source/index.rst                         |  5 +-
 docs/source/lig_cli_options.rst               |  2 +-
 docs/source/lig_prep_advanced.rst             | 47 ++++++++++++++++---
 docs/source/lig_prep_basic.rst                |  6 +++
 meeko/atomtyper.py                            |  8 ++--
 meeko/bondtyper.py                            |  4 +-
 meeko/cli/mk_prepare_ligand.py                | 11 ++---
 meeko/data/params/example_offatom_charge.json |  4 +-
 meeko/openff_xml_parser.py                    |  6 +--
 9 files changed, 67 insertions(+), 26 deletions(-)

diff --git a/docs/source/index.rst b/docs/source/index.rst
index 0bc40b15..8f2a4165 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -26,8 +26,9 @@ integration with external software that also interfaces with RDKit.
 Command line scripts
 --------------------
 
-There are three scripts: ``mk_prepare_ligand.py``, ``mk_prepare_receptor.py``,
-and ``mk_export.py``.
+There are two scripts to write the input files for docking,
+``mk_prepare_ligand.py`` and ``mk_prepare_receptor.py``, and one script to
+convert docking output files  ``mk_export.py``.
 
 Running a docking
 -----------------
diff --git a/docs/source/lig_cli_options.rst b/docs/source/lig_cli_options.rst
index 760bf78a..eff62399 100644
--- a/docs/source/lig_cli_options.rst
+++ b/docs/source/lig_cli_options.rst
@@ -83,7 +83,7 @@ Molecule Preparation Options
 
    (Flag) Allow amide bonds to rotate, making them non-planar (not recommended).
 
-.. option:: -p, --atom_type_smarts <JSON_FILENAME>
+.. option:: -p, --load_atom_params <JSON_FILENAME>
 
    Specify SMARTS-based atom typing in JSON format.
 
diff --git a/docs/source/lig_prep_advanced.rst b/docs/source/lig_prep_advanced.rst
index fceaf1fa..08c64276 100644
--- a/docs/source/lig_prep_advanced.rst
+++ b/docs/source/lig_prep_advanced.rst
@@ -52,19 +52,21 @@ From command line
 ^^^^^^^^^^^^^^^^^
 
 The command line script ``mk_prepare_ligand.py`` can read a JSON file
-using the ``-c`` or ``--config_file`` option. The JSON file corresponding
-to the Python example above is:
+using the ``-c`` or ``--config_file`` option. The contents of a JSON file
+corresponding to the Python example above is:
 
 .. code-block:: json
 
    {"merge_these_atom_types": ["H"], "charge_model": "gasteiger"}
 
+Passing the filename to the script with option ``-c``:
 
 .. code-block:: bash
 
    mk_prepare_ligand.py -i mol.sdf -c my_config.json
 
 Here, ``mol.sdf`` is some file in the current working directory.
+
 Alternatively, each parameter can be passed directly:
 
 .. code-block:: bash
@@ -132,12 +134,45 @@ The full set of atom types can also be specified. This can only be done from
 Python or by passing the equivalent configuration JSON file to ``mk_prepare_ligand.py``.
 The easiest way to do so, is to put all SMARTS in a JSON file. See the default
 file for an example, it is located at ``meeko/data/params/ad4_types.json``.
+The ``IDX`` key can be used as described above. Entries are matched in the
+order they appear in the file, the last SMARTS pattern that matches an atom is
+the one that determines the atom type.
+Then the filename can be passed to option ``-p/--load_atom_params``.
 
 
+Rigidifying bonds
+-----------------
 
+By default, single bonds are made rotatable except bonds in rings and amide bonds.
+Thioamide and amidine bonds are also not rotatable.
+Tertiary amides with non-equivalent substituents on the nitrogen are still made
+rotatable, which often leads to unreasonable geometries, but is necessary to
+visit both amide rotamers during docking.
 
-No rigid macrocycles
---------------------
+Here, we configure Meeko to make single bonds in some conjugated systems rigid,
+as defined byt the SMARTS ``"C=CC=C"``, and rigidify all amide bonds matched
+by ``"[CX3](=O)[NX3]"``, which includes tertiary amides but not thioamides or
+amidines:
+
+.. code-block:: bash
+
+   mk_prepare_ligand.py\
+     --rigidify_bonds_smarts "C=CC=C"\
+     --rigidify_bonds_indices 2 3\
+     --rigidify_bonds_smarts "[CX3](=O)[NX3]"\
+     --rigidify_bonds_indices 1 3\
+     -i mol.sdf
 
-Hydrated docking
-----------------
+The equivalent code in Python to initialize the molecule preparator is: 
+
+.. code-block:: python
+
+   mk_prep = MoleculePreparation(
+       rigidify_bonds_smarts = ["C=CC=C", "[CX3](=O)[NX3]"],
+       rigidify_bonds_indices = [(1, 2), (0, 2)],
+)
+
+The indices are the indices of the atoms in the SMARTS strings. Note that
+we use 0-based indices from the Python API, but 1-based indices from the
+command line script. In a future version of Meeko we may use 0-based indices
+everywhere.
diff --git a/docs/source/lig_prep_basic.rst b/docs/source/lig_prep_basic.rst
index 98de381c..d44eefab 100644
--- a/docs/source/lig_prep_basic.rst
+++ b/docs/source/lig_prep_basic.rst
@@ -24,12 +24,18 @@ Writing a single PDBQT file:
 
     mk_prepare_ligand.py -i molecule.sdf -o molecule.pdbqt
 
+If the ``-o`` option is omitted, the output filename will be the same as the
+input but with ``.pdbqt`` extension. Option ``-`` prints the PDBQT
+string to standard output instead of writting to a file.
+
 Writing multiple PDBQT files from an SD file with multiple molecules:
 
 .. code-block:: bash
 
     mk_prepare_ligand.py -i multi_mol.sdf --multimol_outdir folder_for_pdbqt_files
 
+Optionally, passing ``-z`` or ``--multimol_targz`` will compress the output as
+``.tar.gz`` files with 10000 PDBQT files each.
 
 Python API
 ----------
diff --git a/meeko/atomtyper.py b/meeko/atomtyper.py
index 12409adb..14a15693 100644
--- a/meeko/atomtyper.py
+++ b/meeko/atomtyper.py
@@ -50,7 +50,7 @@ def _type_atoms(molsetup, atom_params):
                     0
                 ]  # by default, the first atom in the smarts gets parameterized
                 if "IDX" in line:
-                    idxs = [i - 1 for i in line["IDX"]]  # convert from 1- to 0-indexing
+                    idxs = [i for i in line["IDX"]]
                 # match SMARTS
                 hits = molsetup.find_pattern(smarts)
                 for atompar in line:
@@ -103,7 +103,7 @@ def _cache_offatoms(molsetup, offatom_params):
                 # atom indexes in smarts string
                 smarts_idxs = [0]
                 if "IDX" in line:
-                    smarts_idxs = [i - 1 for i in line["IDX"]]
+                    smarts_idxs = [i for i in line["IDX"]]
                 for smarts_idx in smarts_idxs:
                     for hit in hits:
                         parent_idx = hit[smarts_idx]
@@ -131,7 +131,7 @@ def _cache_offatoms(molsetup, offatom_params):
                                 # replace SMARTS indexes by the atomic index
                                 elif key in ["z", "x"]:
                                     for i in offatom[key]:
-                                        idx = hit[i - 1]
+                                        idx = hit[i]
                                         tmp[parent_idx][-1]["offatom"][key].append(idx)
                                 # convert degrees to radians
                                 elif key in ["theta", "phi"]:
@@ -217,7 +217,7 @@ def _type_dihedrals(molsetup, dihedral_params):
             hits = molsetup.find_pattern(smarts)
             if len(hits) == 0:
                 continue  # non-rotatable bonds get dihedrals
-            idxs = [i - 1 for i in line["IDX"]]
+            idxs = [i for i in line["IDX"]]
             tid = line["id"] if "id" in line else None
             fourier_series = []
             term_indices = {}
diff --git a/meeko/bondtyper.py b/meeko/bondtyper.py
index 792b0542..08d4bfc9 100644
--- a/meeko/bondtyper.py
+++ b/meeko/bondtyper.py
@@ -24,12 +24,12 @@ def __call__(
         """
 
         amide_bonds = [
-            (x[0], x[1]) for x in setup.find_pattern("[NX3]-[CX3]=[O,N]")
+            (x[0], x[1]) for x in setup.find_pattern("[NX3]-[CX3]=[O,N,S]")
         ]  # includes amidines
 
         # tertiary amides with non-identical substituents will be allowed to rotate
         tertiary_amides = [
-            x for x in setup.find_pattern("[NX3]([!#1])([!#1])-[CX3]=[O,N]")
+            x for x in setup.find_pattern("[NX3]([!#1])([!#1])-[CX3]=[O,N,S]")
         ]
         equivalent_atoms = setup.get_equivalent_atoms()
         num_amides_removed = 0
diff --git a/meeko/cli/mk_prepare_ligand.py b/meeko/cli/mk_prepare_ligand.py
index a9de8d28..2785fa24 100755
--- a/meeko/cli/mk_prepare_ligand.py
+++ b/meeko/cli/mk_prepare_ligand.py
@@ -192,10 +192,11 @@ def cmd_lineparser():
     )
     config_group.add_argument(
         "-p",
-        "--atom_type_smarts",
-        dest="atom_type_smarts_json",
+        "--load_atom_params",
+        nargs="+",
         action="store",
-        help="SMARTS based atom typing (JSON format)",
+        default=["ad4_types"],
+        help="filename with SMARTS defined atom types (JSON format)",
     )
     config_group.add_argument(
         "-aa",
@@ -296,9 +297,7 @@ def cmd_lineparser():
         if key in args.__dict__:
             config[key] = args.__dict__[key]
 
-    if args.atom_type_smarts_json is not None:
-        with open(args.atom_type_smarts_json) as f:
-            config["atom_type_smarts"] = json.load(f)
+    config["load_atom_params"] = args.load_atom_params
 
     if args.add_atom_types_json is not None:
         additional_ats = []
diff --git a/meeko/data/params/example_offatom_charge.json b/meeko/data/params/example_offatom_charge.json
index f906e0e5..a11f07e7 100644
--- a/meeko/data/params/example_offatom_charge.json
+++ b/meeko/data/params/example_offatom_charge.json
@@ -1,7 +1,7 @@
 {"example_offsite_charges": [
-        {"smarts": "[#7X2;v3;!+](=,:[*])[*]", "IDX": [1],
+        {"smarts": "[#7X2;v3;!+](=,:[*])[*]", "IDX": [0],
          "OFFATOMS": [
-            {"z": [2, 3], "phi": 0, "distance": 0.2,
+            {"z": [1, 2], "phi": 0, "distance": 0.2,
             "atype": "OFFCHRG", "pull_charge_fraction": 1.08}
             ]
         }
diff --git a/meeko/openff_xml_parser.py b/meeko/openff_xml_parser.py
index 6a30e7da..5e079d36 100644
--- a/meeko/openff_xml_parser.py
+++ b/meeko/openff_xml_parser.py
@@ -137,7 +137,7 @@ def make_dihedral_entry(attrib_dict_from_xml):
     assert 1 in labels and 2 in labels and 3 in labels and 4 in labels
     dihedral_entry = {
         "smarts": smarts,
-        "IDX": [labels[index] + 1 for index in (1, 2, 3, 4)],
+        "IDX": [labels[index] for index in (1, 2, 3, 4)],
     }
     terms = {}  # keywords found for each term of the fourier series
     expected_keywords = set(["k", "phase", "periodicity", "idivf"])
@@ -185,7 +185,7 @@ def make_vdw_entry(attrib_dict_from_xml):
     assert 1 in labels
     vdw_entry = {
         "smarts": smarts,
-        "IDX": [labels[1] + 1],
+        "IDX": [labels[1]],
         "id": attrib_dict_from_xml["id"],
     }
     assert ("epsilon" in attrib_dict_from_xml) and (
@@ -220,7 +220,7 @@ def assign_atypes(vdw_list, use_openff_id=True, force_uppercase=True):
     for v in vdw_list:
         mol = Chem.MolFromSmarts(v["smarts"])
         atom = mol.GetAtomWithIdx(
-            v["IDX"][0] - 1
+            v["IDX"][0]
         )  # consider only the first if multiple IDX
         element = mini_periodic_table[atom.GetAtomicNum()]
         atomic_numbers.append(atom.GetAtomicNum())

From 469a697fba5bbce6c5fb315d9b099927e40d5b99 Mon Sep 17 00:00:00 2001
From: diogom <diogom@scripps.edu>
Date: Mon, 11 Nov 2024 01:47:00 -0800
Subject: [PATCH 15/18] work on docs

---
 docs/source/cli_export_result.rst  |  92 --------------
 docs/source/cli_rec_prep.rst       | 192 -----------------------------
 docs/source/export_cli_options.rst |  43 +++++++
 docs/source/export_usage.rst       | 117 ++++++++++++++++++
 docs/source/index.rst              |  13 +-
 docs/source/py_export_result.rst   |  37 ------
 docs/source/py_lig_prep.rst        |  32 -----
 docs/source/rec_cli_options.rst    | 130 +++++++++++++++++++
 docs/source/rec_overview.rst       |  73 +++++++++++
 9 files changed, 369 insertions(+), 360 deletions(-)
 delete mode 100644 docs/source/cli_export_result.rst
 create mode 100644 docs/source/export_cli_options.rst
 create mode 100644 docs/source/export_usage.rst
 delete mode 100644 docs/source/py_export_result.rst
 delete mode 100644 docs/source/py_lig_prep.rst
 create mode 100644 docs/source/rec_cli_options.rst
 create mode 100644 docs/source/rec_overview.rst

diff --git a/docs/source/cli_export_result.rst b/docs/source/cli_export_result.rst
deleted file mode 100644
index 3a35d263..00000000
--- a/docs/source/cli_export_result.rst
+++ /dev/null
@@ -1,92 +0,0 @@
-mk_export.py
-============
-
-A command-line script to export docking poses of ligand (and flexible residues) to an SD file (.sdf), and to export the full receptor structures with updated conformations of flexible residues to a PDB file. Currently supports PDBQT files from AutoDock-Vina and DLG files from AutoDock-GPU that have the REMARK lines containing Smiles and index mapping information. 
-
-Basic usage
------------
-
-.. code-block:: bash
-
-    mk_export.py vina_results.pdbqt -s vina_results.sdf
-    mk_export.py autodock-gpu_results.dlg -s autodock-gpu_results.sdf
-
-About
------
-
-Convert docking results to SDF
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-AutoDock-GPU and Vina write docking results in the PDBQT format. The DLG output
-from AutoDock-GPU contains docked poses in PDBQT blocks, plus additional information.
-Meeko generates RDKit molecules from PDBQT using the SMILES
-string in the REMARK lines. The REMARK lines also have the mapping of atom indices
-between SMILES and PDBQT. SD files with docked coordinates are written
-from RDKit molecules. 
-
-Why this matters
-~~~~~~~~~~~~~~~~
-
-Making RDKit molecules from SMILES is safer than guessing bond orders
-from the coordinates, specially because the PDBQT lacks hydrogens bonded
-to carbon. As an example, consider the following conversion, in which
-OpenBabel adds an extra double bond, not because it has a bad algorithm,
-but because this is a nearly impossible task.
-
-.. code-block:: bash
-
-    obabel -:"C1C=CCO1" -o pdbqt --gen3d | obabel -i pdbqt -o smi
-    [C]1=[C][C]=[C]O1
-
-Caveats
-~~~~~~~
-
-If docking does not use explicit Hs, which it often does not, the
-exported positions of hydrogens are calculated from RDKit. This can
-be annoying if a careful forcefield minimization is employed before
-docking, as probably rigorous Hs positions will be replaced by the
-RDKit geometry rules, which are empirical and much simpler than most
-force fields.
-
-Options
--------
-
-Positional Argument (Input)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-.. option:: docking_results_filename
-
-   One or more docking output files in either PDBQT format (from Vina) or DLG format (from AD-GPU).
-
-Output Options
-~~~~~~~~~~~~~~
-
-.. option:: --suffix <suffix>
-
-   Set a suffix for output filenames that are not explicitly specified. The default suffix is ``_docked``. 
-
-.. option:: -s, --write_sdf <output_SDF_filename>
-
-   Specify the output SDF filename. Defaults to the input filename with a suffix defined by ``--suffix``. 
-
-.. option:: -j, --read_json <input_JSON_filename>
-
-   Provide a receptor file generated by ``mk_prepare_receptor`` with the ``-j/--write_json`` option. Currently only effective when used with ``-p, --write_pdb``. 
-
-.. option:: -p, --write_pdb <output_PDB_filename>
-
-   Specify the output PDB filename. Defaults to the input filename with a suffix defined by ``--suffix``. Must be used together with ``-j, --read_json``. 
-
-.. option:: --all_dlg_poses
-
-   (Flag) Write all poses from AutoDock-GPU DLG output files, instead of only the lead of each cluster. Currently only effective for ``-s, --write_sdf``. 
-
-.. option:: -k, --keep_flexres_sdf
-
-   (Flag) Include flexible residues, if any, in the SDF output.
-
-.. option:: -, --redirect_stdout
-
-   (Flag) Instead of writing an SDF file, print it directly to the standard output (STDOUT).
-
-
diff --git a/docs/source/cli_rec_prep.rst b/docs/source/cli_rec_prep.rst
index 10796405..823081fb 100644
--- a/docs/source/cli_rec_prep.rst
+++ b/docs/source/cli_rec_prep.rst
@@ -18,67 +18,6 @@ This is equivalent to:
 
 Read more about the syntax in :ref:`Write flags` and :ref:`Options`. 
 
-About
------
-
-The input structure is matched against templates to
-guarantee chemical correctness and identify problems with the input structures.
-This allows the user to identify and fix problems, resulting in a molecular
-model that is correct with respect to heavy atoms, protonation state,
-connectivity, bond orders, and formal charges.
-
-The matching algorithm uses the connectivity and elements, but not bond orders
-or atom names. Hydrogens are optional. This makes it compatible with input
-files from various sources.
-
-Templates are matched on a per residue basis. Each residue is represented
-as an instance of a PolymerResidue object, which contains:
-
- * an RDKit molecule that represents the actual state
- * a padded RDKit molecule containing a few atoms from the adjacent residues
- * parameters such as partial charges
-
-The positions are set by the input, and the connectivity and formal charges
-are defined by the templates. Heavy atoms must match exactly. If heavy atoms
-are missing or in excess, the templates will fail to match.
-
-Missing hydrogens are added by RDKit, but are not subjected to minimization
-with a force field. Thus, their bond lengths are not super accurate.
-
-Different states of the same residue are stored as different templates,
-for example different protonation states of HIS, N-term, LYN/LYS, etc.
-Residue name is primary key unless user overrides.
-
-Currently not supported: capped residues from charmm-gui.
-
-.. _templates:
-
-Templates
-~~~~~~~~~
-
-The templates contain SMILES strings that are used to create the RDKit
-molecules that constitute every residue in the processed model. In this way,
-the chemistry of the processed model is fully defined by the templates,
-and the only thing that is preserved from the input are the atom positions
-and the connectivity between residues.
-
-The SMILES strings contain all atoms that exist in the final model,
-and none that do not exist. This also applies to hydrogens,
-meaning that the SMILES are expected to have real hydrogens. Note that
-real hydrogens are different from explicit hydrogens. Real hydrogens will be
-represented as an actual atom in an RDKit molecule, while explicit hydrogens
-are a just property of heavy atoms. In the SMILES, real hydrogens are defined
-with square brackets "[H]" and explicit hydrogens without, e.g. "[nH]" to set
-the number of explicit hydrogens on an aromatic nitrogen to one.
-
-Residues that are part of a polymer, which is often all of them, will have
-bonds to adjacent residues. The heavy atoms involved in the bonds will miss
-a real hydrogen and have an implicit (or explicit) one instead. As an
-example, consider modeling an alkyl chain as a polymer, in which the monomer
-is a single carbon atom. Our template SMILES would be "[H]C[H]". The RDKit
-molecule will have three atoms and the carbon will have two implicit hydrogens.
-The implicit hydrogens correspond to bonds to adjacent residues in the
-processed polymer. 
 
 Write flags
 ~~~~~~~~~~~
@@ -136,134 +75,3 @@ The arguments involving assignment of residues to properties:
 use the residue selection lanaguge described above, followed by an equal sign (``=``) as the delimiter and the assigned value, which could be the name of a residue template, the atom index for the blunt end, the wanted altloc ID, or the atom name of the reactive atom. Each residue selection is comibned with the most recent assignment that precedes it, resulting in a further expanded list of residue-assignment pairs. 
 
 For an input like ``"A:5,7=CYX,A:19A,B:17=HID``, this assignment language represents: ``residues (number) 5 in Chain A are set to (template name) CYX`` and ``residue (number) 19 A in Chain A, and residue (number) 17 in Chain B are set to (template name) HID``. 
-
-Options
--------
-
-Input/Output Options
-~~~~~~~~~~~~~~~~~~~~
-
-.. option:: --read_pdb <PDB_FILENAME>
-
-   Read a PDB file using the PDB parser in RDKit.
-
-.. option:: -i, --read_with_prody <MACROMOL_FILENAME>
-
-   Read a PDB or mmCIF file using ProDy (if installed). ProDy can be installed from PyPI or conda-forge.
-
-.. option:: -o, --output_basename <basename>
-
-   Specify a default basename for output files created by `--write` options when no filename is specified.
-
-.. option:: -p, --write_pdbqt <PDBQT_FILENAME> [*]
-
-   Output PDBQT files with `_rigid` or `_flex` suffixes for flexible residues. Defaults to `--output_basename` if no filename is provided.
-
-.. option:: -j, --write_json <JSON_FILENAME> [*]
-
-   Save the receptor's parameterized configuration to JSON format. Defaults to `--output_basename` if unspecified.
-
-.. option:: -g, --write_gpf <GPF_FILENAME> [*]
-
-   Output an AutoGrid input file (GPF). Defaults to `--output_basename` if not specified.
-
-.. option:: -v, --write_vina_box <VINA_BOX_FILENAME> [*]
-
-   Generate a configuration file for Vina with grid box dimensions. Defaults to `--output_basename` if not specified.
-
-.. option:: --write_pdb <PDB_FILENAME> [*]
-
-   Save the prepared receptor in PDB format. Must specify the filename.
-
-Receptor Perception Options
-~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-.. option:: -n, --set_template <template>
-
-   Assign templates to residues, e.g., `A:5,7=CYX,B:17=HID`.
-
-.. option:: -d, --delete_residues <residues>
-
-   Specify residues to delete, e.g., `A:350,B:15,16,17`.
-
-.. option:: -b, --blunt_ends <positions>
-
-   Blunt end definitions, e.g., `A:123,200=2,A:1=0`.
-
-.. option:: --add_templates <JSON_FILENAME> [*]
-
-   Load additional templates from one or more JSON files.
-
-.. option:: -a, --allow_bad_res
-
-   (Flag) Ignore residues with missing atoms instead of raising an error.
-
-.. option:: --default_altloc <location>
-
-   Define a default alternate location for residues, overridden by `--wanted_altloc`.
-
-.. option:: --wanted_altloc <location>
-
-   Specify alternate locations for particular residues, e.g., `:5=B,B:17=A`.
-
-.. option:: --mk_config <JSON_FILENAME>
-
-   Specify a JSON configuration file for receptor preparation.
-
-Grid Box Options
-~~~~~~~~~~~~~~~~
-
-.. option:: --box_size <X Y Z>
-
-   Set the size of the grid box in Angstroms (x, y, z).
-
-.. option:: --box_center <X Y Z>
-
-   Define the center of the grid box in Angstroms (x, y, z).
-
-.. option:: --box_center_off_reactive_res
-
-   (Flag) Shift the grid box center 5 Å along the CA-CB bond from CB. Applicable only when there is one reactive flexible residue.
-
-.. option:: --box_enveloping <FILENAME>
-
-   Adjust the grid box to enclose atoms in the specified file. Supported formats: `.sdf`, `.mol`, `.mol2`, `.pdb`, `.pdbqt`.
-
-.. option:: --padding <value>
-
-   Set padding around atoms specified in `--box_enveloping` (in Å).
-
-Flexible and/or Reactive Options
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-.. option:: -f, --flexres <residues>
-
-   Define flexible residues by chain ID and residue number, e.g., `-f ":42,B:23"`. 
-
-.. option:: -r, --reactive_flexres <residues>
-
-   Define reactive flexible residues by chain ID and residue number, e.g., `-r ":42,B:23"`. Maximum of 8 reactive residues.
-
-.. option:: --reactive_name <residue:atom>
-
-   Specify the reactive atom name for a residue type, e.g., `--reactive_name "TRP:NE1"`. Can be repeated for multiple assignments.
-
-.. option:: -s, --reactive_name_specific <residue:atom>
-
-   Specify the reactive atom for an individual residue by residue ID, e.g., `-s "A:42=NE2"`. The residue becomes reactive.
-
-.. option:: --r_eq_12 <value>
-
-   Set the equilibrium distance (r_eq) for reactive atoms in 1-2 interactions. Default is 1.8 Å.
-
-.. option:: --eps_12 <value>
-
-   Set the epsilon value for reactive atoms in 1-2 interactions. Default is 2.5.
-
-.. option:: --r_eq_13_scaling <factor>
-
-   Scale r_eq for 1-3 interactions across reactive atoms. Default is 0.5.
-
-.. option:: --r_eq_14_scaling <factor>
-
-   Scale r_eq for 1-4 interactions across reactive atoms. Default is 0.5.
diff --git a/docs/source/export_cli_options.rst b/docs/source/export_cli_options.rst
new file mode 100644
index 00000000..e2725adc
--- /dev/null
+++ b/docs/source/export_cli_options.rst
@@ -0,0 +1,43 @@
+
+Options of mk_export.py
+=======================
+
+Positional Argument (Input)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. option:: docking_results_filename
+
+   One or more docking output files in either PDBQT format (from Vina) or DLG format (from AD-GPU).
+
+Output Options
+~~~~~~~~~~~~~~
+
+.. option:: --suffix <suffix>
+
+   Set a suffix for output filenames that are not explicitly specified. The default suffix is ``_docked``. 
+
+.. option:: -s, --write_sdf <output_SDF_filename>
+
+   Specify the output SDF filename. Defaults to the input filename with a suffix defined by ``--suffix``. 
+
+.. option:: -j, --read_json <input_JSON_filename>
+
+   Provide a receptor file generated by ``mk_prepare_receptor`` with the ``-j/--write_json`` option. Currently only effective when used with ``-p, --write_pdb``. 
+
+.. option:: -p, --write_pdb <output_PDB_filename>
+
+   Specify the output PDB filename. Defaults to the input filename with a suffix defined by ``--suffix``. Must be used together with ``-j, --read_json``. 
+
+.. option:: --all_dlg_poses
+
+   (Flag) Write all poses from AutoDock-GPU DLG output files, instead of only the lead of each cluster. Currently only effective for ``-s, --write_sdf``. 
+
+.. option:: -k, --keep_flexres_sdf
+
+   (Flag) Include flexible residues, if any, in the SDF output.
+
+.. option:: -, --redirect_stdout
+
+   (Flag) Instead of writing an SDF file, print it directly to the standard output (STDOUT).
+
+
diff --git a/docs/source/export_usage.rst b/docs/source/export_usage.rst
new file mode 100644
index 00000000..ae882ec1
--- /dev/null
+++ b/docs/source/export_usage.rst
@@ -0,0 +1,117 @@
+Exporting docking results
+=========================
+
+AutoDock writes results in PDBQT format, which lacks bond orders and hydrogens
+bonded to carbon. When other software reads PDBQT files, it needs to infer
+bond orders, which can be impossible to do accurately for some molecules.
+As an example, consider the following conversion, in which OpenBabel adds an
+extra double bond:
+
+.. code-block:: bash
+
+    obabel -:"C1C=CCO1" -o pdbqt --gen3d | obabel -i pdbqt -o smi
+    [C]1=[C][C]=[C]O1
+
+Meeko can export docking results to SDF while preserving bond orders because
+when it writes the PDBQT file for a ligand (as input for docking) it writes
+a SMILES string and index mapping information in remark lines.
+Both AutoDock-GPU and AutoDock-Vina preserve
+the remarks in the output files, and Meeko then uses them to create an accurate
+RDKit molecule and write an SD file.
+
+When docking with flexible sidechains,
+Meeko can write a PDB file for the entire receptor with the modified
+sidechain positions, or include the sidechains as additional fragments in
+the SDF for the ligand. 
+
+Since hydrogens bonded to carbon are excluded by default in AutoDock,
+exported positions of these hydrogens are calculated by RDKit. This can
+be annoying if a careful forcefield minimization is employed before
+docking, as probably rigorous Hs positions will be replaced by the
+RDKit geometry rules, which are empirical and much simpler than most
+force fields.
+
+
+Command line usage
+------------------
+
+The script reads vina output files (``.pdbqt`` extension):
+
+.. code-block:: bash
+
+    mk_export.py vina_results.pdbqt -s vina_results.sdf
+
+And AutoDock-GPU output files (``.dlg`` extension):
+
+.. code-block:: bash
+
+    mk_export.py autodock-gpu_results.dlg -s autodock-gpu_results.sdf
+
+It can also read files zipped with gzip (``.gz`` extension), and multiple
+filenames can be passed, for example using wildcards. In this case, the output
+file names will be the input filenames with an added suffix. The suffix is
+``_docked`` by default but can be modified with the ``--suffix`` flag:
+
+.. code-block:: bash
+
+   mk_export.py docking_results/*.dlg.gz --suffix _docking_result
+
+To write a PDB file with the entire receptor including modified sidechain
+positions, the prepared receptor written by ``mk_prepare_receptor.py`` or by
+the equivalent Python code is required using the ``-j/--read_json`` option,
+and the ``-p/--write_pdb`` option sets the filename of the output PDB file:
+
+.. code-block:: bash
+
+   mk_export.py results.pdbqt -s results.sdf -p results.pdb -j rec.json
+
+To write flexible sidechains as additional fragments to the ligand SDF,
+use the ``-k/--keep_flexres_sdf`` flag.
+
+
+In Python
+---------
+
+This example converts an output file from AutoDock-GPU into a list of
+RDKit molecules, each with multiple conformers corresponding to a different
+docked pose:
+
+.. code-block:: python
+   
+   from meeko import PDBQTMolecule
+   from meeko import RDKitMolCreate
+   
+   fn = "autodock-gpu_results.dlg"
+   pdbqt_mol = PDBQTMolecule.from_file(fn, is_dlg=True, skip_typing=True)
+   rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(
+       pdbqt_mol,
+       only_cluster_leads=True,
+       keep_flexres=False,
+   )
+
+The length of `rdkitmol_list` is one if there are no sidechains and only one
+ligand was docked.
+If multiple ligands and/or sidechains are docked simultaneously, each will be
+an individual RDKit molecule in `rdkitmol_list`.
+In this example, ``keep_flexres=False`` would exclude sidechains from the
+RDKit molecules even if such sidechains existed.
+Note that docking multiple
+ligands simultaneously is only available in Vina, and it differs from docking
+multiple ligands one after the other. Each failed creation of an RDKit molecule
+for a ligand or sidechain results in a `None` in `rdkitmol_list`.
+
+For Vina's output PDBQT files, omit `is_dlg=True`:
+
+.. code-block:: python
+
+   pdbqt_mol = PDBQTMolecule.from_file("vina_results.pdbqt", skip_typing=True)
+
+When using Vina from Python, the output string can be passed directly.
+See [the docs](https://autodock-vina.readthedocs.io/en/latest/docking_python.html)
+for context on the `v` object.
+
+.. code-block:: python
+
+    vina_output_string = v.poses()
+    pdbqt_mol = PDBQTMolecule(vina_output_string, is_dlg=True, skip_typing=True)
+    rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 8f2a4165..4d09009b 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -61,22 +61,24 @@ to run molecular docking and virtual screening.
    Overview <lig_overview>
    Basic usage <lig_prep_basic>
    Advanced usage <lig_prep_advanced>
-   Command line options <lig_cli_options>
+   mk_prepare_ligand.py options <lig_cli_options>
 
 .. toctree::
    :maxdepth: 2
    :hidden:
    :caption: Receptor preparation
 
-   cli_rec_prep
+   Overview <rec_overview>
+   Command line usage <cli_rec_prep>
+   mk_prepare_receptor.py options <rec_cli_options>
 
 .. toctree::
    :maxdepth: 2
    :hidden:
    :caption: Exporting results
 
-   cli_export_result
-   In Python <py_export_result>
+   Usage <export_usage>
+   mk_export.py options <export_cli_options>
 
 .. toctree::
    :maxdepth: 2
@@ -84,6 +86,3 @@ to run molecular docking and virtual screening.
    :caption: Building residue templates
 
    In Python <py_build_temp>
-
-
-
diff --git a/docs/source/py_export_result.rst b/docs/source/py_export_result.rst
deleted file mode 100644
index f6e55da4..00000000
--- a/docs/source/py_export_result.rst
+++ /dev/null
@@ -1,37 +0,0 @@
-Exporting docking results in Python
-===================================
-
-.. code-block:: python
-
-    from meeko import PDBQTMolecule
-    from meeko import RDKitMolCreate
-    
-    fn = "autodock-gpu_results.dlg"
-    pdbqt_mol = PDBQTMolecule.from_file(fn, is_dlg=True, skip_typing=True)
-    rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)
-
-The length of ``rdkitmol_list`` is one if there are no sidechains and only one
-ligand was docked.
-If multiple ligands and/or sidechains are docked simultaneously, each will be
-an individual RDKit molecule in ``rdkitmol_list``.
-Sidechains are truncated at the C-alpha.
-Note that docking multiple
-ligands simultaneously is only available in Vina, and it differs from docking
-multiple ligands one after the other. Each failed creation of an RDKit molecule
-for a ligand or sidechain results in a ``None`` in ``rdkitmol_list``.
-
-For Vina's output PDBQT files, omit ``is_dlg=True``.
-
-.. code-block:: python
-
-    pdbqt_mol = PDBQTMolecule.from_file("vina_results.pdbqt", skip_typing=True)
-    rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)
-
-When using Vina from Python, the output string can be passed directly.
-See `the docs <https://autodock-vina.readthedocs.io/en/latest/docking_python.html>`_ for context on the ``Vina`` object.
-
-.. code-block:: python
-
-    vina_output_string = v.poses()
-    pdbqt_mol = PDBQTMolecule(vina_output_string, is_dlg=True, skip_typing=True)
-    rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)
diff --git a/docs/source/py_lig_prep.rst b/docs/source/py_lig_prep.rst
deleted file mode 100644
index 16c3c367..00000000
--- a/docs/source/py_lig_prep.rst
+++ /dev/null
@@ -1,32 +0,0 @@
-Ligand preparation in Python
-============================
-
-
-Creating a PDBQT string from an RDKit molecule
-----------------------------------------------
-.. code-block:: python
-
-    from meeko import MoleculePreparation
-    from meeko import PDBQTWriterLegacy
-    from rdkit import Chem
-    
-    input_molecule_file = "example/BACE_macrocycle/BACE_4.sdf"
-    
-    # there is one molecule in this SD file, this loop iterates just once
-    for mol in Chem.SDMolSupplier(input_molecule_file, removeHs=False):
-        preparator = MoleculePreparation()
-        mol_setups = preparator.prepare(mol)
-        for setup in mol_setups:
-            setup.show() # optional
-            pdbqt_string = PDBQTWriterLegacy.write_string(setup)
-
-At this point, ``pdbqt_string`` can be written to a file for
-docking with AutoDock-GPU or Vina, or passed directly to Vina within Python
-using ``set_ligand_from_string (pdbqt_string)``. For context, see `the docs on using Vina from Python <https://autodock-vina.readthedocs.io/en/latest/docking_python.html>`_.
-
-One advantage of this approach is that input PDBQT files are not written to the filesystem.
-The PDBQT format is lossy, because it lacks bond orders and non-polar hydrogens,
-making it a poor choice to store molecular information.
-
-Another advantage is to write custom workflows entirely from Python without external
-calls to ``mk_prepare_ligand.py``.
diff --git a/docs/source/rec_cli_options.rst b/docs/source/rec_cli_options.rst
new file mode 100644
index 00000000..20111cff
--- /dev/null
+++ b/docs/source/rec_cli_options.rst
@@ -0,0 +1,130 @@
+Options of mk_prepare_receptor.py
+---------------------------------
+
+Input/Output Options
+~~~~~~~~~~~~~~~~~~~~
+
+.. option:: --read_pdb <PDB_FILENAME>
+
+   Read a PDB file using the PDB parser in RDKit.
+
+.. option:: -i, --read_with_prody <MACROMOL_FILENAME>
+
+   Read a PDB or mmCIF file using ProDy (if installed). ProDy can be installed from PyPI or conda-forge.
+
+.. option:: -o, --output_basename <basename>
+
+   Specify a default basename for output files created by `--write` options when no filename is specified.
+
+.. option:: -p, --write_pdbqt <PDBQT_FILENAME> [*]
+
+   Output PDBQT files with `_rigid` or `_flex` suffixes for flexible residues. Defaults to `--output_basename` if no filename is provided.
+
+.. option:: -j, --write_json <JSON_FILENAME> [*]
+
+   Save the receptor's parameterized configuration to JSON format. Defaults to `--output_basename` if unspecified.
+
+.. option:: -g, --write_gpf <GPF_FILENAME> [*]
+
+   Output an AutoGrid input file (GPF). Defaults to `--output_basename` if not specified.
+
+.. option:: -v, --write_vina_box <VINA_BOX_FILENAME> [*]
+
+   Generate a configuration file for Vina with grid box dimensions. Defaults to `--output_basename` if not specified.
+
+.. option:: --write_pdb <PDB_FILENAME> [*]
+
+   Save the prepared receptor in PDB format. Must specify the filename.
+
+Receptor Perception Options
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. option:: -n, --set_template <template>
+
+   Assign templates to residues, e.g., `A:5,7=CYX,B:17=HID`.
+
+.. option:: -d, --delete_residues <residues>
+
+   Specify residues to delete, e.g., `A:350,B:15,16,17`.
+
+.. option:: -b, --blunt_ends <positions>
+
+   Blunt end definitions, e.g., `A:123,200=2,A:1=0`.
+
+.. option:: --add_templates <JSON_FILENAME> [*]
+
+   Load additional templates from one or more JSON files.
+
+.. option:: -a, --allow_bad_res
+
+   (Flag) Ignore residues with missing atoms instead of raising an error.
+
+.. option:: --default_altloc <location>
+
+   Define a default alternate location for residues, overridden by `--wanted_altloc`.
+
+.. option:: --wanted_altloc <location>
+
+   Specify alternate locations for particular residues, e.g., `:5=B,B:17=A`.
+
+.. option:: --mk_config <JSON_FILENAME>
+
+   Specify a JSON configuration file for receptor preparation.
+
+Grid Box Options
+~~~~~~~~~~~~~~~~
+
+.. option:: --box_size <X Y Z>
+
+   Set the size of the grid box in Angstroms (x, y, z).
+
+.. option:: --box_center <X Y Z>
+
+   Define the center of the grid box in Angstroms (x, y, z).
+
+.. option:: --box_center_off_reactive_res
+
+   (Flag) Shift the grid box center 5 Å along the CA-CB bond from CB. Applicable only when there is one reactive flexible residue.
+
+.. option:: --box_enveloping <FILENAME>
+
+   Adjust the grid box to enclose atoms in the specified file. Supported formats: `.sdf`, `.mol`, `.mol2`, `.pdb`, `.pdbqt`.
+
+.. option:: --padding <value>
+
+   Set padding around atoms specified in `--box_enveloping` (in Å).
+
+Flexible and/or Reactive Options
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. option:: -f, --flexres <residues>
+
+   Define flexible residues by chain ID and residue number, e.g., `-f ":42,B:23"`. 
+
+.. option:: -r, --reactive_flexres <residues>
+
+   Define reactive flexible residues by chain ID and residue number, e.g., `-r ":42,B:23"`. Maximum of 8 reactive residues.
+
+.. option:: --reactive_name <residue:atom>
+
+   Specify the reactive atom name for a residue type, e.g., `--reactive_name "TRP:NE1"`. Can be repeated for multiple assignments.
+
+.. option:: -s, --reactive_name_specific <residue:atom>
+
+   Specify the reactive atom for an individual residue by residue ID, e.g., `-s "A:42=NE2"`. The residue becomes reactive.
+
+.. option:: --r_eq_12 <value>
+
+   Set the equilibrium distance (r_eq) for reactive atoms in 1-2 interactions. Default is 1.8 Å.
+
+.. option:: --eps_12 <value>
+
+   Set the epsilon value for reactive atoms in 1-2 interactions. Default is 2.5.
+
+.. option:: --r_eq_13_scaling <factor>
+
+   Scale r_eq for 1-3 interactions across reactive atoms. Default is 0.5.
+
+.. option:: --r_eq_14_scaling <factor>
+
+   Scale r_eq for 1-4 interactions across reactive atoms. Default is 0.5.
diff --git a/docs/source/rec_overview.rst b/docs/source/rec_overview.rst
new file mode 100644
index 00000000..da9ffde2
--- /dev/null
+++ b/docs/source/rec_overview.rst
@@ -0,0 +1,73 @@
+Overview of receptor preparation
+================================
+
+Polymers and monomers
+---------------------
+
+Receptors are represented as a collection of subunits where each subunit is
+treated as an individual molecule. THe same code path for ligand parameterization
+is also used herein. The Python class for receptors is called
+``Polymer`` and the class for subunits is ``Monomer``.
+
+During most of the v0.6.0 development, the name of the ``Polymer`` class was
+``LinkedRDKitChorizo``. Many of the commit messages, as well as GitHub issues and pull
+requests mention chorizos rather then polymers. The analogy is that each of the
+links in the chorizo is an RDKit molecule, and that these classes often represent protein.
+
+The input structure is matched against templates to
+guarantee chemical correctness and identify problems with the input structures.
+This allows the user to identify and fix problems, resulting in a molecular
+model that is correct with respect to heavy atoms, protonation state,
+connectivity, bond orders, and formal charges.
+
+The matching algorithm uses the connectivity and elements, but not bond orders
+or atom names. Hydrogens are optional. This makes it compatible with input
+files from various sources.
+
+Templates are matched on a per residue basis. Each residue is represented
+as an instance of a PolymerResidue object, which contains:
+
+ * an RDKit molecule that represents the actual state
+ * a padded RDKit molecule containing a few atoms from the adjacent residues
+ * parameters such as partial charges
+
+The positions are set by the input, and the connectivity and formal charges
+are defined by the templates. Heavy atoms must match exactly. If heavy atoms
+are missing or in excess, the templates will fail to match.
+
+Missing hydrogens are added by RDKit, but are not subjected to minimization
+with a force field. Thus, their bond lengths are not super accurate.
+
+Different states of the same residue are stored as different templates,
+for example different protonation states of HIS, N-term, LYN/LYS, etc.
+Residue name is primary key unless user overrides.
+
+Currently not supported: capped residues from charmm-gui.
+
+
+Templates
+---------
+
+The templates contain SMILES strings that are used to create the RDKit
+molecules that constitute every residue in the processed model. In this way,
+the chemistry of the processed model is fully defined by the templates,
+and the only thing that is preserved from the input are the atom positions
+and the connectivity between residues.
+
+The SMILES strings contain all atoms that exist in the final model,
+and none that do not exist. This also applies to hydrogens,
+meaning that the SMILES are expected to have real hydrogens. Note that
+real hydrogens are different from explicit hydrogens. Real hydrogens will be
+represented as an actual atom in an RDKit molecule, while explicit hydrogens
+are a just property of heavy atoms. In the SMILES, real hydrogens are defined
+with square brackets "[H]" and explicit hydrogens without, e.g. "[nH]" to set
+the number of explicit hydrogens on an aromatic nitrogen to one.
+
+Residues that are part of a polymer, which is often all of them, will have
+bonds to adjacent residues. The heavy atoms involved in the bonds will miss
+a real hydrogen and have an implicit (or explicit) one instead. As an
+example, consider modeling an alkyl chain as a polymer, in which the monomer
+is a single carbon atom. Our template SMILES would be "[H]C[H]". The RDKit
+molecule will have three atoms and the carbon will have two implicit hydrogens.
+The implicit hydrogens correspond to bonds to adjacent residues in the
+processed polymer. 

From 61856b0e95c867cc7181cc645d1d47497fc136c4 Mon Sep 17 00:00:00 2001
From: diogom <diogom@scripps.edu>
Date: Mon, 11 Nov 2024 01:53:57 -0800
Subject: [PATCH 16/18] add clarifycation from Amy's 46a2f3a

---
 docs/source/cli_rec_prep.rst | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/docs/source/cli_rec_prep.rst b/docs/source/cli_rec_prep.rst
index 823081fb..84944487 100644
--- a/docs/source/cli_rec_prep.rst
+++ b/docs/source/cli_rec_prep.rst
@@ -75,3 +75,6 @@ The arguments involving assignment of residues to properties:
 use the residue selection lanaguge described above, followed by an equal sign (``=``) as the delimiter and the assigned value, which could be the name of a residue template, the atom index for the blunt end, the wanted altloc ID, or the atom name of the reactive atom. Each residue selection is comibned with the most recent assignment that precedes it, resulting in a further expanded list of residue-assignment pairs. 
 
 For an input like ``"A:5,7=CYX,A:19A,B:17=HID``, this assignment language represents: ``residues (number) 5 in Chain A are set to (template name) CYX`` and ``residue (number) 19 A in Chain A, and residue (number) 17 in Chain B are set to (template name) HID``. 
+
+At present, Meeko always requires **chain ID** and **residue number** to identify a residue (cofactor, ligand, or ion). The only exception occurs when the chain ID in the input file is empty; in this case, the chain ID should be omitted from the selection language, i.e. ``"A:19,:17"`` represents: ``residue (number) 19 in Chain A`` and ``residue (number) 17 with an empty chain ID``. 
+

From 48690089eb1109eb3e5e18de95add856e004cb21 Mon Sep 17 00:00:00 2001
From: diogom <diogom@scripps.edu>
Date: Mon, 11 Nov 2024 02:19:14 -0800
Subject: [PATCH 17/18] update README

---
 README.md | 285 +++++++-----------------------------------------------
 1 file changed, 33 insertions(+), 252 deletions(-)

diff --git a/README.md b/README.md
index 871b5b0b..1d41fc6a 100644
--- a/README.md
+++ b/README.md
@@ -1,286 +1,67 @@
-# Meeko: preparation of small molecules for AutoDock
+# Meeko: interface for AutoDock
 
 [![API stability](https://img.shields.io/badge/stable%20API-no-orange)](https://shields.io/)
 [![PyPI version fury.io](https://img.shields.io/badge/version-0.6.0-green.svg)](https://pypi.python.org/pypi/meeko/)
 [![Documentation Status](https://readthedocs.org/projects/meeko/badge/?version=readthedocs)](https://meeko.readthedocs.io/en/readthedocs/?badge=readthedocs)
 
-Meeko reads an RDKit molecule object and writes a PDBQT string (or file)
-for [AutoDock-Vina](https://github.com/ccsb-scripps/AutoDock-Vina)
-and [AutoDock-GPU](https://github.com/ccsb-scripps/AutoDock-GPU).
-It converts the docking output to RDKit molecules and
-SD files, without loss of bond orders.
+Meeko prepares the input for AutoDock and processes its output.
+It is developed alongside AutoDock-GPU and AutoDock-Vina.
+Meeko parameterizes both small organic molecules (ligands) and proteins
+and nucleic acids (receptors).
 
 Meeko is developed by the [Forli lab](https://forlilab.org/) at the
 [Center for Computational Structural Biology (CCSB)](https://ccsb.scripps.edu)
 at [Scripps Research](https://www.scripps.edu/).
 
 
-## Usage notes
+## Documentation
 
-Meeko does not calculate 3D coordinates or assign protonation states.
-Input molecules must have explicit hydrogens.
+The docs are hosted on [meeko.readthedocs.io](meeko.readthedocs.io)
 
-Sampling of macrocycle conformers ([paper](https://doi.org/10.1017/qrd.2022.18))
-is enabled by default.
 
-SDF format strongly preferred over Mol2.
-See
-[this discussion](https://github.com/rdkit/rdkit/discussions/3647), and
-[this one](https://sourceforge.net/p/rdkit/mailman/message/37668451/),
-[also this](rdkit/rdkit#4061),
-[and this](https://sourceforge.net/p/rdkit/mailman/message/37374678/).
-and RDKit issues
-[1755](https://github.com/rdkit/rdkit/issues/1755) and
-[917](https://github.com/rdkit/rdkit/issues/917). So, what could go wrong?
-For example, reading Mol2 files from ZINC
-[led to incorrect net charge of some molecules.](https://github.com/forlilab/Meeko/issues/63)
+## Reporting bugs
 
+Please check if a similar bug has been reported and, if not, [open an issue](https://github.com/forlilab/Meeko/issues).
 
-## v0.6.0-alpha.3
 
-This release aims to distribute an enhanced `mk_prepare_receptor.py`.
-Some features are still being developed, hence the `-alpha` suffix in the version.
-Reactive docking is not working in v0.6.0-alpha, but should be restored soon.
-Documentation is also work in progress.
+## Installation
 
+Visit the docs for a more complete description. One option is conda or mamba:
 
-## API changes in v0.5
-
-Class `MoleculePreparation` no longer has method `write_pdbqt_string()`.
-Instead, `MoleculePreparation.prepare()` returns a list of `MoleculeSetup` objects
-that must be passed, individually, to `PDBQTWriterLegacy.write_string()`.
-```python
-from meeko import MoleculePreparation
-from meeko import PDBQTWriterLegacy
-
-preparator = MoleculePreparation()
-mol_setups = preparator.prepare(rdkit_molecule_3D_with_Hs)
-for setup in mol_setups:
-    pdbqt_string, is_ok, error_msg = PDBQTWriterLegacy.write_string(setup)
-    if is_ok:
-        print(pdbqt_string, end="")
-```
-
-Argument `keep_nonpolar_hydrogens` is replaced by `merge_these_atom_types`, both in the Python
-interface and for script `mk_prepare_ligand.py`.
-The default is `merge_these_atom_types=("H",)`, which
-merges hydrogens typed `"H"`, keeping the current default behavior.
-To keep all hydrogens, set `merge_these_atom_types` to an empty
-list when initializing `MoleculePreparation`, or pass no atom types
-to `--merge_these_atom_types` from the command line:
-```sh
-mk_prepare_ligand.py -i molecule.sdf --merge_these_atom_types
-``` 
-
-## Dependencies
-
-* Python (>=3.5)
-* Numpy
-* Scipy
-* RDKit (>=2023.09)
-* ProDy (optionally, for covalent docking)
-
-Conda or Miniconda can install the dependencies:
 ```bash
-conda install -c conda-forge numpy scipy rdkit
-pip install prody # optional. pip recommended at http://prody.csb.pitt.edu/downloads/
-```
-
-
-## Python tutorial
-
-
-#### 2. RDKit molecule from docking results
-
-```python
-from meeko import PDBQTMolecule
-from meeko import RDKitMolCreate
-
-fn = "autodock-gpu_results.dlg"
-pdbqt_mol = PDBQTMolecule.from_file(fn, is_dlg=True, skip_typing=True)
-rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)
-```
-The length of `rdkitmol_list` is one if there are no sidechains and only one
-ligand was docked.
-If multiple ligands and/or sidechains are docked simultaneously, each will be
-an individual RDKit molecule in `rdkitmol_list`.
-Sidechains are truncated at the C-alpha.
-Note that docking multiple
-ligands simultaneously is only available in Vina, and it differs from docking
-multiple ligands one after the other. Each failed creation of an RDKit molecule
-for a ligand or sidechain results in a `None` in `rdkitmol_list`.
-
-For Vina's output PDBQT files, omit `is_dlg=True`.
-```python
-pdbqt_mol = PDBQTMolecule.from_file("vina_results.pdbqt", skip_typing=True)
-rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)
-```
-
-When using Vina from Python, the output string can be passed directly.
-See [the docs](https://autodock-vina.readthedocs.io/en/latest/docking_python.html)
-for context on the `v` object.
-```python
-vina_output_string = v.poses()
-pdbqt_mol = PDBQTMolecule(vina_output_string, is_dlg=True, skip_typing=True)
-rdkitmol_list = RDKitMolCreate.from_pdbqt_mol(pdbqt_mol)
+micromamba install meeko
 ```
 
-#### 3. Initializing MoleculePreparation from a dictionary (or JSON)
-
-This is useful for saving and loading configuration files with json.
-```python
-import json
-from meeko import MoleculePreparation
-
-mk_config = {"hydrate": True} # any arguments of MoleculePreparation.__init__
-print(json.dumps(mk_config), file=open('mk_config.json', 'w'))
-with open('mk_config.json') as f:
-    mk_config = json.load(f)
-preparator = MoleculePreparation.from_config(mk_config)
-```
-The command line tool `mk_prepare_ligand.py` can read the json files using
-option `-c` or `--config`.
-
-
-## Possibly useful configurations of MoleculePreparation
-
-Here we create an instance of MoleculePreparation that attaches pseudo
-waters to the ligand ([see paper on hydrated docking](https://pubs.acs.org/doi/abs/10.1021/jm2005145)),
-keeps macrocycles rigid (generally a bad idea),
-and keeps conjugated bonds and amide bonds rigid. 
-By default, most amides are kept rigid, except for tertiary amides with
-different substituents on the nitrogen.
+or from PyPI:
 
-```python
-preparator = MoleculePreparation(
-    hydrate=True,
-    rigid_macrocycles=True,
-    rigidify_bonds_smarts = ["C=CC=C", "[CX3](=O)[NX3]"],
-    rigidify_bonds_indices = [(1, 2), (0, 2)],
-)
-```
-
-The same can be done with the command line script. Note that indices of the
-atoms in the SMARTS are 0-based for the Python API but
-1-based for the command line script:
-```console
-mk_prepare_ligand.py\
-    --hydrate\
-    --rigid_macrocycles\
-    --rigidify_bonds_smarts "C=CC=C"\
-    --rigidify_bonds_indices 2 3\
-    --rigidify_bonds_smarts "[CX3](=O)[NX3]"\
-    --rigidify_bonds_indices 1 3
+```bash
+pip install meeko
 ```
 
-## Docking covalent ligands as flexible sidechains
+## Usage
 
-The input ligand must be the product of the reaction and contain the
-atoms of the flexible sidechain up to (and including) the C-alpha.
-For example, if the ligand is an acrylamide (smiles: `C=CC(=O)N`) reacting
-with a cysteine (sidechain smiles: `CCS`), then the input ligand for
-Meeko must correspond to smiles `CCSCCC(=O)N`.
+Meeko exposes a Python API to enable scripting. Here we share very minimal examples
+using the command line scripts just to give context.
+Please visit the [meeko.readthedocs.io](meeko.readthedocs.io) for more information.
 
-Meeko will align the ligand atoms that match the C-alpha and C-beta of
-the protein sidechain. Options `--tether_smarts` and `--tether_smarts_indices`
-define these atoms. For a cysteine, `--tether_smarts "SCC"` and
-`--tether_smarts_indices 3 2` would work, although it is safer to define
-a more spefic SMARTS to avoid matching the ligand more than once. The first
-index (3 in this example) defines the C-alpha, and the second index defines
-the C-beta. 
-
-For the example in this repository, which is based on PDB entry 7aeh,
-the following options prepare the covalently bound ligand for tethered docking:
-```console
-cd example/covalent_docking
-
-mk_prepare_ligand.py\
-    -i ligand_including_cys_sidechain.sdf\
-    --receptor protein.pdb\
-    --rec_residue ":CYS:6"\
-    --tether_smarts "NC(=O)C(O)(C)SCC"\
-    --tether_smarts_indices 9 8\
-    -o prepared.pdbqt
+Parameterizing a ligand and writing a PDBQT file:
+```bash
+mk_prepare_ligand.py -i molecule.sdf -o molecule.pdbqt
 ```
 
-Prody is required for preparing ligands as flexible sidechains.
-Often, `pip install prody` suffices to install Prody in the currently
-active virtual environment (e.g., conda). For more detailed installation
-instructions visit http://prody.csb.pitt.edu/downloads/.
+Parameterizing a receptor with a flexible sidechain and writing a PDBQT file
+as well as a JSON file that stores the entire receptor datastructure. In this
+example, the `-o` option sets the output base name, `-j` triggers writing the
+.json file, `-p` triggers writting the .pdbqt file, and `-f` makes residue
+42 in chain A flexible.
 
-## Reactive Docking
-
-### 1. Prepare protein with waterkit
-Follow `wk_prepare_receptor.py` instructions and run with `--pdb`.
-The goal of this step is to perform essential fixes to the protein
-(such as missing atoms), to add hydrogens, and to follow the Amber
-naming scheme for atoms and residues, e.g., `HIE` or `HID`
-instead of `HIS`.
-
-### 2. Prepare protein pdbqt
-Here, `wk.pdb` was written by waterkit. The example below will center a gridbox of specified size on the given reactive residue.
-
-```console
-   $ mk_prepare_receptor.py\
-    --pdb wk.pdb\
-    -o receptor.pdbqt\
-    --flexres " :ARG:348"\
-    --reactive_flexres " :SER:308"
-    --reactive_flexres " :SER:308"\
-    --box_center_on_reactive_res\
-    --box_size 40 40 40  # x y z (angstroms)
-```
-A manual box center can be specified with `--box_center`.
-Reactive parameters can also be modified:
-```sh
-  --r_eq_12 R_EQ_12     r_eq for reactive atoms (1-2 interaction)
-  --eps_12 EPS_12       epsilon for reactive atoms (1-2 interaction)
-  --r_eq_13_scaling R_EQ_13_SCALING
-                        r_eq scaling for 1-3 interaction across reactive atoms
-  --r_eq_14_scaling R_EQ_14_SCALING
-                        r_eq scaling for 1-4 interaction across reactive atoms
+```bash
+mk_prepare_receptor.py -i nucleic_acid.cif -o my_receptor -j -p -f A:42
 ```
 
-Receptor preparation can't handle heteroatoms for the moment.
-Also nucleic acids, ions, and post-translational modifications (e.g.
-phosphorilation) are not supported. Only the 20 standard amino acids
-can be parsed, and it is required to have Amber atom names and
-hydrogens. No atoms can be missing.
-
-### 3. Run autogrid
+Finally, converting docking results to SDF for the ligand, and PDB for the
+receptor with updated sidechain positions:
 
-Make affinity maps for the `_rigid.pdbqt` part of the receptor.
-Make affinity maps for the `_rigid.pdbqt` part of the receptor. `mk_prepare_receptor.py` will prepare the GPF for you.
-
-### 4. Write ligand PDBQT
-mk_prepare_ligand.py -i sufex1.sdf --reactive_smarts "S(=O)(=O)F" --reactive_smarts_idx 1 -o sufex1.pdbqt\
-
-### 5. Configure AD-GPU for reactive docking
-
-For reactive docking there are two options that need to be passed to AutoDock-GPU:
-For reactive docking there are an additional option that needs to be passed to AutoDock-GPU:
-    ```console
-    --import_dpf
-    ```
-
-The `--derivtype` option, if needed, was written by `mk_prepare_receptor.py` to a file suffixed with `.derivtype`.
-
-The filename to be passed to `--import_dpf` was written by `mk_prepare_receptor.py`
-and it is suffixed with `reactive_config`.
-```sh
-ADGPU -I *.reactive_config -L sufex1.pdbqt -N sufex1_docked_ -F *_flex.pdbqt -C 1
+```bash
+mk_export.py vina_results.pdbqt -j my_receptor.json -s lig_docked.sdf -p rec_docked.pdb
 ```
-
-## Polymer docs temporary placeholder
-
-Current init is a builder method designed to parse PDB files and topology
-objects. It is not designed to monomerize a polymer that already exists
-as a single RDKit molecule, or to build a single polymer molecule from the
-individual residues (monomers). To do that, one would have to react residues
-with each other. Padding could be adapted.
-
-
-How bonds between residues are handled and blunt residues, and deleting
-bonds "interactively".
-
-Document residue\_chem\_params.

From 635f33807d3d85dca0093399088ee52e4794de0d Mon Sep 17 00:00:00 2001
From: diogom <diogom@scripps.edu>
Date: Mon, 11 Nov 2024 02:25:45 -0800
Subject: [PATCH 18/18] mv templates docs inside rec prep

---
 docs/source/index.rst         | 8 +-------
 docs/source/py_build_temp.rst | 6 +++---
 2 files changed, 4 insertions(+), 10 deletions(-)

diff --git a/docs/source/index.rst b/docs/source/index.rst
index 4d09009b..532cca7d 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -69,6 +69,7 @@ to run molecular docking and virtual screening.
    :caption: Receptor preparation
 
    Overview <rec_overview>
+   Info on templates  <py_build_temp>
    Command line usage <cli_rec_prep>
    mk_prepare_receptor.py options <rec_cli_options>
 
@@ -79,10 +80,3 @@ to run molecular docking and virtual screening.
 
    Usage <export_usage>
    mk_export.py options <export_cli_options>
-
-.. toctree::
-   :maxdepth: 2
-   :hidden:
-   :caption: Building residue templates
-
-   In Python <py_build_temp>
diff --git a/docs/source/py_build_temp.rst b/docs/source/py_build_temp.rst
index 079151ee..91de0428 100644
--- a/docs/source/py_build_temp.rst
+++ b/docs/source/py_build_temp.rst
@@ -1,7 +1,7 @@
 .. _py_build_temp:
 
-Building residue templates in Python
-===================================
+Advanced information about templates
+====================================
 
 The interpretation of the valence (bonds) and formal charge of atoms is an essential step when parsing a PDB/CIF file, and the accuracy of residue mapping is crucial to the creation of a macrobiomolecule system. In Meeko, the input residue names are used as keys and the chemical templates are retrieved accordingly based on :ref:`templates <templates>`. 
 
@@ -278,4 +278,4 @@ And below is the content of ``CRO_templates.json``, which can be loaded by ``mk_
                 "link_labels": {"1": "N-term"}
             }
         }
-    }
\ No newline at end of file
+    }