Ligand Parameterization Basics
==============================

| Zeke A. Piskulich\ :sup:`1`, Patricio Barletta\ :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

- Perform the basic steps of ligand parameterization by hand using Antechamber, ParmChk2, and tleap.
- Look through the files generated by these programs to understand what information is stored where.
- Perform ligand parameterization using RESP charge fitting with Antechamber and Gaussian.

.. end-learning-objectives

Relevant literature
-------------------

Outline
-------

In this tutorial we will be covering the basic use of ligand parameterization with LigandParam, a tool developed by the York group to make ligand parameterization easy to do in a way that matches our best practices.

This tutorial will cover both the old way of doing simple parameterizations by hand, and then show how to do the same thing with LigandParam. 

Tutorial
--------

.. start-tutorial

.. mermaid::
    
  flowchart LR
      Start[Start: ligand.pdb]
      
      subgraph ByHand [Simple AM1-BCC]
          direction LR
          A1[antechamber -c bcc<br/>-> ligand.mol2]
          A2[parmchk2 -f mol2<br/>-> ligand.frcmod]
          A3[tleap -f leap.in<br/>-> ligand.lib]
          A4[python mol2_to_pdb_bfactor.py<br/>-> ligand_charges.pdb]
          A5[VMD visualize_charges.tcl]
      end
  
      subgraph RESP [RESP charge fitting]
      direction LR
      R1[antechamber -fo gcrt -gv<br/>-> ligand.gjf]
      R2[Gaussian: g16 < ligand.gjf > ligand_run.log]
      R3[antechamber -fi gout -c resp<br/>-> ligand_resp.mol2]
      R4[parmchk2 -i ligand_resp.mol2<br/>-> ligand_resp.frcmod]
      R5[tleap -f leap_resp.in <br/>-> ligand_resp.lib]
      R6[python mol2_to_pdb_bfactor.py<br/>-> resp_charges.pdb]
      end
  
      Start --> ByHand
      Start --> RESP
  
      ByHand --> A1 --> A2 --> A3 --> A4 --> A5
      RESP --> R1 --> R2 --> R3 --> R4 --> R5 --> R6
  
      %% Comparison arrow
      A4 -. Compare .-> R6
      CompareNote[(Compare charges / diff pdb<br/>python pdb_bfactor_diff.py<br/>vmd -e visualize_charges.tcl -args diff.pdb)]
      A4 --> CompareNote
      R6 --> CompareNote
  
      %% files summary node
      Files[[Files produced:
      - ligand.mol2, ligand.frcmod, ligand.lib
      - ligand.gjf, ligand_run.log, ligand_resp.mol2, ligand_resp.frcmod, ligand_resp.lib
      - ligand_charges.pdb, resp_charges.pdb, diff.pdb
      ]]
      CompareNote --> Files
  
Parameterizing with AM1-BCC Charges
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

At the very simplest, a ligand parameterization can be done just with Antechamber, ParmChk2, tleap, as follows.

Start by creating a subdirectory called by_hand, and copy the ligand pdb into that subdirectory from your current working directory.

