Direct use of weighted Quantum Chemical Energy for ligands in BUSTER refinement


QM introductory tutorial example:

benzopyran ligand bound to 2i0g

Contents


Background


Starting Files

2i0g.pdb protein/ligand coordinates from the RCSB
2i0g.mtz structure factors
I0G.grade_PDB_ligand.cif grade_PDB_ligand dictionary
grade_PDB_ligand I0G

(A) Preparation: add hydrogen atoms to the ligand(s)

hydrogenate -p 2i0g.pdb \
           -l I0G.grade_PDB_ligand.cif \
           -ligonly -zero \
           -o 2i0g_addhI0G.pdb 
coot -p 2i0g_addhI0G.pdb

(B) Refine with -qm I0G -l I0G.grade_PDB_ligand.cif

time refine                 \
-p 2i0g_addhI0G.pdb         \
-m 2i0g.mtz                 \
-autoncs                    \
-l I0G.grade_PDB_ligand.cif \
-qm I0G                     \
-d partB > partB.log &
firefox partB/01-BUSTER/Cycle-1/LIST.html
first QM calculation of energy/gradient
.....
    QM_HELPER_LOG: Hamiltonian            =            RM1  Number of QM Atoms     =             39
    QM_HELPER_LOG: System Charge          =              0  System Multiplicity    =              1
    QM_HELPER_LOG: Number of Orbitals     =            102  Occupied Orbitals      =             54
    QM_HELPER_LOG: Number Alpha Electrons =             54  Number Beta Electrons  =             54
    QM_HELPER_LOG: Number of Shells       =             60  Number of Primitives   =            360
    QM_HELPER_LOG: Number of 2e Integrals =          24933  Energy Base Line       =     333346.698
    QM_HELPER_LOG: No. of Boundary Atoms  =              0  Matrix Size            =           5253

graph_autobuster_QM -d partB -launch
graph_autobuster_R -d partB -launch
visualise-geometry-coot partB
visualise-geometry-COOT partB

  time refine                 \
  -p partB/refine.pdb         \
  -m 2i0g.mtz                 \
  -autoncs                    \
  -l I0G.grade_PDB_ligand.cif \
  -qm I0G                     \
  -M TLSbasic                 \
  -d partB2 > partB2.log      &


(C) Control refine with dictionary -l I0G.grade_PDB_ligand.cif

time refine                 \
-p 2i0g_addhI0G.pdb         \
-m 2i0g.mtz                 \
-autoncs                    \
-l I0G.grade_PDB_ligand.cif \
-d partC > partC.log &

(D) Analysis comparing the QM and dictionary refine jobs

What original 2i0g.pdb partB -qm partC GRADE dictionary
CPU secs (1) - 1684 1638
Rwork (2) 0.2329 0.2157 0.2173
Rfree (2) 0.2677 0.2524 0.2512
MolProbity remark 40 clashscore 9.96 0.28 0.28
MolProbity remark 40 Ramachandran favoured 95.1% 97.9% 98.4%
RMS bond deviation in Å 0.0117 0.0092 0.0090
RMS angle deviation in degs. 1.35 1.06 1.04
Final X-ray weight - 6.90 6.52
    1. refine.log reported CPU time on a Intel(R) Core(TM)2 CPU 6700 @ 2.66GHz (4 year old basic workstation)
    2. BUSTER values

Let's now look at the I0G ligand, concentrating on the A site (because the B site is really pretty much identical). First let's extract the A site ligand from the various coordinate files (together with the CRYST1 card so the resultant PDB file works with crystallographic tools).

egrep "CRYST1|I0G A" 2i0g.pdb > 2i0g.I0G_A.pdb
egrep "CRYST1|I0G A" partB/refine.pdb > partB.I0G_A.pdb
egrep "CRYST1|I0G A" partC/refine.pdb > partC.I0G_A.pdb

We can assess how well the different ligand conformations fit to the density by eye using COOT with the result that they appear very similar and fit the density about equally well. A correlation coefficient to the BUSTER 2Fo-Fc map can be found by running:

corr -p partC.I0G_A.pdb -m partB/refine.mtz

These confirm that the density fits are equally good.

What original 2i0g.pdb partB -qm partC GRADE dictionary
I0G ligand CC (to partB/refine.mtz) 0.937 0.938 0.937

What about the geometry of the ligand? This can be assessed by the CCDC Mogul program, compared to similar CSD structures. We have found that a useful metric is obtained by running Mogul in default mode and searching for "All fragments", "All bond fragments", sorting by Z and finding the number of bonds with Z above 2, as well as the average Z score (using a spreadsheet such as OpenOffice or gnumeric).

What original 2i0g.pdb partB -qm partC GRADE dictionary
I0G ligand site A Mogul number of bonds with Z score > 2.0 5/24 0/24 0/24
I0G ligand site A Mogul average Z score for bonds 1.3 0.49 0.19


(E) Assessing ligand strain energy and motions for the QM result


   egrep "I0G A|CRYST1" partB/refine.pdb > partB.IOGA.pdb

gelly_refine -p partB.IOGA.pdb \
    -l I0G.grade_PDB_ligand.cif -qm I0G \
   -d partE -o partE_out.pdb \
   -max 1000 -glim 0.01 -jiggle_xyz 0.01 > partE.log &

graph_autobuster_QM -d partE
...
# final raw QM energies   =   -405.27 kJ/mol

graph_autobuster_QM -d partB -qmzero -405.27
...
# final raw QM energies   =   -399.62  -399.61 kJ/mol
# -qmzero specified so report QM energies relative to -405.27
#     big   small   iters      QMener_01      QMener_02
#       5      33     346           5.65           5.66


Back to Ligand QM/force field top index page

Page by Oliver Smart and Ian Tickle from 14 Feb 2011, modified November 2013. Any questions regarding our software or this wiki should be directed to buster-develop@globalphasing.com