Hide keyboard shortcuts

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 

2 

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 

8 

9 

10class SimpleQMMM(Calculator): 

11 """Simple QMMM calculator.""" 

12 

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

14 

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

16 """SimpleQMMM object. 

17 

18 The energy is calculated as:: 

19 

20 _ _ _ 

21 E = E (R ) - E (R ) + E (R ) 

22 QM QM MM QM MM all 

23 

24 parameters: 

25 

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. 

37 

38 """ 

39 self.selection = selection 

40 self.qmcalc = qmcalc 

41 self.mmcalc1 = mmcalc1 

42 self.mmcalc2 = mmcalc2 

43 self.vacuum = vacuum 

44 

45 self.qmatoms = None 

46 self.center = None 

47 

48 Calculator.__init__(self) 

49 

50 def _get_name(self): 

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

52 

53 def initialize_qm(self, atoms): 

54 constraints = atoms.constraints 

55 atoms.constraints = [] 

56 self.qmatoms = atoms[self.selection] 

57 atoms.constraints = constraints 

58 self.qmatoms.pbc = False 

59 if self.vacuum: 

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

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

62 

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

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

65 

66 if self.qmatoms is None: 

67 self.initialize_qm(atoms) 

68 

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

70 if self.vacuum: 

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

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

73 

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

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

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

77 forces = self.mmcalc2.get_forces(atoms) 

78 

79 if self.vacuum: 

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

81 forces[self.selection] += qmforces 

82 

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

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

85 

86 self.results['energy'] = energy 

87 self.results['forces'] = forces 

88 

89 

90class EIQMMM(Calculator, IOContext): 

91 """Explicit interaction QMMM calculator.""" 

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

93 

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

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

96 """EIQMMM object. 

97 

98 The energy is calculated as:: 

99 

100 _ _ _ _ 

101 E = E (R ) + E (R ) + E (R , R ) 

102 QM QM MM MM I QM MM 

103 

104 parameters: 

105 

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

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

108 qmcalc: Calculator object 

109 QM-calculator. 

110 mmcalc: Calculator object 

111 MM-calculator. 

112 interaction: Interaction object 

113 Interaction between QM and MM regions. 

114 vacuum: float or None 

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

116 calculator doesn't need a box. 

117 embedding: Embedding object or None 

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

119 default one. 

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

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

122 

123 """ 

124 

125 self.selection = selection 

126 

127 self.qmcalc = qmcalc 

128 self.mmcalc = mmcalc 

129 self.interaction = interaction 

130 self.vacuum = vacuum 

131 self.embedding = embedding 

132 

133 self.qmatoms = None 

134 self.mmatoms = None 

135 self.mask = None 

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

137 

138 self.output = self.openfile(output) 

139 

140 Calculator.__init__(self) 

141 

142 def _get_name(self): 

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

144 

145 def initialize(self, atoms): 

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

147 self.mask[self.selection] = True 

148 

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 

154 

155 self.qmatoms.pbc = False 

156 

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) 

162 

163 self.qmatoms.calc = self.qmcalc 

164 self.mmatoms.calc = self.mmcalc 

165 

166 if self.embedding is None: 

167 self.embedding = Embedding() 

168 

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

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

171 

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

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

174 

175 if self.qmatoms is None: 

176 self.initialize(atoms) 

177 

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

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

180 

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) 

186 

187 self.embedding.update(shift) 

188 

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

190 self.qmatoms, self.mmatoms, shift) 

191 

192 qmenergy = self.qmatoms.get_potential_energy() 

193 mmenergy = self.mmatoms.get_potential_energy() 

194 energy = ienergy + qmenergy + mmenergy 

195 

196 print('Energies: {0:12.3f} {1:+12.3f} {2:+12.3f} = {3:12.3f}' 

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

198 

199 qmforces = self.qmatoms.get_forces() 

200 mmforces = self.mmatoms.get_forces() 

201 

202 mmforces += self.embedding.get_mm_forces() 

203 

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

205 forces[self.mask] = qmforces + iqmforces 

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

207 

208 self.results['energy'] = energy 

209 self.results['forces'] = forces 

210 

211 

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 

219 

220 

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 

229 

230 def __repr__(self): 

231 return 'Embedding(molecule_size={0})'.format(self.molecule_size) 

232 

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)) 

241 

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 

275 

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 

283 

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

285 # each molecule: 

286 distances = positions[:, 0] - qmcenter 

287 

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

289 offsets = distances - positions[:, 0] 

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

291 

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) 

297 

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

299 

300 positions = wrap_pos.copy() 

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

302 

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] 

307 

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) 

313 

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) 

318 

319 

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) 

335 

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) 

341 

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

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

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

345 

346 return sigma, epsilon 

347 

348 

349class LJInteractionsGeneral: 

350 name = 'LJ-general' 

351 

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. 

356 

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) 

380 

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() 

391 

392 def combine_lj(self): 

393 self.sigma, self.epsilon = combine_lj_lorenz_berthelot( 

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

395 

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

397 epsilon = self.epsilon 

398 sigma = self.sigma 

399 

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] 

417 

418 mask = (mask1, mask2) 

419 e_all = 0 

420 qmforces_all = np.zeros_like(qmatoms.positions) 