.. note::

   This tutorial assumes you have a file called ligand.pdb in your current working directory, this can be any small molecule you want to parameterize. In this tutorial, we will be working with a minimal ligand that has a net -1 charge. 

   .. image:: /_static/files/ModularTutorials/ForceField/ligand_parameterization/minimal.png
    :alt: Minimal Ligand
    :width: 400px

   .. toggle::

        .. code-block:: console

            ATOM      1 N1   LIG A   1       4.400  -0.997   0.000  0.00  0.00           N
            ATOM      2 C2   LIG A   1       3.064  -1.169   0.000  0.00  0.00           C
            ATOM      3 C3   LIG A   1       2.164  -0.103  -0.000  0.00  0.00           C
            ATOM      4 C4   LIG A   1       2.724   1.178  -0.000  0.00  0.00           C
            ATOM      5 N5   LIG A   1       4.054   1.386  -0.000  0.00  0.00           N
            ATOM      6 C6   LIG A   1       4.815   0.279   0.000  0.00  0.00           C
            ATOM      7 C7   LIG A   1       0.679  -0.362  -0.000  0.00  0.00           C
            ATOM      8 O7   LIG A   1       0.294  -1.565  -0.001  0.00  0.00           O
            ATOM      9 N8   LIG A   1      -0.049   0.760  -0.000  0.00  0.00           N
            ATOM     10 C9   LIG A   1      -1.396   0.628  -0.000  0.00  0.00           C
            ATOM     11 S10  LIG A   1      -2.272  -0.950  -0.000  0.00  0.00           S
            ATOM     12 C11  LIG A   1      -3.705   0.008   0.000  0.00  0.00           C
            ATOM     13 N12  LIG A   1      -3.551   1.310   0.000  0.00  0.00           N
            ATOM     14 N13  LIG A   1      -2.221   1.672  -0.000  0.00  0.00           N
            ATOM     15 O14  LIG A   1      -4.950  -0.558   0.001  0.00  0.00           O
            ATOM     16 H2   LIG A   1       2.685  -2.188   0.000  0.00  0.00           H
            ATOM     17 H6   LIG A   1       5.893   0.435   0.000  0.00  0.00           H
            ATOM     18 H4   LIG A   1       2.075   2.050  -0.000  0.00  0.00           H
            ATOM     19 H141 LIG A   1      -5.545   0.220   0.001  0.00  0.00           H
            END
 

.. code-block:: bash

    mkdir by_hand
    cp ligand.pdb by_hand/
    cd by_hand

Assume you are starting with ligand.pdb, first you would call antechamber to create a mol2 file with bcc charges.

.. code-block:: bash

    # Call Antechamber to apply bcc charges to ligand.pdb
    # -i : Input file (ligand.pdb)
    # -o : Output file (ligand.mol2)
    # -fi: Input file format (pdb)
    # -fo: Output file format (mol2)
    # -c : Charge method (bcc)
    # -rn: Residue name (LIG)
    # -at: Atom type (gaff2)
    # -nc: Net charge (the default ligand for this tutorial has a net -1 charge)
    antechamber -i ligand.pdb -o ligand.mol2 -fi pdb -fo mol2 -c bcc -rn LIG -at gaff2 -nc -1

This command creates a mol2 file with BCC-assigned partial charges next to each atom. If this runs correctly you'll see something like this:

.. code-block:: console

    antechamber -i ligand.pdb -o ligand.mol2 -fi pdb -fo mol2 -c bcc -rn LIG -at gaff2 -nc -1
    Info: acdoctor mode is on: check and diagnose problems in the input file.
    Info: The atom type is set to gaff2; the options available to the -at flag are
          gaff, gaff2, amber, bcc, abcg2, and sybyl.

    -- Check Format for pdb File --
       Status: pass
    Info: Bond types are assigned for valence state (1) with penalty (1).

    Running: /home/piskulic/Project/SoftwareDevelopment/install_lbsr_dev/bin/sqm -O -i sqm.in -o sqm.out

At the end, you should now have a mol2 for the ligand. The final mol2 should look similar to what you can find below.

