Parameterizing Ligands and Building Systems with AmberStudio
====================================================================================

| Zeke A. Piskulich\ :sup:`1`,  Patricio Barletta\ :sup:`1`, Ryan Snyder\ :sup:`1`, and Darrin M. York\ :sup:`1`
| :sup:`1`\ Laboratory for Biomolecular Simulation Research, Institute
|   for Quantitative Biomedicine and Department of Chemistry and Chemical
|   Biology, Rutgers University, Piscataway, NJ 08854, USA

Learning objectives
-------------------

.. start-learning-objectives

- Learn what files are required for an AmberStudio ABFE calculation.
- Set up and prepare files for ABFE on Tyk2 inhibitors.

.. end-learning-objectives

Activity
--------

In this activity, you will use AmberStudio to generate ligand parameters and build systems.

.. start-tutorial

Generating Ligand Parameters with AmberStudio
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

To generate ligand parameters with AmberStudio, you need to start with a :term:`PDB <.pdb>`. AmberStudio contains Flows for ligand parameterization that use the Ligand Parameterization tools described in the Ligand Parameterization tutorials located here: :doc:`/ModularTutorials/ForceField/ligand_parameterization/ligandparam`

:download:`Download the example files <https://drive.google.com/file/d/1aZXCEzUK-gYPwfIfCuiyV_xRikRxHf_2/view?usp=drive_link>`

To do this with AmberStudio, you would start with the following directory structure:

.. code-block:: bash

    <top-dir>
        run_bcc_param.py
        binder_ejm31.pdb

The contents of run_bcc_param.py are as follows:

.. code-block:: python

    from pathlib import Path

    from amberstudio.artifacts import ArtifactContainer, BinderLigandPDB
    from amberstudio.worknodes import ParametrizeBinderBCC


    def main() -> None:
        cwd = Path(__file__).resolve().parent
        ligand_pdb = cwd / "binder_ejm31.pdb"

        node = ParametrizeBinderBCC(
            "param_bcc",
            root_dir=cwd,
            resname="LIG",
            atom_type="gaff2",
            charge_model="bcc",
            charge=0,
        )

        input_artifacts = ArtifactContainer("ejm31", [BinderLigandPDB(ligand_pdb)])
        node.run(input_artifacts, cwd=cwd, sysname="ejm31")


    if __name__ == "__main__":
        main()

If you look at the output WorkNode folder (``param_bcc``), you will see all the files from the parameterization. The files you would want are:

- :term:`binder_ejm31.mol2 <.mol2>`
- :term:`binder_ejm31.lib <.lib>`
- :term:`binder_ejm31.frcmod <.frcmod>`
- :term:`binder_ejm31.pdb <.pdb>`

These would be used to build a system.

.. note::

    In this section, we demonstrated how a single WorkNode could be run independently from a Pipeline;
    however, we would rarely want to do this. If we had added this WorkNode to a Pipeline, we wouldn't
    have to specify the ``ArtifactContainer``, and the WorkNode would run over all systems.


Building Systems with AmberStudio
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

:download:`Download the example files <https://drive.google.com/file/d/1rIRvhjFnFP_Qz0TY1MJK_V3SG8ttKkRz/view?usp=drive_link>`

Now, we will demonstrate building a system with AmberStudio.

.. code-block:: bash

    <top-dir>
        run_build.py
        <ejm31>/
            binder_ejm31.mol2
            binder_ejm31.lib
            binder_ejm31.frcmod
            binder_ejm31.pdb
            target_ejm31.pdb

The input files here are:

- The ligand parameter files (:term:`mol2 <.mol2>`, :term:`lib <.lib>`, :term:`frcmod <.frcmod>`) generated from a ligand parameterization (e.g. using the above script)
- binder_ejm31.pdb is a :term:`PDB <.pdb>` of the ligand.
- target_ejm31.pdb is a :term:`PDB <.pdb>` of the protein target.

.. note::

    AmberStudio assumes that the coordinates of **binder_ejm31.pdb** correspond to a docked structure. Presently, AmberStudio does not perform docking.

.. code-block:: python

    from pathlib import Path
    import logging

    from amberstudio import Pipeline
    from amberstudio.flows import FlowBuild
    from amberstudio.schedulers import ReferenceScheduler


    def main() -> None:
        cwd = Path(__file__).resolve().parent

        pipeline = Pipeline(
            "ejm31_build",
            cwd=cwd,
            scheduler=ReferenceScheduler(max_gpus=0),
            logging_level=logging.INFO,
            ignore_checkpoint=True,
        )

        pipeline.append_flow(
            FlowBuild(
                parametrize=False,
                skippable=False,
            )
        )
        pipeline.launch()


    if __name__ == "__main__":
        main()

Note a few things about this script:

- We set up a pipeline this time, which contains the build flow.
- Since we already parameterized the ligands, we set parametrize to False.
- This script brings you from the above files, all the way to a built :term:`parm7 <.parm7>` and :term:`rst7 <.rst7>` pair for both the aqueous (ligand in water) and complex configurations (ligand in protein) of the ligand.

Now let's explore the **ejm31** folder. You will see that there are many different folders that have been added.
Each folder was the working directory of a separate WorkNode in the calculation.

The **leaf_build** folder contains the output systems.

.. note::

    These **leaf** WorkNodes are called "dummy" WorkNodes because they don't usually do anything
    but collect the output from a flow into a single directory. Pipelines do the same when they
    start, they collect all input files into a **Root** WorkNode.

Try opening these with VMD.

.. code-block:: bash

    vmd python complex_ejm31.parm7 complex_ejm31.rst7

You should see the ligand bound into the protein. In the current example, HMR was applied to the system to
use hydrogen-mass repartitioning. No Molecular Dynamics has been run yet.


.. end-tutorial