421 mmforces_all = np.zeros_like(mmatoms.positions) 

422 

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 

430 

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

432 

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 

463 

464 e_all += energy 

465 qmforces_all += qmforces 

466 mmforces_all[m] += mmforces 

467 

468 return e_all, qmforces_all, mmforces_all 

469 

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 

475 

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

477 # each molecule: 

478 distances = positions[:, 0] - qmcenter 

479 

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

481 offsets = distances - positions[:, 0] 

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

483 

484 return positions 

485 

486 

487class LJInteractions: 

488 name = 'LJ' 

489 

490 def __init__(self, parameters): 

491 """Lennard-Jones type explicit interaction. 

492 

493 parameters: dict 

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

495 for that pair. 

496 

497 Example: 

498 

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

500 

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 

508 

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 

530 

531 

532class RescaledCalculator(Calculator): 

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

534 lattice constant and bulk modulus 

535 

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'] 

541 

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) 

549 

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

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

552 

553 # mm_pos = atoms.get_positions() 

554 scaled_atoms = atoms.copy() 

555 

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) 

559 

560 results = {} 

561 

562 if 'energy' in properties: 

563 energy = self.mm_calc.get_potential_energy(scaled_atoms) 

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

565 

566 if 'forces' in properties: 

567 forces = self.mm_calc.get_forces(scaled_atoms) 

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

569 

570 if 'stress' in properties: 

571 stress = self.mm_calc.get_stress(scaled_atoms) 

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

573 

574 self.results = results 

575 

576 

577class ForceConstantCalculator(Calculator): 

578 """ 

579 Compute forces based on provided force-constant matrix 

580 

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

582 of QM method. 

583 """ 

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

585 

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

587 """ 

588 Parameters: 

589 

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) 

606 

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 

615 

616 

617class ForceQMMM(Calculator): 

618 """ 

619 Force-based QM/MM calculator 

620 

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

622 with MM forces: 

623 

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

625 { F^i_MM otherwise 

626 

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'] 

632 

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 

645 

646 Parameters: 

647 

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 """ 

667 

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

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

670 

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 

679 

680 self.qm_buffer_mask = None 

681 

682 Calculator.__init__(self) 

683 

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) 

694 

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

696 

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 

701 

702 def get_qm_cluster(self, atoms): 

703 

704 if self.qm_buffer_mask is None: 

705 self.initialize_qm_buffer_mask(atoms) 

706 

707 qm_cluster = atoms[self.qm_buffer_mask] 

708 del qm_cluster.constraints 

709 

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 

722 

723 if atoms.cell.orthorhombic: 

724 cell_size = np.diagonal(atoms.cell) 

725 else: 

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

727 

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))) 

735 

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)) 

743 

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) 

750 

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

752 qm_cluster.pbc = qm_cluster_pbc 

753 

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

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

756 

757 if 'cell_origin' in qm_cluster.info: 

758 del qm_cluster.info['cell_origin'] 

759 

760 # center the cluster only in non pbc directions 

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

762 

763 return qm_cluster 

764 

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

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

767 

768 qm_cluster = self.get_qm_cluster(atoms) 

769 

770 forces = self.mm_calc.get_forces(atoms) 

771 qm_forces = self.qm_calc.get_forces(qm_cluster) 

772 

773 forces[self.qm_selection_mask] = \ 

774 qm_forces[self.qm_selection_mask[self.qm_buffer_mask]] 

775 

776 if self.zero_mean: 

777 # Target is that: forces.sum(axis=1) == [0., 0., 0.] 

778 forces[:] -= forces.mean(axis=0) 

779 

780 self.results['forces'] = forces 

781 self.results['energy'] = 0.0 

782 

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 

796 

797 region = np.full_like(atoms, "MM") 

798 

799 region[self.qm_selection_mask] = ( 

800 np.full_like(region[self.qm_selection_mask], "QM")) 

801 

802 buffer_only_mask = self.qm_buffer_mask & ~self.qm_selection_mask 

803 

804 region[buffer_only_mask] = np.full_like(region[buffer_only_mask], 

805 "buffer") 

806 

807 if print_mapping: 

808 

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}") 

813 

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}") 

819 

820 return region 

821 

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" 

828 

829 self.qm_buffer_mask = self.qm_selection_mask ^ buffer_mask 

830 

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 

842 

843 region = self.get_region_from_masks(atoms=atoms) 

844 

845 atoms_copy = atoms.copy() 

846 atoms_copy.new_array("region", region) 

847 

848 atoms_copy.calc = self # to keep the calculation results 

849 

850 atoms_copy.write(filename, format='extxyz') 

851 

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 

861 

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. 

868 

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 """ 

874 

875 from ase.io import read 

876 atoms = read(filename, format='extxyz') 

877 

878 if "region" in atoms.arrays: 

879 region = atoms.get_array("region") 

880 else: 

881 raise RuntimeError("Please provide extxyz file with 'region' array") 

882 

883 dummy_qm_mask = np.full_like(atoms, False, dtype=bool) 

884 dummy_qm_mask[0] = True 

885 

886 self = cls(atoms, dummy_qm_mask, 

887 qm_calc, mm_calc, buffer_width=1.0) 

888 

889 self.set_masks_from_region(region) 

890 

891 return self