.. toggle:: 

    .. code-block:: console

        @<TRIPOS>MOLECULE
        LIG
           19    20     1     0     0
        SMALL
        bcc


        @<TRIPOS>ATOM
              1 N1           4.4000    -0.9970     0.0000 nb         1 LIG      -0.723000
              2 C2           3.0640    -1.1690     0.0000 ca         1 LIG       0.447700
              3 C3           2.1640    -0.1030    -0.0000 ca         1 LIG      -0.292200
              4 C4           2.7240     1.1780    -0.0000 ca         1 LIG       0.447700
              5 N5           4.0540     1.3860    -0.0000 nb         1 LIG      -0.723000
              6 C6           4.8150     0.2790     0.0000 ca         1 LIG       0.605900
              7 C7           0.6790    -0.3620    -0.0000 ce         1 LIG       0.538900
              8 O7           0.2940    -1.5650    -0.0010 o          1 LIG      -0.741300
              9 N8          -0.0490     0.7600    -0.0000 nf         1 LIG      -0.562000
             10 C9          -1.3960     0.6280    -0.0000 cd         1 LIG       0.700400
             11 S10         -2.2720    -0.9500    -0.0000 ss         1 LIG      -0.243800
             12 C11         -3.7050     0.0080     0.0000 cd         1 LIG       0.441700
             13 N12         -3.5510     1.3100     0.0000 nc         1 LIG      -0.413000
             14 N13         -2.2210     1.6720    -0.0000 nc         1 LIG      -0.416000
             15 O14         -4.9500    -0.5580     0.0010 oh         1 LIG      -0.611300
             16 H2           2.6850    -2.1880     0.0000 h4         1 LIG       0.052600
             17 H6           5.8930     0.4350     0.0000 h5         1 LIG       0.031100
             18 H4           2.0750     2.0500    -0.0000 h4         1 LIG       0.052600
             19 H141        -5.5450     0.2200     0.0010 ho         1 LIG       0.408000
        @<TRIPOS>BOND
             1     1     2 ar
             2     1     6 ar
             3     2     3 ar
             4     2    16 1
             5     3     4 ar
             6     3     7 1
             7     4     5 ar
             8     4    18 1
             9     5     6 ar
            10     6    17 1
            11     7     8 1
            12     7     9 2
            13     9    10 1
            14    10    11 1
            15    10    14 2
            16    11    12 1
            17    12    13 2
            18    12    15 1
            19    13    14 1
            20    15    19 1
        @<TRIPOS>SUBSTRUCTURE
             1 LIG         1 TEMP              0 ****  ****    0 ROOT


.. hint:: 

    A quick refresher on the mol2 format.

    The mol2 format is a widely used file format for representing molecular structures. It contains information about the atoms, bonds, and other molecular features. The key sections of a mol2 file include:

    - **@<TRIPOS>MOLECULE**: This section provides general information about the molecule, including its name, number of atoms, and other properties.
    - **@<TRIPOS>ATOM**: This section lists all the atoms in the molecule, along with their coordinates, atom types, and partial charges.
    - **@<TRIPOS>BOND**: This section defines the bonds between atoms, including bond types and connectivity.
    - **@<TRIPOS>SUBSTRUCTURE**: This section describes any substructures within the molecule, such as functional groups or rings.

    Now, its important to note that antechamber writes a mol2 file with amber atom types instead of elements.

    This means that when you are looking at the mol2 file, you should pay attention to the atom types and ensure they match the expected amber types.



Now, we need to create a frcmod file from the mol2 file. The frcmod file contains additional force field parameters that are not present in the mol2 file or the GAFF2 forcefield.

.. code-block:: bash

    # Call ParmChk2 to create a frcmod file from the mol2 file
    # -i : Input file (ligand.mol2)
    # -o : Output file (ligand.frcmod)
    # -f : Input file format (mol2)
    # -s : ff parm set (2, which is equivalent to gaff2)

    parmchk2 -i ligand.mol2 -o ligand.frcmod -f mol2 -s 2

If everything runs as expected, you should get no output to your screen and you should see a file ligand.frcmod appear in your filesystem.

The contents of that file can be seen below.

