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 (1998) 7260,
Phys. Stat. Sol. B 217 (2000) 357,
J. Phys. : Condens. Matter. 14 (2002) 3015.
The QM/MM interface in CHARMM has been described in
J. Phys. Chem. B 105 (2001) 569
* Menu:
* Description:: Description of the sccdftb commands.
* Usage:: How to run sccdftb in CHARMM.
* Installation:: How to install sccdftb in CHARMM environment.
* Status:: Status of the interface code.
The SCCDFTB QM potential is initialized with the SCCDFTB command
[SYNTAX SCCDFTB]
SCCDFTB [REMOve] [CHRG] (atom selection) [TEMPerature] [SCFtolerance]
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.
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.
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.
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.
For the example of the system containing O N C and H, sccdftb.dat would
contain:
'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'
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 requested from Marcus Elstner
(elstner@phys.upb.de).
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). Diagonalization routines ------------------------ As default, the library routine dsygv.f (LAPACK) is used for the diagonalization of the hamiltonian matrix. This is called by chmdir/source/scctbint/scctbsrc/ewevge.f. A faster (about factor 2) solution is given by the dsygvd.f routine (but less stable), which is called by ewevge-dsygvd.f: copy ewevge-dsygvd.f to ewevge.f and recompile to invoke this option. Contact Marcus Elstner for more details or questions.
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. Handling of QM/MM electrostatic interaction: currently no cut-off is
used, which will be changed in the future. Ewald summation will
also be implemented.
2. Output from SCCDFTB (such as Mulliken charges) will be better
organized.
3. Dispersion interactions among SCCDFTB atoms will be added.
4. Interface with centroid path-integral simulations and Tsallis
statistics.
5. Interface with BLOCK for free energy simulations.
6. Unrestricted solution for open-shell systems.
7. Time-dependent treatment for electronically excited states.
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory
Modified, updated and generalized by C.L. Brooks, III
The Scripps Research Institute