import numpy as np
from ase.calculators.calculator import Calculator
from ase.data import atomic_numbers
from ase.utils import IOContext
from ase.geometry import get_distances
from ase.cell import Cell
[docs]class SimpleQMMM(Calculator):
"""Simple QMMM calculator."""
implemented_properties = ['energy', 'forces']
def __init__(self, selection, qmcalc, mmcalc1, mmcalc2, vacuum=None):
"""SimpleQMMM object.
The energy is calculated as::
_ _ _
E = E (R ) - E (R ) + E (R )
QM QM MM QM MM all
parameters:
selection: list of int, slice object or list of bool
Selection out of all the atoms that belong to the QM part.
qmcalc: Calculator object
QM-calculator.
mmcalc1: Calculator object
MM-calculator used for QM region.
mmcalc2: Calculator object
MM-calculator used for everything.
vacuum: float or None
Amount of vacuum to add around QM atoms. Use None if QM
calculator doesn't need a box.
"""
self.selection = selection
self.qmcalc = qmcalc
self.mmcalc1 = mmcalc1
self.mmcalc2 = mmcalc2
self.vacuum = vacuum
self.qmatoms = None
self.center = None
self.name = '{0}-{1}+{1}'.format(qmcalc.name, mmcalc1.name)
Calculator.__init__(self)
def initialize_qm(self, atoms):
constraints = atoms.constraints
atoms.constraints = []
self.qmatoms = atoms[self.selection]
atoms.constraints = constraints
self.qmatoms.pbc = False
if self.vacuum:
self.qmatoms.center(vacuum=self.vacuum)
self.center = self.qmatoms.positions.mean(axis=0)
def calculate(self, atoms, properties, system_changes):
Calculator.calculate(self, atoms, properties, system_changes)
if self.qmatoms is None:
self.initialize_qm(atoms)
self.qmatoms.positions = atoms.positions[self.selection]
if self.vacuum:
self.qmatoms.positions += (self.center -
self.qmatoms.positions.mean(axis=0))
energy = self.qmcalc.get_potential_energy(self.qmatoms)
qmforces = self.qmcalc.get_forces(self.qmatoms)
energy += self.mmcalc2.get_potential_energy(atoms)
forces = self.mmcalc2.get_forces(atoms)
if self.vacuum:
qmforces -= qmforces.mean(axis=0)
forces[self.selection] += qmforces
energy -= self.mmcalc1.get_potential_energy(self.qmatoms)
forces[self.selection] -= self.mmcalc1.get_forces(self.qmatoms)
self.results['energy'] = energy
self.results['forces'] = forces
[docs]class EIQMMM(Calculator, IOContext):
"""Explicit interaction QMMM calculator."""
implemented_properties = ['energy', 'forces']
def __init__(self, selection, qmcalc, mmcalc, interaction,
vacuum=None, embedding=None, output=None):
"""EIQMMM object.
The energy is calculated as::
_ _ _ _
E = E (R ) + E (R ) + E (R , R )
QM QM MM MM I QM MM
parameters:
selection: list of int, slice object or list of bool
Selection out of all the atoms that belong to the QM part.
qmcalc: Calculator object
QM-calculator.
mmcalc: Calculator object
MM-calculator.
interaction: Interaction object
Interaction between QM and MM regions.
vacuum: float or None
Amount of vacuum to add around QM atoms. Use None if QM
calculator doesn't need a box.
embedding: Embedding object or None
Specialized embedding object. Use None in order to use the
default one.
output: None, '-', str or file-descriptor.
File for logging information - default is no logging (None).
"""
self.selection = selection
self.qmcalc = qmcalc
self.mmcalc = mmcalc
self.interaction = interaction
self.vacuum = vacuum
self.embedding = embedding
self.qmatoms = None
self.mmatoms = None
self.mask = None
self.center = None # center of QM atoms in QM-box
self.name = '{0}+{1}+{2}'.format(qmcalc.name,
interaction.name,
mmcalc.name)
self.output = self.openfile(output)
Calculator.__init__(self)
def initialize(self, atoms):
self.mask = np.zeros(len(atoms), bool)
self.mask[self.selection] = True
constraints = atoms.constraints
atoms.constraints = [] # avoid slicing of constraints
self.qmatoms = atoms[self.mask]
self.mmatoms = atoms[~self.mask]
atoms.constraints = constraints
self.qmatoms.pbc = False
if self.vacuum:
self.qmatoms.center(vacuum=self.vacuum)
self.center = self.qmatoms.positions.mean(axis=0)
print('Size of QM-cell after centering:',
self.qmatoms.cell.diagonal(), file=self.output)
self.qmatoms.calc = self.qmcalc
self.mmatoms.calc = self.mmcalc
if self.embedding is None:
self.embedding = Embedding()
self.embedding.initialize(self.qmatoms, self.mmatoms)
print('Embedding:', self.embedding, file=self.output)
def calculate(self, atoms, properties, system_changes):
Calculator.calculate(self, atoms, properties, system_changes)
if self.qmatoms is None:
self.initialize(atoms)
self.mmatoms.set_positions(atoms.positions[~self.mask])
self.qmatoms.set_positions(atoms.positions[self.mask])
if self.vacuum:
shift = self.center - self.qmatoms.positions.mean(axis=0)
self.qmatoms.positions += shift
else:
shift = (0, 0, 0)
self.embedding.update(shift)
ienergy, iqmforces, immforces = self.interaction.calculate(
self.qmatoms, self.mmatoms, shift)
qmenergy = self.qmatoms.get_potential_energy()
mmenergy = self.mmatoms.get_potential_energy()
energy = ienergy + qmenergy + mmenergy
print('Energies: {0:12.3f} {1:+12.3f} {2:+12.3f} = {3:12.3f}'
.format(ienergy, qmenergy, mmenergy, energy), file=self.output)
qmforces = self.qmatoms.get_forces()
mmforces = self.mmatoms.get_forces()
mmforces += self.embedding.get_mm_forces()
forces = np.empty((len(atoms), 3))
forces[self.mask] = qmforces + iqmforces
forces[~self.mask] = mmforces + immforces
self.results['energy'] = energy
self.results['forces'] = forces
def wrap(D, cell, pbc):
"""Wrap distances to nearest neighbor (minimum image convention)."""
for i, periodic in enumerate(pbc):
if periodic:
d = D[:, i]
L = cell[i]
d[:] = (d + L / 2) % L - L / 2 # modify D inplace
[docs]class Embedding:
def __init__(self, molecule_size=3, **parameters):
"""Point-charge embedding."""
self.qmatoms = None
self.mmatoms = None
self.molecule_size = molecule_size
self.virtual_molecule_size = None
self.parameters = parameters
def __repr__(self):
return 'Embedding(molecule_size={0})'.format(self.molecule_size)
def initialize(self, qmatoms, mmatoms):
"""Hook up embedding object to QM and MM atoms objects."""
self.qmatoms = qmatoms
self.mmatoms = mmatoms
charges = mmatoms.calc.get_virtual_charges(mmatoms)
self.pcpot = qmatoms.calc.embed(charges, **self.parameters)
self.virtual_molecule_size = (self.molecule_size *
len(charges) // len(mmatoms))
def update(self, shift):
"""Update point-charge positions."""
# Wrap point-charge positions to the MM-cell closest to the
# center of the the QM box, but avoid ripping molecules apart:
qmcenter = self.qmatoms.positions.mean(axis=0)
# if counter ions are used, then molecule_size has more than 1 value
if self.mmatoms.calc.name == 'combinemm':
mask1 = self.mmatoms.calc.mask
mask2 = ~mask1
vmask1 = self.mmatoms.calc.virtual_mask
vmask2 = ~vmask1
apm1 = self.mmatoms.calc.apm1
apm2 = self.mmatoms.calc.apm2
spm1 = self.mmatoms.calc.atoms1.calc.sites_per_mol
spm2 = self.mmatoms.calc.atoms2.calc.sites_per_mol
pos = self.mmatoms.positions
pos1 = pos[mask1].reshape((-1, apm1, 3))
pos2 = pos[mask2].reshape((-1, apm2, 3))
pos = (pos1, pos2)
else:
pos = (self.mmatoms.positions, )
apm1 = self.molecule_size
apm2 = self.molecule_size
# This is only specific to calculators where apm != spm,
# i.e. TIP4P. Non-native MM calcs do not have this attr.
if hasattr(self.mmatoms.calc, 'sites_per_mol'):
spm1 = self.mmatoms.calc.sites_per_mol
spm2 = self.mmatoms.calc.sites_per_mol
else:
spm1 = self.molecule_size
spm2 = spm1
mask1 = np.ones(len(self.mmatoms), dtype=bool)
mask2 = mask1
wrap_pos = np.zeros_like(self.mmatoms.positions)
com_all = []
apm = (apm1, apm2)
mask = (mask1, mask2)
spm = (spm1, spm2)
for p, n, m, vn in zip(pos, apm, mask, spm):
positions = p.reshape((-1, n, 3)) + shift
# Distances from the center of the QM box to the first atom of
# each molecule:
distances = positions[:, 0] - qmcenter
wrap(distances, self.mmatoms.cell.diagonal(), self.mmatoms.pbc)
offsets = distances - positions[:, 0]
positions += offsets[:, np.newaxis] + qmcenter
# Geometric center positions for each mm mol for LR cut
com = np.array([p.mean(axis=0) for p in positions])
# Need per atom for C-code:
com_pv = np.repeat(com, vn, axis=0)
com_all.append(com_pv)
wrap_pos[m] = positions.reshape((-1, 3))
positions = wrap_pos.copy()
positions = self.mmatoms.calc.add_virtual_sites(positions)
if self.mmatoms.calc.name == 'combinemm':
com_pv = np.zeros_like(positions)
for ii, m in enumerate((vmask1, vmask2)):
com_pv[m] = com_all[ii]
# compatibility with gpaw versions w/o LR cut in PointChargePotential
if 'rc2' in self.parameters:
self.pcpot.set_positions(positions, com_pv=com_pv)
else:
self.pcpot.set_positions(positions)
def get_mm_forces(self):
"""Calculate the forces on the MM-atoms from the QM-part."""
f = self.pcpot.get_forces(self.qmatoms.calc)
return self.mmatoms.calc.redistribute_forces(f)
def combine_lj_lorenz_berthelot(sigmaqm, sigmamm,
epsilonqm, epsilonmm):
"""Combine LJ parameters according to the Lorenz-Berthelot rule"""
sigma = []
epsilon = []
# check if input is tuple of vals for more than 1 mm calc, or only for 1.
if type(sigmamm) == tuple:
numcalcs = len(sigmamm)
else:
numcalcs = 1 # if only 1 mm calc, eps and sig are simply np arrays
sigmamm = (sigmamm, )
epsilonmm = (epsilonmm, )
for cc in range(numcalcs):
sigma_c = np.zeros((len(sigmaqm), len(sigmamm[cc])))
epsilon_c = np.zeros_like(sigma_c)
for ii in range(len(sigmaqm)):
sigma_c[ii, :] = (sigmaqm[ii] + sigmamm[cc]) / 2
epsilon_c[ii, :] = (epsilonqm[ii] * epsilonmm[cc])**0.5
sigma.append(sigma_c)
epsilon.append(epsilon_c)
if numcalcs == 1: # retain original, 1 calc function
sigma = np.array(sigma[0])
epsilon = np.array(epsilon[0])
return sigma, epsilon
[docs]class LJInteractionsGeneral:
name = 'LJ-general'
def __init__(self, sigmaqm, epsilonqm, sigmamm, epsilonmm,
qm_molecule_size, mm_molecule_size=3,
rc=np.Inf, width=1.0):
"""General Lennard-Jones type explicit interaction.
sigmaqm: array
Array of sigma-parameters which should have the length of the QM
subsystem
epsilonqm: array
As sigmaqm, but for epsilon-paramaters
sigmamm: Either array (A) or tuple (B)
A (no counterions):
Array of sigma-parameters with the length of the smallests
repeating atoms-group (i.e. molecule) of the MM subsystem
B (counterions):
Tuple: (arr1, arr2), where arr1 is an array of sigmas with
the length of counterions in the MM subsystem, and
arr2 is the array from A.
epsilonmm: array or tuple
As sigmamm but for epsilon-parameters.
qm_molecule_size: int
number of atoms of the smallest repeating atoms-group (i.e.
molecule) in the QM subsystem (often just the number of atoms
in the QM subsystem)
mm_molecule_size: int
as qm_molecule_size but for the MM subsystem. Will be overwritten
if counterions are present in the MM subsystem (via the CombineMM
calculator)
"""
self.sigmaqm = sigmaqm
self.epsilonqm = epsilonqm
self.sigmamm = sigmamm
self.epsilonmm = epsilonmm
self.qms = qm_molecule_size
self.mms = mm_molecule_size
self.rc = rc
self.width = width
self.combine_lj()
def combine_lj(self):
self.sigma, self.epsilon = combine_lj_lorenz_berthelot(
self.sigmaqm, self.sigmamm, self.epsilonqm, self.epsilonmm)
def calculate(self, qmatoms, mmatoms, shift):
epsilon = self.epsilon
sigma = self.sigma
# loop over possible multiple mm calculators
# currently 1 or 2, but could be generalized in the future...
apm1 = self.mms
mask1 = np.ones(len(mmatoms), dtype=bool)
mask2 = mask1
apm = (apm1, )
sigma = (sigma, )
epsilon = (epsilon, )
if hasattr(mmatoms.calc, 'name'):
if mmatoms.calc.name == 'combinemm':
mask1 = mmatoms.calc.mask
mask2 = ~mask1
apm1 = mmatoms.calc.apm1
apm2 = mmatoms.calc.apm2
apm = (apm1, apm2)
sigma = sigma[0] # Was already loopable 2-tuple
epsilon = epsilon[0]
mask = (mask1, mask2)
e_all = 0
qmforces_all = np.zeros_like(qmatoms.positions)
mmforces_all = np.zeros_like(mmatoms.positions)
# zip stops at shortest tuple so we dont double count
# cases of no counter ions.
for n, m, eps, sig in zip(apm, mask, epsilon, sigma):
mmpositions = self.update(qmatoms, mmatoms[m], n, shift)
qmforces = np.zeros_like(qmatoms.positions)
mmforces = np.zeros_like(mmatoms[m].positions)
energy = 0.0
qmpositions = qmatoms.positions.reshape((-1, self.qms, 3))
for q, qmpos in enumerate(qmpositions): # molwise loop
# cutoff from first atom of each mol
R00 = mmpositions[:, 0] - qmpos[0, :]
d002 = (R00**2).sum(1)
d00 = d002**0.5
x1 = d00 > self.rc - self.width
x2 = d00 < self.rc
x12 = np.logical_and(x1, x2)
y = (d00[x12] - self.rc + self.width) / self.width
t = np.zeros(len(d00))
t[x2] = 1.0
t[x12] -= y**2 * (3.0 - 2.0 * y)
dt = np.zeros(len(d00))
dt[x12] -= 6.0 / self.width * y * (1.0 - y)
for qa in range(len(qmpos)):
if ~np.any(eps[qa, :]):
continue
R = mmpositions - qmpos[qa, :]
d2 = (R**2).sum(2)
c6 = (sig[qa, :]**2 / d2)**3
c12 = c6**2
e = 4 * eps[qa, :] * (c12 - c6)
energy += np.dot(e.sum(1), t)
f = t[:, None, None] * (24 * eps[qa, :] *
(2 * c12 - c6) / d2)[:, :, None] * R
f00 = - (e.sum(1) * dt / d00)[:, None] * R00
mmforces += f.reshape((-1, 3))
qmforces[q * self.qms + qa, :] -= f.sum(0).sum(0)
qmforces[q * self.qms, :] -= f00.sum(0)
mmforces[::n, :] += f00
e_all += energy
qmforces_all += qmforces
mmforces_all[m] += mmforces
return e_all, qmforces_all, mmforces_all
def update(self, qmatoms, mmatoms, n, shift):
# Wrap point-charge positions to the MM-cell closest to the
# center of the the QM box, but avoid ripping molecules apart:
qmcenter = qmatoms.cell.diagonal() / 2
positions = mmatoms.positions.reshape((-1, n, 3)) + shift
# Distances from the center of the QM box to the first atom of
# each molecule:
distances = positions[:, 0] - qmcenter
wrap(distances, mmatoms.cell.diagonal(), mmatoms.pbc)
offsets = distances - positions[:, 0]
positions += offsets[:, np.newaxis] + qmcenter
return positions
[docs]class LJInteractions:
name = 'LJ'
def __init__(self, parameters):
"""Lennard-Jones type explicit interaction.
parameters: dict
Mapping from pair of atoms to tuple containing epsilon and sigma
for that pair.
Example:
lj = LJInteractions({('O', 'O'): (eps, sigma)})
"""
self.parameters = {}
for (symbol1, symbol2), (epsilon, sigma) in parameters.items():
Z1 = atomic_numbers[symbol1]
Z2 = atomic_numbers[symbol2]
self.parameters[(Z1, Z2)] = epsilon, sigma
self.parameters[(Z2, Z1)] = epsilon, sigma
def calculate(self, qmatoms, mmatoms, shift):
qmforces = np.zeros_like(qmatoms.positions)
mmforces = np.zeros_like(mmatoms.positions)
species = set(mmatoms.numbers)
energy = 0.0
for R1, Z1, F1 in zip(qmatoms.positions, qmatoms.numbers, qmforces):
for Z2 in species:
if (Z1, Z2) not in self.parameters:
continue
epsilon, sigma = self.parameters[(Z1, Z2)]
mask = (mmatoms.numbers == Z2)
D = mmatoms.positions[mask] + shift - R1
wrap(D, mmatoms.cell.diagonal(), mmatoms.pbc)
d2 = (D**2).sum(1)
c6 = (sigma**2 / d2)**3
c12 = c6**2
energy += 4 * epsilon * (c12 - c6).sum()
f = 24 * epsilon * ((2 * c12 - c6) / d2)[:, np.newaxis] * D
F1 -= f.sum(0)
mmforces[mask] += f
return energy, qmforces, mmforces
class RescaledCalculator(Calculator):
"""Rescales length and energy of a calculators to match given
lattice constant and bulk modulus
Useful for MM calculator used within a :class:`ForceQMMM` model.
See T. D. Swinburne and J. R. Kermode, Phys. Rev. B 96, 144102 (2017)
for a derivation of the scaling constants.
"""
implemented_properties = ['forces', 'energy', 'stress']
def __init__(self, mm_calc,
qm_lattice_constant, qm_bulk_modulus,
mm_lattice_constant, mm_bulk_modulus):
Calculator.__init__(self)
self.mm_calc = mm_calc
self.alpha = qm_lattice_constant / mm_lattice_constant
self.beta = mm_bulk_modulus / qm_bulk_modulus / (self.alpha**3)
def calculate(self, atoms, properties, system_changes):
Calculator.calculate(self, atoms, properties, system_changes)
# mm_pos = atoms.get_positions()
scaled_atoms = atoms.copy()
# scaled_atoms.positions = mm_pos/self.alpha
mm_cell = atoms.get_cell()
scaled_atoms.set_cell(mm_cell / self.alpha, scale_atoms=True)
results = {}
if 'energy' in properties:
energy = self.mm_calc.get_potential_energy(scaled_atoms)
results['energy'] = energy / self.beta
if 'forces' in properties:
forces = self.mm_calc.get_forces(scaled_atoms)
results['forces'] = forces / (self.beta * self.alpha)
if 'stress' in properties:
stress = self.mm_calc.get_stress(scaled_atoms)
results['stress'] = stress / (self.beta * self.alpha**3)
self.results = results
class ForceConstantCalculator(Calculator):
"""
Compute forces based on provided force-constant matrix
Useful with `ForceQMMM` to do harmonic QM/MM using force constants
of QM method.
"""
implemented_properties = ['forces', 'energy']
def __init__(self, D, ref, f0):
"""
Parameters:
D: matrix or sparse matrix, shape `(3*len(ref), 3*len(ref))`
Force constant matrix.
Sign convention is `D_ij = d^2E/(dx_i dx_j), so
`force = -D.dot(displacement)`
ref: ase.atoms.Atoms
Atoms object for reference configuration
f0: array, shape `(len(ref), 3)`
Value of forces at reference configuration
"""
assert D.shape[0] == D.shape[1]
assert D.shape[0] // 3 == len(ref)
self.D = D
self.ref = ref
self.f0 = f0
self.size = len(ref)
Calculator.__init__(self)
def calculate(self, atoms, properties, system_changes):
Calculator.calculate(self, atoms, properties, system_changes)
u = atoms.positions - self.ref.positions
f = -self.D.dot(u.reshape(3 * self.size))
forces = np.zeros((len(atoms), 3))
forces[:, :] = f.reshape(self.size, 3)
self.results['forces'] = forces + self.f0
self.results['energy'] = 0.0
[docs]class ForceQMMM(Calculator):
"""
Force-based QM/MM calculator
QM forces are computed using a buffer region and then mixed abruptly
with MM forces:
F^i_QMMM = { F^i_QM if i in QM region
{ F^i_MM otherwise
cf. N. Bernstein, J. R. Kermode, and G. Csanyi,
Rep. Prog. Phys. 72, 026501 (2009)
and T. D. Swinburne and J. R. Kermode, Phys. Rev. B 96, 144102 (2017).
"""
implemented_properties = ['forces', 'energy']
def __init__(self,
atoms,
qm_selection_mask,
qm_calc,
mm_calc,
buffer_width,
vacuum=5.,
zero_mean=True,
qm_cell_round_off=3,
qm_radius=None):
"""
ForceQMMM calculator
Parameters:
qm_selection_mask: list of ints, slice object or bool list/array
Selection out of atoms that belong to the QM region.
qm_calc: Calculator object
QM-calculator.
mm_calc: Calculator object
MM-calculator (should be scaled, see :class:`RescaledCalculator`)
Can use `ForceConstantCalculator` based on QM force constants, if
available.
vacuum: float or None
Amount of vacuum to add around QM atoms.
zero_mean: bool
If True, add a correction to zero the mean force in each direction
qm_cell_round_off: float
Tolerance value in Angstrom to round the qm cluster cell
qm_radius: 3x1 array of floats qm_radius for [x, y, z]
3d qm radius for calculation of qm cluster cell. default is None
and the radius is estimated from maximum distance between the atoms
in qm region.
"""
if len(atoms[qm_selection_mask]) == 0:
raise ValueError("no QM atoms selected!")
self.qm_selection_mask = qm_selection_mask
self.qm_calc = qm_calc
self.mm_calc = mm_calc
self.vacuum = vacuum
self.buffer_width = buffer_width
self.zero_mean = zero_mean
self.qm_cell_round_off = qm_cell_round_off
self.qm_radius = qm_radius
self.qm_buffer_mask = None
Calculator.__init__(self)
def initialize_qm_buffer_mask(self, atoms):
"""
Initialises system to perform qm calculation
"""
# calculate the distances between all atoms and qm atoms
# qm_distance_matrix is a [N_QM_atoms x N_atoms] matrix
_, qm_distance_matrix = get_distances(
atoms.positions[self.qm_selection_mask],
atoms.positions,
atoms.cell, atoms.pbc)
self.qm_buffer_mask = np.zeros(len(atoms), dtype=bool)
# every r_qm is a matrix of distances
# from an atom in qm region and all atoms with size [N_atoms]
for r_qm in qm_distance_matrix:
self.qm_buffer_mask[r_qm < self.buffer_width] = True
def get_qm_cluster(self, atoms):
if self.qm_buffer_mask is None:
self.initialize_qm_buffer_mask(atoms)
qm_cluster = atoms[self.qm_buffer_mask]
del qm_cluster.constraints
round_cell = False
if self.qm_radius is None:
round_cell = True
# get all distances between qm atoms.
# Treat all X, Y and Z directions independently
# only distance between qm atoms is calculated
# in order to estimate qm radius in thee directions
R_qm, _ = get_distances(atoms.positions[self.qm_selection_mask],
cell=atoms.cell, pbc=atoms.pbc)
# estimate qm radius in three directions as 1/2
# of max distance between qm atoms
self.qm_radius = np.amax(np.amax(R_qm, axis=1), axis=0) * 0.5
if atoms.cell.orthorhombic:
cell_size = np.diagonal(atoms.cell)
else:
raise RuntimeError("NON-orthorhombic cell is not supported!")
# check if qm_cluster should be left periodic
# in periodic directions of the cell (cell[i] < qm_radius + buffer
# otherwise change to non pbc
# and make a cluster in a vacuum configuration
qm_cluster_pbc = (atoms.pbc &
(cell_size <
2.0 * (self.qm_radius + self.buffer_width)))
# start with the original orthorhombic cell
qm_cluster_cell = cell_size.copy()
# create a cluster in a vacuum cell in non periodic directions
qm_cluster_cell[~qm_cluster_pbc] = (
2.0 * (self.qm_radius[~qm_cluster_pbc] +
self.buffer_width +
self.vacuum))
if round_cell:
# round the qm cell to the required tolerance
qm_cluster_cell[~qm_cluster_pbc] = (np.round(
(qm_cluster_cell[~qm_cluster_pbc]) /
self.qm_cell_round_off) *
self.qm_cell_round_off)
qm_cluster.set_cell(Cell(np.diag(qm_cluster_cell)))
qm_cluster.pbc = qm_cluster_pbc
qm_shift = (0.5 * qm_cluster.cell.diagonal() -
qm_cluster.positions.mean(axis=0))
if 'cell_origin' in qm_cluster.info:
del qm_cluster.info['cell_origin']
# center the cluster only in non pbc directions
qm_cluster.positions[:, ~qm_cluster_pbc] += qm_shift[~qm_cluster_pbc]
return qm_cluster
def calculate(self, atoms, properties, system_changes):
Calculator.calculate(self, atoms, properties, system_changes)
qm_cluster = self.get_qm_cluster(atoms)
forces = self.mm_calc.get_forces(atoms)
qm_forces = self.qm_calc.get_forces(qm_cluster)
forces[self.qm_selection_mask] = \
qm_forces[self.qm_selection_mask[self.qm_buffer_mask]]
if self.zero_mean:
# Target is that: forces.sum(axis=1) == [0., 0., 0.]
forces[:] -= forces.mean(axis=0)
self.results['forces'] = forces
self.results['energy'] = 0.0
def get_region_from_masks(self, atoms=None, print_mapping=False):
"""
creates region array from the masks of the calculators. The tags in
the array are:
QM - qm atoms
buffer - buffer atoms
MM - atoms treated with mm calculator
"""
if atoms is None:
if self.atoms is None:
raise ValueError('Calculator has no atoms')
else:
atoms = self.atoms
region = np.full_like(atoms, "MM")
region[self.qm_selection_mask] = (
np.full_like(region[self.qm_selection_mask], "QM"))
buffer_only_mask = self.qm_buffer_mask & ~self.qm_selection_mask
region[buffer_only_mask] = np.full_like(region[buffer_only_mask],
"buffer")
if print_mapping:
print(f"Mapping of {len(region):5d} atoms in total:")
for region_id in np.unique(region):
n_at = np.count_nonzero(region == region_id)
print(f"{n_at:16d} {region_id}")
qm_atoms = atoms[self.qm_selection_mask]
symbol_counts = qm_atoms.symbols.formula.count()
print("QM atoms types:")
for symbol, count in symbol_counts.items():
print(f"{count:16d} {symbol}")
return region
def set_masks_from_region(self, region):
"""
Sets masks from provided region array
"""
self.qm_selection_mask = region == "QM"
buffer_mask = region == "buffer"
self.qm_buffer_mask = self.qm_selection_mask ^ buffer_mask
def export_extxyz(self, atoms=None, filename="qmmm_atoms.xyz"):
"""
exports the atoms to extended xyz file with additional "region"
array keeping the mapping between QM, buffer and MM parts of
the simulation
"""
if atoms is None:
if self.atoms is None:
raise ValueError('Calculator has no atoms')
else:
atoms = self.atoms
region = self.get_region_from_masks(atoms=atoms)
atoms_copy = atoms.copy()
atoms_copy.new_array("region", region)
atoms_copy.calc = self # to keep the calculation results
atoms_copy.write(filename, format='extxyz')
@classmethod
def import_extxyz(cls, filename, qm_calc, mm_calc):
"""
A static method to import the the mapping from an estxyz file saved by
export_extxyz() function
Parameters
----------
filename: string
filename with saved configuration
qm_calc: Calculator object
QM-calculator.
mm_calc: Calculator object
MM-calculator (should be scaled, see :class:`RescaledCalculator`)
Can use `ForceConstantCalculator` based on QM force constants, if
available.
Returns
-------
New object of ForceQMMM calculator with qm_selection_mask and
qm_buffer_mask set according to the region array in the saved file
"""
from ase.io import read
atoms = read(filename, format='extxyz')
if "region" in atoms.arrays:
region = atoms.get_array("region")
else:
raise RuntimeError("Please provide extxyz file with 'region' array")
dummy_qm_mask = np.full_like(atoms, False, dtype=bool)
dummy_qm_mask[0] = True
self = cls(atoms, dummy_qm_mask,
qm_calc, mm_calc, buffer_width=1.0)
self.set_masks_from_region(region)
return self