.. toggle::

    .. code-block:: console

        Remark line goes here
        MASS

        BOND

        ANGLE
        ca-ce-o    79.900     121.180   same as ce-ce-o , penalty score=  2.2
        cd-nf-ce   84.130     112.590   same as cd-ne-ce, penalty score=  1.0
        nf-ce-o    89.540     125.450   same as ne-cf-o , penalty score= 11.9
        oh-cd-ss   78.580     119.440   same as os-cd-ss, penalty score=  2.0

        DIHE
        ca-ca-ce-o    4    2.800       180.000           2.000      same as X -c2-ca-X , penalty score=237.0
        ca-ca-ce-nf   4    4.000       180.000          -2.000      same as ca-ce-ce-n2
        ca-ca-ce-nf   1    0.750       180.000          -3.000      same as ca-ce-ce-n2
        ca-ca-ce-nf   1    0.820       180.000           1.000      same as ca-ce-ce-n2, penalty score=638.0
        ca-ce-nf-cd   1    0.400         0.000           3.000      same as c2-ce-ne-c2, penalty score=312.0
        ss-cd-nf-ce   2    1.600       180.000           2.000      same as X -cf-nf-X , penalty score=136.0
        nc-cd-nf-ce   2    1.600       180.000           2.000      same as X -cf-nf-X , penalty score=136.0
        o -ce-nf-cd   2    1.600       180.000           2.000      same as X -ce-ne-X , penalty score=138.0
        nf-cd-ss-cd   2    3.200       180.000           2.000      same as n2-c2-ss-ca, penalty score=406.0
        nc-cd-ss-cd   2    3.200       180.000           2.000      same as n2-c2-ss-ca, penalty score=406.0
        oh-cd-ss-cd   2    2.200       180.000           2.000      same as X -c2-ss-X , penalty score=232.0
        ss-cd-oh-ho   2    2.100       180.000           2.000      same as X -c2-oh-X , penalty score=232.0
        nc-cd-oh-ho   2    2.100       180.000           2.000      same as X -c2-oh-X , penalty score=232.0

        IMPROPER
        ca-h4-ca-nb        10.5          180.0         2.0          Same as X -X -cc-X , penalty score= 41.9 (use general term))
        ca-ca-ca-ce         1.1          180.0         2.0          Same as c2-ca-ca-ca, penalty score= 61.9)
        h5-nb-ca-nb        10.5          180.0         2.0          Same as X -X -cc-X , penalty score= 41.9 (use general term))
        ca-nf-ce-o         10.5          180.0         2.0          Same as X -X -cc-X , penalty score= 37.0 (use general term))
        nc-nf-cd-ss        10.5          180.0         2.0          Using general improper torsional angle  X- X-cd- X, penalty score=  9.0)
        nc-oh-cd-ss        10.5          180.0         2.0          Using general improper torsional angle  X- X-cd- X, penalty score=  9.0)

        NONBON


You'll see that there are a number of different sections. The MASS and BOND sections are empty, because those are all described by the GAFF2 forcefield that we are using. 

The ANGLE, DIHE, and IMPROPER sections all contain entries that were not found in the standard GAFF2 forcefield, so parmchk2 has generated parameters for them based on similar parameters found in the forcefield. The quality of those dihedral parameters is assessed by the penalty score, with lower scores indicating better confidence in the parameters. For terms that have high penalty scores, you might want to consider using electronic sturcture or machine-learning based approaches to improve their parameters.

.. note::

    These are not the only parameters that describe bonded/angular/dihedral interactions - the standard forcefields have files shipped with amber that define their parameters (located in amber/data/leap/gaff2.dat for gaff2).

Now that we have base intramolecular force-field parameters, the last thing we need is a lib file. These files store those partial charges that we provided earlier so that if you are building a system from a docked structure
you aren't forced to re-generate a mol2 with the correct coordinates and can use a pdb instead.

Now we need to create the lib file, to do that we will use a program called tleap.

Write a quick file, called tleap.in, that contains the following commands.

.. code-block:: bash
    
    source leaprc.gaff2
    loadamberparams ligand.frcmod
    LIG = loadmol2 ligand.mol2

    saveOff LIG ligand.lib
    quit

