Enzyme modelling from a NAMD setup
In this tutorial, we will perform a QM/MM calculation on a cytochrome P450 enzyme complexed with camphor in the Compound I state. A CHARMM forcefield will again be used, but this time starting from input files generated by NAMD.
1. NAMD setup
We start with NAMD format structure input files which can be found in samples/p450_namd/namd_setup/. These are start.pdb containing the coordinates, and start.psf, the protein structure file.
The ChemShell script namd_setup.chm will prepare the initial system ready for a ChemShell QM/MM calculation. The PDB and PSF files are read in with the following commands:
read_xplor_psf atom_types= types atom_charges= charges atom_masses= mass file= start.psf read_pdb coords= charmm.c file= start.pdb
The command read_xplor_psf imports PSF files in the XPLOR format used by NAMD/VMD. The atom types, charges and masses are saved into the variable names provided, and in our case written to disk in the file save.chm for later use.
2. Quartet state QM/MM energy calculation
Once the PDB and PSF information has been imported a ChemShell QM/MM calculation can be performed on the P450 system. The quartet state quartet/quartet.chm should be run first.
To begin with the previously prepared setup information is read in, together with the locations of the CHARMM forcefield topology and parameter files in the toppar directory.
# Read in charges, atom types, etc. These files are created by namd_setup.chm
source ../namd_setup/save.chm
source ../namd_setup/selection.chm
# Define topology and parameter files
set top { ../toppar/top_all22_heme.inp ../toppar/camphor.rtf }
set par ../toppar/par_mod.inp
# Read starting coordinates from PDB file
read_pdb coords=start.c file=$file_name.pdb
The QM region is then defined (by default just the haem but camphor can be added for a more involved calculation) together with DFT options. In this case we use the B3LYP functional with 6-31G on all atoms except Fe, where the LANL2DZ ECP basis is used. For the quartet state a UHF calculation with spin multiplicity of 4 is specified.
set qmatoms { 4576 - 4579, 5148 - 5177 }
# QM calculation options
set qmflags { hamiltonian=dft functional=b3lyp \
basisspec = { { 6-31g { c n o h s } } { lanl2dz fe } } \
mult=4 scftype=uhf \
maxcyc=200 }
The QM/MM calculation is then ready to run. In the example provided we use ORCA as the QM code; other codes can be used by changing the theory= option.
energy coords=start.c \
theory=hybrid : [ list \
coupling=shift \
atom_charges= $charges \
qm_theory=orca : $qmflags \
qm_region= $qmatoms \
mm_theory=dl_poly : [ list \
list_option=none \
exact_srf=yes \
use_pairlist=no \
cutoff=1000 \
scale14 = { 1.0 1.0 } \
atom_types= $types \
use_charmm_psf=yes \
charmm_psf_file= $file_name.psf \
charmm_parameter_file= $par \
charmm_mass_file= $top ] ]
Here we are using electrostatic embedding (specified by coupling=shift). This means that atom_charges should be specified in the hybrid options section, not in mm_theory, so that the hybrid driver can delete charges on the QM region atoms as appropriate.
The atom types and the CHARMM forcefield information are loaded in much the same way as in the previous tutorial, and the other DL_POLY options are essentially identical.
As the QM region is relatively large, the calculation will take approximately 30 minutes to run on a single processor. If you can run ORCA in an MPI parallel environment, ChemShell can take advantage of this by specifying the nproc= option to run ORCA in parallel.
In addition to the usual information provided by the ChemShell output, it can also be helpful to look at the output of the QM calculation, here in orca1.out. ORCA performs a Mulliken population analysis which we can use to check that the correct electronic state has been found.
-------------------------------------------- MULLIKEN ATOMIC CHARGES AND SPIN POPULATIONS -------------------------------------------- 0 C : -0.649175 -0.017855 1 H : 0.172804 0.016469 2 H : 0.215423 0.004285 3 S : -0.078789 0.355088 4 O : -0.405010 0.876042 5 Fe: 0.491707 1.157048 ...
Here you can see that the spin density on the Fe and O atoms adds up to roughly 2, corresponding to the two single occupied Fe-O π* orbitals, and there is some unpaired electron density on the sulphur atom, as expected.
2. Doublet state QM/MM energy calculation
The doublet calculation doublet/doublet.chm differs only in the QM flags set. As well as changing the spin multiplicity of the calculation, it also demonstrates how ORCA can be restarted from a wavefunction guess resulting from the quartet calculation, in order to aid convergence of the SCF procedure.
set qmflags { hamiltonian=dft functional=b3lyp \
basisspec = { { 6-31g { c n o h s } } { lanl2dz fe } } \
mult=2 scftype=uhf restart=yes moinp=../quartet/orca1.gbw \
maxcyc=200 }
In this case the Mulliken analysis in orca1.out should show the opposite spin density on the sulphur:
-------------------------------------------- MULLIKEN ATOMIC CHARGES AND SPIN POPULATIONS -------------------------------------------- 0 C : -0.649446 0.021974 1 H : 0.173016 -0.017432 2 H : 0.215289 -0.005462 3 S : -0.077582 -0.379189 4 O : -0.409642 0.844262 5 Fe: 0.492735 1.265113 ...