Hands-on deMon2k: charge transfer for DNA oxidation



Date: Janvier 2021

N. Gillet, E. Dumont



some_text >
This short, Covid-19 adapted tutoral, aims at illustrating some key notions with QM/MM methods. We will be working with the deMon2K code, which is interfaced with gromacs for QM/MM calculations. You will need to download the following tar file : demon2k-143D.tgz.

We are considering the formation of guanine radical cations, after one-electron reduction. This reaction is ubiquitous in DNA. Depending on the coupling with other residues in DNA (double helix or other topologies such as G-quadruplexes), the ease of formation and stability of such radical cations will depend on which nucleobase is targeted.


A survey of QM/MM methods can be downloaded here : qmmm-edumont-012021.pdf.


1. A first series of questions on the intrinsic system: DFT and C-DFT

The input for a simple DFT optimization is relatively straightforward (more than with other interfaces...) --- see Input.in in the directory G3G4_gas_opt : 

CHARGE 0
MULTIPLICITY 1
BASIS (AUG-CC-PVTZ)   ! basis set 
AUXIS (GEN-A2)	      ! auxiliary basis set 
VXCTYPE AUXIS PBE50   ! PBE50 functional (50% HFX respectively) 
SCFTYPE MAX=500 UKS   ! unrestricted Kohn-Sham 
DISPERSION            ! unrestricted orbitals are crucial here 
OPTIMIZATION          ! Not only the single point energy 
GEOMETRY CARTESIAN ANGSTROM
 C       0.000000     0.000000     0.000000
 O       0.000000     0.000000     1.210205
 N3      1.248155     0.000000    -0.545314
 C       1.457077    -0.071271    -1.903699
 N       2.700820    -0.033121    -2.315889
 N       0.444628    -0.135963    -2.814955
 C      -0.840005    -0.109306    -2.246578
 C      -1.115429     0.004695    -0.926042
 N      -2.482609     0.057508    -0.672784
 C      -2.982821     0.026876    -1.879113
 N11    -2.091876    -0.247238    -2.854179
 H       2.018819    -0.129715     0.093444
 H       3.029149     0.098121    -3.261638
 H       3.384345    -0.056976    -1.572649
 H      -4.026280    -0.018678    -2.151500
 H      -2.292374    -0.417990    -3.818878
 N      -0.543664     3.212475    -2.028536
 C      -1.522060     3.139225    -2.887617
 N19    -2.791118     3.411752    -2.370013
 C      -2.635910     3.421707    -0.987790
 C      -1.252653     3.389497    -0.800836
 N      -3.595757     3.556044     0.007900
 C      -3.023543     3.740982     1.182237
 N      -3.868228     3.570127     2.161453
 N25    -1.669804     3.611943     1.500200
 C      -0.701583     3.425177     0.546833
 O       0.472135     3.290805     0.991306
 H      -1.327566     3.626097     2.449783
 H      -4.843110     3.559984     1.899279
 H      -3.590714     3.333756     3.103674
 H      -1.341013     3.005184    -3.944016
 H      -3.643995     3.536219    -2.877073
#
CONSTANTS
N19 XYZ
N11 XYZ
N3 XYZ
N25 XYZ
Another related directory G3G4_gas_cdft calls for constrained DFT (C-DFT), this time on the radical cation of a guanine dimer (hence with a delocalization of the charge and spin character). 
CHARGE 1
MULTIPLICITY 2
BASIS (AUG-CC-PVTZ)
AUXIS (GEN-A2)
VXCTYPE AUXIS PBE50
SCFTYPE MAX=500 UKS
CNSDFT TOL=0.0001 HIRSH		! C-DFT with Hirshfeld charges 
0.0 DONOR  1-16			! one guanine as donor 
1.0 ACCEPTOR  17-32		! one guanine as acceptor 
ETMOD
TOPOLOGY FINE
GRID FINE
DISPERSION
POPULATION HIRSH		! Hirshfeld population analysis 
GEOMETRY CARTESIAN ANGSTROM

Questions:
  1. To which family belongs the PBE50 functional ?
  2. Comment on the quality of the basis set. What is the reason to resort to an auxiliary basis set here ?
  3. Based on the outputs provided, estimate the energy for the one-electron oxidation of GG. This energy could be compared to the ionization potential of guanine as a mononer (link to the NIST database).
  4. Which constraint is applied along the optimizatios here ? Do you see why it is applied ?



2. A QM/MM study for a solvated system.


We now want to situate the very same reaction, this time in the macromolecular environment as guanines G3 and G4 are embedded in a G-quadruplex. The structure chosen corresponds to the PDB ID 143D. It is not directly used, since one needs to include a box of water molecules and counterions, but usually treated beforehand with a short classical MD simulation from which one or several representative trajectories are extracted. The structure can be seen by downloading the pdb file : 143d_hieragglo_all_rep.c0.pdb.

The input to specify a very similar calculation now looks like: 
CHARGE 0
MULTIPLICITY 1
BASIS (AUG-CC-PVTZ)
AUXIS (GEN-A2)
VXCTYPE AUXIS PBE50	
SCFTYPE MAX=500 UKS	! 500 cycles of Kohn-Sham 
DISPERSION
QM/MM FF=AMBER-FF99	! for the force field 
 N1 1161
 C2 1165
 H3 1175
 N4 1164
 C5 1163
 C6 1169
 O7 1174
 N8 1168
 H9 1170
 C10 1167
 N11 1171
 H12 1172
 H13 1173
 N14 1166
 C15 1162
FORCEFIELD FF=AMBER-FF99
GEOMETRY CARTESIAN ANGSTROMS GFILE
RESTRAINTS KJ
   NaW0 ATOM 100.  2
   OW0 ATOM 100.  2
   HW0 ATOM 100.  2
Finally, the QM/MM for the corresponding radical cation is given in the directory G3G4_environment_CDFT.

Questions:
  1. Describe the system under investigation (PDB can be opened with VMD, molden, or another software you are familiar with). Which partition scheme has been chosen here ? Provide the numbers of atoms respecitvely described at the QM and MM levels (to be read from the output file, the QM part consists in the guanines G3 and G4).
  2. How is the subsystem (MM part) described here ?
  3. Estimate the difference of energy for one-electron oxidation of GG with the G-quadruplex and compare it to the value for the GG system.
  4. Compute and comment on the charge and spin densities distribution for the radical cation GG given in the output file.
  5. Same question for the QM-only energy (which can be read directly from the ouput). Provide a comparison with the gas phase system.
  6. How would you account for charge tranfer beyond the QM part ?



Last update : Jan. 2021