Now, we can run tleap on this file as:

.. code-block:: bash

    tleap -f leap.in

When you run this, you will see output on your screen that looks something like this.

.. toggle::

    .. code-block:: console

        -I: Adding /home/piskulic/Project/SoftwareDevelopment/install_lbsr_dev/dat/leap/prep to search path.
        -I: Adding /home/piskulic/Project/SoftwareDevelopment/install_lbsr_dev/dat/leap/lib to search path.
        -I: Adding /home/piskulic/Project/SoftwareDevelopment/install_lbsr_dev/dat/leap/parm to search path.
        -I: Adding /home/piskulic/Project/SoftwareDevelopment/install_lbsr_dev/dat/leap/cmd to search path.
        -f: Source tleap.in.

        Welcome to LEaP!
        (no leaprc in search path)
        Sourcing: ./tleap.in
        ----- Source: /home/piskulic/Project/SoftwareDevelopment/install_lbsr_dev/dat/leap/cmd/leaprc.gaff2
        ----- Source of /home/piskulic/Project/SoftwareDevelopment/install_lbsr_dev/dat/leap/cmd/leaprc.gaff2 done
        Log file: ./leap.log
        Loading parameters: /home/piskulic/Project/SoftwareDevelopment/install_lbsr_dev/dat/leap/parm/gaff2.dat
        Reading title:
        AMBER General Force Field for organic molecules (Version 2.2.30, Oct 2025)
        Loading parameters: ./ligand.frcmod
        Reading force field modification type file (frcmod)
        Reading title:
        Remark line goes here
        Loading Mol2 file: ./ligand.mol2
        Reading MOLECULE named LIG
         Creating ligand.lib
        Building topology.
        Building atom parameters.
        	Quit

        Exiting LEaP: Errors = 0; Warnings = 0; Notes = 0.

Everytime you parameterize with tleap, you should check the Errors and Warnings section to see whether any issues were reported. This can help you catch potential problems early and ensure that your parameterization is as accurate as possible.

Now - lets visualize the charges using a simple python script and VMD. Note - you may want to work on this on your local machine, due to occaisional challenges with X forwarding.

.. toggle::

    .. literalinclude:: /_static/files/ModularTutorials/ForceField/ligand_parameterization/mol2_to_pdb_bfactor.py
        :language: python
        

This python script can be run on your mol2 file to produce a pdb file with B-factors set to the partial charges of the atoms. To run it, do as follows:

.. code-block:: bash

    python mol2_to_pdb_bfactor.py ligand.mol2 ligand_charges.pdb

You'll see an output file, ligand_charges.pdb. 

Now you can visualize the charges using VMD or any other molecular visualization tool.

You can either do this by hand, OR, you can use the provided VMD tcl script below.

.. toggle::

    .. literalinclude:: /_static/files/ModularTutorials/ForceField/ligand_parameterization/visualize_charges.tcl
        :language: tcl


.. code-block:: bash

    vmd -e visualize_charges.tcl -args ligand_charges.pdb -1.0 1.0 


If things worked out - you should see something like this.

.. image:: /_static/files/ModularTutorials/ForceField/ligand_parameterization/byhand.png
    :alt: Ligand Charges Visualization
    :width: 400px

Here, we are defining the color range for the visualization, with -1.0 being the minimum and 1.0 being the maximum. 



Now this is an easy parameterization, and this already involves running three different programs: antechamber, parmchk2, and tleap. Now imagine you add in the by-hand adjustments described above, and also add in RESP charge fitting.
It is pretty clear that this process can get complicated quickly, it is easy to make mistakes, and keeping track of all the different files and parameters can be challenging.


Parameterizing with RESP Charges
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. warning::

    This section requires access to the Gaussian Electronic Structure Program. If you do not have access to Gaussian, you can skip this section and move on to the LigandParam tutorial.

    On {{DISPLAY_NAME}}, Gaussian is available through the following command:

    .. code-block:: console

        {{GAUSSIAN_LOAD}}


