QM/MM Tutorial for Charmm / Turbomole

by Ville Kaila – Jan 2, 2014
see also hands-on Tutorial that I gave at CSC, Espoo, Finland in March 2018

Group web page: http://www.compbio.ch.tum.de
E-mail: ville.kaila .AT. ch.tum.de

NOTE!
I can unfortunately take no responsibility for the correctness of this tutorial. If you have suggestions or comments, please e-mail me. This tutorial is created for helping my group members getting started with setting up QM/MM calculations. However, I’m very happy if it’s useful for others interested in QM/MM simulations.

This document helps you prepare a QM/MM simulation using CHARMM and TURBOMOLE. See below for instructions for using Q-Chem/CHARMM.

You will need:
– CHARMM compiled for Turbomole. On AG Kaila cluster:
module load charmm/openmpi/c38b1_QT
– Turbomole
– turbo.py where all paths have been specified. turbo.py have been developed by Christopher Rowley’s team. For AG Kaila team members: our version is available on AG Kaila cluster e.g. at /home/qchem_scratch/turbo.py

0. Creating a structure

Relax a molecular model using classical MD and split your structure into pdb, which each contain a separate segments. You can do this with e.g. NAMD

The pdb formatting is:

ATOM 1 N LEU 1 -25.450 29.290 -19.372 1.00 22.23 PROA N
…
END
TER

CHARMM is sometimes sensitive with the END and TER stuff.

1. Create a CHARMM input file (mycharm.inp):

*QM/MM calculation of protein/s system: protein/s system, Lipids, water box
*

! CHARMM should begin with two “*”-lines
! an exclamation mark, indicates a comment
! stop tells charmm to exit

BOMBLEV -2
! error level that tells the program how severe errors will lead stop the execution of the program
! The lower the value, the more severe mistakes can pass
! For example BOMBLEV -5 allows severe mistakes to go through – ok for debugging, but not recommended for production run

! read in the topology and parameters from directory toppar36

set dir toppar36

set a @dir/top_all36_prot.rtf
set b @dir/par_all36_prot.prm
set c @dir/top_all36_lipid.rtf
set d @dir/par_all36_lipid.prm
set e @dir/toppar_water_ions.str

! read the protein topology
read rtf card name @a

! read the protein parameters
read para card flex name @b

! append the lipids topology
read rtf card append name @c

! append the nucleic acid parameters
read para card flex append name @d

! append water and ion top/par information
stream @e

! read in parameters for link atoms – see Appendix
OPEN UNIT 1 READ CARD NAME “qqh.rtf”
READ RTF CARD UNIT 1 append
CLOSE UNIT 1

OPEN UNIT 1 READ CARD NAME “qqh.prm”
READ PARA CARD UNIT 1 append flex
CLOSE UNIT 1

! read in structure of segment PROA
open read card unit 10 name “PROA.pdb”
read sequence pdb unit 10
generate PROA setup warn first NTER last CTER

open read card unit 10 name “PROA.pdb”
read coor pdb unit 10 resid
! change protonation state of Asp-93 to AspH
patch ASPP PROA 93 setup warn
autogenerate angles dihedrals

! read in some ions
open read card unit 10 name “sod.pdb”
read sequence pdb unit 10
generate IO1 setup warn first none last none
open read card unit 10 name “sod.pdb”
read coor pdb unit 10 resid

! read in some water molecules
open read card unit 10 name “water1.pdb”
read sequence pdb unit 10
generate WAT1 setup warn noangle nodihedral
open read card unit 10 name “water1.pdb”
read coor pdb unit 10 resid

! create missing hydrogens etc. with HBUILD – if necessary remove the “!”
!ic param
!ic build
!prnlev 0
!hbuild
!prnlev 5

! you can also read in final relaxed coordinates from previous runs
open read unit 22 card name “tmp.crd”
read coor unit 22 card

! call TURBOMOLE
! create directories turbo_input; data;turbo_output; add turbo.py to your path
envi qturboinpath “exact_path_to/data”
envi qturbooutpath “exact_path_to/turbodir”
envi qturboexe “exact_path_to/turbo.py”

! add link atoms between CA and CB of residues 160 and 161 of segname PROA
addlinkatom QQH1 PROA 160 CB PROA 160 CA
addlinkatom QQH2 PROA 161 CB PROA 161 CA

! define QM region
! syntax atom SEGNAME RESID ATOMNAME
! note “.or.” command and continuation of line with “-”

! syntax: atom + [SEGNAME] + [RESID] + [ATOMNAME]
define qm sele –
atom PROA 160 CB –
.or. atom PROA 160 HB –
.or. atom PROA 160 CG2 –
.or. atom PROA 160 HG21 –
.or. atom PROA 160 HG22 –
.or. atom PROA 160 HG23 –
! until  your qm region has been specified
show end

! specify the direction of the link atoms by charmm lonepair command
! order qm atom – mm atom
lonepair colinear scaled -.7261 –
sele type QQH1 show end sele atom PROA 160 CB show end sele atom PROA 160 CA
lonepair colinear scaled -.7261 –
sele type QQH2 show end sele atom PROA 161 CB show end sele atom PROA 161 CA

! scale the mass io QHH* atoms to 1.008 (H)
scalar mass show select type qq* show end
scalar mass set 1.008000 select type qq* show end
scalar mass show select type qq* show end

! For preparing a turbomole input you need to create the qm system and prepare an input for TURBOMOLE (see appendix)
!
print coord sele qm .or. atom * * QQH* end
stop

! after you have created a turbomole input, you should remove the “stop” command above and continue

! tell CHARMM how you want to calculate non-bonded interactions
! here with a switching function
NBOND –
atom cdiel eps 1.0 – !Electrostatics
vatom vswitch – ! VdW
elec switch – ! Elec
ctonnb 100.0 ctofnb 120.0 cutnb 160.0

faster off

! call turbomole; remove all classical interaction within qm region
qturbo remove sele qm end

! make a constrained minimization with

!open trajectory file for saving the minimization trajectory
open write unit 41 file name “optimized.dcd”

! select atoms for a reaction coordinate
SET ATOM1 PROA 160 O1
SET ATOM2 PROA 161 H1
SET ATOM3 PROA 161 O2

! create a variable V with a value of 1.30 Å
SET V 1.30

! enable all energy terms
SKIP NONE
! create a reaction coordinate “r0= d12-d23”
!that is influenced by an external potential of 0.5*k* (r-r0)*2
! here k = 100 kcal/mol/Å
RESDistance RESET KVAL 100.0 RVAL @v –
1.0 @atom1 @atom2 -1.0 @atom2 @atom3

! QM/MM minimizations will typically not converge if you haven’t fixed some parts of your full molecular system – here only the QM region can move; you can also create a 15 Å sphere around your qm region that can fully relax and fix or constrain the rest.
cons fix sele all .and. .not. qm end

! perform 500 steps of Adopted Basis Newton-Raphson !minimization or until tolerance in gradients and energies have reached 0.006 kcal/mol save into unit 41 (specified above), save every step

mini abnr nstep 500 nprint 1 tolg 0.0006 tole 0.0006 IUNCrd 41 NSAVC 1

! write out final coordinates
open write card unit 10 name “mini.pdb”
write coor pdb unit 10
* minimized
*

open write card unit 10 name “mini.crd”
write coor unit 10 card
* minimized
*

open write unit 10 card name “mini.psf”
write psf unit 10 card
* minimized
*

stop
! remove “stop” if you want to continue to MD

! MD part starts here
! open trajectory file, and other files for writing
open write unit 42 file name “./trajectory.dcd”
open write unit 31 card name “./production.dat”
open unit 20 write unform name “./velocities.trj”

! perform 1000 fs = 1ps of dynamics
set STEP 1000

SCALAR FBETA SET 8.0 select all end

