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@par10.mgsl.dcrt.nih.gov,milan@kihp6.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).
* Menu:
* Description:: Description of the qchem commands.
* Usage:: How to run Q-Chem in CHARMM.
* Installation:: How to install Q-Chem in CHARMM environment.
* Functionality:: Functionality of the interface code.
* Status:: Status of the interface code.
The Q-Chem QM potential is initialized with the QCHEM command.
[SYNTAX QCHEm]
QCHEm [REMOve] [EXGRoup] NOGUess] (atom selection)
REMOve: Classical energies within QM atoms are removed.
EXGRoup: QM/MM Electrostatics for link host groups removed.
NOGUess: Obtains initial orbital guess from previous calculation.
Default is to recalculate initial orbitals each time.
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.
=======================================================================
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 2004) - Parallel version is fully functional - Replica/Path and Nudged Elastic Band Methods function in a highly parallel and parallel/parallel fashion (parallel/parallel is not implemented at the source code level, accomplished via CHARMM scripting). - 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 ==============================================================================
1. QM/MM optimizations (analytic gradients) using Q-Chem can be performed using the following methods. - HF* (RHF, UHF, ROHF) - DFT* (RHF, UHF, ROHF) - MP2 (RHF, UHF, ROHF) - 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) 3. The remainder of Q-Chem's analytic derivative and energy point methods will be made available in future releases. ==============================================================================
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. ==============================================================================
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory
Modified, updated and generalized by C.L. Brooks, III
The Scripps Research Institute