comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-28 22:25:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

...

Browse source
This commit is contained in:
Eric Bylaska 2005-02-09 23:47:37 +00:00
parent 3ba73ec4ec
commit 7f4f74d396
1 changed files with 228 additions and 62 deletions
Download patch file Download diff file Expand all files Collapse all files

290
doc/user/qmmm.tex
View file

@ -1,13 +1,13 @@
% $Id: qmmm.tex,v 1.10 2003-04-11 14:10:05 d3j191 Exp $
% $Id: qmmm.tex,v 1.11 2005-02-09 23:47:37 bylaska Exp $
\label{sec:qmmm}
Combined or hybrid Quantum Mechanics and Molecular Mechanics (QM/MM)
Combined or hybrid Quantum Mechanics and Molecular Mechanics (QMMM)
is a simulation methodology that is about 15 years old but in all the
literature there are cautions that calibration computations must be
done to validate the model for each particular chemical system
studied. This is not a black box style computation and the NWChem
users are advised that without calibration QM/MM may not give the
users are advised that without calibration QMMM may not give the
appropriate results\footnote{c.f., Singh and Kollman, J. Comp. Chem.
{\bf 7}, 718 (1986); M.~J.~Field, P.~A.~Bash and M.~Karplus, J.
Comp. Chem. {\bf 11}, 700, (1990); J. Gao, ``Methods and
@ -15,100 +15,266 @@ appropriate results\footnote{c.f., Singh and Kollman, J. Comp. Chem.
Potentials.'' In {\it Reviews in Computational Chemistry};
K.~B.~Lipkowitz, D.~B.~Boyd, Eds.; VCH Publishers: New York;
Vol. 7, pp 119-185 (1995); and M. A. Thompson and G. K. Schenter, J.
Phys. Chem {\bf 99} 6374 (1995) }.
Phys. Chem {\bf 99} 6374 (1995) }. Since both quantum-mechanical and classical
molecular mechanics are involved in the calculation good working knowledge of the two methods
is required to ensure meaningful results.
The QM/MM module in NWChem is driven by the molecular dynamics module.
This module currently works for any QM method that has
analytic gradients\footnote{The QM/MM method will work with numerical
gradients available in NWChem, but it is expected that the
performance will not allow any substantive simulations}. The input
for this requires the definition of chemical system via the same
interface that is used by the MD module (c.f. Section
\ref{sec:nwmd}). The extensions to this interface include the
definition of ``Quantum'' atoms and ``Link'' where appropriate. The
QM information must be present in the traditional NWChem input deck
except for the geometry\footnote{Any geometry information in the
traditional form will be ignored}. The geometrical information will
be constructed automatically by nwmd. For dynamics and free energy
simulations the input is again identical to that for nwmd with
limitations on the kinds of simulations that can be done.
Link atoms are defined by the molecular dynamics module, based on
the specification of the quantum atoms in the prepare module. For
the link atoms, basis sets \verb+X_L+ need to be defined in the
basis set input block.
The QM/MM module is invoked with the task directive where the
``theory'' is QMMM. The recognized operations on the QM/MM theory
directive are energy, optimize, and dynamics.
The QMMM module is invoked with the following task directive.
\begin{verbatim}
task qmmm (energy | optimize | dynamics)
task qmmm <string qmtheory> <string operation> [numerical] [ignore]
\end{verbatim}
where {\it qmtheory} specifies quantum method for the calculation of the quantum region. It is expected that
most of QMMM simulations will be performed with
with HF or DFT theories, but any other QM theory supported by NWChem should also work.
Currently the supported operations for QMMM runs
are energy, optimize, saddle, dynamics, numerical hessian, and numerical
frequencies.
Unlike pure
quantum mechanical calculations the information about the chemical system
for QMMM simulations is contained not in the geometry block but in the externally prepared topology and restart files.
These files have to be present prior to any QMMM simulation.
The input file for QMMM simulation can be divided into three major parts -- specification of the molecular
mechanics parameters for the classical region, specification of the quantum mechanical method for the quantum region,
and the parameters of the interaction between quantum and classical methods.
All this discussed in detail in the sections below.
\section{Preparation of the restart and topology files}
Generated by the prepare
module (see section \ref{sec:prepare}) restart and topology files contain
information about the classical force field as well as the
coordinates of quantum (qm) and molecular
mechanics (mm) regions.
In a typical setting this "preparation stage"
will be run separately from main QMMM simulation. This will require a
properly formatted
PDB file for the system. In more complex cases (e.g.non-standard residues or nucleotides) additional fragment and parameter
files might have to be provided by the user. The definition of the quantum region
in the input for the prepare module is specified by either {\it modify atom} directive (see Section \ref{sec:prepare}):
\begin{verbatim}
modify atom <integer isgm>:<string atomname> quantum
\end{verbatim}
Tasks \verb+gradient+, \verb+saddle+, \verb+frequencies+ and
\verb+thermodynamics+ are currently not available in the QM/MM mode.
or {\it modify segment} directive
\begin{verbatim}
modify segment <integer isgm> quantum
\end{verbatim}
Here {\it isgm} and {\it atomname} refer to the residue number and atom name record
as given in the PDB file.
It is important to note that
that the leading blanks
in atom name record should be indicated with underscores.
Per PDB format quidelines the atom name record starts at column 13. If, for example,
the atom name record "OW" starts
in the 14th column in PDB file, it will appear
as "\_OW" in the modify atom directive in the prepare block. In the current implementation
only solute atoms can be declared as quantum. If part of the solvent has to be treated quantum mechanically
then it has to redeclared to be solute.
In addition to modify commands the prepare input block should
also contain {\it update lists} and {\it ignore} directives. There are other options
that can be used in the input block for the prepare module ( e.g. solvating the structure, etc ),
those discussed
in more details in Section \ref{sec:prepare}.
The successful run of the prepare module will result in generation of
topology and restart files. Similar to classical MD, both files are required for QMMM simulations
and have to be placed in the same directory
as the input file. Here is an example input file that will generate QMMM restart and topology files for the ethanol molecule
\begin{verbatim}
title "Prepare QMMM calculation of ethanol"
start etl
prepare
#--name of the pdb file
source etl0.pdb
#--generate new topology and sequence file
new_top new_seq
#--generate new restart file
new_rst
#--define quantum region (note the use of underscore)
modify atom 1:_C1 quantum
modify atom 1:2H1 quantum
modify atom 1:3H1 quantum
modify atom 1:4H1 quantum
#
update lists
ignore
end
task prepare
\end{verbatim}
These are contents of etl0.pdb file used in the above input file.
\begin{verbatim}
ATOM 1 O etl 1 1.201 -0.271 -0.000 1.00 0.00 O
ATOM 2 H etl 1 1.995 0.329 -0.000 1.00 0.00 H
ATOM 3 C1 etl 1 -1.180 -0.393 0.000 1.00 0.00 C
ATOM 4 2H1 etl 1 -2.128 0.155 -0.000 1.00 0.00 H
ATOM 5 3H1 etl 1 -1.130 -1.030 0.887 1.00 0.00 H
ATOM 6 4H1 etl 1 -1.130 -1.030 -0.887 1.00 0.00 H
ATOM 7 C2 etl 1 0.006 0.573 0.000 1.00 0.00 C
ATOM 8 2H2 etl 1 -0.042 1.220 0.890 1.00 0.00 H
ATOM 9 3H2 etl 1 -0.042 1.220 -0.890 1.00 0.00 H
END
\end{verbatim}
Running the input shown above will produce (among other things) the topology file (etl.top) and the restart file
(etl{\_}md.rst). The naming of the topology file follows after the rtdb name specified in the start directive in the input (i.e. "start etl"),
while the "{\_}md" suffix in the restart file name is specific to the way prepare module works in this particular case. If necessary, this
particular naming scheme can be altered using {\it system} keyword in the prepare input block (for more details see Section \ref{sec:prepare}).
\section{Molecular Mechanics Parameters}
The molecular mechanics parameters are given in the form of standard MD input block as
used by the MD module (c.f. Section \ref{sec:nwmd}). This input block
is required for QMMM simulations. It specifies the
restart and topology file that will be used in the calculation.
It also contains information relevant to the calculation
of the classical region
(e.g. cutoff distances, constraints, optimization and dynamics parameters, etc)
in the system. In this input block one can also set fixed atom constraints on both classical and quantum atoms. Continuing with our
example for ethanol molecule here is a simple input block that may be used for this system.
\begin{verbatim}
md
# this specifies that etl_md.rst will be used as a restart file
# and etl.top will be a topology file
system etl_md
# if we ever wanted to fix C1 atom
fix solute 1 _C1
end
\end{verbatim}
The QM/MM input consists of the standard NWChem input block:
\section{Quantum Mechanical Parameters}
The parameters defining calculation of the QM region (including basis sets)
must be present in the traditional NWChem input format
except for the geometry block. The geometrical information will
be constructed automatically by QMMM using information from the MD module.
\section{QMMM interface parameters}
The QMMM interface parameters define the interaction between classical and quantum regions.
The input follows standard NWChem format:
\begin{verbatim}
qmmm
...
[ eref <double precision default 0.0d0>]
[ bqzone <double precision default 9.0d0>]
[ bq_exclude <(none||all||linkbond||linkbond_H) default none>]
[ bq_update <(static||dynamic) >]
[ link_atoms <(hydrogen||halogen) default halogen>]
[ link_ecp <(auto||user) default auto>]
[ optimization <(all||mm|qm) default qm>]
end
\end{verbatim}
The \verb+qmmm+ has the following the additional sub-directive that the user
may specify for the particular simulation. These options currently are:
Detailed explanation of the subdirectives in the QMMM input block is given below:
\begin{itemize}
\item
\begin{verbatim}
eatoms < real eatoms>
eref <double precision default 0.0d0>
\end{verbatim}
There is one compound input directive that must exist for the QM/MM
simulation to proceed. This sets the relative zero of energy for the
QM component of the system. It is not incorrect to leave this value as
zero but the energetics of the QM system will likely over shadow the
MM component of the system. Properties based on energy fluctuations
of the system will be overly sensitive to the energy of the QM
component of the system. The zero of energy for the MM system is by
This directive sets the relative zero of energy for the
QM component of the system. The need for this directive
arizes from different definitions of zero energy for QM and MM methods.
Most QM methods define the zero of energy for the system as
vacuum. The zero of energy for the MM system is by
definition of most parameterized force fields the separated atom
energy. The zero of energy for QM systems by definition of most QM
methods is the vacuum. The {\it a priori} determination of the
separated atom energy for a particular QM method is not well defined
and thus leads to a number of assumptions or guess work depending upon
the particular QM method being utilized. Therefor, the determination
of the QM separated atom energy (``eatoms'') is left to the user.
There is no default for this and the input {\bf must} be present for a
QM/MM simulation.
energy. Therefore in many cases the energetics of the QM system
will likely overshadow the
MM component of the system. This imbalance can be corrected by
suitably chosen value of {\it eref}
\item
\begin{verbatim}
qatoms < real qatoms>
bqzone <double precision default 9.0d0>
\end{verbatim}
This input directive is used to specify the total charge on the set
of quantum atoms.
This directive defines the radius of the zone (in angstroms) around the quantum region
where classical residues/segments
will be allowed to interact with quantum region electrostatically. It should be noted
that classical atoms interacting with quantum region via bonded interactions are always
included in the bqzone (this is true even if bqzone is set to $0.0$). In addition, even if one atom
of a given charged group is in the bqzone (residues are typically treated as one charged group) then the whole
group will be included in the bqzone.
\item
\begin{verbatim}
link hydrogen
bq_exclude <(none||all||linkbond||linkbond_H) default none>
\end{verbatim}
This directive specified to use hydrogens as link atoms. The default
is to use the next group VII atom.
This directive operates in conjunction with bqzone keyword offering additional level of control
on electrostatic interactions between quantum and classical regions. Default value {\it none} will
leave bqzone unchanged, {\it all} will remove all atoms from the bqzone resulting in
no electrostatic interactions between
classical and quantum regions,
{\it linkbond} will result in the removal of all atoms that are connected to a quantum region
by at most two bonds,
{\it linkbond{\_}H} is similar to {\it linkbond} but will only remove hydrogen atoms.
If necccessary, further modifications of the electrostatic interactions can be achieved
using the prepare module{\ref{sec:prepare}}.
\item
\begin{verbatim}
nobq ( hydrogen | all | none )
bq_update <(static||dynamic)
\end{verbatim}
This directive controls whether the bqzone stays fixed ({\it static}) or constantly
updated ({\it dynamic}) during the calculation. In most cases direct specification of this keyword
will not be necessary as its value will be set automatically based on the nature of the calculation.
\item
\begin{verbatim}
link_atoms <(hydrogen||halogen) default halogen>
\end{verbatim}
This directive specifies that hydrogen atoms, all atoms, or no atoms,
respectively, bonded to link atoms will carry a zero charge in the
QM part of the calculation.
This directive controls the treatment of bonds crossing the boundary between quantum and classical regions.
The use of {\it hydrogen } keyword will trigger truncation of such bonds with hydrogen link atoms. The position of the hydrogen
atom will be calculated from the coordinates of the quantum and classical atom of the truncated bond using the
following expression
\begin{displaymath}
\mathbf{R}_{hlink} = (1-g)\mathbf{R}_{quant} + g*\mathbf{R}_{class}
\end{displaymath}
where $g$ is the scale factor set at $0.709$
Setting link{\_}atoms to {\it halogen } will result in the modification of the {\it \underline{quantum}}
atom of the truncated bond to
to the fluoride atom. This fluoride atom will typically carry an effective core potential (ECP) basis set as specified
in {\it link{\_}ecp} directive.
\item
\begin{verbatim}
link_ecp <(auto||user) default auto>
\end{verbatim}
This directive specifies ECP basis set on fluoride link atoms. If set to {\it auto }
the ECP basis set given by Zhang, Lee, Yang for 6-31G* basis.\footnote{Y. Zhang, T. Lee, and W. Yang, J. Chem. Phys. 110, 46 (1999)}
will be used. Strictly speaking, this implies the use of 6-31G* spherical basis as the main basis set.
If other choices are desired then keyword {\it user } should be used and ECP basis set should be entered separatelly
following the format given in
section \ref{sec:ecp}.
The name tag for fluoride link atoms is F{\_}L.
\item
\begin{verbatim}
optimization <(all||mm|qm) default qm>
\end{verbatim}
This directive specifies which region will be optimized during qmmm optimization. If set to {\it all} both qm and mm
region will be optimized simultaneously using mm optimization module, if {\it mm} only mm portion will be optimized
using mm optimization module, finally if {\it qm} is specified then only qm region will be optimized using driver
module.
\end{itemize}
All other parameters that control the QM/MM simulation are set via the
input to nwmd (see chapter \ref{sec:nwmd}).