ase_qmmm_manyqm

Introduction

This an general interface to run QM/MM calculations with a ase-QM and ase-MM calculator. Currently QM can be FHI-aims and MM can be Gromacs.

QM and MM region can be covalently bonded. In principle you can cut wherever you like, but other cutting schemes are probably better than others. The standard recommendation is to cut only the C-C bonds.

There can be many QM regions embedded in the MM system.

Ase_qmmm_manyqm Calculator

QM/MM interface with QM=FHI-aims, MM=gromacs

QM could be any QM calculator, but you need to read in QM-atom charges from the QM program (in method ‘get_qm_charges’).

One can have many QM regions, each with a different calculator. There can be only one MM calculator, which is calculating the whole system.

Non-bonded interactions:

Generally energies and forces are treated by:
  • Within the same QM-QM: by QM calculator
  • MM-MM: by MM calculator
  • QM-MM: by MM using MM vdw parameters and QM charges.
  • Different QM-different QM: by MM using QM and MM charges and MM-vdw parameters

The Hirschfeld charges (or other atomic charges) on QM atoms are calculated by QM in a H terminated cluster in vacuum. The charge of QM atom next to MM atom (edge-QM-atom) and its H neighbors are set as in the classical force field.

The extra(missing) charge results from:
  • linkH atoms
  • The edge-QM atoms, and their singly bonded neighbors (typically -H or =O). These have their original MM charges.
  • From the fact that the charge of the QM fraction is not usually an integer when using the original MM charges.
  • The extra/missing charge is added equally to all QM atoms (not being linkH and not being edge-QM-atom or its singly bonded neighbor) so that the total charge of the MM-fragment involving QM atoms will be the same as in the original MM-description.

