Tutorial 2: Ligand Preparation
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 Tutorial 1: A Standard Protein. 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:
1HETATM 2707 C10 92V A 203 5.681 6.839 11.647 1.00 16.14 C C
2HETATM 2708 C13 92V A 203 6.936 8.149 9.504 1.00 24.30 C C
3HETATM 2709 C21 92V A 203 2.291 0.347 17.510 1.00 41.07 C C
4HETATM 2710 C22 92V A 203 1.970 -0.529 16.333 1.00 44.07 C C
5HETATM 2711 C01 92V A 203 4.544 2.998 15.933 1.00 18.17 C C
6HETATM 2712 C02 92V A 203 4.885 3.156 14.589 1.00 16.22 C C
7HETATM 2713 C05 92V A 203 4.759 4.228 16.562 1.00 22.02 C C
8HETATM 2714 C06 92V A 203 5.755 5.073 13.263 1.00 16.65 C C
9HETATM 2715 C07 92V A 203 6.887 4.558 12.658 1.00 15.45 C C
10HETATM 2716 C08 92V A 203 7.414 5.162 11.535 1.00 18.98 C C
11HETATM 2717 C09 92V A 203 6.817 6.303 11.035 1.00 16.08 C C
12HETATM 2718 C11 92V A 203 5.140 6.222 12.771 1.00 16.49 C C
13HETATM 2719 C14 92V A 203 4.812 2.130 13.492 1.00 18.17 C C
14HETATM 2720 C18 92V A 203 4.042 1.727 16.526 1.00 18.18 C C
15HETATM 2721 C23 92V A 203 2.091 -2.011 16.414 1.00 34.51 C C
16HETATM 2722 F15 92V A 203 4.352 2.697 12.388 1.00 17.16 F F
17HETATM 2723 F16 92V A 203 4.003 1.128 13.799 1.00 16.57 F F
18HETATM 2724 F17 92V A 203 6.023 1.607 13.258 1.00 18.36 F F
19HETATM 2725 N03 92V A 203 5.280 4.422 14.439 1.00 18.78 N N
20HETATM 2726 N04 92V A 203 5.233 5.093 15.624 1.00 21.26 N N
21HETATM 2727 N20 92V A 203 2.726 1.642 17.070 1.00 21.82 N N
22HETATM 2728 O12 92V A 203 7.352 6.895 9.910 1.00 21.21 O O
23HETATM 2729 O19 92V A 203 4.735 0.718 16.505 1.00 20.16 O O
24HETATM 2730 S24 92V A 203 0.818 -2.821 17.310 1.00 24.90 S S
25HETATM 2731 H10192V A 203 5.263 7.634 11.281 1.00 19.37 H H
26HETATM 2732 H13292V A 203 6.090 8.071 9.022 1.00 29.16 H H
27HETATM 2733 H13192V A 203 7.609 8.541 8.918 1.00 29.16 H H
28HETATM 2734 H13392V A 203 6.814 8.719 10.285 1.00 29.16 H H
29HETATM 2735 H21292V A 203 3.000 -0.065 18.031 1.00 49.28 H H
30HETATM 2736 H21192V A 203 1.497 0.429 18.063 1.00 49.28 H H
31HETATM 2737 H22192V A 203 2.543 -0.240 15.602 1.00 52.88 H H
32HETATM 2738 H22292V A 203 1.054 -0.333 16.076 1.00 52.88 H H
33HETATM 2739 H05192V A 203 4.620 4.423 17.500 1.00 26.43 H H
34HETATM 2740 H07192V A 203 7.307 3.761 13.019 1.00 18.54 H H
35HETATM 2741 H08192V A 203 8.202 4.795 11.107 1.00 22.78 H H
36HETATM 2742 H11192V A 203 4.359 6.591 13.210 1.00 19.79 H H
37HETATM 2743 H20192V A 203 2.159 2.350 17.059 1.00 26.18 H H
38HETATM 2744 H23292V A 203 2.941 -2.225 16.833 1.00 41.41 H H
39HETATM 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:
$ pdbtop ligand -i 92v.pdb -o 92v-1
The output are shown below:
$ 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:
1...
2BONDS
3CG2O1 NG2S1 370.00000000 1.34500000
4CG2O1 OG2D1 620.00000000 1.23000000
5CG2O6 FGA3 300.00000000 1.32359548 ! Guess
6CG2O6 FGA3 300.00000000 1.32390861 ! Guess
7CG2O6 FGA3 300.00000000 1.33970370 ! Guess
8CG321 HGA2 309.00000000 1.11100000
9CG321 CG321 222.50000000 1.53000000
10CG321 CG323 300.00000000 1.48913599 ! Guess
11...
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:
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:
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:
1basis
2 def2-svp
3end
4
5scf
6 charge 0 # The total charge.
7 spin2p1 1
8end
9
10task
11 energy b3lyp
12end
13
14mol
15 C 5.68100 6.83900 11.64700
16 C 6.93600 8.14900 9.50400
17 C 2.29100 0.34700 17.51000
18 C 1.97000 -0.52900 16.33300
19 C 4.54400 2.99800 15.93300
20 C 4.88500 3.15600 14.58900
21 C 4.75900 4.22800 16.56200
22 C 5.75500 5.07300 13.26300
23 C 6.88700 4.55800 12.65800
24 C 7.41400 5.16200 11.53500
25 C 6.81700 6.30300 11.03500
26 C 5.14000 6.22200 12.77100
27 C 4.81200 2.13000 13.49200
28 C 4.04200 1.72700 16.52600
29 C 2.09100 -2.01100 16.41400
30 F 4.35200 2.69700 12.38800
31 F 4.00300 1.12800 13.79900
32 F 6.02300 1.60700 13.25800
33 N 5.28000 4.42200 14.43900
34 N 5.23300 5.09300 15.62400
35 N 2.72600 1.64200 17.07000
36 O 7.35200 6.89500 9.91000
37 O 4.73500 0.71800 16.50500
38 S 0.81800 -2.82100 17.31000
39 H 5.26300 7.63400 11.28100
40 H 6.09000 8.07100 9.02200
41 H 7.60900 8.54100 8.91800
42 H 6.81400 8.71900 10.28500
43 H 3.00000 -0.06500 18.03100
44 H 1.49700 0.42900 18.06300
45 H 2.54300 -0.24000 15.60200
46 H 1.05400 -0.33300 16.07600
47 H 4.62000 4.42300 17.50000
48 H 7.30700 3.76100 13.01900
49 H 8.20200 4.79500 11.10700
50 H 4.35900 6.59100 13.21000
51 H 2.15900 2.35000 17.05900
52 H 2.94100 -2.22500 16.83300
53 H 2.09800 -2.36500 15.50800
54 H 1.07077 -2.74618 18.59320
55end
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:
$ 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:
$ Multiwfn qbics.mwfn
7
18
1
y
Then, charages will saved in qbics.chg. Like this:
1...
2O 7.352000 6.895000 9.910000 -0.4150163963
3O 4.735000 0.718000 16.505000 -0.4750785142
4S 0.818000 -2.821000 17.310000 -0.4133866976
5H 5.263000 7.634000 11.281000 0.1341722006
6...
7H 2.098000 -2.365000 15.508000 0.0597608406
8H 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:
1...
2O 7.352000 6.895000 9.910000 -0.4150163963
3O 4.735000 0.718000 16.505000 -0.4750785142
4S 0.818000 -2.821000 17.310000 -0.2102434918
5H 5.263000 7.634000 11.281000 0.1341722006
6...
7H 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:
1 39 !NATOM
2 1 A 203 92V C10 CG2R61 -0.1884038891 12.0096 0
3 2 A 203 92V C13 CG331 0.3591542151 12.0096 0
4 3 A 203 92V C21 CG321 0.3627540705 12.0096 0
5 4 A 203 92V C22 CG321 -0.2214684281 12.0096 0
6 5 A 203 92V C01 CG25C1 -0.2977155670 12.0096 0
7 6 A 203 92V C02 CG2R51 -0.1510380349 12.0096 0
8 7 A 203 92V C05 CG2R51 0.1852921737 12.0096 0
9 8 A 203 92V C06 CG2R61 -0.1440665182 12.0096 0
10 9 A 203 92V C07 CG2R61 0.0533632909 12.0096 0
11 10 A 203 92V C08 CG2R61 -0.3666550095 12.0096 0
12 11 A 203 92V C09 CG2R61 0.4065859289 12.0096 0
13 12 A 203 92V C11 CG2R61 -0.1640553820 12.0096 0
14 13 A 203 92V C14 CG2O6 0.5702169293 12.0096 0
15 14 A 203 92V C18 CG2O1 0.6124231757 12.0096 0
16 15 A 203 92V C23 CG323 0.1112031570 12.0096 0
17 16 A 203 92V F15 FGA3 -0.1940642847 18.9984 0
18 17 A 203 92V F16 FGA3 -0.1435770277 18.9984 0
19 18 A 203 92V F17 FGA3 -0.1875808418 18.9984 0
20 19 A 203 92V N03 NG2R51 0.4922276791 14.0064 0
21 20 A 203 92V N04 NG2R50 -0.5104595741 14.0064 0
22 21 A 203 92V N20 NG2S1 -0.6536415473 14.0064 0
23 22 A 203 92V O12 OG301 -0.4150163963 15.9990 0
24 23 A 203 92V O19 OG2D1 -0.4750785142 15.9990 0
25 24 A 203 92V S24 SG302 -0.2102434918 32.0590 0
26 25 A 203 92V H101 HGR61 0.1341722006 1.0078 0
27 26 A 203 92V H132 HGA3 -0.0264696064 1.0078 0
28 27 A 203 92V H131 HGA3 -0.0264696064 1.0078 0
29 28 A 203 92V H133 HGA3 -0.0264696064 1.0078 0
30 29 A 203 92V H212 HGA2 -0.0259570279 1.0078 0
31 30 A 203 92V H211 HGA2 -0.0259570279 1.0078 0
32 31 A 203 92V H221 HGA2 0.0769967516 1.0078 0
33 32 A 203 92V H222 HGA2 0.0769967516 1.0078 0
34 33 A 203 92V H051 HGR62 0.1072192902 1.0078 0
35 34 A 203 92V H071 HGR61 0.1059494334 1.0078 0
36 35 A 203 92V H081 HGR61 0.1842935479 1.0078 0
37 36 A 203 92V H111 HGR61 0.1619565257 1.0078 0
38 37 A 203 92V H201 HGP1 0.3340605792 1.0078 0
39 38 A 203 92V H232 HGA2 0.0597608406 1.0078 0
40 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.