Coverage for /builds/kinetik161/ase/ase/calculators/qmmm.py: 90.97%

1from typing import Sequence 

2 

3import numpy as np 

4 

5from ase.calculators.calculator import Calculator 

6from ase.cell import Cell 

7from ase.data import atomic_numbers 

8from ase.geometry import get_distances 

9from ase.utils import IOContext 

10 

11 

12class SimpleQMMM(Calculator): 

13 """Simple QMMM calculator.""" 

14 

15 implemented_properties = ['energy', 'forces'] 

16 

17 def __init__(self, selection, qmcalc, mmcalc1, mmcalc2, vacuum=None): 

18 """SimpleQMMM object. 

19 

20 The energy is calculated as:: 

21 

22 _ _ _ 

23 E = E (R ) - E (R ) + E (R ) 

24 QM QM MM QM MM all 

25 

26 parameters: 

27 

28 selection: list of int, slice object or list of bool 

29 Selection out of all the atoms that belong to the QM part. 

30 qmcalc: Calculator object 

31 QM-calculator. 

32 mmcalc1: Calculator object 

33 MM-calculator used for QM region. 

34 mmcalc2: Calculator object 

35 MM-calculator used for everything. 

36 vacuum: float or None 

37 Amount of vacuum to add around QM atoms. Use None if QM 

38 calculator doesn't need a box. 

39 

40 """ 

41 self.selection = selection 

42 self.qmcalc = qmcalc 

43 self.mmcalc1 = mmcalc1 

44 self.mmcalc2 = mmcalc2 

45 self.vacuum = vacuum 

46 

47 self.qmatoms = None 

48 self.center = None 

49 

50 Calculator.__init__(self) 

51 

52 def _get_name(self): 

53 return f'{self.qmcalc.name}-{self.mmcalc1.name}+{self.mmcalc1.name}' 

54 

55 def initialize_qm(self, atoms): 

56 constraints = atoms.constraints 

57 atoms.constraints = [] 

58 self.qmatoms = atoms[self.selection] 

59 atoms.constraints = constraints 

60 self.qmatoms.pbc = False 

61 if self.vacuum: 

62 self.qmatoms.center(vacuum=self.vacuum) 

63 self.center = self.qmatoms.positions.mean(axis=0) 

64 

65 def calculate(self, atoms, properties, system_changes): 

66 Calculator.calculate(self, atoms, properties, system_changes) 

67 

68 if self.qmatoms is None: 

69 self.initialize_qm(atoms) 

70 

71 self.qmatoms.positions = atoms.positions[self.selection] 

72 if self.vacuum: 

73 self.qmatoms.positions += (self.center - 

74 self.qmatoms.positions.mean(axis=0)) 

75 

76 energy = self.qmcalc.get_potential_energy(self.qmatoms) 

77 qmforces = self.qmcalc.get_forces(self.qmatoms) 

78 energy += self.mmcalc2.get_potential_energy(atoms) 

79 forces = self.mmcalc2.get_forces(atoms) 

80 

81 if self.vacuum: 

82 qmforces -= qmforces.mean(axis=0) 

83 forces[self.selection] += qmforces 

84 

85 energy -= self.mmcalc1.get_potential_energy(self.qmatoms) 

86 forces[self.selection] -= self.mmcalc1.get_forces(self.qmatoms) 

87 

88 self.results['energy'] = energy 

89 self.results['forces'] = forces 

90 

91 

92class EIQMMM(Calculator, IOContext): 

93 """Explicit interaction QMMM calculator.""" 

94 implemented_properties = ['energy', 'forces'] 

95 

96 def __init__(self, selection, qmcalc, mmcalc, interaction, 

97 vacuum=None, embedding=None, output=None): 

