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.
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.
- 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.
- 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 firstname.lastname@example.org
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!
||The number of QM regions|
||List of members of a Class defining a ase-QM calculator for each QM region|
||A member of a Class defining a ase-MM calculator|
||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|
- 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
- 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 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()
- 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 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()