Combined Quantum Mechanical and Molecular Mechanics Method
Based on GAMESS in CHARMM
by Milan Hodoscek
(milan@par10.mgsl.dcrt.nih.gov,milan@kihp6.ki.si)
Ab initio program GAMESS (General Atomic and Molecular
Electronic Structure System) is connected to CHARMM program in a QM/MM
method. This method is extension of the QUANTUM code which is
described in J. Comp. Chem., Vol. 11, No. 6, 700-733 (1990).
* Menu:
* Description:: Description of the gamess commands.
* Using:: How to run GAMESS in CHARMM.
* Replica path:: How to run GAMESS/CHARMM with REPLICA/PATH.
* Installation:: How to install GAMESS in CHARMM environment.
* Status:: Status of the interface code.
* Functionality:: Functionality of the interface code.
* Implementation:: Implementation.
The GAMESS QM potential is initialized with the GAMEss command.
[SYNTAX GAMEss]
GAMEss [REMOve] [EXGRoup] [QINPut] [BLURred] [NOGUess] [FMO] (atom selection)
REMOve: Classical energies within QM atoms are removed.
EXGRoup: QM/MM Electrostatics for link host groups removed.
QINPut: Charges are taken from PSF for the QM atoms. Charges
may be non integer numbers. Use this with the REMOve!
NOGUess: Obtains initial orbital guess from previous calculation.
Default is to recalculate initial orbitals each time.
FMO: Enable Fragment MO method with CHARMM
BLURred: MM charges are scaled by a gaussian function (equivalent to ECP)
Width of the gaussian function is specified in WMAIN array
(usually by SCALar command)
The value for charge is taken from PSF. Some values of WMAIN have
special meaning:
WMAIN.GT.999.0 ignore this atom from the QM/MM interaction
WMAIN.EQ. 0.0 treat this atom as point charge in the QM/MM potential
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. Also note that QQH type has
to be added in the RTF file (see test/c25test/gmstst.inp).
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.
=======================================================================
In order to run GAMESS and CHARMM on parallel machines I/O of
GAMESS and CHARMM was separated. This is now true even for scalar
runs. CHARMM input scripts are the same as before except the addition of
ENVIronment commands and GAMEss command itself. GAMESS commands are in a
separate file which is pointed to by INPUT environment variable.
Names of the files for GAMESS are specefied with environment
variables as follows:
use ENVIronment command inside CHARMM
envi INPUT "test.gms" ! quotes needed for lowercase names
envi OUTPUT "test.out"
envi PUNCH "scratch/test.dat"
envi DICTNRY "scratch/test.F10"
envi WORK15 "scratch/test.F15"
envi DASORT "scratch/test.F20"
or use (t)csh
setenv INPUT test.gms
setenv OUTPUT test.out
setenv PUNCH scratch/test.dat
setenv DICTNRY scratch/test.F10
setenv WORK15 scratch/test.F15
setenv DASORT scratch/test.F20
or ksh,sh,bash
export INPUT = test.gms
export OUTPUT = test.out
export PUNCH = scratch/test.dat
export DICTNRY = scratch/test.F10
export WORK15 = scratch/test.F15
export DASORT = scratch/test.F20
For complete information about GAMESS input see INPUT.DOC file in
GAMESS distribution. (NIH: ~milan/gamess/hp/INPUT.DOC)
Example:
--------
GAMESS commands have to be in a separate file. Example for the GAMESS input
follows:
-----------------------------------------------------------------------------
$CONTRL COORD=CART NOSYM=1 NPRINT=-5 ! This is rarely changed
SCFTYP=RHF ICHARG=0 ! This usually has to be changed
RunTyp=Gradient ! Normally forces are needed
! RunTyp=Energy ! If only energy is needed
$END
$SYSTEM MEMORY=1000000 ! memory allocation
TIMLIM=100000 $END
$BASIS GBASIS=N31 NGAUSS=6 $END
$SCF DIRSCF=.True. $END ! DIRSCF=.true. recommended
! if there are convergence problems
! try SOSCF=.FALSE.
$DATA
! This can be empty
$END
------------------------------------------------------------------------------
The above is for 6-31G calculation of any neutral molecule. $DATA
section may be left empty or filled with basis set information in the case
when it cannot be specified by the $BASIS keyword.
[NOTE: For more examples look at test/c25test/gmstst.inp]
==========================================================
Replica/Path method (parallel/parallel setup)
---------------------------------------------
Running GAMESS/CHARMM interface with Replica/Path method needs
few additional steps:
- GAMESS/CHARMM must be compiled with the parallel
functionality. Make sure that the GENCOMM keyword is in
pref.dat. (Run CHARMM interactively and type pref).
- The number of processes must be equal to number of replicas
multiplied by an integer (1,2,3...). This ensures that each
replica is an independent process. If the factor is more
than 1, it means each replica will run itself in parallel
(parallel/parallel).
- GAMESS control file (the one assigned to the INPUT environment
variable) must be linked the number of replica times. Each
symbolic link must have _<int> appended to the original
name:
ln -s test.gms test.gms_1
ln -s test.gms test.gms_2, etc
the number of links must be greater or equal to the number
of replicas
- The path to the above link must be absolute. This depends on
the way CHARMM is run in parallel. For example for MPICH
library on must use the following command:
charmm -p4wd /data/rpath/reaction -p4pg 20cpus < inp > out
The /data/rpath/reaction must be the same on all the
processes, either exact copies or NFS mounted.
- The gamess output files have also _<int> appended to their names.
Installation
------------
Look at the GAMESS home page for instructions how to obtain the code.
Installation itself cannot be automated yet so one has to
follow this procedure (if there are any problems ask milan@cmm.ki.si):
1. Put all the source (*.src and *.c) files in source/gamint/gamess
2. Follow instructions in the begining of the gamess.src file, to
change:
C
C ----- CHARMM INTERFACE -----
C TO USE GAMESS FROM INSIDE OF CHARMM, YOU MUST
C 1. INITIALIZE KCHRMM JUST BELOW TO 1
C 2. CHANGE "PROGRAM GAMESS" ABOVE TO "SUBROUTINE GAMESS"
C 3. CHANGE THE "STOP" STATEMENT BELOW TO "RETURN"
C 4. DELETE DUMMY SUBROUTINES -CHGMIU- AND -CHMDAT- BELOW
C 5. CHANGE -MXCHRM- FROM 1 TO 25120 IN ALL PARAMETER DEFINITIONS
C FOUND IN GAMESS,GRD1,INPUTB,INPUTC,INT1 MODULES
3. install.com <machine-type> <size-keyword> Q <other-options>
The compile scripts are available for the following platforms:
T3E, T3D, IBMRS, IBM/SP, SUN, SGI, HP-UX, Convex SPP, DEC alpha,
PC/Linux with 3 compilers (Absoft, f2c, g77)
==========================================================
GAMESS/CHARMM interface status (July 1996) - Parallel version is fully functional - I/O including standard input and output are separated for GAMESS. - All CHARMM testcases are still OK when CHARMM is compiled with GAMESS inside. - GAMESS, CADPAC and QUANTUM keywords cannot coexist in pref.dat - MNDO, AM1, PM3 hamiltonians work only in pure QM calculations. No QM/MM energies and derivatives with this wavefunctions yet. - GAMESS recognizes atoms by their masses as specified in the RTF file
The following methods work with the GAMESS/CHARMM
(from GAMESS INTRO.DOC file)
I. A wide range of quantum chemical computations are
possible using GAMESS in the CHARMM MM field, which
1. Calculates RHF, UHF, ROHF, GVB, or MCSCF self-
consistent field molecular wavefunctions.
2. Calculates CI or MP2 corrections to the energy
of these SCF functions.
3. Calculates analytic energy gradients for all SCF
wavefunctions, plus closed shell MP2 or CI.
4. Optimizes molecular geometries using the energy
gradient, in terms of Cartesian or internal coords.
5. Searches for potential energy surface saddle points.
6. Computes the energy hessian, and thus normal modes,
vibrational frequencies, and IR intensities.
7. Traces the intrinsic reaction path from a saddle
point to reactants or products.
8. Traces gradient extremal curves, which may lead from
one stationary point such as a minimum to another,
which might be a saddle point.
9. Follows the dynamic reaction coordinate, a classical
mechanics trajectory on the potential energy surface.
10. Computes radiative transition probabilities.
11. Evaluates spin-orbit coupled wavefunctions.
12. Applies finite electric fields, extracting the
molecule's linear polarizability, and first and
second order hyperpolarizabilities.
13. Evaluates analytic frequency dependent non-linear
optical polarizability properties, for RHF functions.
14. Obtains localized orbitals by the Foster-Boys,
Edmiston-Ruedenberg, or Pipek-Mezey methods, with
optional SCF or MP2 energy analysis of the LMOs.
15. Calculates the following molecular properties:
a. dipole, quadrupole, and octupole moments
b. electrostatic potential
c. electric field and electric field gradients
d. electron density and spin density
e. Mulliken and Lowdin population analysis
f. virial theorem and energy components
g. Stone's distributed multipole analysis
16. Models solvent effects by
a. effective fragment potentials (EFP)
b. polarizable continuum model (PCM)
c. self-consistent reaction field (SCRF)
II. A quick summary of the current program capabilities
is given below.
SCFTYP= RHF ROHF UHF GVB MCSCF
--- ---- --- --- -----
Energy CDP CDP CDP CDP CDP
analytic gradient CDP CDP CDP CDP CDP
numerical Hessian CDP CDP CDP CDP CDP
analytic Hessian CDP CDP - CDP -
MP2 energy CDP CDP CDP - C
MP2 gradient CD - - - -
CI energy CDP CDP - CDP CDP
CI gradient CD - - - -
MOPAC energy yes yes yes yes -
MOPAC gradient yes yes yes - -
C= conventional storage of AO integrals on disk
D= direct evaluation of AO integrals
P= parallel execution
III. The methods listed above which don't have analytic gradients
are not available for CHARMM minimizations and dynamic
calculations.
IV. The following are available only in the pure QM calculations:
1. Calculates semi-empirical MNDO, AM1, or PM3
RHF, UHF, or ROHF wavefunctions.
Implementation
--------------
This is for version 22 NOV 2004 R1 of GAMESS or later:
The new DDI (Distributed Data Interface) library is used by parallel
GAMESS. This implementation in GAMESS uses 2 threads per CPU, and
client-server model which CANNOT be made compatible with CHARMM
parallel library. The new library was made from scratch and put in
source/gamint/ddi.src, which replaces the ddi.src file in GAMESS
distribution.
The following files need to be removed from gamint/gamess directory:
vector.src
The following files are modified from original GAMESS:
gamint/gamess/gamess.src
gamint/gamess/grd1.src
gamint/gamess/grd2a.src
gamint/gamess/guess.src
gamint/gamess/inputa.src
gamint/gamess/inputb.src
gamint/gamess/inputc.src
gamint/gamess/int1.src
gamint/gamess/iolib.src
gamint/gamess/mccas.src
gamint/gamess/prpel.src
gamint/gamess/prppop.src
gamint/gamess/prplib.src
gamint/gamess/rhfuhf.src
gamint/gamess/scflib.src
gamint/gamess/unport.c
gamint/gamess/zunix.c
Some changes for CHARMM interface into upstream version of the GAMESS
still didn't made it so contact milan@cmm.ki.si for more info.
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory
Modified, updated and generalized by C.L. Brooks, III
The Scripps Research Institute