Tutorial 2: Ligand Preparation =============================== .. contents:: :local: In this tutorial, we will show you how to use pdbtop to prepare a ligand. For any nonstandard ligand, you can follow the procedure in this tutorial. .. hint:: If you use the procedure in this tutorial, please cite the following papers: - Zhang, J.* `Atom Typing Using Graph Representation Learning: How Do Models Learn Chemistry? `_ *J. Chem. Phys.* **2022**, *156*, 204108. - Zhang, J.*; Tian, L. `Efficient Evaluation of Electrostatic Potential with Computerized Optimized Code. `_ *Phys. Chem. Chem. Phys.* **2021**, *23*, 20323-20328. - Tian, L.* Chen, F. `Multiwfn: A multifunctional wavefunction analyzer. `_ *J. Comput. Chem.* **2012**, *33*, 580-592. Get Structure -------------------------- We want to study the ligand ``92V`` in ``5VBM.pdb`` mentioned in :doc:`tutorial-1`. In original ``5vbm.pdb``, there are 2 ligands ``GDP`` and ``92V``. In this study, we only consider ``92V``. So, we need to extract the ligand ``92V`` from the original file and save it to a new file, say ``92v.pdb``: .. code-block:: bash :caption: 92v.pdb :linenos: HETATM 2707 C10 92V A 203 5.681 6.839 11.647 1.00 16.14 C C HETATM 2708 C13 92V A 203 6.936 8.149 9.504 1.00 24.30 C C HETATM 2709 C21 92V A 203 2.291 0.347 17.510 1.00 41.07 C C HETATM 2710 C22 92V A 203 1.970 -0.529 16.333 1.00 44.07 C C HETATM 2711 C01 92V A 203 4.544 2.998 15.933 1.00 18.17 C C HETATM 2712 C02 92V A 203 4.885 3.156 14.589 1.00 16.22 C C HETATM 2713 C05 92V A 203 4.759 4.228 16.562 1.00 22.02 C C HETATM 2714 C06 92V A 203 5.755 5.073 13.263 1.00 16.65 C C HETATM 2715 C07 92V A 203 6.887 4.558 12.658 1.00 15.45 C C HETATM 2716 C08 92V A 203 7.414 5.162 11.535 1.00 18.98 C C HETATM 2717 C09 92V A 203 6.817 6.303 11.035 1.00 16.08 C C HETATM 2718 C11 92V A 203 5.140 6.222 12.771 1.00 16.49 C C HETATM 2719 C14 92V A 203 4.812 2.130 13.492 1.00 18.17 C C HETATM 2720 C18 92V A 203 4.042 1.727 16.526 1.00 18.18 C C HETATM 2721 C23 92V A 203 2.091 -2.011 16.414 1.00 34.51 C C HETATM 2722 F15 92V A 203 4.352 2.697 12.388 1.00 17.16 F F HETATM 2723 F16 92V A 203 4.003 1.128 13.799 1.00 16.57 F F HETATM 2724 F17 92V A 203 6.023 1.607 13.258 1.00 18.36 F F HETATM 2725 N03 92V A 203 5.280 4.422 14.439 1.00 18.78 N N HETATM 2726 N04 92V A 203 5.233 5.093 15.624 1.00 21.26 N N HETATM 2727 N20 92V A 203 2.726 1.642 17.070 1.00 21.82 N N HETATM 2728 O12 92V A 203 7.352 6.895 9.910 1.00 21.21 O O HETATM 2729 O19 92V A 203 4.735 0.718 16.505 1.00 20.16 O O HETATM 2730 S24 92V A 203 0.818 -2.821 17.310 1.00 24.90 S S HETATM 2731 H10192V A 203 5.263 7.634 11.281 1.00 19.37 H H HETATM 2732 H13292V A 203 6.090 8.071 9.022 1.00 29.16 H H HETATM 2733 H13192V A 203 7.609 8.541 8.918 1.00 29.16 H H HETATM 2734 H13392V A 203 6.814 8.719 10.285 1.00 29.16 H H HETATM 2735 H21292V A 203 3.000 -0.065 18.031 1.00 49.28 H H HETATM 2736 H21192V A 203 1.497 0.429 18.063 1.00 49.28 H H HETATM 2737 H22192V A 203 2.543 -0.240 15.602 1.00 52.88 H H HETATM 2738 H22292V A 203 1.054 -0.333 16.076 1.00 52.88 H H HETATM 2739 H05192V A 203 4.620 4.423 17.500 1.00 26.43 H H HETATM 2740 H07192V A 203 7.307 3.761 13.019 1.00 18.54 H H HETATM 2741 H08192V A 203 8.202 4.795 11.107 1.00 22.78 H H HETATM 2742 H11192V A 203 4.359 6.591 13.210 1.00 19.79 H H HETATM 2743 H20192V A 203 2.159 2.350 17.059 1.00 26.18 H H HETATM 2744 H23292V A 203 2.941 -2.225 16.833 1.00 41.41 H H HETATM 2745 H23192V A 203 2.098 -2.365 15.508 1.00 41.41 H H Generate Topology -------------------------- There is not preset topology for this ligand. Since hydrogen atoms are already there, we do not need to add other ones. Just guess topology with the following command: .. code-block:: bash $ pdbtop ligand -i 92v.pdb -o 92v-1 The output are shown below: .. code-block:: bash $ pdbtop ligand -i 92v.pdb -o 92v-1 Guess topology and parameters ... Warning: Impropers cannot be generated automatically. You have to add them manually if necessary. Guess topology and parameters done. Write PDB: 92v-1.pdb Write PSF: 92v-1.psf Write naive parameters: 92v-1-naive.prm Total charge: -2.39000 Now we have the topology for the ligand in the file ``92v-1.psf``. #. pdbtop can generate bonds based on the distance between atoms using covalent bond radii. However, it cannot generate the impropers automatically. You have to add them manually if necessary. pdbtop will give you a warning if there are any impropers in the ligand. #. pdbtop will type the atoms automatically with a graph neural network. This will **NOT** guarantee the 100% correctness of the typing. You should check the typing carefully. Also, if there are some atoms that rarely appeared in biological systems (like Pd), pdbtop will definitely fail. You must manually type them. #. The force field parameters are generated by searching for the most similar atom types in CHARMM36 force field. In ``92v-1-navie.prm``, you will see this: .. code-block:: bash :caption: 92v-1-naive.prm :linenos: ... BONDS CG2O1 NG2S1 370.00000000 1.34500000 CG2O1 OG2D1 620.00000000 1.23000000 CG2O6 FGA3 300.00000000 1.32359548 ! Guess CG2O6 FGA3 300.00000000 1.32390861 ! Guess CG2O6 FGA3 300.00000000 1.33970370 ! Guess CG321 HGA2 309.00000000 1.11100000 CG321 CG321 222.50000000 1.53000000 CG321 CG323 300.00000000 1.48913599 ! Guess ... The parameters with ``! Guess`` are the ones that are guessed by pdbtop. You should check them carefully. An important point is that the total charge is ``-2.39``. This is ridiculous. We must do a calculation to get reasonable charges. Generate Charges -------------------------- The ligand ``92V`` is shown below: .. image:: _static/figs/p9.png Since it is covalently bonded to the protein, there is a dangling bond at the sulfur atom. We need to saturate it with a hydrogen atom. You have to add it manually with `Qbics Mol* `_. The new ligand is shown below: .. image:: _static/figs/p10.png Then we perform a standard quantum chemical calcualtion at B3LYP/def2-SVP level of theory (or other method that is suitable for your system). We can use `Qbics `_ to do this, with the following input: .. code-block:: bash :caption: qbics.inp :linenos: basis def2-svp end scf charge 0 # The total charge. spin2p1 1 end task energy b3lyp end mol C 5.68100 6.83900 11.64700 C 6.93600 8.14900 9.50400 C 2.29100 0.34700 17.51000 C 1.97000 -0.52900 16.33300 C 4.54400 2.99800 15.93300 C 4.88500 3.15600 14.58900 C 4.75900 4.22800 16.56200 C 5.75500 5.07300 13.26300 C 6.88700 4.55800 12.65800 C 7.41400 5.16200 11.53500 C 6.81700 6.30300 11.03500 C 5.14000 6.22200 12.77100 C 4.81200 2.13000 13.49200 C 4.04200 1.72700 16.52600 C 2.09100 -2.01100 16.41400 F 4.35200 2.69700 12.38800 F 4.00300 1.12800 13.79900 F 6.02300 1.60700 13.25800 N 5.28000 4.42200 14.43900 N 5.23300 5.09300 15.62400 N 2.72600 1.64200 17.07000 O 7.35200 6.89500 9.91000 O 4.73500 0.71800 16.50500 S 0.81800 -2.82100 17.31000 H 5.26300 7.63400 11.28100 H 6.09000 8.07100 9.02200 H 7.60900 8.54100 8.91800 H 6.81400 8.71900 10.28500 H 3.00000 -0.06500 18.03100 H 1.49700 0.42900 18.06300 H 2.54300 -0.24000 15.60200 H 1.05400 -0.33300 16.07600 H 4.62000 4.42300 17.50000 H 7.30700 3.76100 13.01900 H 8.20200 4.79500 11.10700 H 4.35900 6.59100 13.21000 H 2.15900 2.35000 17.05900 H 2.94100 -2.22500 16.83300 H 2.09800 -2.36500 15.50800 H 1.07077 -2.74618 18.59320 end Of course, if the ligand is charged, say with -2, you should change ``charge 0`` to ``charge -2``. Run the calculation with the following command: .. code-block:: bash $ qbics qbics.inp -n 32 > qbics.out # -n is the number of threads After calculation, there will be a file called ``qbics.mwfn``. Then, we can use `Multiwfn `_ to calculate RESP charges: .. code-block:: bash $ Multiwfn qbics.mwfn 7 18 1 y Then, charages will saved in ``qbics.chg``. Like this: .. code-block:: bash :caption: qbics.chg :linenos: ... O 7.352000 6.895000 9.910000 -0.4150163963 O 4.735000 0.718000 16.505000 -0.4750785142 S 0.818000 -2.821000 17.310000 -0.4133866976 H 5.263000 7.634000 11.281000 0.1341722006 ... H 2.098000 -2.365000 15.508000 0.0597608406 H 1.070770 -2.746180 18.593200 0.2031432058 The last column is the charge. Since the sulfur atom in ``92V`` does not have a hydrogen atom, we can add the charge of that hydrogen atom to the sulfur atom. The charge of the hydrogen atom is ``0.2031432058``. So, ``-0.4133866976`` + ``0.2031432058`` = ``-0.2102434918`` is the new charge: .. code-block:: bash :caption: qbics-2.chg :linenos: ... O 7.352000 6.895000 9.910000 -0.4150163963 O 4.735000 0.718000 16.505000 -0.4750785142 S 0.818000 -2.821000 17.310000 -0.2102434918 H 5.263000 7.634000 11.281000 0.1341722006 ... H 2.098000 -2.365000 15.508000 0.0597608406 We replace the charge in ``ATOM`` section in ``92v-1.psf`` file, to get a new psf ``92v-2.psf``: .. code-block:: bash :caption: 92v-2.psf :linenos: 39 !NATOM 1 A 203 92V C10 CG2R61 -0.1884038891 12.0096 0 2 A 203 92V C13 CG331 0.3591542151 12.0096 0 3 A 203 92V C21 CG321 0.3627540705 12.0096 0 4 A 203 92V C22 CG321 -0.2214684281 12.0096 0 5 A 203 92V C01 CG25C1 -0.2977155670 12.0096 0 6 A 203 92V C02 CG2R51 -0.1510380349 12.0096 0 7 A 203 92V C05 CG2R51 0.1852921737 12.0096 0 8 A 203 92V C06 CG2R61 -0.1440665182 12.0096 0 9 A 203 92V C07 CG2R61 0.0533632909 12.0096 0 10 A 203 92V C08 CG2R61 -0.3666550095 12.0096 0 11 A 203 92V C09 CG2R61 0.4065859289 12.0096 0 12 A 203 92V C11 CG2R61 -0.1640553820 12.0096 0 13 A 203 92V C14 CG2O6 0.5702169293 12.0096 0 14 A 203 92V C18 CG2O1 0.6124231757 12.0096 0 15 A 203 92V C23 CG323 0.1112031570 12.0096 0 16 A 203 92V F15 FGA3 -0.1940642847 18.9984 0 17 A 203 92V F16 FGA3 -0.1435770277 18.9984 0 18 A 203 92V F17 FGA3 -0.1875808418 18.9984 0 19 A 203 92V N03 NG2R51 0.4922276791 14.0064 0 20 A 203 92V N04 NG2R50 -0.5104595741 14.0064 0 21 A 203 92V N20 NG2S1 -0.6536415473 14.0064 0 22 A 203 92V O12 OG301 -0.4150163963 15.9990 0 23 A 203 92V O19 OG2D1 -0.4750785142 15.9990 0 24 A 203 92V S24 SG302 -0.2102434918 32.0590 0 25 A 203 92V H101 HGR61 0.1341722006 1.0078 0 26 A 203 92V H132 HGA3 -0.0264696064 1.0078 0 27 A 203 92V H131 HGA3 -0.0264696064 1.0078 0 28 A 203 92V H133 HGA3 -0.0264696064 1.0078 0 29 A 203 92V H212 HGA2 -0.0259570279 1.0078 0 30 A 203 92V H211 HGA2 -0.0259570279 1.0078 0 31 A 203 92V H221 HGA2 0.0769967516 1.0078 0 32 A 203 92V H222 HGA2 0.0769967516 1.0078 0 33 A 203 92V H051 HGR62 0.1072192902 1.0078 0 34 A 203 92V H071 HGR61 0.1059494334 1.0078 0 35 A 203 92V H081 HGR61 0.1842935479 1.0078 0 36 A 203 92V H111 HGR61 0.1619565257 1.0078 0 37 A 203 92V H201 HGP1 0.3340605792 1.0078 0 38 A 203 92V H232 HGA2 0.0597608406 1.0078 0 39 A 203 92V H231 HGA2 0.0597608406 1.0078 0 Now, the file ``92v-1.pdb`` and ``92v-2.psf`` is ready for use.