.. warning::

    Gaussian runs take up quite a few resources and should not be run on a login node. Please request an interactive node if you are going to follow along with this tutorial.

    .. code-block:: console

        {{INTERACTIVE_REQUEST}}



One of the problems with AM1-BCC charges is that they are not always very accurate, especially for molecules with unusual electronic structures. A better way to get partial charges is to use RESP charge fitting, which fits charges to reproduce the electrostatic potential around the molecule as calculated by a quantum chemistry program. This is a more involved process, but can yield much better results.

Lets go through the steps to do this. We will start from the same pdb file that we used above.

Lets set up a working directory for RESP charge fitting. Again,ligand.pdb should be located in your working directory.

.. code-block:: bash

    mkdir -p resp_charge_fitting
    cd resp_charge_fitting
    cp ../ligand.pdb .

First, we need to generate a Gaussian input file to calculate the electrostatic potential. We can do this with antechamber as follows:

.. code-block:: bash

    # Antechamber to generate Gaussian input file
    # -i input file
    # -fi input file format
    # -o output file
    # -fo output file format
    # -gv Gaussian adds keyword to generate Gaussian input file
    # -ge Gaussian output file
    # -nc net charge
    # -at atom type
    antechamber -i ligand.pdb -fi pdb -o ligand.gjf -fo gcrt -gv 1 -ge ligand.log -nc -1 -at gaff2


The options here are the same as before, with the addition of the gaussian options , which tells antechamber to generate a Gaussian input file for calculating the electrostatic potential.

.. note:: Details on Gaussian Input File

    .. toggle:: 

        The generated Gaussian input file (ligand.gjf) looks like this.

        .. code-block:: text

            --Link1--
            %chk=molecule
            #HF/6-31G* SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6) opt
            # iop(6/50=1)

            remark line goes here

            -1   1
                N    4.4000000000       -0.9970000000        0.0000000000
                C    3.0640000000       -1.1690000000        0.0000000000
                C    2.1640000000       -0.1030000000       -0.0000000000
                C    2.7240000000        1.1780000000       -0.0000000000
                N    4.0540000000        1.3860000000       -0.0000000000
                C    4.8150000000        0.2790000000        0.0000000000
                C    0.6790000000       -0.3620000000       -0.0000000000
                O    0.2940000000       -1.5650000000       -0.0010000000
                N   -0.0490000000        0.7600000000       -0.0000000000
                C   -1.3960000000        0.6280000000       -0.0000000000
                S   -2.2720000000       -0.9500000000       -0.0000000000
                C   -3.7050000000        0.0080000000        0.0000000000
                N   -3.5510000000        1.3100000000        0.0000000000
                N   -2.2210000000        1.6720000000       -0.0000000000
                O   -4.9500000000       -0.5580000000        0.0010000000
                H    2.6850000000       -2.1880000000        0.0000000000
                H    5.8930000000        0.4350000000        0.0000000000
                H    2.0750000000        2.0500000000       -0.0000000000
                H   -5.5450000000        0.2200000000        0.0010000000

            ligand.log

            ligand.log


        Lets go through these sections:

        - The first section contains the Gaussian route section, which specifies the level of theory and basis set to be used (HF/6-31G*), as well as other options for the calculation.
            - SCF=tight: Tight convergence criteria for the self-consistent field (SCF) procedure.
            - Test: This is a keyword that tells Gaussian to perform a test calculation.
            - Pop=MK: This specifies that the electrostatic potential should be calculated using the Merz-Kollman scheme.
            - iop(6/33=2): This option specifies that the electrostatic potential should be calculated on a grid of points around the molecule.
            - iop(6/42=6): This option specifies that the electrostatic potential should be calculated using a grid spacing of 0.2 angstroms.
            - opt: This tells Gaussian to optimize the geometry of the molecule before calculating the electrostatic potential.
            - iop(6/50=1): This option specifies that the electrostatic potential should be calculated using a single point calculation.
        - The second section is a comment line, which can be used to describe the calculation. "remark line goes here"
        - The third section contains the charge and multiplicity of the molecule. In this case, the molecule has a net charge of -1 and a multiplicity of 1 (singlet state).
        - The fourth section contains the Cartesian coordinates of the atoms in the molecule.