dynamics lang leap start –
timestep 0.001 – ! timestep in picoseconds
nstep @STEP – ! number of steps and energy evaluations
nprint 1 – ! step frequency for writing in kunit and pr energy on unit 6
iunwri -1 – ! unit to write restart file
iunrea -1 – ! unit to read restart file
iuncrd 42 – ! unit to write coordinates (unformatted)
iunvel 20 – ! unit to write velocities out to
kunit 31 – ! unit to write energy and temp inforamtion out to file
iprfrq 20 – ! frequency for avg and rms energy
firstt 310. – ! initial T for velocity assignments
tstruc 310. – ! initial T for velocity assignments
finalt 310. – ! final T
tbath 310. – ! temperature of heat bath
iasors 1 – ! asgn (NOT scale) veldur heat/equil (.ne. 0 = asgn; .eq. 0 = scale)
iasvel 1 – ! use gaussian distrib of vel (.gt. 0 =gaussian ; .lt. 0 = uniform)
nsavc 1 – ! freq for writing coordinates
nsavv 1 – ! freq for writing velocities
isvfrq 100 – ! frequency for writing restart files
ihtfrq 0 – ! frequency for heating steps during the dynamics
ieqfrq 0 – ! frequency for scaling velocities during heating
inbfrq 1 – ! lists updated when necessary (heuristic test)
ihbfrq 0 – ! frequency for updating hydrogen bond list
ntrfrq 0 – ! step freq for stopping rotation and translation
echeck 100000000.0 – ! max variation of energy step-to-step
ichecw 0 ! check to see if avg. temp. lies w/temp range (0 = do not check)

! write out final coordinates
open write card unit 10 name “dynamics.pdb”
write coor pdb unit 10
* dynamics run
*

open write card unit 10 name “dynamics.crd”
write coor unit 10 card
* dynamics run
*

open write unit 10 card name “dynamics.psf”
write psf unit 10 card
* dynamics run
*

stop

Appendix

You need to tell CHARMM about link atoms

qqh.prm

* —-
* Built parameters for QQH link atoms
* by user AGKAILA
* —-
*

!******************************
! QQH link atoms
!******************************
ATOMS
MASS 101 QQH1 1.00800 !
MASS 102 QQH2 1.00800 !
MASS 103 QQH3 1.00800 !
MASS 104 QQH4 1.00800 !
MASS 105 QQH5 1.00800 !
MASS 106 QQH6 1.00800 !
MASS 107 QQH7 1.00800 !
MASS 108 QQH8 1.00800 !
MASS 109 QQH9 1.00800 !

BONDS
QQH1 CT1 0.000 1.0000 !
QQH1 CT2 0.000 1.0000 !
QQH1 CT2A 0.000 1.0000 !
QQH1 C 0.000 1.0000 !
QQH1 NH1 0.000 1.0000 !
ANGLES
QQH1 CT2 CT1 0.000 1.0000 !
QQH1 CT2 C 0.000 1.0000 !
QQH1 CT2 NH1 0.000 1.0000 !
QQH1 CT2A CT1 0.000 1.0000 !
QQH1 C CT1 0.000 1.0000 !
QQH1 NH1 CT1 0.000 1.0000 !
DIHEDRAL

NONBONDED nbxmod 5 atom cdiel shift vatom vdistance vswitch –
cutnb 14.0 ctofnb 12.0 ctonnb 10.0 eps 1.0 e14fac 1.0 wmin 1.5

!**************************************************************
! QQH
!**************************************************************
QQH1 0.000000 0.000000 0.000000 ! Link Atom
QQH2 0.000000 0.000000 0.000000 ! Link Atom
QQH3 0.000000 0.000000 0.000000 ! Link Atom
QQH4 0.000000 0.000000 0.000000 ! Link Atom
QQH5 0.000000 0.000000 0.000000 ! Link Atom
QQH6 0.000000 0.000000 0.000000 ! Link Atom
QQH7 0.000000 0.000000 0.000000 ! Link Atom
QQH8 0.000000 0.000000 0.000000 ! Link Atom
QQH9 0.000000 0.000000 0.000000 ! Link Atom
!
END