98 """EIQMMM object. 

99 

100 The energy is calculated as:: 

101 

102 _ _ _ _ 

103 E = E (R ) + E (R ) + E (R , R ) 

104 QM QM MM MM I QM MM 

105 

106 parameters: 

107 

108 selection: list of int, slice object or list of bool 

109 Selection out of all the atoms that belong to the QM part. 

110 qmcalc: Calculator object 

111 QM-calculator. 

112 mmcalc: Calculator object 

113 MM-calculator. 

114 interaction: Interaction object 

115 Interaction between QM and MM regions. 

116 vacuum: float or None 

117 Amount of vacuum to add around QM atoms. Use None if QM 

118 calculator doesn't need a box. 

119 embedding: Embedding object or None 

120 Specialized embedding object. Use None in order to use the 

121 default one. 

122 output: None, '-', str or file-descriptor. 

123 File for logging information - default is no logging (None). 

124 

125 """ 

126 

127 self.selection = selection 

128 

129 self.qmcalc = qmcalc 

130 self.mmcalc = mmcalc 

131 self.interaction = interaction 

132 self.vacuum = vacuum 

133 self.embedding = embedding 

134 

135 self.qmatoms = None 

136 self.mmatoms = None 

137 self.mask = None 

138 self.center = None # center of QM atoms in QM-box 

139 

140 self.output = self.openfile(output) 

141 

142 Calculator.__init__(self) 

143 

144 def _get_name(self): 

145 return f'{self.qmcalc.name}+{self.interaction.name}+{self.mmcalc.name}' 

146 

147 def initialize(self, atoms): 

148 self.mask = np.zeros(len(atoms), bool) 

149 self.mask[self.selection] = True 

150 

151 constraints = atoms.constraints 

152 atoms.constraints = [] # avoid slicing of constraints 

153 self.qmatoms = atoms[self.mask] 

154 self.mmatoms = atoms[~self.mask] 

155 atoms.constraints = constraints 

156 

157 self.qmatoms.pbc = False 

158 

159 if self.vacuum: 

160 self.qmatoms.center(vacuum=self.vacuum) 

161 self.center = self.qmatoms.positions.mean(axis=0) 

162 print('Size of QM-cell after centering:', 

163 self.qmatoms.cell.diagonal(), file=self.output) 

164 

165 self.qmatoms.calc = self.qmcalc 

166 self.mmatoms.calc = self.mmcalc 

167 

168 if self.embedding is None: 

169 self.embedding = Embedding() 

170 

171 self.embedding.initialize(self.qmatoms, self.mmatoms) 

172 print('Embedding:', self.embedding, file=self.output) 

173 

174 def calculate(self, atoms, properties, system_changes): 

175 Calculator.calculate(self, atoms, properties, system_changes) 

176 

177 if self.qmatoms is None: 

178 self.initialize(atoms) 

179 

180 self.mmatoms.set_positions(atoms.positions[~self.mask]) 

181 self.qmatoms.set_positions(atoms.positions[self.mask]) 

182 

183 if self.vacuum: 

184 shift = self.center - self.qmatoms.positions.mean(axis=0) 

185 self.qmatoms.positions += shift 

186 else: 

187 shift = (0, 0, 0) 

188 

189 self.embedding.update(shift) 

190 

191 ienergy, iqmforces, immforces = self.interaction.calculate( 

192 self.qmatoms, self.mmatoms, shift) 

193 

194 qmenergy = self.qmatoms.get_potential_energy() 

195 mmenergy = self.mmatoms.get_potential_energy() 

196 energy = ienergy + qmenergy + mmenergy 

197 

198 print('Energies: {:12.3f} {:+12.3f} {:+12.3f} = {:12.3f}' 

199 .format(ienergy, qmenergy, mmenergy, energy), file=self.output) 

200 

201 qmforces = self.qmatoms.get_forces() 

202 mmforces = self.mmatoms.get_forces() 

203 

204 mmforces += self.embedding.get_mm_forces() 

205 

206 forces = np.empty((len(atoms), 3)) 

207 forces[self.mask] = qmforces + iqmforces 

208 forces[~self.mask] = mmforces + immforces 

209 

210 self.results['energy'] = energy 

211 self.results['forces'] = forces 

212 

213 

214def wrap(D, cell, pbc): 

215 """Wrap distances to nearest neighbor (minimum image convention).""" 

216 for i, periodic in enumerate(pbc): 

217 if periodic: 

218 d = D[:, i] 

219 L = cell[i] 

220 d[:] = (d + L / 2) % L - L / 2 # modify D inplace 

