Coverage for /builds/kinetik161/ase/ase/vibrations/vibrations.py: 88.93%

1"""A class for computing vibrational modes"""

3import sys

4from collections import namedtuple

5from math import log, pi, sqrt

6from pathlib import Path

8import numpy as np

10import ase.io

11import ase.units as units

12from ase.parallel import paropen, world

13from ase.utils.filecache import get_json_cache

15from .data import VibrationsData

18class AtomicDisplacements:

19 def _disp(self, a, i, step):

20 if isinstance(i, str): # XXX Simplify by removing this.

21 i = 'xyz'.index(i)

22 return Displacement(a, i, np.sign(step), abs(step), self)

24 def _eq_disp(self):

25 return self._disp(0, 0, 0)

27 @property

28 def ndof(self):

29 return 3 * len(self.indices)

32class Displacement(namedtuple('Displacement', ['a', 'i', 'sign', 'ndisp',

33 'vib'])):

34 @property

35 def name(self):

36 if self.sign == 0:

37 return 'eq'

39 axisname = 'xyz'[self.i]

40 dispname = self.ndisp * ' +-'[self.sign]

41 return f'{self.a}{axisname}{dispname}'

43 @property

44 def _cached(self):

45 return self.vib.cache[self.name]

47 def forces(self):

48 return self._cached['forces'].copy()

50 @property

51 def step(self):

52 return self.ndisp * self.sign * self.vib.delta

54 # XXX dipole only valid for infrared

55 def dipole(self):

56 return self._cached['dipole'].copy()

58 # XXX below stuff only valid for TDDFT excitation stuff

59 def save_ov_nn(self, ov_nn):

60 np.save(Path(self.vib.exname) / (self.name + '.ov'), ov_nn)

62 def load_ov_nn(self):

63 return np.load(Path(self.vib.exname) / (self.name + '.ov.npy'))

65 @property

66 def _exname(self):

67 return Path(self.vib.exname) / f'ex.{self.name}{self.vib.exext}'

69 def calculate_and_save_static_polarizability(self, atoms):

70 exobj = self.vib._new_exobj()

71 excitation_data = exobj(atoms)

72 np.savetxt(self._exname, excitation_data)

74 def load_static_polarizability(self):

75 return np.loadtxt(self._exname)

77 def read_exobj(self):

78 # XXX each exobj should allow for self._exname as Path

79 return self.vib.read_exobj(str(self._exname))

81 def calculate_and_save_exlist(self, atoms):

82 # exo = self.vib._new_exobj()

83 excalc = self.vib._new_exobj()

84 exlist = excalc.calculate(atoms)

85 # XXX each exobj should allow for self._exname as Path

86 exlist.write(str(self._exname))

89class Vibrations(AtomicDisplacements):

90 """Class for calculating vibrational modes using finite difference.

92 The vibrational modes are calculated from a finite difference

93 approximation of the Hessian matrix.

95 The *summary()*, *get_energies()* and *get_frequencies()* methods all take

96 an optional *method* keyword. Use method='Frederiksen' to use the method

97 described in:

99 T. Frederiksen, M. Paulsson, M. Brandbyge, A. P. Jauho:

100 "Inelastic transport theory from first-principles: methodology and

101 applications for nanoscale devices", Phys. Rev. B 75, 205413 (2007)

102

103 atoms: Atoms object

104 The atoms to work on.

105 indices: list of int

106 List of indices of atoms to vibrate. Default behavior is

107 to vibrate all atoms.

108 name: str

109 Name to use for files.

110 delta: float

111 Magnitude of displacements.

112 nfree: int

113 Number of displacements per atom and cartesian coordinate, 2 and 4 are

114 supported. Default is 2 which will displace each atom +delta and

115 -delta for each cartesian coordinate.

116

117 Example:

118

119 >>> from ase import Atoms

120 >>> from ase.calculators.emt import EMT

121 >>> from ase.optimize import BFGS

122 >>> from ase.vibrations import Vibrations

123 >>> n2 = Atoms('N2', [(0, 0, 0), (0, 0, 1.1)],

124 ... calculator=EMT())

125 >>> BFGS(n2).run(fmax=0.01)

126 BFGS: 0 16:01:21 0.440339 3.2518

127 BFGS: 1 16:01:21 0.271928 0.8211

128 BFGS: 2 16:01:21 0.263278 0.1994

129 BFGS: 3 16:01:21 0.262777 0.0088

130 >>> vib = Vibrations(n2)

131 >>> vib.run()

132 >>> vib.summary()

133 ---------------------

134 # meV cm^-1

135 ---------------------

136 0 0.0 0.0

137 1 0.0 0.0

138 2 0.0 0.0

139 3 1.4 11.5

140 4 1.4 11.5

141 5 152.7 1231.3

142 ---------------------

143 Zero-point energy: 0.078 eV

144 >>> vib.write_mode(-1) # write last mode to trajectory file

145

146 """

147

148 def __init__(self, atoms, indices=None, name='vib', delta=0.01, nfree=2):

149 assert nfree in [2, 4]

150 self.atoms = atoms

151 self.calc = atoms.calc

152 if indices is None:

153 indices = range(len(atoms))

154 if len(indices) != len(set(indices)):

155 raise ValueError(

156 'one (or more) indices included more than once')

157 self.indices = np.asarray(indices)

158

159 self.delta = delta

160 self.nfree = nfree

161 self.H = None

162 self.ir = None

163 self._vibrations = None

164

165 self.cache = get_json_cache(name)

166

167 @property

168 def name(self):

169 return str(self.cache.directory)

170

171 def run(self):

172 """Run the vibration calculations.

173

174 This will calculate the forces for 6 displacements per atom +/-x,

175 +/-y, +/-z. Only those calculations that are not already done will be

176 started. Be aware that an interrupted calculation may produce an empty

177 file (ending with .json), which must be deleted before restarting the

178 job. Otherwise the forces will not be calculated for that

179 displacement.

180

181 Note that the calculations for the different displacements can be done

182 simultaneously by several independent processes. This feature relies

183 on the existence of files and the subsequent creation of the file in

184 case it is not found.

185

186 If the program you want to use does not have a calculator in ASE, use

187 ``iterdisplace`` to get all displaced structures and calculate the

188 forces on your own.

189 """

190

191 if not self.cache.writable:

192 raise RuntimeError(

193 'Cannot run calculation. '

194 'Cache must be removed or split in order '

195 'to have only one sort of data structure at a time.')

196

197 self._check_old_pickles()

198

199 for disp, atoms in self.iterdisplace(inplace=True):

200 with self.cache.lock(disp.name) as handle:

201 if handle is None:

202 continue

203

204 result = self.calculate(atoms, disp)

205

206 if world.rank == 0:

207 handle.save(result)

208

209 def _check_old_pickles(self):

210 from pathlib import Path

211 eq_pickle_path = Path(f'{self.name}.eq.pckl')

212 pickle2json_instructions = f"""\

213Found old pickle files such as {eq_pickle_path}. \

214Please remove them and recalculate or run \

215"python -m ase.vibrations.pickle2json --help"."""

216 if len(self.cache) == 0 and eq_pickle_path.exists():

217 raise RuntimeError(pickle2json_instructions)

218

219 def iterdisplace(self, inplace=False):

220 """Yield name and atoms object for initial and displaced structures.

221

222 Use this to export the structures for each single-point calculation

223 to an external program instead of using ``run()``. Then save the

224 calculated gradients to <name>.json and continue using this instance.

225 """

226 # XXX change of type of disp

227 atoms = self.atoms if inplace else self.atoms.copy()

228 displacements = self.displacements()

229 eq_disp = next(displacements)

230 assert eq_disp.name == 'eq'

231 yield eq_disp, atoms

232

233 for disp in displacements:

234 if not inplace:

235 atoms = self.atoms.copy()

236 pos0 = atoms.positions[disp.a, disp.i]

237 atoms.positions[disp.a, disp.i] += disp.step

238 yield disp, atoms

239

240 if inplace:

241 atoms.positions[disp.a, disp.i] = pos0

242

243 def iterimages(self):

244 """Yield initial and displaced structures."""

245 for name, atoms in self.iterdisplace():

246 yield atoms

247

248 def _iter_ai(self):

249 for a in self.indices:

250 for i in range(3):

251 yield a, i

252

253 def displacements(self):

254 yield self._eq_disp()

255

256 for a, i in self._iter_ai():

257 for sign in [-1, 1]:

258 for ndisp in range(1, self.nfree // 2 + 1):

259 yield self._disp(a, i, sign * ndisp)

260

261 def calculate(self, atoms, disp):

262 results = {}

263 results['forces'] = self.calc.get_forces(atoms)

264

265 if self.ir:

266 results['dipole'] = self.calc.get_dipole_moment(atoms)

267

268 return results

269

270 def clean(self, empty_files=False, combined=True):

271 """Remove json-files.

272

273 Use empty_files=True to remove only empty files and

274 combined=False to not remove the combined file.

275

276 """

277

278 if world.rank != 0:

279 return 0

280

281 if empty_files:

282 return self.cache.strip_empties() # XXX Fails on combined cache

283

284 nfiles = self.cache.filecount()

285 self.cache.clear()

286 return nfiles

287

288 def combine(self):

289 """Combine json-files to one file ending with '.all.json'.

290

291 The other json-files will be removed in order to have only one sort

292 of data structure at a time.

293

294 """

295 nelements_before = self.cache.filecount()

296 self.cache = self.cache.combine()

297 return nelements_before

298

299 def split(self):

300 """Split combined json-file.

301

302 The combined json-file will be removed in order to have only one

303 sort of data structure at a time.

304

305 """

306 count = self.cache.filecount()

307 self.cache = self.cache.split()

308 return count

309

310 def read(self, method='standard', direction='central'):

311 self.method = method.lower()

312 self.direction = direction.lower()

313 assert self.method in ['standard', 'frederiksen']

314 assert self.direction in ['central', 'forward', 'backward']

315

316 n = 3 * len(self.indices)

317 H = np.empty((n, n))

318 r = 0

319

320 eq_disp = self._eq_disp()

321

322 if direction != 'central':

323 feq = eq_disp.forces()

324

325 for a, i in self._iter_ai():

326 disp_minus = self._disp(a, i, -1)

327 disp_plus = self._disp(a, i, 1)

328

329 fminus = disp_minus.forces()

330 fplus = disp_plus.forces()

331 if self.method == 'frederiksen':

332 fminus[a] -= fminus.sum(0)

333 fplus[a] -= fplus.sum(0)

334 if self.nfree == 4:

335 fminusminus = self._disp(a, i, -2).forces()

336 fplusplus = self._disp(a, i, 2).forces()

337 if self.method == 'frederiksen':

338 fminusminus[a] -= fminusminus.sum(0)

339 fplusplus[a] -= fplusplus.sum(0)

340 if self.direction == 'central':

341 if self.nfree == 2:

342 H[r] = .5 * (fminus - fplus)[self.indices].ravel()

343 else:

344 assert self.nfree == 4

345 H[r] = H[r] = (-fminusminus +

346 8 * fminus -

347 8 * fplus +

348 fplusplus)[self.indices].ravel() / 12.0

349 elif self.direction == 'forward':

350 H[r] = (feq - fplus)[self.indices].ravel()

351 else:

352 assert self.direction == 'backward'

353 H[r] = (fminus - feq)[self.indices].ravel()

354 H[r] /= 2 * self.delta

355 r += 1

356 H += H.copy().T

357 self.H = H

358 masses = self.atoms.get_masses()

359 if any(masses[self.indices] == 0):

360 raise RuntimeError('Zero mass encountered in one or more of '

361 'the vibrated atoms. Use Atoms.set_masses()'

362 ' to set all masses to non-zero values.')

363

364 self.im = np.repeat(masses[self.indices]**-0.5, 3)

365 self._vibrations = self.get_vibrations(read_cache=False,

366 method=self.method,

367 direction=self.direction)

368

369 omega2, modes = np.linalg.eigh(self.im[:, None] * H * self.im)

370 self.modes = modes.T.copy()

371

372 # Conversion factor:

373 s = units._hbar * 1e10 / sqrt(units._e * units._amu)

374 self.hnu = s * omega2.astype(complex)**0.5

375

376 def get_vibrations(self, method='standard', direction='central',

377 read_cache=True, **kw):

378 """Get vibrations as VibrationsData object

379

380 If read() has not yet been called, this will be called to assemble data

381 from the outputs of run(). Most of the arguments to this function are

382 options to be passed to read() in this case.

383

384 Args:

385 method (str): Calculation method passed to read()

386 direction (str): Finite-difference scheme passed to read()

387 read_cache (bool): The VibrationsData object will be cached for

388 quick access. Set False to force regeneration of the cache with

389 the current atoms/Hessian/indices data.

390 **kw: Any remaining keyword arguments are passed to read()

391

392 Returns:

393 VibrationsData

394

395 """

396 if read_cache and (self._vibrations is not None):

397 return self._vibrations

398

399 else:

400 if (self.H is None or method.lower() != self.method or

401 direction.lower() != self.direction):

402 self.read(method, direction, **kw)

403

404 return VibrationsData.from_2d(self.atoms, self.H,

405 indices=self.indices)

406

407 def get_energies(self, method='standard', direction='central', **kw):

408 """Get vibration energies in eV."""

409 return self.get_vibrations(method=method,

410 direction=direction, **kw).get_energies()

411

412 def get_frequencies(self, method='standard', direction='central'):

413 """Get vibration frequencies in cm^-1."""

414 return self.get_vibrations(method=method,

415 direction=direction).get_frequencies()

416

417 def summary(self, method='standard', direction='central', freq=None,

418 log=sys.stdout):

419 if freq is not None:

420 energies = freq * units.invcm

421 else:

422 energies = self.get_energies(method=method, direction=direction)

423

424 summary_lines = VibrationsData._tabulate_from_energies(energies)

425 log_text = '\n'.join(summary_lines) + '\n'

426

427 if isinstance(log, str):

428 with paropen(log, 'a') as log_file:

429 log_file.write(log_text)

430 else:

431 log.write(log_text)

432

433 def get_zero_point_energy(self, freq=None):

434 if freq:

435 raise NotImplementedError()

436 return self.get_vibrations().get_zero_point_energy()

437

438 def get_mode(self, n):

439 """Get mode number ."""

440 return self.get_vibrations().get_modes(all_atoms=True)[n]

441

442 def write_mode(self, n=None, kT=units.kB * 300, nimages=30):

443 """Write mode number n to trajectory file. If n is not specified,

444 writes all non-zero modes."""

445 if n is None:

446 for index, energy in enumerate(self.get_energies()):

447 if abs(energy) > 1e-5:

448 self.write_mode(n=index, kT=kT, nimages=nimages)

449 return

450

451 else:

452 n %= len(self.get_energies())

453

454 with ase.io.Trajectory('%s.%d.traj' % (self.name, n), 'w') as traj:

455 for image in (self.get_vibrations()

456 .iter_animated_mode(n,

457 temperature=kT, frames=nimages)):

458 traj.write(image)

459

460 def show_as_force(self, n, scale=0.2, show=True):

461 return self.get_vibrations().show_as_force(n, scale=scale, show=show)

462

463 def write_jmol(self):

464 """Writes file for viewing of the modes with jmol."""

465

466 with open(self.name + '.xyz', 'w') as fd:

467 self._write_jmol(fd)

468

469 def _write_jmol(self, fd):

470 symbols = self.atoms.get_chemical_symbols()

471 freq = self.get_frequencies()

472 for n in range(3 * len(self.indices)):

473 fd.write('%6d\n' % len(self.atoms))

474

475 if freq[n].imag != 0:

476 c = 'i'

477 freq[n] = freq[n].imag

478

479 else:

480 freq[n] = freq[n].real

481 c = ' '

482

483 fd.write('Mode #%d, f = %.1f%s cm^-1'

484 % (n, float(freq[n].real), c))

485

486 if self.ir:

487 fd.write(', I = %.4f (D/Å)^2 amu^-1.\n' % self.intensities[n])

488 else:

489 fd.write('.\n')

490

491 mode = self.get_mode(n)

492 for i, pos in enumerate(self.atoms.positions):

493 fd.write('%2s %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n' %

494 (symbols[i], pos[0], pos[1], pos[2],

495 mode[i, 0], mode[i, 1], mode[i, 2]))

496

497 def fold(self, frequencies, intensities,

498 start=800.0, end=4000.0, npts=None, width=4.0,

499 type='Gaussian', normalize=False):

500 """Fold frequencies and intensities within the given range

501 and folding method (Gaussian/Lorentzian).

502 The energy unit is cm^-1.

503 normalize=True ensures the integral over the peaks to give the

504 intensity.

505 """

506

507 lctype = type.lower()

508 assert lctype in ['gaussian', 'lorentzian']

509 if not npts:

510 npts = int((end - start) / width * 10 + 1)

511 prefactor = 1

512 if lctype == 'lorentzian':

513 intensities = intensities * width * pi / 2.

514 if normalize:

515 prefactor = 2. / width / pi

516 else:

517 sigma = width / 2. / sqrt(2. * log(2.))

518 if normalize:

519 prefactor = 1. / sigma / sqrt(2 * pi)

520

521 # Make array with spectrum data

522 spectrum = np.empty(npts)

523 energies = np.linspace(start, end, npts)

524 for i, energy in enumerate(energies):

525 energies[i] = energy

526 if lctype == 'lorentzian':

527 spectrum[i] = (intensities * 0.5 * width / pi /

528 ((frequencies - energy)**2 +

529 0.25 * width**2)).sum()

530 else:

531 spectrum[i] = (intensities *

532 np.exp(-(frequencies - energy)**2 /

533 2. / sigma**2)).sum()

534 return [energies, prefactor * spectrum]

535

536 def write_dos(self, out='vib-dos.dat', start=800, end=4000,

537 npts=None, width=10,

538 type='Gaussian', method='standard', direction='central'):

539 """Write out the vibrational density of states to file.

540

541 First column is the wavenumber in cm^-1, the second column the

542 folded vibrational density of states.

543 Start and end points, and width of the Gaussian/Lorentzian

544 should be given in cm^-1."""

545 frequencies = self.get_frequencies(method, direction).real

546 intensities = np.ones(len(frequencies))

547 energies, spectrum = self.fold(frequencies, intensities,

548 start, end, npts, width, type)

549

550 # Write out spectrum in file.

551 outdata = np.empty([len(energies), 2])

552 outdata.T[0] = energies

553 outdata.T[1] = spectrum

554

555 with open(out, 'w') as fd:

556 fd.write(f'# {type.title()} folded, width={width:g} cm^-1\n')

557 fd.write('# [cm^-1] arbitrary\n')

558 for row in outdata:

559 fd.write('%.3f %15.5e\n' %

560 (row[0], row[1]))