Combined Quantum Mechanical and Molecular Mechanics Method
Based on SCCDFTB in CHARMM
by Qiang Cui and Marcus Elstner
(cui@chem.wisc.edu, elstner@phys.upb.de)
The approximate Density Functional program SCCDFTB (Self-
consistent charge Density-Functional Tight-Binding) is interfaced with
CHARMM program in a QM/MM method.
This method is described in
Phys. Rev. B, 58,7260 (1998)
Phys. Stat. Sol. B, 217, 357 (2000)
J. Phys. : Condens. Matter., 14, 3015 (2002)
Recent reviews are:
WIREs Comput. Mol. Sci., 4, 49-61 (2014)
PCCP, 16, 14368-14377 (2014)
The QM/MM interface in CHARMM has been described in
J. Phys. Chem. B 105 (2001) 569
The GHO-SCC-DFTB/MM boundary treatment has been described
in J. Phys. Chem. A 108 (2004) 5454.
A recent review of SCC-DFTB/MM can be found in
J. Phys. Chem. B. 110, 6458-6469 (2006). Recent extensions can be found in
J. Phys. Chem. A. 111, 10861-10873 (2007).
J. Chem. Theory Comput. 7, 931-948 (2011).
The extension of the SCC-DFTB method to work with the Replica Path
and the Nudged Elastic Band methods has been described in the following
paper and should be cited when applied:
H. L. Woodcock, M. Hodoscek, and B. R. Brooks Exploring SCC-DFTB Paths
for Mapping QM/MM Reaction Mechanisms J. Phys. Chem. A; 2007; 111(26)
5720-5728.
* Menu:
* Description:: Description of the sccdftb commands.
* Usage:: How to run sccdftb in CHARMM.
* Installation:: How to install sccdftb in CHARMM environment.
* FEP:: Free energy perturbations with SCC-DFTB/MM
* Electrostatics:: Electrostatics in SCC-DFTB/MM
* GLNK:: GLNK Commands
* Status:: Status of the interface code.
The SCCDFTB QM potential is initialized with the SCCDFTB command
[SYNTAX SCCDFTB]
SCCDFTB [REMOve] [CHRG] (atom selection) [GLNK atom-selection]
[TEMPerature] [SCFtolerance]
[CUTF] [EWAD EOPT KAPPA kappa KMAX kmax KSQMAX ksqmax]
[PMEW FFTX fftx FFTY ffty FFTZ fftz ORDEr order]
[UPDT 0]
[MULL]
[DISP]
[HBON GAUS]
[DHUB DHGA]
[D3RD]
[SCFC]
[QDIP UDIP unit]
[SHES] [RHES]
[CDKO]
[UNPE unpe]
[LDEP]
[DMEThod int]
[MIXE int] [GENS int] [SYMM real]
[IGUESS]
[IDAMP]
[PLUS]
REMOve: Classical energies within QM atoms are removed.
CHRG: Net charge in the QM subsystem.
The atoms in selection will be treated as QM atoms.
GLNK atom-selection: contains a list of atoms that are selected as
GHO boundary atoms.
TEMPerature: Specifies the electronic temperature (Fermi distribution).
Can be used to accelerate or achieve SCF convergence
(default =0.0).
SCFtolerance: Convergence criteria for the SCF cycle. As default
a value of 1.d-7 is used.
CUTF: a flag that turns on cut-off, will use the same scheme used for
MM interactions (allows atom-based fshift and shift). This option
is recommended for QM-MM interactions in EWAD calculation.
EWAD: a flag that turns on ewald summation for QM-QM and QM-MM
electrostatic interactions
EOPT: performs an internal optimization for kappa and kmax as well as
real-space sum. NOT recommended - very inefficient. incompatible
with CUTF
KAPPA, KMAX, KSQMAX: parameters used for QM-QM and QM-MM ewald
contributions.
PMEW: a flag that turns on particle mesh ewald (PME) summation for
QM-QM and QM-MM electrostatic interactions
FFTX, FFTY, FFTZ: Integer numbers of Fast Fourier Transform grid points
for the charge mesh (default 32)
ORDEr: specifies the order of the B-spline interpolation, e.g. cubic is
order 4, fifth degree is ORDEr 6 (default). The ORDEr must be
an even number and at least 4.
UPDT: Whether to update the box during a MD run. Default is 0 (not update!)
Do NOT forget to specify this to 1 when running CPT calculations
(although additional tests should be done for CPT in terms of virial
calculations!)
MULL: Transfer the mulliken charges to the CG array such that they can be
printed out by transferring the CG array to any vectors (e.g.,
scalar wmain = charge; print coor)
DISP: Dispersion interactions among QM atoms can be calculated using an
empirical formular (Elstner et al. J. Chem. Phys. 114, 5149, 2001).
One needs to specify a set of parameters in the DISPERSION.INP file
(see the above ref for details).
TWOBody THREebody: Add the D3(BJ) two-body or three-body dispersion correction
due to Grimme (Grimme et al. J. Chem. Phys. 132 (2010), 154104 and Grimme
et al. J. Comput. Chem. 32 (2011), 1456-1465). Parameters are set
automatically with the DFTB method (I.e. DFTB2/DFTB3/CPE correction, etc).
HBON: The short-range behavior of XH gamma function is modified with a damping
function. This significantly enhances the hydrogen bonding interactions
(Elstner, Cui, unpublished). The damping exponent needs to be specified
in the sccdftb.dat file. However, with the modified gamma, the repulsive
potential has to be adjusted accordingly. This can be done in an
empirical fashion by including a Gaussian (constrained to operate in
a range by a switching function) in the relevant repulsive potential;
the parameters for the Gaussian and the switching function need to be
specified in the spl file and turned on using the GAUS keyword.
DHUB: --> obsolete keyword (but still fully functional): use D3RD instead.
To improve proton affinities, which depend much on the charges in the
protonated DHGA and deprotonated molecules, the SCC-DFTB is expanded
to the 3rd order. Currently only on-site terms have been included,
which were observed to have a major impact on the calculated PAs for
many molecules. The relevant parameters are the derivative of the
Hubbard parameters (related to chemical hardness), which need to be
specified in the sccdftb.dat file. The DHGA keyword includes further
flexibility in the behavior of the Hubbard parameters as a function of
charge (Elstner, Cui, JPC-A, 111, 10861-10873 (2007)).
D3RD: Full 3rd order extension of SCC-DFTB. For invoking DFTB3
(JCTC 2011,7,931) D3RD should be used in combination with HBON.
The necessary parameters (Hubbard derivatives) need to be
specified in the sccdftb.dat file (see below).
SCFC: instead of using energy difference between iterations to determine SCF
convergence, use difference in Mulliken charge between iterations. The
threshold is 0.05*sqrt(SCFtolerance). This option can be useful when
very tight convergence is desired (see, e.g., SI of JCP, 127, 234504
(2007))
CPE0, CPEQ: Use the charge-independent DFTB3/CPE model with the parameters
from Christensen et al. (2015) http://arxiv.org/abs/1507.00370 - Note:
The parameters are optimized for simultaneous use of the D3(BJ)
three-body dispersion correction (i.e., the THREebody keyword).
Note: the CPE methods requires CHARMM to be compiled with the +DFTBMKL
option - this is described under the MIXE keyword in this document.
QDIP: indicate that QM dipole moment will be written to file of unit # UDIP.
Can be useful for IR calculations using ACF of QM dipoles.
SHES: Save "dd1.dat" file in finite differences frequency calculations
(DIAG FINITE and REDU FIX FINITE)
RHES: Read in "dd1.dat" file (from parent directory of the calculation)
CDKO: The Charge-Dependent Klopman-Ohno flag activates a short-range
damping of the 1/r electrostatic interaction between QM (SCC-DFTB)
and MM atoms taking into account charge-penetration effects. It
is particularly useful for describing highly charged species in
a QM/MM setup such as phosphate hydrolysis reactions in
biological systems and was investigated in combination with
an optimized set of element type dependent QM Van der Waals
parameters (Hou, Cui, JCTC, 8, 4293-4304 (2012)). The relevant parameters
need to be specified in a separate input file named ko_para.inp
(see below for details) that needs to be located in the directory
where CHARMM is executed (similar to sccdftb.dat). Note, when
using in combination with LDEP the s-orbital Hubbard (and deriv)
is taken as the QM Hubbard for the QM/MM interaction.
UNPE: Followed by a floating-point number this flag specifies the
number of unpaired electrons and invokes the inclusion of
spin-polarization effects in a collinear description (see e.g.
JPC-A 2007,111,5622. and refs therein).
LDEP: L-dependence of the Hubbard parameters. This option should only be
used with DFTB parameters that are designed for l-dependence! If
D3RD is switched on also the Hubbard derivative parameters for s-,
p-, and d-orbitals are necessary and need to be specified in the
sccdftb.dat (see below). An example for this option is provided
in test/c38test/sccdftb_ldep.inp.
DMET: Needs an integer(0/1/2/3) as input to specify the diagonalizer for solving
the secular equations in DFTB.
0 : Legacy diagonalizer
1 : LAPACK DSYGV, requires +DFTBMKL
2 : LAPACK DSYGVD, requires +DFTBMKL - (recommended)
Options 1 and 2 requires CHARMM to be linked to an implementation of LAPACK,
and to be compiled with the +DFTBMKL option. If +DFTBMKL is used, BLAS routines
are used to calculate the density matrix and calculate the Mulliken charges and
DFTB gradient. DMET 2 is usually the fastest, and will normally speed up the
DFTB calculation substantially.
MIXE: Needs an integer(0/1/2/3) as input to specify the charge-mixing scheme
to be used in the SCF iterations.
0 : Simple mixing scheme
1 : Anderson mixing scheme (recommended)
2 : Broyden mixing scheme (default)
3 : DIIS mixing scheme
For the Anderson and DIIS mixing schemes, the LAPACK subroutine DGESV
is required. Hence the CHARMM executable has to be linked to LAPACK
(Intel MKL is recommended). For the Anderson and DIIS implementation to
work, the following need to be done:
(a) Modify the last line of Makefile_(machine type) for linking to LAPACK.
(b) Add +DFTBMKL to the installation command line such that DFTBMKL gets
included in the generated pref.dat. An example can be:
./install.com gnu large T IFORT +DFTBMKL
If MIXE is specified as '1' or '3' in the input without having compiled CHARMM
with LAPACK, the Broyden mixer will be used instead of Anderson/DIIS.
As an example, when compiling with the ifort compiler, the following is
one possible option to be APPENDED to the last line of build/UNX/Makefile_em64t
(where ${MKLROOT} has been pre-defined appropriately):
-Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_lp64.a \
${MKLROOT}/lib/intel64/libmkl_core.a ${MKLROOT}/lib/intel64/libmkl_sequential.a \
-Wl,--end-group -lpthread -lm -ldl
(the above is ONE line only.)
In situations where Intel MKL is unavailable, simply compile with the OS
default version of BLAS and LAPACK by linking CHARMM with:
-llapack -lblas
Note that Intel MKL is usually the fastest LAPACK implementation for CHARMM.
GENS: Number of iterations (including the current one) whose input and output
charges are mixed in the Anderson mixing scheme. Generally, 5-6 should be
sufficient. Default: 8
SYMM: A real number, in order to avoid linear dependency near the convergence
region in the Anderson mixing scheme. Default: 0.01 (recommended)
IGUESS: Order of DXL-BOMD algorithm for charge propagation from previous
dynamics. The input value and real order correspondence is [ 105: 1,
107: 3, 109: 5, 111: 7, 113: 9, 115: 11, 117: 13].
The details of the algorithm can be found in J. Chem. Phys. (135,
044122, 2011). The default value is 0, which means getting charge from
the previous geometry converged values.
IDAMP: Number of damping steps in DXL-BOMD. The number of steps before turning
on DXL-BOMD algorithm.
PLUS: Use the libDFTB+ back-end to carry out the SCC-DFTB calculations. This
option requires CHARMM to be compiled with the +DFTBPLUS option.
Instructions for this can be found in the installation section in this
document.
In the SCCDFTB program the atomtypes are represented by consecutive
numbers. The definition of SCCDFTB atom numbers has to be accomplished
before invoking the SCCDFTB command. The numbers are stored in WMAIN.
If the QM system e.g contains only O, N, C and H atoms,
the the numbering can be executed as follows:
scalar WMAIN set 1.0 sele type O* SHOW end
scalar WMAIN set 2.0 sele type N* SHOW end
scalar WMAIN set 3.0 sele type C* SHOW end
scalar WMAIN set 4.0 sele type H* SHOW end
Now, the O atoms are represented by 1.0, the N atoms by 2.0 etc.
Link atom may be added between an QM and MM atoms with the
following command:
ADDLinkatom link-atom-name QM-atom-spec MM-atom-spec
link-atom-name ::= a four character descriptor starting with QQ.
atom-spec::= {residue-number atom-name}
{ segid resid atom-name }
{ BYNUm atom-number }
When using link atoms to break a bond between QM and MM
regions bond and angle parameters have to be added to parameter file
or better use READ PARAm APPEnd command.
If define is used for selection of QM region put it after all
ADDLink commands so the numbers of atoms in the selections are not
changed. Link atoms are always selected as QM atoms.
Currently, three different link atom schemes are implemented.
SLA(default), EXGR and DIV. For a detailed comparison between them, see
J. Phys. Chem. B 109, 9082-9095 (2005). Briefly, SLA should be avoided
if there is major charge change during the QM/MM calculations (e.g.,
deprotonation). EXGR in general works well, but can be problematic if the
QM region interacts directly with the mainchain NH (in the same group as
CA and therefore can be excluded from interacting with QM).
Example:
SCCDFTB DIV remove CHRG 2 SELE resn @m END TEMP 0.00 SCFT 0.00000001
SCCDFTB input files
-------------------
SCCDFTB needs to read in the parameter files, which have
a two-body character. Therefore, the interaction parmeters
for all pairs of atoms have to be read in.
These files are named like oo.spl, on.spl, oc.spl, no.spl etc.,
where oo.spl contains the two-center integrals for the O-O interaction,
on.spl the two-center integrals for the O-N interaction etc.
DFTB needs these parameters for the O-N and N-O interaction,
similarily for all other pairwise interactions.
The file sccdftb.dat contains the paths to these parameters, as:
'potential:atom-1-atom-1'
'potential:atom-1-atom-2'
'potential:atom-1-atom-3'
... \\
'potential:atom-1-atom-N'
'potential:atom-2-atom-1'
... \\
'potential:atom-2-atom-N'
'potential:atom-N-atom-1'
... \\
'potential:atom-N-atom-N'
where atom-1 is the atom defined by 1.0, as described above,
atom-2 defined in WMAIN by 2.0 etc.
If open-shell treatment via UNPE is specified, SCC-DFTB needs to read
in the spin-polarization constants. These are listed directly after the
parameter files section, one line for each atom type as defined in
WMAIN and following the same order as above (first line for atom-1,
second line for atom-2, and so on). Each line needs to look like:
spin-'dummy' Wss Wsp Wps Wpp Wsd Wpd Wdd Wds Wdp
where for readability 'dummy' can be the name of the atom type (e.g. c
for carbon). Values for the constants Wss,... can be found in the
Dissertation thesis of Christof Koehler, 2003, Unviversity Paderborn,
page 114f, which is also available online via
http://digital.ub.uni-paderborn.de/hs/content/titleinfo/3163
When third order DFTB is invoked via D3RD (or the now obsolete keyword
DHUB) the Hubbard derivatives Ud need to be specified following the
parameter file section, or the spin-section if the UNPE keyword is used:
'dummy' Ud-atom1
'dummy' Ud-atom2
...
For values of Ud see JCTC 2011,7,931.
When D3RD and LDEP is specified the Hubbard derivatives section looks
like:
'dummy' Ud-atom1-d-orbital Ud-atom1-p-orbital Ud-atom1-s-orbital
'dummy' Ud-atom2-d-orbital Ud-atom2-p-orbital Ud-atom2-s-orbital
...
If DHGA is specified (again, this is an obsolete option and was only
used in combination with DHUB), one line containing three number follow
in the form:
v0_hgau alp_hgau q0_hgau
For details see JCTC 2008,4,2067.
For HBON one additional parameter 'zeta' is necessary which is set as a
single number in a single line following all previous sections (if
invoked via corresponding keyword).
An example of a system containing O N C and H and having specified the
keywords UNPE 1.0, D3RD, and HBON, sccdftb.dat would look like:
'PATH/oo.spl'
'PATH/on.spl'
'PATH/oc.spl'
'PATH/oh.spl'
'PATH/no.spl'
'PATH/nn.spl'
'PATH/nc.spl'
'PATH/nh.spl'
'PATH/co.spl'
'PATH/cn.spl'
'PATH/cc.spl'
'PATH/ch.spl'
'PATH/ho.spl'
'PATH/hn.spl'
'PATH/hc.spl'
'PATH/hh.spl'
spin-o -0.035 -0.030 -0.030 -0.028 0.0 0.0 0.0 0.0 0.0
spin-n -0.033 -0.027 -0.027 -0.026 0.0 0.0 0.0 0.0 0.0
spin-c -0.031 -0.025 -0.025 -0.023 0.0 0.0 0.0 0.0 0.0
spin-h -0.072 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
'o' -0.19
'n' -0.13
'c' -0.23
'h' -0.16
4.2
where PATH specifies the path to the directory where the data files
are located. Be careful, an error in the sequence or a wrong assingnment
of parameters to atoms (coordinates) will make results meaningless.
Parameter files can be downloaded at www.dftb.org or
requested from Marcus Elstner (marcus.elstner@kit.edu).
Another input file is necessary if the CDKO keyword is applied.
It contains two parameters for each QM atom type and one atom-size
describing parameter for each element appearing in the MM part. The
file is called ko_para.inp and needs to look like:
nQMatomtypes nMMelements
WMAIN:atom-1 a:atom-1 b:atom-1
WMAIN:atom-2 a:atom-2 b:atom-2
... ...
WMAIN:nQMatomtypes a:atom-nQMatomtypes b:atom-nQMatomtypes
<blank line>
name:MMelement1 U:MMelement1
name:MMelement2 U:MMelement2
... ...
name:nMMelements U:nMMelements
For the sccdftb.dat example above, the number of QM atomtypes is
nQMatomtypes = 4 (O,N,C, and H), if the MM part is composed solely out
of water the number of MM elements nMMelements = 2, where
name:MMelement1 is O and name:MMelement2 is H. a, b, and U are the
corresponding parameters published for phosphate
hydrolysis reactions (Hou, Cui, 8, 4293-4304 (2012)).
An example for this scheme is provided in test/c37test/sccdftb_cdko.inp
For additional discussion on the CDKO parameters, see, JPC, B 118, 11007-11027 (2014).
SCCDFTB output files (currently disabled)
-----------------------------------------
SPE.DAT : contains the Kohn-Sham energies with occupations numbers.
CHR.DAT : contains the atomic (Mulliken) charges of the atoms
(first row) and for the orbitals (s, px,py,pz,dxx.. ) in the
following columns.
REST.DAT: contains dipolemoment (D), calculated from the Mulliken
charges (not a reliable estimate of Dipolemoment in general!)
Installation of SCCDFTB
-----------------------
The source code of SCCDFTB ist distributed with CHARMM.
To compile the SCCDFTB method as the quantum part:
./install machine size T
T invokes the SCCDFTB
The parameter files have to be reqeusted and stored in a directory,
which can be reached by 'PATH' (see up).
To enable the DFTB+ interface, three steps must be taken:
1) Obtain and compile the DFTB+ library - instructions
must be found elsewhere.
2) Modify the Makefile_${em64t|gnu} for proper linking, etc.:
2a) Add the directory where libdftb+ is compiled to the $INCLUDE
environment variable in the makefile:
LIBDFTBDIR = /path/to/dftbplus_source/api/charmm/_obj_x86_64-linux-ifort
INCLUDE = -I$(LIBDFTBDIR)
2b) Modify the link-line to properly link to libdftb+, e.g. for Makefile_em64t
$(LD) -o charmm.exe $(LIB)/*.o $(LIBDFTBDIR)/_dir_extlib_xmlf90/*.o $(LIBDFTBDIR)/*.o \
$(LIBS) $(LIBS) $(LIBS) $(ADDLIB) $(QLIB) $(MSGLIB) -I$(LIBDFTBDIR)/ \
-I$(LIBDFTBDIR)/_dir_extlib_xmlf90 $(MATHLIBS)
where $(MATHLIBS) is the link options to BLAS/LAPACK (usually MKL).
3) use the "+DFTBPLUS keyword" to compile, e.g.:
./install machine size T +DFTBPLUS
4) Use the PLUS keyword in the SCCDFTB group to run using the DFTB+ engine.
Free energy perturbations with SCC-DFTB/MM
The code currently allows dual-topology based SCC-DFTB/MM free energy
perturbation calculations; since all scaling related to the QM component
of the free energy derivative is done inside SCC-DFTB, the FEP
calculations do not have to use BLOCK.
As discussed in JPC, 107, 8643 (2003), a practical problem of using
FEP with QM/MM potentials is that the structure of the QM region
undergoes significant distortions at end-points if one scales the entire
QM molecule; there is no such problem if one chooses to scale only
QM/MM interactions, but that requires calculation of new terms. The
general solution is to add harmonic constraints on the QM part -
either only at the end-points and then re-weight the calculated free
energy derivatives - or, more elegantly, add harmonic constraints
as "chaperones" throughout the "alchemy" simulation and compute
corrections based on local configuration integrals. See W. Yang et al.
J. Chem. Phys. 2004.
For the special case where the two end-states have very similar chemical
structures - such as in redox, metal-exchange and pKa applications,
which we believe are scenarios where QM/MM treatment is useful, a simple
dual-topology-single-coordinate (DTSC) approach has been introduced. As
the name implies, one uses only one set of coordinates for the two
states (e.g., reduced and oxidized states). Due to the fact that the
free energy is path-independent, such an approach is formally exact. In
practical applications, error might arise due to SHAKE - i.e. X-H
distances are assumed to the same in the two states - which usually has
negligible effects.
At each configuration (hence single-coordinate) along the trajectory,
two electronic structure calculations are carried out (dual topology)
and the free energy derivative with respect to the coupling parameter
is evaluated and averaged on the fly.
With minor modifications, the algorithm also works for pKa prediction
for a specific group in large molecules. For more details, refer to the
following publications:
M. Formaneck, G. Li, X. Zhang, Q. Cui, J. Thero. Comput. Chem.
1, 53-68 (2002)
G. Li, X. Zhang, Q. Cui, J. Phys. Chem. B (2003) 107, 8643
G. Li, Q. Cui, J. Phys. Chem. B, (2003) 107, 14521
NOTE BENE: It MUST be used with "FAST OFF" because only generic
atom-atom codes have been modified so far (made default).
Due to the fact that ALL QM related components are handled within SCC
(including GSBP and eWald, see next section), FEP (such as pKa)
calculations can be used with both eWald and GSBP - provided that it is
the QM part that undergoes "alchemical" mutation.
The code has NOT been extensively tested in which both QM and MM
undergo changes.
Two examples are given to illustrate computational details; the
first one deals with redox potential calculations for FAD in Cholesterol
oxidase, and the second one concerns pKa calculations of ethanethiol
(Ch3CH2SH) in water. The test files are scc_fep_dtsc.inp and
scc_pka1/2.inp.
Example(1)
----------
For redox potential calculations, the following set-up is used,
......
SCCDFTB LAMDa [REST] STOP [CUTF] OUTPut int -
[REMOve] (atom selection 1) [CHRG] [TEMPerature] [SCFtolerance] -
[UNPE] -
INIT @lam PASS int STEP int TIAV int -
[CHRG] [TEMPerature] [SCFtolerance] -
[UNPE] -
(atom selection 2) -
(atom selection 3)
LAMD: invoke the TI method to perform free energy calculations
REST: Restart option for accumulating statistics concerning <DU/DL>
i.e., necessary values will be read in from dynamics restart
file
STOP: employ the dual-topology-single-coordinate approach
CUTF: invoke cutoff for QM/MM electrostatic interactions
OUTP: unit number for storing the free energy derivative <DU/DL>
INIT: the current lamda value
PASS: numbers of MD steps to be skipped when accumulating <DU/DL>
STEP: the frequency of collecting statistics for <DU/DL>
TIAV: the frequency of computing the average of DU/DL.
atom selection 1: Reactant+Product to set up MM list for QM atoms
atom selection 2: Reactant state
atom selection 3: Product state
--------------------------------------------------------------------
For pKa calculations, two free energy simulations are in
principle required; in the first step, the protonated state is mutated
into the ionized state as the acidic proton is mutated into a dummy atom
in the second step, the dummy atom is transferred into the gas phase.
Test calculations indicate that the contribution from the 2nd step is
likely to be small.
Example(2a)
----------
In the first step, BLOCK is used together with SCCDFTB
......
BLOCK 3
SCCDFTB STOP PKAC ISTP 1
CALL 2 SELE qm1 END
CALL 3 SELE qmh END
CALL 1 SELE .not. (qm1 .or. qmh) END
COEF 1 1 1.0
COEF 1 2 1.0
COEF 1 3 1.0
COEF 2 2 0.0
COEF 2 3 @lam
COEF 3 3 0.0
END
SCCDFTB PKAC ISTP 1 HYGN int [CUTF] OUTPut int -
[REMOve] (atom selection 1) [CHRG] [TEMPerature] [SCFtolerance] -
[UNPE] -
INIT @lam PASS int STEP int TIAV int -
[CHRG] [TEMPerature] [SCFtolerance] -
[UNPE] -
atom selection 2 -
atom selection 3
......
In the BLOCK section :
the SCCD keyword is used to set up coefficent matrix for calculating
bonded contribution involving the dummy atom to <DU/DL>.
ISTP 1: the first step in pKa calculations
qm1 is the ionized state (e.g., CH3CH2S-);
qmh is the acidic proton.
In the SCCDFTB section:
PKAC : invoke pKa calculation
ISTP 1: the first step in pKa calculations
HYGN : atomic index (number) of the acidic proton in the psf
atom selection 1: protonated state (CHRG: protonated state)
atom selection 2: protonated state (CHRG: deprotonated state)
atom selection 3: deprotonated state
Example(2b)
----------
In the second step for pKa calculations, the dummy atom is transferred
into vacuum,
......
calc 1mlam 1.0-@lam
BLOCK 3
SCCDFTB STOP PKAC ISTP 2
CALL 2 SELE qm1 END
CALL 3 SELE qmh END
CALL 1 SELE .not. (qm1 .or. qmh) END
COEF 1 1 1.0
COEF 1 2 1.0
COEF 1 3 @1mlam
COEF 2 2 0.0
COEF 2 3 0.0 bond 1.0 angl 1.0 dihe 1.0
COEF 3 3 0.0
END
SCCDFTB PKAC ISTP 2 HYGN int [CUTF] OUTPut int -
[REMOve] [CHRG] (atom selection 1) [TEMPerature] [SCFtolerance] -
INIT @lam PASS int STEP int TIAV int
......
Note that with BLOCK, the coefficient matrix is different in the
second step: we are only scaling the non-bond (vdW) interaction between
the environment and the dummy atom (1 and 3). The bonded terms between
the QM and the dummy atom (2 and 3) is kept (coefficient as 1.0) and
will be taken out with local configuration integrals.
In SCC-DFTB, atom selection 1: deprotonated state
Since 2004, electrostatics in SCC-DFTB/MM simulations can be
treated in several ways for both spherical and periodic conditions:
i). As for other QM packages, the default is no cut-off for QM/MM
electrostatic interactions. This is NOT recommended when cut-off
is used for MM; the imbalance will cause over-polarization of the media
(e.g.,see discussion in classical simulations by Woods, J. Chem. Phys.
103, 6177, 1995). A useful option is to use extended electrostatics for
MM.
ii). Cut-off is introduced for QM/MM electrostatics, similar to MM
interactions; i.e., the same scaling factors are the same as those for
MM interactions. Simply add "CUTF" to the SCC-DFTB command line.
Currently only supports energy/force-shifts based on atoms
iii). For spherical boundary conditions, the
GSBP approach can now be used with SCC-DFTB. The current implementation
takes GSBP contributions into the SCF iteration, although for a large
inner region, this may not be necessary. Further tests are being carried
out. The code will be extended to other boundary conditions and QM
methods in the future. If one uses sorting (i.e., truncate size of
basis in GSBP), make sure a SCC-DFTB/MM
energy calculation with MULL (save Mulliken charge) is carried out
before issuing GSBP, since the Mulliken charges are used to estimate
contributions from various basis functions to the QM related terms.
See test cases for examples.
The GSBP and PB reference calculations are made consistent in terms of
reference state and the boundaries are updated in the PB calculations.
Simple test on a simple sodium ion led to correct answer with different EPSP
values.
The implementation is described in: J. Chem. Phys. 123, 014905 (2005).
iv). For PBC simulations, one can use either cut-off or eWald sum for
SCC-DFTB/MM interactions. Also PME has been implemented for the QM/MM
interactions but no extensive performance checks have been carried out
yet. The current QM/MM implementation allows in principle all cell shapes.
For eWald, one can either let the code optimize the exponent to get the
best balance between real space sum and the reciprocal space sum (EOPT)
or one can specify a set of parameters (Kappa, KMAX, KSQMAX).
The real space sum is done till convergence is met with EOPT or
without CUTF (so more expensive); EOPT is incompatible with CUTF.
With cutoff (CUTF), the real sum is limited to atoms within the cutoff-
which is recommended (much more efficient). In any case, one should
carefully test kappa, KMAX to ensure the convergence of energy and,
more importantly, force from SCC-DFTB/MM calculations.
A sample command line would be:
......
SCCDFTB remove CHRG 2 SELE resn @m END TEMP 0.00 SCFT 0.00000001 EWAD -
CUTF Kappa 0.45 KMAX 6 KSQMAX 100
......
With the 2009 implementation, QM/MM-eWald without CUTF is about 4-6
times slower than a QM/MM calculation without eWald (i.e., with only cutoff).
With real-space cutoff (CUTF), QM/MM-eWald is only slightly slower than QM/MM
without eWald. In 2012, QM/MM-PME has been implemented for SCC-DFTB by the
Gao group. For an example, see test/c37test/sccdftb_pme.inp.
The eWald implementation is described in: J. Phys. Chem. B 109, 17715-17733
(2005). The PME implementation follows the same approach as established for
SQUANTM as described in J. Chem. Theory Comput, 1, 2-13 (2005).
An important point for PBC simulations is that all image must be used
with "UPDAte IMAL ". This is because symmetry operations have not been
considered in the SCC-DFTB/MM code - which obviously needs to be fixed
in the future.
The eWald implementation is described in: J. Phys. Chem. B 109, 17715-17733
(2005)
v). SCC-DFTB has been implemented to work with Poisson-Boltzmann, which
allows to calculate solvation free energy and optimize reaction path with
implicit solvent. Several keywords are added to the PBEQ module:
PSTL(0.01 kcal/mol): Energy convergence criterion of the iterative
QM-PB calculation
MXPS(5000): Maximum # of SCC and PB solver iterations.
IGAS: flag to initiate calculation from gas phase in every step,
otherwise from last step (this is default for solvation free
energy calculation).
CHDR: instead of a fixed set of atomic radii, use charge-depdent radii,
which can be useful for charged species (Hou, Zhu and QC,
unpublished). In this case, a file named radius.inp, which
contains optimized parameters for the charge-dependence of atomic
radius, must be included in the local directory, similar to
sccdftb.dat.
For details please refer test cases: test_sccpb.inp
Description of the GLNK Command
[GLNK atom-selection]
atom-selection: contains a list of atoms that are boundary atoms.
Restrictions: see the correponding entry for GLNK in quantum.doc
Description: see the correponding entry for GLNK in quantum.doc
Limitations: The present implementation allows up to 5 QM-boundary
atoms. To improve the geometry for the QM/MM boundary bond, an
empirical correction (Ecor) term is added. Currently, Ecor parameters
are only available for cases where the QM/MM partition cuts a C-C,
a C-O, or a C-S bond. For other cases, no empirical corrections
will be included. Unrestricted GHO-SCC-DFTB for open-shell system
is not implemented.
Reference: Reference made to the following paper, which contains
a more thorough description and discussion of test cases, is appreciated.
Jingzhi Pu, Jiali Gao, and Donald G. Truhlar,
J. Phys. Chem. A 108, 5454-5463 (1998). "Combining Self-Consistent-Charge
Density-Functional Tight-Binding (SCC-DFTB) with Molecular Mechanics by
the Generalized Hybrid Orbital (GHO) Method."
The current implementation has analytical first derivative and thus
allows energy minimizations, reaction path search (e.g., travel) and
molecular dynamics simulations; SCC-DFTB/MM also works with Monte Carlo.
Replica can also be used, which makes it possible to use replica path
and related approaches (such the nudged elastic band) for determining
reaction path with the SCC-DFTB/MM potential; along the same line,
path integral simulations can be carried out as well, although only for
equilibrium properties at this stage.
Several aspects of the code will be improved in the near future,
and new functionalities will be added:
1. Interface with centroid path-integral simulations and Tsallis
statistics.
2. More flexible interface with BLOCK for general free energy
simulations.
3. Better methods for open-shell systems; constrained density functional
theories.
4. Time-dependent treatment for electronically excited states; non-adiabtic MD.
5. Integration with polarizable force field models (Drude).
CHARMM Documentation / Rick_Venable@nih.gov