221 

222 

223class Embedding: 

224 def __init__(self, molecule_size=3, **parameters): 

225 """Point-charge embedding.""" 

226 self.qmatoms = None 

227 self.mmatoms = None 

228 self.molecule_size = molecule_size 

229 self.virtual_molecule_size = None 

230 self.parameters = parameters 

231 

232 def __repr__(self): 

233 return f'Embedding(molecule_size={self.molecule_size})' 

234 

235 def initialize(self, qmatoms, mmatoms): 

236 """Hook up embedding object to QM and MM atoms objects.""" 

237 self.qmatoms = qmatoms 

238 self.mmatoms = mmatoms 

239 charges = mmatoms.calc.get_virtual_charges(mmatoms) 

240 self.pcpot = qmatoms.calc.embed(charges, **self.parameters) 

241 self.virtual_molecule_size = (self.molecule_size * 

242 len(charges) // len(mmatoms)) 

243 

244 def update(self, shift): 

245 """Update point-charge positions.""" 

246 # Wrap point-charge positions to the MM-cell closest to the 

247 # center of the the QM box, but avoid ripping molecules apart: 

248 qmcenter = self.qmatoms.positions.mean(axis=0) 

249 # if counter ions are used, then molecule_size has more than 1 value 

250 if self.mmatoms.calc.name == 'combinemm': 

251 mask1 = self.mmatoms.calc.mask 

252 mask2 = ~mask1 

253 vmask1 = self.mmatoms.calc.virtual_mask 

254 vmask2 = ~vmask1 

255 apm1 = self.mmatoms.calc.apm1 

256 apm2 = self.mmatoms.calc.apm2 

257 spm1 = self.mmatoms.calc.atoms1.calc.sites_per_mol 

258 spm2 = self.mmatoms.calc.atoms2.calc.sites_per_mol 

259 pos = self.mmatoms.positions 

260 pos1 = pos[mask1].reshape((-1, apm1, 3)) 

261 pos2 = pos[mask2].reshape((-1, apm2, 3)) 

262 pos = (pos1, pos2) 

263 else: 

264 pos = (self.mmatoms.positions, ) 

265 apm1 = self.molecule_size 

266 apm2 = self.molecule_size 

267 # This is only specific to calculators where apm != spm, 

268 # i.e. TIP4P. Non-native MM calcs do not have this attr. 

269 if hasattr(self.mmatoms.calc, 'sites_per_mol'): 

270 spm1 = self.mmatoms.calc.sites_per_mol 

271 spm2 = self.mmatoms.calc.sites_per_mol 

272 else: 

273 spm1 = self.molecule_size 

274 spm2 = spm1 

275 mask1 = np.ones(len(self.mmatoms), dtype=bool) 

276 mask2 = mask1 

277 

278 wrap_pos = np.zeros_like(self.mmatoms.positions) 

279 com_all = [] 

280 apm = (apm1, apm2) 

281 mask = (mask1, mask2) 

282 spm = (spm1, spm2) 

283 for p, n, m, vn in zip(pos, apm, mask, spm): 

284 positions = p.reshape((-1, n, 3)) + shift 

285 

286 # Distances from the center of the QM box to the first atom of 

287 # each molecule: 

288 distances = positions[:, 0] - qmcenter 

289 

290 wrap(distances, self.mmatoms.cell.diagonal(), self.mmatoms.pbc) 

291 offsets = distances - positions[:, 0] 

292 positions += offsets[:, np.newaxis] + qmcenter 

293 

294 # Geometric center positions for each mm mol for LR cut 

295 com = np.array([p.mean(axis=0) for p in positions]) 

296 # Need per atom for C-code: 

297 com_pv = np.repeat(com, vn, axis=0) 

298 com_all.append(com_pv) 

299 

300 wrap_pos[m] = positions.reshape((-1, 3)) 

301 

302 positions = wrap_pos.copy() 

303 positions = self.mmatoms.calc.add_virtual_sites(positions) 

304 

305 if self.mmatoms.calc.name == 'combinemm': 

306 com_pv = np.zeros_like(positions) 

307 for ii, m in enumerate((vmask1, vmask2)): 

308 com_pv[m] = com_all[ii] 

309 

310 # compatibility with gpaw versions w/o LR cut in PointChargePotential 

311 if 'rc2' in self.parameters: 

312 self.pcpot.set_positions(positions, com_pv=com_pv) 

313 else: 

314 self.pcpot.set_positions(positions) 

315 

316 def get_mm_forces(self): 

317 """Calculate the forces on the MM-atoms from the QM-part.""" 

318 f = self.pcpot.get_forces(self.qmatoms.calc) 

319 return self.mmatoms.calc.redistribute_forces(f) 

320 

321 

322def combine_lj_lorenz_berthelot(sigmaqm, sigmamm, 

323 epsilonqm, epsilonmm): 

324 """Combine LJ parameters according to the Lorenz-Berthelot rule""" 

325 sigma = [] 

326 epsilon = [] 

327 # check if input is tuple of vals for more than 1 mm calc, or only for 1. 

328 if isinstance(sigmamm, Sequence): 

329 numcalcs = len(sigmamm) 

330 else: 

331 numcalcs = 1 # if only 1 mm calc, eps and sig are simply np arrays 

332 sigmamm = (sigmamm, ) 

333 epsilonmm = (epsilonmm, ) 

334 for cc in range(numcalcs): 

335 sigma_c = np.zeros((len(sigmaqm), len(sigmamm[cc]))) 

336 epsilon_c = np.zeros_like(sigma_c) 

337 

338 for ii in range(len(sigmaqm)): 

339 sigma_c[ii, :] = (sigmaqm[ii] + sigmamm[cc]) / 2 

340 epsilon_c[ii, :] = (epsilonqm[ii] * epsilonmm[cc])**0.5 

341 sigma.append(sigma_c) 

342 epsilon.append(epsilon_c) 

343 

344 if numcalcs == 1: # retain original, 1 calc function 

345 sigma = np.array(sigma[0]) 

346 epsilon = np.array(epsilon[0]) 

347 

348 return sigma, epsilon 

349 

350 

351class LJInteractionsGeneral: 

352 name = 'LJ-general' 

353 

354 def __init__(self, sigmaqm, epsilonqm, sigmamm, epsilonmm, 

355 qm_molecule_size, mm_molecule_size=3, 

356 rc=np.Inf, width=1.0): 

357 """General Lennard-Jones type explicit interaction. 

358 

359 sigmaqm: array 

360 Array of sigma-parameters which should have the length of the QM 

361 subsystem 

362 epsilonqm: array 

363 As sigmaqm, but for epsilon-paramaters 

364 sigmamm: Either array (A) or tuple (B) 

365 A (no counterions): 

366 Array of sigma-parameters with the length of the smallests 

367 repeating atoms-group (i.e. molecule) of the MM subsystem 

368 B (counterions): 

369 Tuple: (arr1, arr2), where arr1 is an array of sigmas with 

370 the length of counterions in the MM subsystem, and 

371 arr2 is the array from A. 

372 epsilonmm: array or tuple 

373 As sigmamm but for epsilon-parameters. 

374 qm_molecule_size: int 

375 number of atoms of the smallest repeating atoms-group (i.e. 

376 molecule) in the QM subsystem (often just the number of atoms 

377 in the QM subsystem) 

378 mm_molecule_size: int 

379 as qm_molecule_size but for the MM subsystem. Will be overwritten 

380 if counterions are present in the MM subsystem (via the CombineMM 

381 calculator) 

382 

383 """ 

384 self.sigmaqm = sigmaqm 

385 self.epsilonqm = epsilonqm 

386 self.sigmamm = sigmamm 

387 self.epsilonmm = epsilonmm 

388 self.qms = qm_molecule_size 

389 self.mms = mm_molecule_size 

390 self.rc = rc 

391 self.width = width 

392 self.combine_lj() 

393 

394 def combine_lj(self): 

395 self.sigma, self.epsilon = combine_lj_lorenz_berthelot( 

396 self.sigmaqm, self.sigmamm, self.epsilonqm, self.epsilonmm) 

397 

398 def calculate(self, qmatoms, mmatoms, shift): 

399 epsilon = self.epsilon 

400 sigma = self.sigma 

401 

402 # loop over possible multiple mm calculators 

403 # currently 1 or 2, but could be generalized in the future... 

404 apm1 = self.mms 

405 mask1 = np.ones(len(mmatoms), dtype=bool) 

406 mask2 = mask1 

407 apm = (apm1, ) 

408 sigma = (sigma, ) 

409 epsilon = (epsilon, ) 

410 if hasattr(mmatoms.calc, 'name'): 

411 if mmatoms.calc.name == 'combinemm': 

412 mask1 = mmatoms.calc.mask 

413 mask2 = ~mask1 

414 apm1 = mmatoms.calc.apm1 

415 apm2 = mmatoms.calc.apm2 

416 apm = (apm1, apm2) 

417 sigma = sigma[0] # Was already loopable 2-tuple 

418 epsilon = epsilon[0] 

419 

420 mask = (mask1, mask2) 

421 e_all = 0 

422 qmforces_all = np.zeros_like(qmatoms.positions) 

423 mmforces_all = np.zeros_like(mmatoms.positions) 

424 

425 # zip stops at shortest tuple so we dont double count 

426 # cases of no counter ions. 

427 for n, m, eps, sig in zip(apm, mask, epsilon, sigma): 

428 mmpositions = self.update(qmatoms, mmatoms[m], n, shift) 

429 qmforces = np.zeros_like(qmatoms.positions) 

430 mmforces = np.zeros_like(mmatoms[m].positions) 

431 energy = 0.0 

432 

433 qmpositions = qmatoms.positions.reshape((-1, self.qms, 3)) 

434 

435 for q, qmpos in enumerate(qmpositions): # molwise loop 

436 # cutoff from first atom of each mol 

437 R00 = mmpositions[:, 0] - qmpos[0, :] 

438 d002 = (R00**2).sum(1) 

439 d00 = d002**0.5 

440 x1 = d00 > self.rc - self.width 

441 x2 = d00 < self.rc 

442 x12 = np.logical_and(x1, x2) 

443 y = (d00[x12] - self.rc + self.width) / self.width 

444 t = np.zeros(len(d00)) 

445 t[x2] = 1.0 

446 t[x12] -= y**2 * (3.0 - 2.0 * y) 

447 dt = np.zeros(len(d00)) 

448 dt[x12] -= 6.0 / self.width * y * (1.0 - y) 

449 for qa in range(len(qmpos)): 

450 if ~np.any(eps[qa, :]): 

451 continue 

452 R = mmpositions - qmpos[qa, :] 

453 d2 = (R**2).sum(2) 

454 c6 = (sig[qa, :]**2 / d2)**3 

455 c12 = c6**2 

456 e = 4 * eps[qa, :] * (c12 - c6) 

457 energy += np.dot(e.sum(1), t) 

458 f = t[:, None, None] * (24 * eps[qa, :] * 

459 (2 * c12 - c6) / d2)[:, :, None] * R 

460 f00 = - (e.sum(1) * dt / d00)[:, None] * R00 

461 mmforces += f.reshape((-1, 3)) 

462 qmforces[q * self.qms + qa, :] -= f.sum(0).sum(0) 

463 qmforces[q * self.qms, :] -= f00.sum(0) 

464 mmforces[::n, :] += f00 

465 

466 e_all += energy 

467 qmforces_all += qmforces 

468 mmforces_all[m] += mmforces 

469 

470 return e_all, qmforces_all, mmforces_all 

471 

472 def update(self, qmatoms, mmatoms, n, shift): 

473 # Wrap point-charge positions to the MM-cell closest to the 

474 # center of the the QM box, but avoid ripping molecules apart: 

475 qmcenter = qmatoms.cell.diagonal() / 2 

476 positions = mmatoms.positions.reshape((-1, n, 3)) + shift 

477 

478 # Distances from the center of the QM box to the first atom of 

479 # each molecule: 

480 distances = positions[:, 0] - qmcenter 

481 

482 wrap(distances, mmatoms.cell.diagonal(), mmatoms.pbc) 

483 offsets = distances - positions[:, 0] 

484 positions += offsets[:, np.newaxis] + qmcenter 

485 

486 return positions 

487 

488 

489class LJInteractions: 

490 name = 'LJ' 

491 

492 def __init__(self, parameters): 

493 """Lennard-Jones type explicit interaction. 

494 

495 parameters: dict 

496 Mapping from pair of atoms to tuple containing epsilon and sigma 

497 for that pair. 

498 

499 Example: 

500 

501 lj = LJInteractions({('O', 'O'): (eps, sigma)}) 

502 

503 """ 

504 self.parameters = {} 

505 for (symbol1, symbol2), (epsilon, sigma) in parameters.items(): 

506 Z1 = atomic_numbers[symbol1] 

507 Z2 = atomic_numbers[symbol2] 

508 self.parameters[(Z1, Z2)] = epsilon, sigma 

509 self.parameters[(Z2, Z1)] = epsilon, sigma 

510 

511 def calculate(self, qmatoms, mmatoms, shift): 

512 qmforces = np.zeros_like(qmatoms.positions) 

513 mmforces = np.zeros_like(mmatoms.positions) 

514 species = set(mmatoms.numbers) 

515 energy = 0.0 

516 for R1, Z1, F1 in zip(qmatoms.positions, qmatoms.numbers, qmforces): 

517 for Z2 in species: 

518 if (Z1, Z2) not in self.parameters: 

519 continue 

520 epsilon, sigma = self.parameters[(Z1, Z2)] 

521 mask = (mmatoms.numbers == Z2) 

522 D = mmatoms.positions[mask] + shift - R1 

523 wrap(D, mmatoms.cell.diagonal(), mmatoms.pbc) 

524 d2 = (D**2).sum(1) 

525 c6 = (sigma**2 / d2)**3 

526 c12 = c6**2 

527 energy += 4 * epsilon * (c12 - c6).sum() 

528 f = 24 * epsilon * ((2 * c12 - c6) / d2)[:, np.newaxis] * D 

529 F1 -= f.sum(0) 

530 mmforces[mask] += f 

531 return energy, qmforces, mmforces 

532 

533 

534class RescaledCalculator(Calculator): 

535 """Rescales length and energy of a calculators to match given 

536 lattice constant and bulk modulus 

537 

538 Useful for MM calculator used within a :class:`ForceQMMM` model. 

539 See T. D. Swinburne and J. R. Kermode, Phys. Rev. B 96, 144102 (2017) 

540 for a derivation of the scaling constants. 

541 """ 

542 implemented_properties = ['forces', 'energy', 'stress'] 

543 

544 def __init__(self, mm_calc, 

545 qm_lattice_constant, qm_bulk_modulus, 

546 mm_lattice_constant, mm_bulk_modulus): 

547 Calculator.__init__(self) 

548 self.mm_calc = mm_calc 

549 self.alpha = qm_lattice_constant / mm_lattice_constant 

550 self.beta = mm_bulk_modulus / qm_bulk_modulus / (self.alpha**3) 

551 

552 def calculate(self, atoms, properties, system_changes): 

553 Calculator.calculate(self, atoms, properties, system_changes) 

554 

555 # mm_pos = atoms.get_positions() 

556 scaled_atoms = atoms.copy() 

557 

558 # scaled_atoms.positions = mm_pos/self.alpha 

559 mm_cell = atoms.get_cell() 

560 scaled_atoms.set_cell(mm_cell / self.alpha, scale_atoms=True) 

561 

562 results = {} 

563 

564 if 'energy' in properties: 

565 energy = self.mm_calc.get_potential_energy(scaled_atoms) 

566 results['energy'] = energy / self.beta 

567 

568 if 'forces' in properties: 

569 forces = self.mm_calc.get_forces(scaled_atoms) 

570 results['forces'] = forces / (self.beta * self.alpha) 

571 

572 if 'stress' in properties: 

573 stress = self.mm_calc.get_stress(scaled_atoms) 

574 results['stress'] = stress / (self.beta * self.alpha**3) 

575 

576 self.results = results 

577 

578 

579class ForceConstantCalculator(Calculator): 

580 """ 

581 Compute forces based on provided force-constant matrix 

582 

583 Useful with `ForceQMMM` to do harmonic QM/MM using force constants 

584 of QM method. 

585 """ 

586 implemented_properties = ['forces', 'energy'] 

587 

588 def __init__(self, D, ref, f0): 

589 """ 

590 Parameters: 

591 

592 D: matrix or sparse matrix, shape `(3*len(ref), 3*len(ref))` 

593 Force constant matrix. 

594 Sign convention is `D_ij = d^2E/(dx_i dx_j), so 

595 `force = -D.dot(displacement)` 

596 ref: ase.atoms.Atoms 

597 Atoms object for reference configuration 

598 f0: array, shape `(len(ref), 3)` 

599 Value of forces at reference configuration 

600 """ 

601 assert D.shape[0] == D.shape[1] 

602 assert D.shape[0] // 3 == len(ref) 

603 self.D = D 

604 self.ref = ref 

605 self.f0 = f0 

606 self.size = len(ref) 

607 Calculator.__init__(self) 

608 

609 def calculate(self, atoms, properties, system_changes): 

610 Calculator.calculate(self, atoms, properties, system_changes) 

611 u = atoms.positions - self.ref.positions 

612 f = -self.D.dot(u.reshape(3 * self.size)) 

613 forces = np.zeros((len(atoms), 3)) 

614 forces[:, :] = f.reshape(self.size, 3) 

615 self.results['forces'] = forces + self.f0 

616 self.results['energy'] = 0.0 

617 

618 

619class ForceQMMM(Calculator): 

620 """ 

621 Force-based QM/MM calculator 

622 

623 QM forces are computed using a buffer region and then mixed abruptly 

624 with MM forces: 

625 

626 F^i_QMMM = { F^i_QM if i in QM region 

627 { F^i_MM otherwise 

628 

629 cf. N. Bernstein, J. R. Kermode, and G. Csanyi, 

630 Rep. Prog. Phys. 72, 026501 (2009) 

631 and T. D. Swinburne and J. R. Kermode, Phys. Rev. B 96, 144102 (2017). 

632 """ 

633 implemented_properties = ['forces', 'energy'] 

634 

635 def __init__(self, 

636 atoms, 

637 qm_selection_mask, 

638 qm_calc, 

639 mm_calc, 

640 buffer_width, 

641 vacuum=5., 

642 zero_mean=True, 

643 qm_cell_round_off=3, 

644 qm_radius=None): 

645 """ 

646 ForceQMMM calculator 

647 

648 Parameters: 

649 

650 qm_selection_mask: list of ints, slice object or bool list/array 

651 Selection out of atoms that belong to the QM region. 

652 qm_calc: Calculator object 

653 QM-calculator. 

654 mm_calc: Calculator object 

655 MM-calculator (should be scaled, see :class:`RescaledCalculator`) 

656 Can use `ForceConstantCalculator` based on QM force constants, if 

657 available. 

658 vacuum: float or None 

659 Amount of vacuum to add around QM atoms. 

660 zero_mean: bool 

661 If True, add a correction to zero the mean force in each direction 

662 qm_cell_round_off: float 

663 Tolerance value in Angstrom to round the qm cluster cell 

664 qm_radius: 3x1 array of floats qm_radius for [x, y, z] 

665 3d qm radius for calculation of qm cluster cell. default is None 

666 and the radius is estimated from maximum distance between the atoms 

667 in qm region. 

668 """ 

669 

670 if len(atoms[qm_selection_mask]) == 0: 

671 raise ValueError("no QM atoms selected!") 

672 

673 self.qm_selection_mask = qm_selection_mask 

674 self.qm_calc = qm_calc 

675 self.mm_calc = mm_calc 

676 self.vacuum = vacuum 

677 self.buffer_width = buffer_width 

678 self.zero_mean = zero_mean 

679 self.qm_cell_round_off = qm_cell_round_off 

680 self.qm_radius = qm_radius 

681 

682 self.qm_buffer_mask = None 

683 

684 Calculator.__init__(self) 

685 

686 def initialize_qm_buffer_mask(self, atoms): 

687 """ 

688 Initialises system to perform qm calculation 

689 """ 

690 # calculate the distances between all atoms and qm atoms 

691 # qm_distance_matrix is a [N_QM_atoms x N_atoms] matrix 

692 _, qm_distance_matrix = get_distances( 

693 atoms.positions[self.qm_selection_mask], 

694 atoms.positions, 

695 atoms.cell, atoms.pbc) 

696 

697 self.qm_buffer_mask = np.zeros(len(atoms), dtype=bool) 

698 

699 # every r_qm is a matrix of distances 

700 # from an atom in qm region and all atoms with size [N_atoms] 

701 for r_qm in qm_distance_matrix: 

702 self.qm_buffer_mask[r_qm < self.buffer_width] = True 

703 

704 def get_qm_cluster(self, atoms): 

705 

706 if self.qm_buffer_mask is None: 

707 self.initialize_qm_buffer_mask(atoms) 

708 

709 qm_cluster = atoms[self.qm_buffer_mask] 

710 del qm_cluster.constraints 

711 

712 round_cell = False 

713 if self.qm_radius is None: 

714 round_cell = True 

715 # get all distances between qm atoms. 

716 # Treat all X, Y and Z directions independently 

717 # only distance between qm atoms is calculated 

718 # in order to estimate qm radius in thee directions 

719 R_qm, _ = get_distances(atoms.positions[self.qm_selection_mask], 

720 cell=atoms.cell, pbc=atoms.pbc) 

721 # estimate qm radius in three directions as 1/2 

722 # of max distance between qm atoms 

723 self.qm_radius = np.amax(np.amax(R_qm, axis=1), axis=0) * 0.5 

724 

725 if atoms.cell.orthorhombic: 

726 cell_size = np.diagonal(atoms.cell) 

727 else: 

728 raise RuntimeError("NON-orthorhombic cell is not supported!") 

729 

730 # check if qm_cluster should be left periodic 

731 # in periodic directions of the cell (cell[i] < qm_radius + buffer 

732 # otherwise change to non pbc 

733 # and make a cluster in a vacuum configuration 

734 qm_cluster_pbc = (atoms.pbc & 

735 (cell_size < 

736 2.0 * (self.qm_radius + self.buffer_width))) 

737 

738 # start with the original orthorhombic cell 

739 qm_cluster_cell = cell_size.copy() 

740 # create a cluster in a vacuum cell in non periodic directions 

741 qm_cluster_cell[~qm_cluster_pbc] = ( 

742 2.0 * (self.qm_radius[~qm_cluster_pbc] + 

743 self.buffer_width + 

744 self.vacuum)) 

745 

746 if round_cell: 

747 # round the qm cell to the required tolerance 

748 qm_cluster_cell[~qm_cluster_pbc] = (np.round( 

749 (qm_cluster_cell[~qm_cluster_pbc]) / 

750 self.qm_cell_round_off) * 

751 self.qm_cell_round_off) 

752 

753 qm_cluster.set_cell(Cell(np.diag(qm_cluster_cell))) 

754 qm_cluster.pbc = qm_cluster_pbc 

755 

756 qm_shift = (0.5 * qm_cluster.cell.diagonal() - 

757 qm_cluster.positions.mean(axis=0)) 

758 

759 if 'cell_origin' in qm_cluster.info: 

760 del qm_cluster.info['cell_origin'] 

761 

762 # center the cluster only in non pbc directions 

763 qm_cluster.positions[:, ~qm_cluster_pbc] += qm_shift[~qm_cluster_pbc] 

764 

765 return qm_cluster 

766 

767 def calculate(self, atoms, properties, system_changes): 

768 Calculator.calculate(self, atoms, properties, system_changes) 

769 

770 qm_cluster = self.get_qm_cluster(atoms) 

771 

772 forces = self.mm_calc.get_forces(atoms) 

773 qm_forces = self.qm_calc.get_forces(qm_cluster) 

774 

775 forces[self.qm_selection_mask] = \ 

776 qm_forces[self.qm_selection_mask[self.qm_buffer_mask]] 

777 

778 if self.zero_mean: