The molecular dynamics simulation package Gromacs contains a QM/MM implementation. The quantum chemical calculations need to be performed by an external program, and interface with several quantum chemistry packages are included: Gamess UK, Gaussian, Mopac and Orca. Please feel invited to read the description on the Gromacs website.

DFTB is a parametrized quantum chemical method based on density-functional theory. Its newest version DFTB3 was developed in 2011 and has been successful in describing organic and biomolecular systems while being 100× to 1000× faster than usual DFT methods. For the description of DFTB, refer to the original papers on DFTB2 and DFTB3.

We have implemented the DFTB method in Gromacs version 5. This means that all of the calculations run within Gromacs and no external program is required for QM/MM simulations. The implementation is described in a paper that is in press now:
T. Kubař, K. Welke and G. Groenhof: New QM/MM Implementation of the DFTB3 method in the Gromacs Package, J Comput Chem DOI:10.1002/jcc.24029 (2015), link.

The DFTB method cannot describe the London dispersion interaction, however this shortcoming can be overcome by application of a post-SCF empirical contribution to the energy. Currently, the D3 dispersion correction is implemented; when using, refer to the respective papers (D3D3-BJ).

 

How to obtain and install the program

First, get Gromacs version 5.0 from the Gromacs website. link

Also, we use Gromacs in connection with Plumed, so you will need to get Plumed version 2.1 or 2.1.x as well. link

Install Plumed, but do not patch the Gromacs source yet.

Get the actual DFTB implementation in the form of a patch to Gromacs, new bugfix release from 30 Oct 2015: link
Please let us know via e-mail who and where you are when you download. Your data will not be used to send any unsolicited mail and will not be given to third parties. Thank you.

Patch the Gromacs source with the file just downloaded:

tar xvzf gromacs-5.0-dftb-v6-plumed.patch.tgz
patch -p0 < gromacs-5.0-dftb-v6a-plumed.patch

Manually correct the references to the libplumed.so or libplumed.a file in the Gromacs main directory, files Plumed.cmake and Plumed.inc (two occurrences in each file), to match your installation.

Compile Gromacs. Choose 'dftb' for GMX_QMMM_PROGRAM, and disable CUDA as well as all of the parallelization options (MPI, THREAD_MPI and OPENMP). It may be a good idea to compile in double precision (GMX_DOUBLE=ON). There may be issues when you link against MKL for linear algebra libraries; it you have solved these, please let us know.

 

How to obtain the parameter sets

The DFTB method involves several parameter sets; refer to the original publications to learn more. The parameter files are available free of charge after prior registration on www.dftb.org. For organic and biomolecular systems, we recommend using the DFTB3 method and the 3OB parameter set; note that the currently availabe chemical elements are H, C, N, O, P, S, F, Cl, Br, I, Na, K and Ca.

The parameter files deposited on the website cannot be read in directly by a C program like Gromacs. Thus, it is necessary to convert them into a suitable format. You may download a script to facilitate this task; note that the chemical elements need to be specified in this bash script.

In case you wish to employ the D3 dispersion correction in the DFTB calculations, download the necessary parameter files and uncompress them in the same directory where the DFTB parameter files reside.

 

How to run simulations

Since the original Gromacs QM/MM interface is used, you need to get familiar with it first: go to Gromacs website. Note that the topology has to be modified for a QM/MM simulation. Note that the QM zone has been entirely within one molecule (in the topology sense) in all applications so far.

Let us summarize the QM/MM-related options in the run input file (.mdp):

QMMM = yes  
QMMM-grps = QMsystem ; needs to be specified in the index file
QMMMscheme = normal  
QMmethod = RHF ; required but ignored
QMbasis = STO-3G ; required but ignored
QMcharge = 0 ; an integer
MMchargescalefactor = 1. ; or smaller

 

For a QM/MM run with DFTB, there are several additional options:

QMdftbsccmode = 3 ; for DFTB3, or 2 for DFTB2 (a.k.a. SCC-DFTB)
QMdftb-telec = 10. ; "electronic temperature" for the Fermi distribution
QMdftb-slko-path = /path/to/skf/ ; path to the parameter files including the trailing slash
QMdftb-slko-separator = ; possibly the character between element names
QMdftb-slko-lowercase   = yes ; or no, in the file names ("Ca" vs. "ca")
QMdftb-slko-suffix = -c.spl ; of the file names
QMdftb-partial-pme = 1 ; accelerated PME calculations - recommended (0 to switch off)
QMdftb-dispersion = 1 ; for D3, or 0 for no dispersion
QMdftb-cdko = 0 ; currently not implemented
QMdftb-mmhub-inf = 1 ; currently not implemented

 

Also, note the following comments on some general settings while doing QM/MM:

  • Cutoff-scheme=group is recommended as we do not know for sure how Verlet would co-operate with QM/MM.
  • In some simulations, wrong results were obtained in NPT simulations without the nstpcouple variable assigned explicitly. We recommend either setting nstpcouple=1 in the run input file, or resorting to NVT.

 

Some extended functionaly is also available and may be requested by setting an appropriate environment variable prior to issuing 'mdrun'. These are:

export GMX_DFTB_SLKO_PATH="/path/..." (overrides .mdp settings)
export GMX_DFTB_SLKO_LOWERCASE=y  (overrides .mdp settings)
export GMX_DFTB_SLKO_SEPARATOR="-" (overrides .mdp settings)
export GMX_DFTB_SLKO_SUFFIX=".skf" (overrides .mdp settings)
export GMX_DFTB_TELEC=300.  (overrides .mdp settings)
export GMX_DFTB_QM_COORD=n (coordinates of QM atoms will be written out every n steps)
export GMX_DFTB_MM_COORD=n (coordinates of MM atoms will be written out every n steps)
export GMX_DFTB_CUTOFF=1 (switched-force cut-off for QM/MM electrostatics instead of PME)
export GMX_DFTB_RFIELD=1 (reaction-field QM/MM electrostatics instead of PME)
export GMX_DFTB_SHIFT=1 (shifted-force cut-off to treat QM/MM electrostatics instead of PME)
export GMX_DFTB_SURFACE_CORRECTION=1 (dipole correction in PME instead of tin-foil boundary conditions)

 

Apart from the usual output, mdrun will write an additional file with Mulliken atom charges named qm_dftb_charges.xvg.

 

Version history

v6a: bugfix release (30 Oct 2015): link

v6: original release for Gromacs 5.0 + Plumed 2.1/2.1.x (26 Mar 2015): link