Coverage for /builds/debichem-team/python-ase/ase/calculators/qmmm.py : 90.91%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import numpy as np
3from ase.calculators.calculator import Calculator
4from ase.data import atomic_numbers
5from ase.utils import IOContext
6from ase.geometry import get_distances
7from ase.cell import Cell
10class SimpleQMMM(Calculator):
11 """Simple QMMM calculator."""
13 implemented_properties = ['energy', 'forces']
15 def __init__(self, selection, qmcalc, mmcalc1, mmcalc2, vacuum=None):
16 """SimpleQMMM object.
18 The energy is calculated as::
20 _ _ _
21 E = E (R ) - E (R ) + E (R )
22 QM QM MM QM MM all
24 parameters:
26 selection: list of int, slice object or list of bool
27 Selection out of all the atoms that belong to the QM part.
28 qmcalc: Calculator object
29 QM-calculator.
30 mmcalc1: Calculator object
31 MM-calculator used for QM region.
32 mmcalc2: Calculator object
33 MM-calculator used for everything.
34 vacuum: float or None
35 Amount of vacuum to add around QM atoms. Use None if QM
36 calculator doesn't need a box.
38 """
39 self.selection = selection
40 self.qmcalc = qmcalc
41 self.mmcalc1 = mmcalc1
42 self.mmcalc2 = mmcalc2
43 self.vacuum = vacuum
45 self.qmatoms = None
46 self.center = None
48 self.name = '{0}-{1}+{1}'.format(qmcalc.name, mmcalc1.name)
50 Calculator.__init__(self)
52 def initialize_qm(self, atoms):
53 constraints = atoms.constraints
54 atoms.constraints = []
55 self.qmatoms = atoms[self.selection]
56 atoms.constraints = constraints
57 self.qmatoms.pbc = False
58 if self.vacuum:
59 self.qmatoms.center(vacuum=self.vacuum)
60 self.center = self.qmatoms.positions.mean(axis=0)
62 def calculate(self, atoms, properties, system_changes):
63 Calculator.calculate(self, atoms, properties, system_changes)
65 if self.qmatoms is None:
66 self.initialize_qm(atoms)
68 self.qmatoms.positions = atoms.positions[self.selection]
69 if self.vacuum:
70 self.qmatoms.positions += (self.center -
71 self.qmatoms.positions.mean(axis=0))
73 energy = self.qmcalc.get_potential_energy(self.qmatoms)
74 qmforces = self.qmcalc.get_forces(self.qmatoms)
75 energy += self.mmcalc2.get_potential_energy(atoms)
76 forces = self.mmcalc2.get_forces(atoms)
78 if self.vacuum:
79 qmforces -= qmforces.mean(axis=0)
80 forces[self.selection] += qmforces
82 energy -= self.mmcalc1.get_potential_energy(self.qmatoms)
83 forces[self.selection] -= self.mmcalc1.get_forces(self.qmatoms)
85 self.results['energy'] = energy
86 self.results['forces'] = forces
89class EIQMMM(Calculator, IOContext):
90 """Explicit interaction QMMM calculator."""
91 implemented_properties = ['energy', 'forces']
93 def __init__(self, selection, qmcalc, mmcalc, interaction,
94 vacuum=None, embedding=None, output=None):
95 """EIQMMM object.
97 The energy is calculated as::
99 _ _ _ _
100 E = E (R ) + E (R ) + E (R , R )
101 QM QM MM MM I QM MM
103 parameters:
105 selection: list of int, slice object or list of bool
106 Selection out of all the atoms that belong to the QM part.
107 qmcalc: Calculator object
108 QM-calculator.
109 mmcalc: Calculator object
110 MM-calculator.
111 interaction: Interaction object
112 Interaction between QM and MM regions.
113 vacuum: float or None
114 Amount of vacuum to add around QM atoms. Use None if QM
115 calculator doesn't need a box.
116 embedding: Embedding object or None
117 Specialized embedding object. Use None in order to use the
118 default one.
119 output: None, '-', str or file-descriptor.
120 File for logging information - default is no logging (None).
122 """
124 self.selection = selection
126 self.qmcalc = qmcalc
127 self.mmcalc = mmcalc
128 self.interaction = interaction
129 self.vacuum = vacuum
130 self.embedding = embedding
132 self.qmatoms = None
133 self.mmatoms = None
134 self.mask = None
135 self.center = None # center of QM atoms in QM-box
137 self.name = '{0}+{1}+{2}'.format(qmcalc.name,
138 interaction.name,
139 mmcalc.name)
141 self.output = self.openfile(output)
143 Calculator.__init__(self)
145 def initialize(self, atoms):
146 self.mask = np.zeros(len(atoms), bool)
147 self.mask[self.selection] = True
149 constraints = atoms.constraints
150 atoms.constraints = [] # avoid slicing of constraints
151 self.qmatoms = atoms[self.mask]
152 self.mmatoms = atoms[~self.mask]
153 atoms.constraints = constraints
155 self.qmatoms.pbc = False
157 if self.vacuum:
158 self.qmatoms.center(vacuum=self.vacuum)
159 self.center = self.qmatoms.positions.mean(axis=0)
160 print('Size of QM-cell after centering:',
161 self.qmatoms.cell.diagonal(), file=self.output)
163 self.qmatoms.calc = self.qmcalc
164 self.mmatoms.calc = self.mmcalc
166 if self.embedding is None:
167 self.embedding = Embedding()
169 self.embedding.initialize(self.qmatoms, self.mmatoms)
170 print('Embedding:', self.embedding, file=self.output)
172 def calculate(self, atoms, properties, system_changes):
173 Calculator.calculate(self, atoms, properties, system_changes)
175 if self.qmatoms is None:
176 self.initialize(atoms)
178 self.mmatoms.set_positions(atoms.positions[~self.mask])
179 self.qmatoms.set_positions(atoms.positions[self.mask])
181 if self.vacuum:
182 shift = self.center - self.qmatoms.positions.mean(axis=0)
183 self.qmatoms.positions += shift
184 else:
185 shift = (0, 0, 0)
187 self.embedding.update(shift)
189 ienergy, iqmforces, immforces = self.interaction.calculate(
190 self.qmatoms, self.mmatoms, shift)
192 qmenergy = self.qmatoms.get_potential_energy()
193 mmenergy = self.mmatoms.get_potential_energy()
194 energy = ienergy + qmenergy + mmenergy
196 print('Energies: {0:12.3f} {1:+12.3f} {2:+12.3f} = {3:12.3f}'
197 .format(ienergy, qmenergy, mmenergy, energy), file=self.output)
199 qmforces = self.qmatoms.get_forces()
200 mmforces = self.mmatoms.get_forces()
202 mmforces += self.embedding.get_mm_forces()
204 forces = np.empty((len(atoms), 3))
205 forces[self.mask] = qmforces + iqmforces
206 forces[~self.mask] = mmforces + immforces
208 self.results['energy'] = energy
209 self.results['forces'] = forces
212def wrap(D, cell, pbc):
213 """Wrap distances to nearest neighbor (minimum image convention)."""
214 for i, periodic in enumerate(pbc):
215 if periodic:
216 d = D[:, i]
217 L = cell[i]
218 d[:] = (d + L / 2) % L - L / 2 # modify D inplace
221class Embedding:
222 def __init__(self, molecule_size=3, **parameters):
223 """Point-charge embedding."""
224 self.qmatoms = None
225 self.mmatoms = None
226 self.molecule_size = molecule_size
227 self.virtual_molecule_size = None
228 self.parameters = parameters
230 def __repr__(self):
231 return 'Embedding(molecule_size={0})'.format(self.molecule_size)
233 def initialize(self, qmatoms, mmatoms):
234 """Hook up embedding object to QM and MM atoms objects."""
235 self.qmatoms = qmatoms
236 self.mmatoms = mmatoms
237 charges = mmatoms.calc.get_virtual_charges(mmatoms)
238 self.pcpot = qmatoms.calc.embed(charges, **self.parameters)
239 self.virtual_molecule_size = (self.molecule_size *
240 len(charges) // len(mmatoms))
242 def update(self, shift):
243 """Update point-charge positions."""
244 # Wrap point-charge positions to the MM-cell closest to the
245 # center of the the QM box, but avoid ripping molecules apart:
246 qmcenter = self.qmatoms.positions.mean(axis=0)
247 # if counter ions are used, then molecule_size has more than 1 value
248 if self.mmatoms.calc.name == 'combinemm':
249 mask1 = self.mmatoms.calc.mask
250 mask2 = ~mask1
251 vmask1 = self.mmatoms.calc.virtual_mask
252 vmask2 = ~vmask1
253 apm1 = self.mmatoms.calc.apm1
254 apm2 = self.mmatoms.calc.apm2
255 spm1 = self.mmatoms.calc.atoms1.calc.sites_per_mol
256 spm2 = self.mmatoms.calc.atoms2.calc.sites_per_mol
257 pos = self.mmatoms.positions
258 pos1 = pos[mask1].reshape((-1, apm1, 3))
259 pos2 = pos[mask2].reshape((-1, apm2, 3))
260 pos = (pos1, pos2)
261 else:
262 pos = (self.mmatoms.positions, )
263 apm1 = self.molecule_size
264 apm2 = self.molecule_size
265 # This is only specific to calculators where apm != spm,
266 # i.e. TIP4P. Non-native MM calcs do not have this attr.
267 if hasattr(self.mmatoms.calc, 'sites_per_mol'):
268 spm1 = self.mmatoms.calc.sites_per_mol
269 spm2 = self.mmatoms.calc.sites_per_mol
270 else:
271 spm1 = self.molecule_size
272 spm2 = spm1
273 mask1 = np.ones(len(self.mmatoms), dtype=bool)
274 mask2 = mask1
276 wrap_pos = np.zeros_like(self.mmatoms.positions)
277 com_all = []
278 apm = (apm1, apm2)
279 mask = (mask1, mask2)
280 spm = (spm1, spm2)
281 for p, n, m, vn in zip(pos, apm, mask, spm):
282 positions = p.reshape((-1, n, 3)) + shift
284 # Distances from the center of the QM box to the first atom of
285 # each molecule:
286 distances = positions[:, 0] - qmcenter
288 wrap(distances, self.mmatoms.cell.diagonal(), self.mmatoms.pbc)
289 offsets = distances - positions[:, 0]
290 positions += offsets[:, np.newaxis] + qmcenter
292 # Geometric center positions for each mm mol for LR cut
293 com = np.array([p.mean(axis=0) for p in positions])
294 # Need per atom for C-code:
295 com_pv = np.repeat(com, vn, axis=0)
296 com_all.append(com_pv)
298 wrap_pos[m] = positions.reshape((-1, 3))
300 positions = wrap_pos.copy()
301 positions = self.mmatoms.calc.add_virtual_sites(positions)
303 if self.mmatoms.calc.name == 'combinemm':
304 com_pv = np.zeros_like(positions)
305 for ii, m in enumerate((vmask1, vmask2)):
306 com_pv[m] = com_all[ii]
308 # compatibility with gpaw versions w/o LR cut in PointChargePotential
309 if 'rc2' in self.parameters:
310 self.pcpot.set_positions(positions, com_pv=com_pv)
311 else:
312 self.pcpot.set_positions(positions)
314 def get_mm_forces(self):
315 """Calculate the forces on the MM-atoms from the QM-part."""
316 f = self.pcpot.get_forces(self.qmatoms.calc)
317 return self.mmatoms.calc.redistribute_forces(f)
320def combine_lj_lorenz_berthelot(sigmaqm, sigmamm,
321 epsilonqm, epsilonmm):
322 """Combine LJ parameters according to the Lorenz-Berthelot rule"""
323 sigma = []
324 epsilon = []
325 # check if input is tuple of vals for more than 1 mm calc, or only for 1.
326 if type(sigmamm) == tuple:
327 numcalcs = len(sigmamm)
328 else:
329 numcalcs = 1 # if only 1 mm calc, eps and sig are simply np arrays
330 sigmamm = (sigmamm, )
331 epsilonmm = (epsilonmm, )
332 for cc in range(numcalcs):
333 sigma_c = np.zeros((len(sigmaqm), len(sigmamm[cc])))
334 epsilon_c = np.zeros_like(sigma_c)
336 for ii in range(len(sigmaqm)):
337 sigma_c[ii, :] = (sigmaqm[ii] + sigmamm[cc]) / 2
338 epsilon_c[ii, :] = (epsilonqm[ii] * epsilonmm[cc])**0.5
339 sigma.append(sigma_c)
340 epsilon.append(epsilon_c)
342 if numcalcs == 1: # retain original, 1 calc function
343 sigma = np.array(sigma[0])
344 epsilon = np.array(epsilon[0])
346 return sigma, epsilon
349class LJInteractionsGeneral:
350 name = 'LJ-general'
352 def __init__(self, sigmaqm, epsilonqm, sigmamm, epsilonmm,
353 qm_molecule_size, mm_molecule_size=3,
354 rc=np.Inf, width=1.0):
355 """General Lennard-Jones type explicit interaction.
357 sigmaqm: array
358 Array of sigma-parameters which should have the length of the QM
359 subsystem
360 epsilonqm: array
361 As sigmaqm, but for epsilon-paramaters
362 sigmamm: Either array (A) or tuple (B)
363 A (no counterions):
364 Array of sigma-parameters with the length of the smallests
365 repeating atoms-group (i.e. molecule) of the MM subsystem
366 B (counterions):
367 Tuple: (arr1, arr2), where arr1 is an array of sigmas with
368 the length of counterions in the MM subsystem, and
369 arr2 is the array from A.
370 epsilonmm: array or tuple
371 As sigmamm but for epsilon-parameters.
372 qm_molecule_size: int
373 number of atoms of the smallest repeating atoms-group (i.e.
374 molecule) in the QM subsystem (often just the number of atoms
375 in the QM subsystem)
376 mm_molecule_size: int
377 as qm_molecule_size but for the MM subsystem. Will be overwritten
378 if counterions are present in the MM subsystem (via the CombineMM
379 calculator)
381 """
382 self.sigmaqm = sigmaqm
383 self.epsilonqm = epsilonqm
384 self.sigmamm = sigmamm
385 self.epsilonmm = epsilonmm
386 self.qms = qm_molecule_size
387 self.mms = mm_molecule_size
388 self.rc = rc
389 self.width = width
390 self.combine_lj()
392 def combine_lj(self):
393 self.sigma, self.epsilon = combine_lj_lorenz_berthelot(
394 self.sigmaqm, self.sigmamm, self.epsilonqm, self.epsilonmm)
396 def calculate(self, qmatoms, mmatoms, shift):
397 epsilon = self.epsilon
398 sigma = self.sigma
400 # loop over possible multiple mm calculators
401 # currently 1 or 2, but could be generalized in the future...
402 apm1 = self.mms
403 mask1 = np.ones(len(mmatoms), dtype=bool)
404 mask2 = mask1
405 apm = (apm1, )
406 sigma = (sigma, )
407 epsilon = (epsilon, )
408 if hasattr(mmatoms.calc, 'name'):
409 if mmatoms.calc.name == 'combinemm':
410 mask1 = mmatoms.calc.mask
411 mask2 = ~mask1
412 apm1 = mmatoms.calc.apm1
413 apm2 = mmatoms.calc.apm2
414 apm = (apm1, apm2)
415 sigma = sigma[0] # Was already loopable 2-tuple
416 epsilon = epsilon[0]
418 mask = (mask1, mask2)
419 e_all = 0
420 qmforces_all = np.zeros_like(qmatoms.positions)
421 mmforces_all = np.zeros_like(mmatoms.positions)
423 # zip stops at shortest tuple so we dont double count
424 # cases of no counter ions.
425 for n, m, eps, sig in zip(apm, mask, epsilon, sigma):
426 mmpositions = self.update(qmatoms, mmatoms[m], n, shift)
427 qmforces = np.zeros_like(qmatoms.positions)
428 mmforces = np.zeros_like(mmatoms[m].positions)
429 energy = 0.0
431 qmpositions = qmatoms.positions.reshape((-1, self.qms, 3))
433 for q, qmpos in enumerate(qmpositions): # molwise loop
434 # cutoff from first atom of each mol
435 R00 = mmpositions[:, 0] - qmpos[0, :]
436 d002 = (R00**2).sum(1)
437 d00 = d002**0.5
438 x1 = d00 > self.rc - self.width
439 x2 = d00 < self.rc
440 x12 = np.logical_and(x1, x2)
441 y = (d00[x12] - self.rc + self.width) / self.width
442 t = np.zeros(len(d00))
443 t[x2] = 1.0
444 t[x12] -= y**2 * (3.0 - 2.0 * y)
445 dt = np.zeros(len(d00))
446 dt[x12] -= 6.0 / self.width * y * (1.0 - y)
447 for qa in range(len(qmpos)):
448 if ~np.any(eps[qa, :]):
449 continue
450 R = mmpositions - qmpos[qa, :]
451 d2 = (R**2).sum(2)
452 c6 = (sig[qa, :]**2 / d2)**3
453 c12 = c6**2
454 e = 4 * eps[qa, :] * (c12 - c6)
455 energy += np.dot(e.sum(1), t)
456 f = t[:, None, None] * (24 * eps[qa, :] *
457 (2 * c12 - c6) / d2)[:, :, None] * R
458 f00 = - (e.sum(1) * dt / d00)[:, None] * R00
459 mmforces += f.reshape((-1, 3))
460 qmforces[q * self.qms + qa, :] -= f.sum(0).sum(0)
461 qmforces[q * self.qms, :] -= f00.sum(0)
462 mmforces[::n, :] += f00
464 e_all += energy
465 qmforces_all += qmforces
466 mmforces_all[m] += mmforces
468 return e_all, qmforces_all, mmforces_all
470 def update(self, qmatoms, mmatoms, n, shift):
471 # Wrap point-charge positions to the MM-cell closest to the
472 # center of the the QM box, but avoid ripping molecules apart:
473 qmcenter = qmatoms.cell.diagonal() / 2
474 positions = mmatoms.positions.reshape((-1, n, 3)) + shift
476 # Distances from the center of the QM box to the first atom of
477 # each molecule:
478 distances = positions[:, 0] - qmcenter
480 wrap(distances, mmatoms.cell.diagonal(), mmatoms.pbc)
481 offsets = distances - positions[:, 0]
482 positions += offsets[:, np.newaxis] + qmcenter
484 return positions
487class LJInteractions:
488 name = 'LJ'
490 def __init__(self, parameters):
491 """Lennard-Jones type explicit interaction.
493 parameters: dict
494 Mapping from pair of atoms to tuple containing epsilon and sigma
495 for that pair.
497 Example:
499 lj = LJInteractions({('O', 'O'): (eps, sigma)})
501 """
502 self.parameters = {}
503 for (symbol1, symbol2), (epsilon, sigma) in parameters.items():
504 Z1 = atomic_numbers[symbol1]
505 Z2 = atomic_numbers[symbol2]
506 self.parameters[(Z1, Z2)] = epsilon, sigma
507 self.parameters[(Z2, Z1)] = epsilon, sigma
509 def calculate(self, qmatoms, mmatoms, shift):
510 qmforces = np.zeros_like(qmatoms.positions)
511 mmforces = np.zeros_like(mmatoms.positions)
512 species = set(mmatoms.numbers)
513 energy = 0.0
514 for R1, Z1, F1 in zip(qmatoms.positions, qmatoms.numbers, qmforces):
515 for Z2 in species:
516 if (Z1, Z2) not in self.parameters:
517 continue
518 epsilon, sigma = self.parameters[(Z1, Z2)]
519 mask = (mmatoms.numbers == Z2)
520 D = mmatoms.positions[mask] + shift - R1
521 wrap(D, mmatoms.cell.diagonal(), mmatoms.pbc)
522 d2 = (D**2).sum(1)
523 c6 = (sigma**2 / d2)**3
524 c12 = c6**2
525 energy += 4 * epsilon * (c12 - c6).sum()
526 f = 24 * epsilon * ((2 * c12 - c6) / d2)[:, np.newaxis] * D
527 F1 -= f.sum(0)
528 mmforces[mask] += f
529 return energy, qmforces, mmforces
532class RescaledCalculator(Calculator):
533 """Rescales length and energy of a calculators to match given
534 lattice constant and bulk modulus
536 Useful for MM calculator used within a :class:`ForceQMMM` model.
537 See T. D. Swinburne and J. R. Kermode, Phys. Rev. B 96, 144102 (2017)
538 for a derivation of the scaling constants.
539 """
540 implemented_properties = ['forces', 'energy', 'stress']
542 def __init__(self, mm_calc,
543 qm_lattice_constant, qm_bulk_modulus,
544 mm_lattice_constant, mm_bulk_modulus):
545 Calculator.__init__(self)
546 self.mm_calc = mm_calc
547 self.alpha = qm_lattice_constant / mm_lattice_constant
548 self.beta = mm_bulk_modulus / qm_bulk_modulus / (self.alpha**3)
550 def calculate(self, atoms, properties, system_changes):
551 Calculator.calculate(self, atoms, properties, system_changes)
553 # mm_pos = atoms.get_positions()
554 scaled_atoms = atoms.copy()
556 # scaled_atoms.positions = mm_pos/self.alpha
557 mm_cell = atoms.get_cell()
558 scaled_atoms.set_cell(mm_cell / self.alpha, scale_atoms=True)
560 results = {}
562 if 'energy' in properties:
563 energy = self.mm_calc.get_potential_energy(scaled_atoms)
564 results['energy'] = energy / self.beta
566 if 'forces' in properties:
567 forces = self.mm_calc.get_forces(scaled_atoms)
568 results['forces'] = forces / (self.beta * self.alpha)
570 if 'stress' in properties:
571 stress = self.mm_calc.get_stress(scaled_atoms)
572 results['stress'] = stress / (self.beta * self.alpha**3)
574 self.results = results
577class ForceConstantCalculator(Calculator):
578 """
579 Compute forces based on provided force-constant matrix
581 Useful with `ForceQMMM` to do harmonic QM/MM using force constants
582 of QM method.
583 """
584 implemented_properties = ['forces', 'energy']
586 def __init__(self, D, ref, f0):
587 """
588 Parameters:
590 D: matrix or sparse matrix, shape `(3*len(ref), 3*len(ref))`
591 Force constant matrix.
592 Sign convention is `D_ij = d^2E/(dx_i dx_j), so
593 `force = -D.dot(displacement)`
594 ref: ase.atoms.Atoms
595 Atoms object for reference configuration
596 f0: array, shape `(len(ref), 3)`
597 Value of forces at reference configuration
598 """
599 assert D.shape[0] == D.shape[1]
600 assert D.shape[0] // 3 == len(ref)
601 self.D = D
602 self.ref = ref
603 self.f0 = f0
604 self.size = len(ref)
605 Calculator.__init__(self)
607 def calculate(self, atoms, properties, system_changes):
608 Calculator.calculate(self, atoms, properties, system_changes)
609 u = atoms.positions - self.ref.positions
610 f = -self.D.dot(u.reshape(3 * self.size))
611 forces = np.zeros((len(atoms), 3))
612 forces[:, :] = f.reshape(self.size, 3)
613 self.results['forces'] = forces + self.f0
614 self.results['energy'] = 0.0
617class ForceQMMM(Calculator):
618 """
619 Force-based QM/MM calculator
621 QM forces are computed using a buffer region and then mixed abruptly
622 with MM forces:
624 F^i_QMMM = { F^i_QM if i in QM region
625 { F^i_MM otherwise
627 cf. N. Bernstein, J. R. Kermode, and G. Csanyi,
628 Rep. Prog. Phys. 72, 026501 (2009)
629 and T. D. Swinburne and J. R. Kermode, Phys. Rev. B 96, 144102 (2017).
630 """
631 implemented_properties = ['forces', 'energy']
633 def __init__(self,
634 atoms,
635 qm_selection_mask,
636 qm_calc,
637 mm_calc,
638 buffer_width,
639 vacuum=5.,
640 zero_mean=True,
641 qm_cell_round_off=3,
642 qm_radius=None):
643 """
644 ForceQMMM calculator
646 Parameters:
648 qm_selection_mask: list of ints, slice object or bool list/array
649 Selection out of atoms that belong to the QM region.
650 qm_calc: Calculator object
651 QM-calculator.
652 mm_calc: Calculator object
653 MM-calculator (should be scaled, see :class:`RescaledCalculator`)
654 Can use `ForceConstantCalculator` based on QM force constants, if
655 available.
656 vacuum: float or None
657 Amount of vacuum to add around QM atoms.
658 zero_mean: bool
659 If True, add a correction to zero the mean force in each direction
660 qm_cell_round_off: float
661 Tolerance value in Angstrom to round the qm cluster cell
662 qm_radius: 3x1 array of floats qm_radius for [x, y, z]
663 3d qm radius for calculation of qm cluster cell. default is None
664 and the radius is estimated from maximum distance between the atoms
665 in qm region.
666 """
668 if len(atoms[qm_selection_mask]) == 0:
669 raise ValueError("no QM atoms selected!")
671 self.qm_selection_mask = qm_selection_mask
672 self.qm_calc = qm_calc
673 self.mm_calc = mm_calc
674 self.vacuum = vacuum
675 self.buffer_width = buffer_width
676 self.zero_mean = zero_mean
677 self.qm_cell_round_off = qm_cell_round_off
678 self.qm_radius = qm_radius
680 self.qm_buffer_mask = None
682 Calculator.__init__(self)
684 def initialize_qm_buffer_mask(self, atoms):
685 """
686 Initialises system to perform qm calculation
687 """
688 # calculate the distances between all atoms and qm atoms
689 # qm_distance_matrix is a [N_QM_atoms x N_atoms] matrix
690 _, qm_distance_matrix = get_distances(
691 atoms.positions[self.qm_selection_mask],
692 atoms.positions,
693 atoms.cell, atoms.pbc)
695 self.qm_buffer_mask = np.zeros(len(atoms), dtype=bool)
697 # every r_qm is a matrix of distances
698 # from an atom in qm region and all atoms with size [N_atoms]
699 for r_qm in qm_distance_matrix:
700 self.qm_buffer_mask[r_qm < self.buffer_width] = True
702 def get_qm_cluster(self, atoms):
704 if self.qm_buffer_mask is None:
705 self.initialize_qm_buffer_mask(atoms)
707 qm_cluster = atoms[self.qm_buffer_mask]
708 del qm_cluster.constraints
710 round_cell = False
711 if self.qm_radius is None:
712 round_cell = True
713 # get all distances between qm atoms.
714 # Treat all X, Y and Z directions independently
715 # only distance between qm atoms is calculated
716 # in order to estimate qm radius in thee directions
717 R_qm, _ = get_distances(atoms.positions[self.qm_selection_mask],
718 cell=atoms.cell, pbc=atoms.pbc)
719 # estimate qm radius in three directions as 1/2
720 # of max distance between qm atoms
721 self.qm_radius = np.amax(np.amax(R_qm, axis=1), axis=0) * 0.5
723 if atoms.cell.orthorhombic:
724 cell_size = np.diagonal(atoms.cell)
725 else:
726 raise RuntimeError("NON-orthorhombic cell is not supported!")
728 # check if qm_cluster should be left periodic
729 # in periodic directions of the cell (cell[i] < qm_radius + buffer
730 # otherwise change to non pbc
731 # and make a cluster in a vacuum configuration
732 qm_cluster_pbc = (atoms.pbc &
733 (cell_size <
734 2.0 * (self.qm_radius + self.buffer_width)))
736 # start with the original orthorhombic cell
737 qm_cluster_cell = cell_size.copy()
738 # create a cluster in a vacuum cell in non periodic directions
739 qm_cluster_cell[~qm_cluster_pbc] = (
740 2.0 * (self.qm_radius[~qm_cluster_pbc] +
741 self.buffer_width +
742 self.vacuum))
744 if round_cell:
745 # round the qm cell to the required tolerance
746 qm_cluster_cell[~qm_cluster_pbc] = (np.round(
747 (qm_cluster_cell[~qm_cluster_pbc]) /
748 self.qm_cell_round_off) *
749 self.qm_cell_round_off)
751 qm_cluster.set_cell(Cell(np.diag(qm_cluster_cell)))
752 qm_cluster.pbc = qm_cluster_pbc
754 qm_shift = (0.5 * qm_cluster.cell.diagonal() -
755 qm_cluster.positions.mean(axis=0))
757 if 'cell_origin' in qm_cluster.info:
758 del qm_cluster.info['cell_origin']
760 # center the cluster only in non pbc directions
761 qm_cluster.positions[:, ~qm_cluster_pbc] += qm_shift[~qm_cluster_pbc]
763 return qm_cluster
765 def calculate(self, atoms, properties, system_changes):
766 Calculator.calculate(self, atoms, properties, system_changes)
768 qm_cluster = self.get_qm_cluster(atoms)
770 forces = self.mm_calc.get_forces(atoms)
771 qm_forces = self.qm_calc.get_forces(qm_cluster)
773 forces[self.qm_selection_mask] = \
774 qm_forces[self.qm_selection_mask[self.qm_buffer_mask]]
776 if self.zero_mean:
777 # Target is that: forces.sum(axis=1) == [0., 0., 0.]
778 forces[:] -= forces.mean(axis=0)
780 self.results['forces'] = forces
781 self.results['energy'] = 0.0
783 def get_region_from_masks(self, atoms=None, print_mapping=False):
784 """
785 creates region array from the masks of the calculators. The tags in
786 the array are:
787 QM - qm atoms
788 buffer - buffer atoms
789 MM - atoms treated with mm calculator
790 """
791 if atoms is None:
792 if self.atoms is None:
793 raise ValueError('Calculator has no atoms')
794 else:
795 atoms = self.atoms
797 region = np.full_like(atoms, "MM")
799 region[self.qm_selection_mask] = (
800 np.full_like(region[self.qm_selection_mask], "QM"))
802 buffer_only_mask = self.qm_buffer_mask & ~self.qm_selection_mask
804 region[buffer_only_mask] = np.full_like(region[buffer_only_mask],
805 "buffer")
807 if print_mapping:
809 print(f"Mapping of {len(region):5d} atoms in total:")
810 for region_id in np.unique(region):
811 n_at = np.count_nonzero(region == region_id)
812 print(f"{n_at:16d} {region_id}")
814 qm_atoms = atoms[self.qm_selection_mask]
815 symbol_counts = qm_atoms.symbols.formula.count()
816 print("QM atoms types:")
817 for symbol, count in symbol_counts.items():
818 print(f"{count:16d} {symbol}")
820 return region
822 def set_masks_from_region(self, region):
823 """
824 Sets masks from provided region array
825 """
826 self.qm_selection_mask = region == "QM"
827 buffer_mask = region == "buffer"
829 self.qm_buffer_mask = self.qm_selection_mask ^ buffer_mask
831 def export_extxyz(self, atoms=None, filename="qmmm_atoms.xyz"):
832 """
833 exports the atoms to extended xyz file with additional "region"
834 array keeping the mapping between QM, buffer and MM parts of
835 the simulation
836 """
837 if atoms is None:
838 if self.atoms is None:
839 raise ValueError('Calculator has no atoms')
840 else:
841 atoms = self.atoms
843 region = self.get_region_from_masks(atoms=atoms)
845 atoms_copy = atoms.copy()
846 atoms_copy.new_array("region", region)
848 atoms_copy.calc = self # to keep the calculation results
850 atoms_copy.write(filename, format='extxyz')
852 @classmethod
853 def import_extxyz(cls, filename, qm_calc, mm_calc):
854 """
855 A static method to import the the mapping from an estxyz file saved by
856 export_extxyz() function
857 Parameters
858 ----------
859 filename: string
860 filename with saved configuration
862 qm_calc: Calculator object
863 QM-calculator.
864 mm_calc: Calculator object
865 MM-calculator (should be scaled, see :class:`RescaledCalculator`)
866 Can use `ForceConstantCalculator` based on QM force constants, if
867 available.
869 Returns
870 -------
871 New object of ForceQMMM calculator with qm_selection_mask and
872 qm_buffer_mask set according to the region array in the saved file
873 """
875 from ase.io import read
876 atoms = read(filename, format='extxyz')
878 if "region" in atoms.arrays:
879 region = atoms.get_array("region")
880 else:
881 raise RuntimeError("Please provide extxyz file with 'region' array")
883 dummy_qm_mask = np.full_like(atoms, False, dtype=bool)
884 dummy_qm_mask[0] = True
886 self = cls(atoms, dummy_qm_mask,
887 qm_calc, mm_calc, buffer_width=1.0)
889 self.set_masks_from_region(region)
891 return self