Vdw interactions are calculated by MM-gromacs for MM and MM-QM interactions. The QM-QM vdw interaction s could be done by the FHI-aims if desired (by modifying the input for QM-FHI-aims input accordingly.

Bonded interactions:

E = E_qm(QM-H) + E_mm(ALL ATOMS), where
  • E_qm(QM-H): qm energy of H terminated QM cluster(s)
  • E_mm(ALL ATOMS): MM energy of all atoms, except for terms in which all MM-interacting atoms are in the same QM region.

Forces do not act on link atoms but they are positioned by scaling. Forces on link atoms are given to their QM and MM neighbors by chain rule. (see Top Curr Chem (2007) 268: 173–290, especially pages 192-194). At the beginning of a run, the optimal edge-qm-atom-linkH bond length(s) is (are) calculated by QM in ‘get_eq_qm_atom_link_h_distances’ or they are read from a file (this is determined by the argument ‘link_info’).

The covalent bonds between QM and MM are found automatically based on ase-covalent radii in neighbor search. Link atoms are positioned automatically (according to .

Questions & Comments markus.kaukonen@iki.fi

Interesting applications could be cases when we need two or more QM regions. For instance two redox centers in a protein, cathode and anode of a fuel cell … you name it!

Arguments

keyword type default value description
nqm_regions int   The number of QM regions
qm_calculators list   List of members of a Class defining a ase-QM calculator for each QM region
mm_calculator Calc   A member of a Class defining a ase-MM calculator
link_info str 'byQM' Can be either ‘byQM’: the edge_qm_atom-link_h_atom distances are calculated by QM or ‘byFile’:the edge_qm_atom-link_h_atom distances are read from a file

Example

  1. Prepare classical input with Gromacs
  • Get THE INDEPENDENT STRUCTURE OF THE ANTITRYPTIC REACTIVE SITE LOOP OF BOWMAN-BIRK INHIBITOR AND SUNFLOWER TRYPSIN INHIBITOR-1 (pdb code 1GM2) from pdb data bank http://www.rcsb.org/pdb/home/home.do, name it to 1GM2.pdb

  • In file 1GM2.pdb take only MODEL1 replace ‘CYS ‘ by ‘CYS2’ in order to get deprotonated CYS-CYS bridge.

  • Generate gromacs coordinates and topology

    >>> pdb2gmx -ff oplsaa -f 1GM2.pdb -o 1GM2.gro -p 1GM2.top -water tip3p -ignh
    
  • Generate the simulation box (cubic is not the most efficient one…)

    >>> editconf -bt cubic -f 1GM2.gro -o 1GM2_box.gro -c -d 0.2
    
  • Solvate the protein in water box

    >>> genbox -cp 1GM2_box.gro -cs spc216.gro -o 1GM2_sol.gro -p 1GM2.top
    
  • Generate index file (needed later also for defining QM atoms)

    >>> make_ndx -f 1GM2_sol.gro -o index.ndx
    >>> select 'q' <ENTER>
    
  • add to first lines to index.ndx (we define here 3 QM regions, indexing from 1):

    [qm_ss]
    18 19 20 21 139 140 141 142
    
    [qm_thr]
    28 29 30 31 32 33 34 35
    
    [qm_ser]
    64 65 66 67 68
    
  1. Relax MM system by fixing QM atoms in their (crystallographic?) positions (see documentation for Gromacs MM calculator).
"""An example for using gromacs calculator in ase.

Atom positions are relaxed.
If qm atoms are found in index file, they are kept fixed in the relaxation.
QM atoms are defined in index file (numbering from 1)
in sets containing QM or Qm or qm in their name.

A sample call::

   ./gromacs_example_mm_relax.py 1GM2_sol.gro
"""

from ase.calculators.gromacs import Gromacs

import sys

infile_name = sys.argv[1]

CALC_MM_RELAX = Gromacs(
    init_structure_file=infile_name,
    structure_file='gromacs_mm-relax.g96',
    force_field='oplsaa',
    water_model='tip3p',
    base_filename='gromacs_mm-relax',
    doing_qmmm=False, freeze_qm=True,
    index_filename='index.ndx',
    extra_mdrun_parameters=' -nt 1 ',
    define='-DFLEXIBLE',
    integrator='cg',
    nsteps='10000',
    nstfout='10',
    nstlog='10',
    nstenergy='10',
    nstlist='10',
    ns_type='grid',
    pbc='xyz',
    rlist='1.15',
    coulombtype='PME-Switch',
    rcoulomb='0.8',
    vdwtype='shift',
    rvdw='0.8',
    rvdw_switch='0.75',
    DispCorr='Ener')
CALC_MM_RELAX.generate_topology_and_g96file()
CALC_MM_RELAX.generate_gromacs_run_file()
CALC_MM_RELAX.run()
  1. QM/MM relaxation (all atoms are free)
""" demo run for ase_qmmm_manyqm calculator """

#  ./test_ase_qmmm_manyqm.py gromacs_mm-relax.g96

from ase.calculators.gromacs import Gromacs
from ase.calculators.aims import Aims
from ase.calculators.ase_qmmm_manyqm import AseQmmmManyqm
from ase.optimize import BFGS

import sys
from ase.io import read

RUN_COMMAND = '/home/mka/bin/aims.071711_6.serial.x'
SPECIES_DIR = '/home/mka/Programs/fhi-aims.071711_6/species_defaults/light/'

LOG_FILE = open("ase-qm-mm-output.log","w")
sys.stdout = LOG_FILE

infile_name = sys.argv[1]

CALC_QM1 = Aims(charge = 0,
                xc = 'pbe',
                sc_accuracy_etot = 1e-5,
                sc_accuracy_eev = 1e-2,
                sc_accuracy_rho = 1e-5,
                sc_accuracy_forces = 1e-3,
                species_dir = SPECIES_DIR,
                run_command = RUN_COMMAND)
CALC_QM1.set(output = 'hirshfeld')

CALC_QM2 = Aims(charge = 0,
                xc = 'pbe',
                sc_accuracy_etot = 1e-5,
                sc_accuracy_eev = 1e-2,
                sc_accuracy_rho = 1e-5,
                sc_accuracy_forces = 1e-3,
                species_dir = SPECIES_DIR,
                run_command = RUN_COMMAND)
CALC_QM2.set(output = 'hirshfeld')

CALC_QM3 = Aims(charge = 0,
                xc = 'pbe',
                sc_accuracy_etot = 1e-5,
                sc_accuracy_eev = 1e-2,
                sc_accuracy_rho = 1e-5,
                sc_accuracy_forces = 1e-3,
                species_dir = SPECIES_DIR,
                run_command = RUN_COMMAND)
CALC_QM3.set(output = 'hirshfeld')

CALC_MM = Gromacs(
    init_structure_file = infile_name,
    structure_file = 'gromacs_qm.g96', \
    force_field = 'oplsaa',
    water_model = 'tip3p',
    base_filename = 'gromacs_qm',
    doing_qmmm = True, freeze_qm = False,
    index_filename = 'index.ndx',
    define = '-DFLEXIBLE',
    integrator = 'md',
    nsteps = '0',
    nstfout = '1',
    nstlog = '1',
    nstenergy = '1',
    nstlist = '1',
    ns_type = 'grid',
    pbc = 'xyz',
    rlist = '1.15',
    coulombtype = 'PME-Switch',
    rcoulomb = '0.8',
    vdwtype = 'shift',
    rvdw = '0.8',
    rvdw_switch = '0.75',
    DispCorr = 'Ener')
CALC_MM.generate_topology_and_g96file()
CALC_MM.generate_gromacs_run_file()

CALC_QMMM = AseQmmmManyqm(nqm_regions = 3,
                          qm_calculators = [CALC_QM1, CALC_QM2, CALC_QM3],
                          mm_calculator = CALC_MM,
                          link_info = 'byQM')
#                         link_info = 'byFILE')

SYSTEM = read('gromacs_qm.g96')
SYSTEM.set_calculator(CALC_QMMM)
DYN = BFGS(SYSTEM)
DYN.run(fmax = 0.05)

print('exiting fine')
LOG_FILE.close()