file qqh.rtf

* —- ——————————
* Built RTF for QMMM link atoms
* by user apgamiz
* ———————————–
*
22 0

!******************************
! QQH link atoms
!******************************
MASS 101 QQH1 1.00800 !
MASS 102 QQH2 1.00800 !
MASS 103 QQH3 1.00800 !
MASS 104 QQH4 1.00800 !
MASS 105 QQH5 1.00800 !
MASS 106 QQH6 1.00800 !
MASS 107 QQH7 1.00800 !
MASS 108 QQH8 1.00800 !
MASS 109 QQH9 1.00800 !

END

2. Prepare a turbomole run file with “define” in your data folder

for a proper general TURBOMOLE tutorial see: http://www.turbomole-gmbh.com/turbomole-manuals.html

coord.xyz contains the qm region

x2t coor.xyz > coord
run define

IF YOU WANT TO READ DEFAULT-DATA FROM ANOTHER control-TYPE FILE,
THEN ENTER ITS LOCATION/NAME OR OTHERWISE HIT >return<.
hit enter

INPUT TITLE OR
ENTER & TO REPEAT DEFINITION OF DEFAULT INPUT FILE

hit enter

IF YOU APPEND A QUESTION MARK TO ANY COMMAND AN EXPLANATION
OF THAT COMMAND MAY BE GIVEN

a coord
IF YOU APPEND A QUESTION MARK TO ANY COMMAND AN EXPLANATION
OF THAT COMMAND MAY BE GIVEN
*

IF YOU DO NOT WANT TO USE INTERNAL COORDINATES ENTER no
no
! internal coordinates will mess up the qm system positionitioning relative to mm point charges

ATOMIC ATTRIBUTE DEFINITION MENU ( #atoms=198 #bas=198 #ecp=2 )

b

ENTER A SET OF ATOMS TO WHICH YOU WANT TO ASSIGN BASIS SETS
( ATOMIC SET : all none )
TO OUTPUT ATOMIC SET SYNTAX ENTER A QUESTION MARK ?
E.G. all dz (DZ BASIS FOR ALL ATOMS)
all sto-3g hondo (SCALED STO-3G)
1,2,4-6 dzp (DZP BASIS FOR ATOMS 1,2,4,5,6)
“c” tz (TZ BASIS FOR ALL CARBON ATOMS)
ANY DISPLAY COMMAND dis MAY BE ENTERED OR YOU MAY
HIT >return< TO QUIT (GOING BACK TO MAIN MENU)

all def2-SVP

*

CHOOSE COMMAND
infsao : OUTPUT SAO INFORMATION
atb : Switch for writing MOs in ASCII or binary format
eht : PROVIDE MOS && OCCUPATION NUMBERS FROM EXTENDED HUECKEL GUESS
use : SUPPLY MO INFORMATION USING DATA FROM
man : MANUAL SPECIFICATION OF OCCUPATION NUMBERS
hcore : HAMILTON CORE GUESS FOR MOS
flip : FLIP SPIN OF A SELECTED ATOM
& : MOVE BACK TO THE ATOMIC ATTRIBUTES MENU
THE COMMANDS use OR eht OR * OR q(uit) TERMINATE THIS MENU !!!
FOR EXPLANATIONS APPEND A QUESTION MARK (?) TO ANY COMMAND

eht

ENTER THE MOLECULAR CHARGE (DEFAULT=0)
+2
! put in here the total charge of you qm region

DO YOU ACCEPT THIS OCCUPATION ? DEFAULT=y

hit enter

GENERAL MENU : SELECT YOUR TOPIC
scf : SELECT NON-DEFAULT SCF PARAMETER
mp2 : OPTIONS AND DATA GROUPS FOR rimp2 and mpgrad
cc : OPTIONS AND DATA GROUPS FOR ricc2
ex : EXCITED STATE AND RESPONSE OPTIONS
prop : SELECT TOOLS FOR SCF-ORBITAL ANALYSIS
…

dft
STATUS OF DFT_OPTIONS:
DFT is NOT used
functional b-p
gridsize m3

ENTER DFT-OPTION TO BE MODIFIED

func : TO CHANGE TYPE OF FUNCTIONAL
grid : TO CHANGE GRIDSIZE
on: TO SWITCH ON DFT
Just , q or ‘*’ terminate this menu.

on

func b-p

*

STATUS OF RI-OPTIONS:
RI IS NOT USED
Memory for RI: 500 Mb
Filename for auxbasis: auxbasis

ENTER RI-OPTION TO BE MODIFIED
m: CHANGE MEMORY FOR RI
f: CHANGE FILENAME
jbas: ASSIGN AUXILIARY RI-J BASIS SETS
on: TO SWITCH ON RI
Use , q, end, or ‘*’ to leave this menu

ri
m 1024
*

marij

dsp
on

in control file change:

$scfiterlimit 500
$thize 0.10000000E+99
$scfdamp start=1.700 step=0.050 min=0.050

remove: $scfdump

I recommend running a single point for the qm system to produce converged mos before starting the QM/MM runs. This will speed up the calculations.

3. In the control file add the following lines:

in $drvopt section add:
point charges

at the end of the control file add:
$point_charges file=point_charges
$point_charge_gradients file=pc_gradient

this will calculate a QM force on MM atoms, add MM charges around QM system, and calculate the gradients form MM on QM

4. Our modifications to turbo.py
original script can be found on https://github.com/RowleyGroup/charmm-turbomole-examples

# for parallel jobs, if mpi is true the MPI executables are used
# otherwise, the SMP executables are used
# On our cluster, MPI Turbomole binaries don’t work, so go SMP instead:
mpi=False

# see our run script below
nprocs=int(os.environ[‘SLURM_NPROCS’])
#nprocs = 1

# set the environment variable OMP_NUM_THREADS to nprocs for SMP jobs
#os.environ[‘OMP_NUM_THREADS’]=str(nprocs)

# the paths of the turbomole executables should be specified here
mpiexec=’/cm/shared/apps/mpiexec/0.84_427/bin/mpiexec’
turbomole=’/cm/shared/apps/TURBOMOLE66/bin/x86_64-unknown-linux-gnu’

ridft_serial=turbomole+’/ridft’
rdgrad_serial=turbomole+’/rdgrad’
dscf_serial=turbomole+’/dscf’
grad_serial=turbomole+’/grad’
ricc2_serial=turbomole+’/ricc2′
escf_serial=turbomole+’/escf’
egrad_serial=turbomole+’/egrad’

ridft_smp=turbomole+’_smp/ridft_smp’
rdgrad_smp=turbomole+’_smp/rdgrad_smp’
dscf_smp=turbomole+’_smp/dscf_smp’
grad_smp=turbomole+’_smp/grad_smp’
ricc2_smp=turbomole+’_smp/ricc2_omp’
escf_smp=turbomole+’_smp/escf_smp’
egrad_smp=turbomole+’_smp/egrad_smp’

ridft_mpi=turbomole+’_mpi/ridft_mpi’
rdgrad_mpi=turbomole+’_mpi/rdgrad_mpi’
dscf_mpi=turbomole+’_mpi/dscf_mpi’
grad_mpi=turbomole+’_mpi/grad_mpi’
ricc2_mpi=turbomole+’_mpi/ricc2_mpi’
escf_mpi=turbomole+’_mpi/escf’
egrad_mpi=turbomole+’_mpi/egrad’

Our run-script:

#!/bin/bash -l
#SBATCH -N 1 # Number of nodes, HAS TO BE 1
#SBATCH -n 4 # Number of CPUs, max 48
#SBATCH -J qmmm
#SBATCH -p medium
#SBATCH -e jobfile.err
#SBATCH -o jobfile.out
#SBATCH –mem-per-cpu=2000
#### start of jobscript

SDIR=`pwd`
ulimit -s unlimited

echo “SDIR is :” $SDIR
echo “SLURM_JOB_ID is:” $SLURM_JOB_ID
echo “Node used :” $SLURM_NODELIST
echo “Cores used:” $SLURM_NPROCS

export HOSTS_FILE=$SDIR/turbomole.machines
rm -f $HOSTS_FILE
srun hostname > $HOSTS_FILE

module load turbomole
module load charmm/openmpi/c38b1_QT

# Make a tempdir on /local disk:
export TURBOTMPDIR=/local/$USER/$SLURM_JOB_ID
mkdir -p $TURBOTMPDIR

# Define name of CHARMM input file:
CHARMMFILE=mycharm.inp

# Change TM path settings in the CHARMM input file automatically
# qturboexe -> SDIR/turbo.py
# qturbooutpath -> SDIR/TURBODIR
# qturboinpath -> SDIR/DATA
sed -i “s:^envi qturboexe.*:envi qturboexe \”${SDIR}/turbo.py\”:” $CHARMMFILE
sed -i “s:^envi qturbooutpath.*:envi qturbooutpath \”${SDIR}/turbodir\”:” $CHARMMFILE
sed -i “s:^envi qturboinpath.*:envi qturboinpath \”${SDIR}/data\”:” $CHARMMFILE

# Run CHARMM:
charmm < $CHARMMFILE > turbo-charmm.out

rm -rf $TURBOTMPDIR

5. Running with Q-Chem and CHARMM

Instead of turbomole envi commands, replace with

envi qchemcnt “qchem.inp”
envi qcheminp “q1.inp”
envi qchemexe “qchem”
envi qchemout “q1.out”

B. instead of “qturbo remove sele qm end” use for 16 CPUs:
qchem remove exgroup PARAllel 16 noguess QCLJ sele qm end

qchem.inp can contain for example, the following lines:
$rem
METHOD BP86
basis VDZ**
qm_mm true
jobtype force
DFT_D EMPIRICAL_GRIMME3
symmetry off
sym_ignore true
print_input false
qmmm_print true
THRESH 11
SCF_CONVERGENCE 6
SCF_ALGORITHM DIIS
MAX_SCF_CYCLES 1000
$end

$molecule
0 1
$end

CHARMM will automatically create all files, but your run script probably needs additionally some commands that call Q-Chem. For example:

#!/bin/bash -l
#SBATCH -N 1 # Number of nodes
#SBATCH -n 48 #
#SBATCH -J QMMM
#SBATCH –ntasks-per-node=48 # Tasks per node
#SBATCH -p medium
#SBATCH -e job.err%J #
#SBATCH -o job.out%J
#SBATCH –mem-per-cpu=1000
#### start of jobscript
SDIR=`pwd`

export QCMACHINEFILE=$SDIR/machines
rm -f $QCMACHINEFILE
srun hostname > $SDIR/machines

module load gcc/4.7.2
module load openmpi/gcc/64/1.8.1
module load charmm/openmpi/c38b1_QC
module list

export QC=/cm/shared/apps/QCHEM/MPICH/4.2/
source $QC/qcenv.sh

ulimit -s unlimited

export QCAUX=$QC/qcaux
export QCSCRATCH=/home/qchem_scratch/$USER/$SLURM_JOB_ID.1
export QCLOCALSCR=/local/$USER/$SLURM_JOB_ID.2
mkdir -p $QCLOCALSCR
mkdir -p $QCSCRATCH
# parallel setting
export QCRSH=ssh
export QCMPI=mpich
export QCMPIRUN=/cm/shared/apps/QCHEM/MPICH/4.2/bin/mpi/mpirun_qchem
export QCMACHINEFILE=$SDIR/machines

module list
module load mpich/ge/gcc/64/3.1.2

echo ‘mpicc version is:’
/usr/bin/which mpicc
echo ‘mpiexec version is:’
/usr/bin/which mpiexec

charmm < mycharm.inp > mycharm.out