Combined Quantum Mechanical and Molecular Mechanics Method
Based on Q-Chem in CHARMM
H. Lee Woodcock
(hlwood@nih.gov)
based on the GAMESS(US) interface from Milan Hodoscek
(milan@cmm.ki.si)
and
the GAMESS(UK) interface from Paul Sherwood
(p.sherwood@dl.ac.uk)
Ab initio program Q-Chem is connected to CHARMM program in a
QM/MM method. This method is based on the interface to the GAMESS (US
version), the latter being an extension of the QUANTUM code which is
described in J. Comp. Chem., Vol. 11, No. 6, 700-733 (1990).
The QM/MM interface between Q-Chem and CHARMM is described in the
following work which should be cited when used...
H. Lee Woodcock, M. Hodosceck, A. T. B. Gilbert, P. M. W. Gill, H. F. Schaefer,
B. R. Brooks; Interfacing CHARMM and Q-Chem to perform QM/MM and QM/MM reaction
pathway calculations. J. Comp. Chem.; 2007; 28 (9); 1485-1502.
* Menu:
* Description:: Description of the qchem commands.
* Usage:: How to run Q-Chem in CHARMM.
* Installation:: How to install Q-Chem in CHARMM environment.
* Status:: Status of the interface code.
* Functionality:: Functionality of the interface code.
* RPath:: Replica Path Command
* Pert:: ab inition QM/MM free energy perturbation
* Normal Mode Analysis:: Full QM/MM Normal Mode Anal. through VIBRAN
* Microiterations:: QM/MM Microiteration optimizations
* MMQM:: Write internal / external Q-Chem input file
* PCM:: Specifications for QM/MM/PCM with Q-Chem
The Q-Chem QM potential is initialized with the QCHEM command.
[SYNTAX QCHEm]
QCHEm [REMOve] [EXGRoup] [DIV] [NOGUess] [BLURred [RECAll INT]] [COORdinates]
[QCLJ] [PARAllel [INT]] [PCM] [[NREStart [INT]] [SINPut] [SOUTput]
[SGRAdient] [RGRAdient] [SHESsian] [RHESsian] [CHARge] [MICRo] [SAVE]
[RESTart] [RESEt] [see below for more options...] (atom selection)
REMOve: Classical energies within QM atoms are removed.
EXGRoup: QM/MM Electrostatics for link host groups removed.
DIV: Charge on MM link host atom divided equally among other MM atoms in
the same group.
NOGUess: Obtains initial orbital guess from previous calculation.
Default is to recalculate initial orbitals each time.
BLURred: MM charges are treated as a gaussian function (equivalent to ECP)
width of the gaussian function is specified by default in WMAIN
array (usually by SCALar command). The value for charge is taken
from PSF. Some values of WMAIN have special meaning:
WMAIN.LE. 0.0 treat this atom as point charge in the QM/MM potential
WMAIN.GE.5000.0 treat this atom as an infinitely diffuse Gaussian
RECAll: Use the RECAll array (as specified in scalar.doc) to set BLUR
widths instead of the main WMAIN array. This is necessary when
using Gaussian BLUR MM charges with the Replica Path or NEB
methods as these make use of the WMAIN array. See QM-MM_DGMM.inp
in the test directory for an example.
COORdinates: This keyword will activate CHARMM to obtain an updated geometry
from Q-Chem as the calculation proceeds. This can be particularly
useful as Q-Chem can perform QM optimiztions using delocalized
internal coordinates in the presence of a fix field of point charges.
This can significantly speed QM/MM minimizations and can be used in
an iterative approach. Note: to use this the JOBTYPE in the Q-Chem
control file should be set to OPT (i.e. JOBTYPE = OPT).
QCLJ: Activates Q-Chem to use CHARMM's Lennard-Jones parameters when
performing QM calculations in a fixed field of point charges. This
can be particularlly useful as the QM region can be overly attracted
to bare point charges.
PARAllel: Allows the user to specify how many processors they wany the Q-Chem
calculation to utilize. Previously, Q-Chem would use the same number
of processors as CHARMM was using, however, in most cases the Q-Chem
calculation will be much more expensive so having 1 CHARMM process
and 4 Q-Chem processes is more efficient. Note: This currently only
works with parallel versions of CHARMM although it can be extended
to work with serial versions.
PCM: Turns on the use of QM based continuium solvent methods in Q-Chem.
The implicit solvent methods in Q-Chem support both QM/CPCM and
QM/MM/CPCM. For full details on the CPCM method see: J. Chem. Phys.
133, 244111 (2010). See below for more detail...
NREStart: To prevent calculations / simulations from abruptly terminating due
to SCF parallel communication failures the NREStart keyword was added
that allows the user to specify the number of times to retry a
particular Q-Chem QM/MM energy and force calculation before terminating
the overall CHARMM process. This keyword takes an integer as its
argument (e.g. NREStart 3).
SINP / SOUT: The SINP and SOUT keywords activate Q-Chem to save input and output
file, respectively, for each step of a minimization or simulation. Input
and output files are saved into the following directories created by
CHARMM: saved_inputs and saved_outputs, respectively, with the step
number of the minimization or simulation appended to the end of the file
(e.g. q1.inp_22).
SGRAdient / RGRAdient: The SGRAdient and RGRAdient keywords activates the Q-Chem/
CHARMM interface to save an individual gradient and read that gradient
back in at a later point in the calculation. This can be particularly
useful when employing methods that require manipulation done to the
gradient at the same time they are acting on the Hessian. Specifically,
this option was added to facilitate QM/MM Mobile Block Hessian
calculations. See the following paper for full details of the QM/MM MBH
method: JCTC, DOI: 10.1021/ct100473f, 2011.
SHESsian: Save Hessian computed via the QM/MM Normal Mode Procedure. The
Hessian will be saved as an ascii file named: hessian.dat. Typically
VIBRAN recomputes the Hessian each time it is needed; for QM or QM/MM
calculations this is inefficient and thus saving the Hessian becomes
very important. Particularly, this is used when employing the VSA
method (see vibran.doc).
RHESsian: Read a previously saved Hessian (hessian.dat) from a file (see SHES).
NOLInkatoms: This option forces all atoms in a Q-Chem/CHARMM calculation to be
treated as QM atoms. This is not intended for use in standard
energy or force (i.e., minimization, dynamics, etc.). Rather, this
option is only intended for use when doing QM/MM normal mode
analysis via ONIOM-type MSCALE calculations. And even then should
only be used on the QM subsystem(s) that is
CHARge: Read QM charges from an file (charges.dat). A charges.dat file is
created by Q-Chem by seting the REM keyword QMMM_CHARGES = TRUE.
This file contains the Muliken charges for the QM region in the
same order that is specified in the PSF file.
MICRo: Turns on the QM/MM Micro-iteration scheme of Kastner et al. J. Chem.
Theory Comput., 3 (3), 1064-1072, 2007. Currently, this is best used
in conjunction with a loop where CHARMM alternates between MM and
QM/MM micro and macro cycles. Also, the CHARge keyword should be
used to set the charges on QM region during the MM cycles. This is a
new feature that requires further testing so be careful!
WQIN This option will write a Q-Chem input file, but will not execute it
(W=Write,Q=QChem, IN=INput file --> WQIN).
OMP This open tell Q-Chem to use the multi-threaded parallel version. This
is the recommended option as performance is improved over the
distrubted parallel version.
MIXed Option for setting up a mixed basis set calculation calculation (i.e.,
using different basis sets for different parts of the QM region). This
should enable significant saving when large QM regions are needed, but
some atoms can be treated more approximately. For more information see
the Q-Chem 4.2 manual (Section 7.5).
BAS1 Use this in conjunction with MIXed (see above). This should be followed
the name of the basis set you want to assign to your first atom
selection (which you must "define" as "basis1").
BAS2 Use this in conjunction with MIXed (see above). This should be followed
the name of the basis set you want to assign to your second atom
selection (which you must "define" as "basis2").
SCRAtch Allows the seletion of a user defined scratch space for Q-Chem jobs. It
should be followed by the path you want to use.
EWALd This activated CHARMM to write out a Q-Chem specific parameter file and
sets up / runs a QM/MM EWALD calculation via Q-Chem. Currently, only single
point energy calculations are supported. Please see the Q-Chem manual for
more details about this procedure and rules on setting QMALpha and MMalpha
(see below). To fully run utilize this command the following command will
also need to be executed after setting the options to "qchem"...
write para qchem card used name usedpara.prm
this is command pull the currently used force field parameters from CHARMM
and uses these in the QM/MM EWALD calculation. This is very helpful as there
is no guarantee that force fields used for the Q-Chem EWALD calculation will
be up to date with the most current CHARMM force fields.
QMALpha The QM alpha(kappa) value that gets passed to Q-Chem during QM/MM EWALD.
MMALpha The MM alpha(kappa) value that gets passed to Q-Chem during QM/MM EWALD.
MESS Activates the QM/MM MESS procedure.
NROOts The number of roots Q-Chem will solve for in the QM/MM MESS procedure.
RESDi Activates Q-Chem to use CHARMM's RESDistance information to perform
restrained QM calculations with a fixed field of point charges. Note,
this is best used in conjunction with the COORd command and setting
jobtype=opt in the Q-Chem control file. This functionality requires
CHARMM to have a restraint set as a linear combination of distances.
This option has been tested with standard QM/MM calculations and the
Replica Path functionality in CHARMM.
CONS This instructs the RESDi command (above) that CHARMM is only passing
restraint information between 2 atoms rather than a linear combination
of distances.
RESEt: Resets all QM/MM options to their initial defaults. This is needed
for the QM/MM Micro-iteration approach to alternate between MM and
QM/MM stages. Note: after using this a new "QCHEM" command must be
issued!
SAVE: Activates CHARMM to save the converged SCF orbitals from a given
energy calculation. Note: This ideally should be called once using
a specific Q-Chem control file that contains specialized SCF
convergence options. This option shoudl then be followed by a new
"QCHEm" call that specifies "RESTart" and performs the actual
QM/MM minimization. See below for example...
RESTart: This tells CHARMM to restart a QM/MM calculation using previously
saved orbitals. The "QCHEm" command that uses this as an option
should be proceeded by a "QCHEm SAVE" command to preform an
initial calculation saving the orbitals. See below for example...
Example:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
envi qchemcnt "qcnt1.inp" ! File that contains special SCF
envi qcheminp "q1.inp" ! convergence options
envi qchemexe "qchem"
envi qchemout "q1.out"
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
QCHEm SAVE REMOve SELEct RESId 1 SHOW END
ENERgy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
envi qchemcnt "qcnt2.inp" ! Regular Q-Chem control file
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
QCHEm RESTart NOGUess REMOve SELEct RESId 1 SHOW END
MINI ABNR NSTEp 10 NPRInt 1
=======================================================================
The atoms in selection will be treated as QM atoms.
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.
If you see the following error in your output script:
FNIDEL> Cannot find element type for number....
That means you either have wrong order in the ADDLink command or the atom
that should be MM is in the QM selection.
=======================================================================
CHARMM input scripts are the same as before except the addition of ENVIronment
commands and the QCHEm command itself. Q-Chem commands are in a separate file
call qchem.inp, (or with an alternative name indicated by the "QCHEMCNT"
environment variable). The Q-Chem input file has the same structure as it
would have for a normal Q-Chem run, except that the specification of the
geometry, in the molecule section, is omitted. Note: the charge and
multiplicity are still included in the molecule section.
Names of the files for Q-Chem are specefied with environment
variables as follows. These four ENVIronment variables must be set!
use ENVIronment command inside CHARMM
ENVI qchemcnt "qchem.inp"
ENVI qcheminp "q1.inp"
ENVI qchemexe "qchem"
ENVI qchemout "qchem.out"
or use the following for (t)csh
setenv qchemcnt qchem.inp
setenv qcheminp q1.inp
setenv qchemexe qchem
setenv qchemout qchem.out
or use the following for ksh,sh,bash
export qchemcnt=qchem.inp
export qcheminp=q1.inp
export qchemexe=qchem
export qchemout=qchem.out
1. The QCHEMCNT variable specifies the main Q-Chem input file which contains
the $rem section, $molecule section (without geometry), $comment section,
ect..,
2. The QCHEMINP variable is the final input file that will get passed to
Q-Chem. CHARMM actually writes this file and adds the correct geometry and
any external/point charges (e.g. MM atoms) to an $external_charges section.
3. The QCHEMEXE is the location of the qchem script. Specify the entire path
unless $QC/bin is included in your default path.
4. The QCHEMOUT file specifies the Q-Chem output file. This file get
overwritten for each optimization/time step. In the future, there will be a
mechanism to save old output files.
Q-Chem input file parameters
----------------------------
The following $rem variables must be specified in the QCHEMCNT file in order
to perform CHARMM QM/MM or pure QM calculations.
qm_mm true
jobtype force
symmetry off
sym_ignore true
print_input false
qmmm_print true
1. qm_mm = true: Turns QM/MM on in Q-Chem
2. jobtype = force: Needed to do QM/MM optimizations. Set to "SP" if QM/MM
energy is desired.
3. symmetry = off: Turn off symmetry
4. sym_ignore = true: Prevents Q-Chem from reorienting molecule
5. print_input = false: Use this if you have a large molecule and do not want
1000s of atoms echoed back to the output file.
6. qmmm_print = true: Reduces some of the print out during QM/MM calculations.
This prevents external charges from being printed out if
there are more than 50 of them.
Sample QCHEMCNT file (qchem.inp):
---------------------------------
$comment
Input file comes from CHARMM
$end
$rem
exchange HF
basis 6-31G*
qm_mm true
jobtype force
symmetry off
sym_ignore true
print_input false
qmmm_print true
$end
$molecule
0 1
$end
-----------------------------------------------------------------------------
The above is for 6-31G calculation of any neutral molecule.
[NOTE: For another example look at test/cquantumtest/alanine_qchem.inp]
==============================================================================
One of the main benefits of using Q-Chem to do QM/MM calculations with CHARMM is the ease of which you can get up and running jobs. All you have to do is compile CHARMM in the following way.... install.com <machine-type> <CHARMM size> QC <other CHARMM options> This will compile the serial version of CHARMM to run with a serial version of Q-Chem. To compile a parallel version of CHARMM to run with a parallel or serial version of Q-Chem you could use the following script.... ----------------------------------------------------------------------------- #!/bin/csh # Compile Parallel CHARMM with Q-Chem support # USE STANDARD MPI (i.e. MPICH) setenv MPI /base/mpi/directory setenv MPI_LIB $MPI/lib setenv MPI_LIB $MPI/include # SET THE PATH TO MPIF77 set path=($MPI/bin $path) install.com <machine-type> <CHARMM size> M QC MPICH <other CHARMM options> ----------------------------------------------------------------------------- ==============================================================================
Q-Chem/CHARMM interface status (July 2007) - Parallel version is fully functional - Replica/Path and Nudged Elastic Band Methods function in a highly parallel and parallel/parallel fashion. - I/O including standard input and output are separated for Q-Chem. - All CHARMM testcases are still OK when CHARMM is compiled with Q-Chem inside. - QCHEM, GAMESS, GAMESSUK, CADPAC and QUANTUM keywords cannot coexist in pref.dat - Q-Chem recognizes atoms by their masses as specified in the RTF file - Delocalized Gaussian Blurred MM charges have been implemented for both energies and analytic gradients - Full (i.e. no restraints/constraints) QM/MM 2nd derivatives (i.e. Hessians) are available. ==============================================================================
1. QM/MM optimizations (analytic gradients) using Q-Chem can be performed using the following methods. - HF* (RHF, UHF, ROHF) - DFT* (RHF, UHF, ROHF) - RIMP2 (RHF, UHF) - MOS-MP2 (RHF, UHF) - SOS-MP2 (RHF, UHF) - SCS-MP2 (RHF, UHF) - MP2 (RHF, UHF) - CCSD (RHF, UHF) * Analytic derivatives run in parallel. 2. QM/MM single point energies using Q-Chem can be performed using the following methods (in addition to the above). - Local MP2 (RHF, UHF) - CCSD(T) (RHF, UHF) 3. Additional analytic derivative and energy point methods will be made available in future releases. To request support for methods please contact H. Lee Woodcock (hlwoodr_at_nih_dot_gov) and/or post request to the CHARMM forums. ==============================================================================
1. Additional ENVIronment variable: To do QM/MM Replica/Path or Nudged Elastic
Band calculations with CHARMM and Q-Chem you must define one extra variable.
ENVI QCHEMPWD "/path/to/working/rpath/directory"
2. After defining this above ENVIronment variable all that is left to do is
add the "rpath" keyword to the QCHEm call. For example...
QCHEm RPATh REMOve select qm_region end
This will create nrep directories in /path/to/working/rpath/directory and each
point of the pathway will be computed in a different directory.
Note: you must be running a parallel version of CHARMM with the same number of
processors as you have replicas (i.e. pathway points).
==============================================================================
To run ab initio QM/MM free energy perturbation you need to specify additional environment variables in the QM/MM setup... 1. sainp: state A control file (same as QCHEMCNT; specific for state A) 2. sbinp: state B control file (same as QCHEMCNT; specific for state B) 3. stateainp: auto generated Q-Chem input file for state A 4. statebinp: auto generated Q-Chem input file for state B 5. stateaout: specify Q-Chem output for state A QM calculation 6. statebout: specify Q-Chem output for state B QM calculation Example... envi qchemexe "qchem" ! Command to call quantum program envi qchemcnt "data/qchem_pert.inp" ! Non Pert Control file envi qcheminp "q1.inp" ! Non Pert Quantum input file envi qchemout "qchem.out" ! Non Pert Quantum output file envi sainp "data/s0.inp" ! State 0 control file envi sbinp "data/s1.inp" ! State 1 control file envi stateainp "state0.inp" ! State 0 quantum input file envi statebinp "state1.inp" ! State 1 quantum input file envi stateaout "state0.out" ! State 0 quantum output file envi statebout "state1.out" ! State 1 quantum output file See test/cquantumtest/qmmm_pert.inp for a complete example. Please see pert.doc for a complete description of running free energy perturbation in CHARMM. ==============================================================================
To run full QM/MM Normal Mode Analysis (i.e. QM/MM 2nd derivatives, Hessians) you need to run QM/MM with the VIBRan module (see vibran.doc) of CHARMM. To perform this calculation just run QCHEm has usual... Example: QCHEm REMOve SELEct RESId 1 SHOW END Then invoke the VIBRan module... VIBRan DIAG END In addition, you must add the following line to the QCHEMCNT file (the file that controls the the REM variables passed to Q-Chem). QMMM_FULL_HESSIAN TRUE Please see the "QM-MM_Normal_Modes.inp" testcase in the test directory for the full example. Currently, this only works with standard point charge QM/MM models (i.e. not Gaussian blurred charges), but this will be extended in the future. Additionally, we have introduced a new keyword to facilitate QM/MM MSCALE (i.e., ONIOM-type) normal mode analysis: NOLInkatoms (see above). This ensures that if some QM atoms in a subsystem are labeled as link atoms, they will be treated correctly (i.e, as a real QM atom) by the subsystem calculation. For energy or force calculations this is not a problem, but for normal mode calculations this can cause problems due to link atom projection methods used. ==============================================================================
Here is an example of how to typical microiteration setup may work: ! Verify current charges ! ---------------------- scalar charge show select segid MAIN end qchem noguess remove sele segid MAIN show end energy ! All QM Charges should 0.0 here b/c they ! zeroed out as part of the QCHEM call ! --------------------------------------- scalar charge show select segid MAIN end ! Reset Q-Chem counters qchem reset select segid MAIN end ! Run a QM job and get charges ! ---------------------------- qchem micro charge remove sele segid MAIN show end energy scalar charge show select segid MAIN end ! Fix QM regions and run MM minimization ! -------------------------------------- cons fix select segid MAIN end mini abnr nstep 100 cons fix select none end scalar charge show select segid MAIN end ! MM energy with cons fix removed ! ------------------------------- energy scalar charge show select segid MAIN end ! Reset Q-Chem counters ! --------------------- qchem reset select segid MAIN end ! Compute QM energy again ! ----------------------- qchem remove sele segid MAIN show end energy ==============================================================================
Write internal (to use with CHARMM) or external (to use as stand alone input) input files for Q-Chem: As part of running Q-Chem with CHARMM a Q-Chem control file is needed. This can now be created on the fly using the [MMQM] functionality. open write unit 3 card name qchem.inp MMQM select qmregion end unit 3 CNTL $REM EXCHANGE HF BASIS STO-3G QM_MM TRUE JOBTYPE FORCE SYMMETRY OFF SYM_IGNORE TRUE PRINT_INPUT TRUE QMMM_PRINT FALSE $END $MOLECULE 0 1 $END END close unit 3 An alternative use of this functionality is to generate input files for use with Q-Chem as an independent of CHARMM. These can be created as follows: open write unit 4 card name qchem.inp MMQM select qmregion end unit 4 $REM EXCHANGE HF BASIS STO-3G QM_MM TRUE JOBTYPE FORCE SYMMETRY OFF SYM_IGNORE TRUE PRINT_INPUT TRUE QMMM_PRINT FALSE $END $MOLECULE 0 1 QCHEM_MOLECULE $END END ==============================================================================
The [PCM] keyword activates the Q-Chem/CHARMM interface to expect either a QM/PCM or QM/MM/PCM calculation. The use of QM continuim solvent methods in Q-Chem requires special input; see below for an example with full details being found in the Q-Chem manual. $rem exchange hf basis sto-3g qm_mm true jobtype force symmetry off sym_ignore true print_input true qmmm_print false solvent_method pcm $end $molecule 0 1 $end $pcm_solvent dielectric 1.0001 $end $pcm Theory CPCM Method SWIG HeavyPoints 50 HPoints 50 Radii FF vdwScale 1.0 $end ==============================================================================
CHARMM Documentation / Rick_Venable@nih.gov