Now, there are a few modifications you can make if you are interested in doing so.

1. You can add %mem=10GB to in crease the memory allocated for the Gaussian job.
2. You can add %nproc=2 to specify that you want to use multiple processors
3. You can add p after the # to indicate that you want to use the parallel version of Gaussian (right before HF/6-31G*, as #p HF/6-31G* ...)

Next, we need to run Gaussian on this input file to generate the electrostatic potential. This can be done with the following command:

.. code-block:: bash

    g16 < ligand.gjf > ligand_run.log


This will run Gaussian and output the results to ligand.log. Once this is done, we can use antechamber again to extract the electrostatic potential and fit RESP charges.


.. warning::

    This will take some time! On 10 cores, this took 42 minutes on our test system.

Now we will take ligand_run.log and use it to generate the RESP charges.

.. code-block:: bash

    antechamber -i ligand_run.log -o ligand_resp.mol2 -fi gout -fo mol2 -c resp -rn LIG -at gaff2 -nc -1


Now lets compare the charges generated by RESP to those generated by AM1-BCC. You can do this by looking at the mol2 files generated in each case.

Finally, we can generate the frcmod and lib files as before, using the RESP-charged mol2 file.

.. code-block:: bash

    parmchk2 -i ligand_resp.mol2 -o ligand_resp.frcmod -f mol2 -s 2

Then create tleap.in again.

.. code-block:: text

    # tleap input file
    source leaprc.gaff2
    loadamberparams ligand_resp.frcmod
    LIG = loadmol2 ligand_resp.mol2

    saveOff LIG ligand_resp.lib
    quit

And run tleap. 

.. code-block:: bash

    tleap -f leap.in

Like before, we can visualize the ligand in VMD using the same procedure as above. 


.. code-block:: bash

    python mol2_to_pdb_bfactor.py ligand_resp.mol2 resp_charges.pdb
    vmd -e visualize_charges.tcl -args resp_charges.pdb -1.0 1.0


You can see that the RESP charges have the nitrogens a bit more negatively charged than the AM1-BCC charges. 

.. image:: /_static/files/ModularTutorials/ForceField/ligand_parameterization/resp.png
    :alt: Ligand Charges Visualization
    :width: 400px

To make this clearer, we can visualize the difference in charges between the two methods. 

To do this, you can use this python code on the two pdb files generated with AM1-BCC and RESP.

Save this file as pdb_bfactor_diff.py.

.. toggle::

    .. literalinclude:: /_static/files/ModularTutorials/ForceField/ligand_parameterization/pdb_bfactor_diff.py
        :language: python

You can run this script with the following command.

.. code-block:: bash

    python pdb_bfactor_diff.py ../ligand_charges.pdb resp_charges.pdb diff.pdb
     vmd -e visualize_charges.tcl -args diff.pdb -0.2 0.2

Note that we've reduced the range for the color scale since we are now examining the difference 


.. image:: /_static/files/ModularTutorials/ForceField/ligand_parameterization/diff.png
    :alt: Ligand Charges Visualization
    :width: 400px


.. note::

    Because RESP calculates charges based on a finite grid around the molecule, the charges that are calculated can be sensitive to the orientation of the molecule. It is often a good idea to run RESP multiple times with different orientations and average the results to get more reliable charges. However, this requires customized programs or scripts to do so and is out of the scope of this tutorial. See the ligandparam tutorial for an automated way to do this more complicated procedure.


.. end-tutorial