Coverage for /builds/debichem-team/python-ase/ase/vibrations/vibrations.py : 88.93%
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
1"""A class for computing vibrational modes"""
3from math import pi, sqrt, log
4import sys
6import numpy as np
7from pathlib import Path
9import ase.units as units
10import ase.io
11from ase.parallel import world, paropen
13from ase.utils.filecache import get_json_cache
14from .data import VibrationsData
16from collections import namedtuple
19class AtomicDisplacements:
20 def _disp(self, a, i, step):
21 if isinstance(i, str): # XXX Simplify by removing this.
22 i = 'xyz'.index(i)
23 return Displacement(a, i, np.sign(step), abs(step), self)
25 def _eq_disp(self):
26 return self._disp(0, 0, 0)
28 @property
29 def ndof(self):
30 return 3 * len(self.indices)
33class Displacement(namedtuple('Displacement', ['a', 'i', 'sign', 'ndisp',
34 'vib'])):
35 @property
36 def name(self):
37 if self.sign == 0:
38 return 'eq'
40 axisname = 'xyz'[self.i]
41 dispname = self.ndisp * ' +-'[self.sign]
42 return f'{self.a}{axisname}{dispname}'
44 @property
45 def _cached(self):
46 return self.vib.cache[self.name]
48 def forces(self):
49 return self._cached['forces'].copy()
51 @property
52 def step(self):
53 return self.ndisp * self.sign * self.vib.delta
55 # XXX dipole only valid for infrared
56 def dipole(self):
57 return self._cached['dipole'].copy()
59 # XXX below stuff only valid for TDDFT excitation stuff
60 def save_ov_nn(self, ov_nn):
61 np.save(self.name + '.ov', ov_nn)
63 def load_ov_nn(self):
64 return np.load(self.name + '.ov.npy')
66 @property
67 def _exname(self):
68 return Path(self.vib.exname) / f'ex.{self.name}{self.vib.exext}'
70 def calculate_and_save_static_polarizability(self, atoms):
71 exobj = self.vib._new_exobj()
72 excitation_data = exobj.calculate(atoms)
73 np.savetxt(self._exname, excitation_data)
75 def load_static_polarizability(self):
76 return np.loadtxt(self._exname)
78 def read_exobj(self):
79 # XXX each exobj should allow for self._exname as Path
80 return self.vib.read_exobj(str(self._exname))
82 def calculate_and_save_exlist(self, atoms):
83 # exo = self.vib._new_exobj()
84 excalc = self.vib._new_exobj()
85 exlist = excalc.calculate(atoms)
86 # XXX each exobj should allow for self._exname as Path
87 exlist.write(str(self._exname))
90class Vibrations(AtomicDisplacements):
91 """Class for calculating vibrational modes using finite difference.
93 The vibrational modes are calculated from a finite difference
94 approximation of the Hessian matrix.
96 The *summary()*, *get_energies()* and *get_frequencies()* methods all take
97 an optional *method* keyword. Use method='Frederiksen' to use the method
98 described in:
100 T. Frederiksen, M. Paulsson, M. Brandbyge, A. P. Jauho:
101 "Inelastic transport theory from first-principles: methodology and
102 applications for nanoscale devices", Phys. Rev. B 75, 205413 (2007)
104 atoms: Atoms object
105 The atoms to work on.
106 indices: list of int
107 List of indices of atoms to vibrate. Default behavior is
108 to vibrate all atoms.
109 name: str
110 Name to use for files.
111 delta: float
112 Magnitude of displacements.
113 nfree: int
114 Number of displacements per atom and cartesian coordinate, 2 and 4 are
115 supported. Default is 2 which will displace each atom +delta and
116 -delta for each cartesian coordinate.
118 Example:
120 >>> from ase import Atoms
121 >>> from ase.calculators.emt import EMT
122 >>> from ase.optimize import BFGS
123 >>> from ase.vibrations import Vibrations
124 >>> n2 = Atoms('N2', [(0, 0, 0), (0, 0, 1.1)],
125 ... calculator=EMT())
126 >>> BFGS(n2).run(fmax=0.01)
127 BFGS: 0 16:01:21 0.440339 3.2518
128 BFGS: 1 16:01:21 0.271928 0.8211
129 BFGS: 2 16:01:21 0.263278 0.1994
130 BFGS: 3 16:01:21 0.262777 0.0088
131 >>> vib = Vibrations(n2)
132 >>> vib.run()
133 >>> vib.summary()
134 ---------------------
135 # meV cm^-1
136 ---------------------
137 0 0.0 0.0
138 1 0.0 0.0
139 2 0.0 0.0
140 3 1.4 11.5
141 4 1.4 11.5
142 5 152.7 1231.3
143 ---------------------
144 Zero-point energy: 0.078 eV
145 >>> vib.write_mode(-1) # write last mode to trajectory file
147 """
149 def __init__(self, atoms, indices=None, name='vib', delta=0.01, nfree=2):
150 assert nfree in [2, 4]
151 self.atoms = atoms
152 self.calc = atoms.calc
153 if indices is None:
154 indices = range(len(atoms))
155 if len(indices) != len(set(indices)):
156 raise ValueError(
157 'one (or more) indices included more than once')
158 self.indices = np.asarray(indices)
160 self.delta = delta
161 self.nfree = nfree
162 self.H = None
163 self.ir = None
164 self._vibrations = None
166 self.cache = get_json_cache(name)
168 @property
169 def name(self):
170 return str(self.cache.directory)
172 def run(self):
173 """Run the vibration calculations.
175 This will calculate the forces for 6 displacements per atom +/-x,
176 +/-y, +/-z. Only those calculations that are not already done will be
177 started. Be aware that an interrupted calculation may produce an empty
178 file (ending with .json), which must be deleted before restarting the
179 job. Otherwise the forces will not be calculated for that
180 displacement.
182 Note that the calculations for the different displacements can be done
183 simultaneously by several independent processes. This feature relies
184 on the existence of files and the subsequent creation of the file in
185 case it is not found.
187 If the program you want to use does not have a calculator in ASE, use
188 ``iterdisplace`` to get all displaced structures and calculate the
189 forces on your own.
190 """
192 if not self.cache.writable:
193 raise RuntimeError(
194 'Cannot run calculation. '
195 'Cache must be removed or split in order '
196 'to have only one sort of data structure at a time.')
198 self._check_old_pickles()
200 for disp, atoms in self.iterdisplace(inplace=True):
201 with self.cache.lock(disp.name) as handle:
202 if handle is None:
203 continue
205 result = self.calculate(atoms, disp)
207 if world.rank == 0:
208 handle.save(result)
210 def _check_old_pickles(self):
211 from pathlib import Path
212 eq_pickle_path = Path(f'{self.name}.eq.pckl')
213 pickle2json_instructions = f"""\
214Found old pickle files such as {eq_pickle_path}. \
215Please remove them and recalculate or run \
216"python -m ase.vibrations.pickle2json --help"."""
217 if len(self.cache) == 0 and eq_pickle_path.exists():
218 raise RuntimeError(pickle2json_instructions)
220 def iterdisplace(self, inplace=False):
221 """Yield name and atoms object for initial and displaced structures.
223 Use this to export the structures for each single-point calculation
224 to an external program instead of using ``run()``. Then save the
225 calculated gradients to <name>.json and continue using this instance.
226 """
227 # XXX change of type of disp
228 atoms = self.atoms if inplace else self.atoms.copy()
229 displacements = self.displacements()
230 eq_disp = next(displacements)
231 assert eq_disp.name == 'eq'
232 yield eq_disp, atoms
234 for disp in displacements:
235 if not inplace:
236 atoms = self.atoms.copy()
237 pos0 = atoms.positions[disp.a, disp.i]
238 atoms.positions[disp.a, disp.i] += disp.step
239 yield disp, atoms
241 if inplace:
242 atoms.positions[disp.a, disp.i] = pos0
244 def iterimages(self):
245 """Yield initial and displaced structures."""
246 for name, atoms in self.iterdisplace():
247 yield atoms
249 def _iter_ai(self):
250 for a in self.indices:
251 for i in range(3):
252 yield a, i
254 def displacements(self):
255 yield self._eq_disp()
257 for a, i in self._iter_ai():
258 for sign in [-1, 1]:
259 for ndisp in range(1, self.nfree // 2 + 1):
260 yield self._disp(a, i, sign * ndisp)
262 def calculate(self, atoms, disp):
263 results = {}
264 results['forces'] = self.calc.get_forces(atoms)
266 if self.ir:
267 results['dipole'] = self.calc.get_dipole_moment(atoms)
269 return results
271 def clean(self, empty_files=False, combined=True):
272 """Remove json-files.
274 Use empty_files=True to remove only empty files and
275 combined=False to not remove the combined file.
277 """
279 if world.rank != 0:
280 return 0
282 if empty_files:
283 return self.cache.strip_empties() # XXX Fails on combined cache
285 nfiles = self.cache.filecount()
286 self.cache.clear()
287 return nfiles
289 def combine(self):
290 """Combine json-files to one file ending with '.all.json'.
292 The other json-files will be removed in order to have only one sort
293 of data structure at a time.
295 """
296 nelements_before = self.cache.filecount()
297 self.cache = self.cache.combine()
298 return nelements_before
300 def split(self):
301 """Split combined json-file.
303 The combined json-file will be removed in order to have only one
304 sort of data structure at a time.
306 """
307 count = self.cache.filecount()
308 self.cache = self.cache.split()
309 return count
311 def read(self, method='standard', direction='central'):
312 self.method = method.lower()
313 self.direction = direction.lower()
314 assert self.method in ['standard', 'frederiksen']
315 assert self.direction in ['central', 'forward', 'backward']
317 n = 3 * len(self.indices)
318 H = np.empty((n, n))
319 r = 0
321 eq_disp = self._eq_disp()
323 if direction != 'central':
324 feq = eq_disp.forces()
326 for a, i in self._iter_ai():
327 disp_minus = self._disp(a, i, -1)
328 disp_plus = self._disp(a, i, 1)
330 fminus = disp_minus.forces()
331 fplus = disp_plus.forces()
332 if self.method == 'frederiksen':
333 fminus[a] -= fminus.sum(0)
334 fplus[a] -= fplus.sum(0)
335 if self.nfree == 4:
336 fminusminus = self._disp(a, i, -2).forces()
337 fplusplus = self._disp(a, i, 2).forces()
338 if self.method == 'frederiksen':
339 fminusminus[a] -= fminusminus.sum(0)
340 fplusplus[a] -= fplusplus.sum(0)
341 if self.direction == 'central':
342 if self.nfree == 2:
343 H[r] = .5 * (fminus - fplus)[self.indices].ravel()
344 else:
345 assert self.nfree == 4
346 H[r] = H[r] = (-fminusminus +
347 8 * fminus -
348 8 * fplus +
349 fplusplus)[self.indices].ravel() / 12.0
350 elif self.direction == 'forward':
351 H[r] = (feq - fplus)[self.indices].ravel()
352 else:
353 assert self.direction == 'backward'
354 H[r] = (fminus - feq)[self.indices].ravel()
355 H[r] /= 2 * self.delta
356 r += 1
357 H += H.copy().T
358 self.H = H
359 masses = self.atoms.get_masses()
360 if any(masses[self.indices] == 0):
361 raise RuntimeError('Zero mass encountered in one or more of '
362 'the vibrated atoms. Use Atoms.set_masses()'
363 ' to set all masses to non-zero values.')
365 self.im = np.repeat(masses[self.indices]**-0.5, 3)
366 self._vibrations = self.get_vibrations(read_cache=False)
368 omega2, modes = np.linalg.eigh(self.im[:, None] * H * self.im)
369 self.modes = modes.T.copy()
371 # Conversion factor:
372 s = units._hbar * 1e10 / sqrt(units._e * units._amu)
373 self.hnu = s * omega2.astype(complex)**0.5
375 def get_vibrations(self, method='standard', direction='central',
376 read_cache=True, **kw):
377 """Get vibrations as VibrationsData object
379 If read() has not yet been called, this will be called to assemble data
380 from the outputs of run(). Most of the arguments to this function are
381 options to be passed to read() in this case.
383 Args:
384 method (str): Calculation method passed to read()
385 direction (str): Finite-difference scheme passed to read()
386 read_cache (bool): The VibrationsData object will be cached for
387 quick access. Set False to force regeneration of the cache with
388 the current atoms/Hessian/indices data.
389 **kw: Any remaining keyword arguments are passed to read()
391 Returns:
392 VibrationsData
394 """
395 if read_cache and (self._vibrations is not None):
396 return self._vibrations
398 else:
399 if (self.H is None or method.lower() != self.method or
400 direction.lower() != self.direction):
401 self.read(method, direction, **kw)
403 return VibrationsData.from_2d(self.atoms, self.H,
404 indices=self.indices)
406 def get_energies(self, method='standard', direction='central', **kw):
407 """Get vibration energies in eV."""
408 return self.get_vibrations(method=method,
409 direction=direction, **kw).get_energies()
411 def get_frequencies(self, method='standard', direction='central'):
412 """Get vibration frequencies in cm^-1."""
413 return self.get_vibrations(method=method,
414 direction=direction).get_frequencies()
416 def summary(self, method='standard', direction='central', freq=None,
417 log=sys.stdout):
418 if freq is not None:
419 energies = freq * units.invcm
420 else:
421 energies = self.get_energies(method=method, direction=direction)
423 summary_lines = VibrationsData._tabulate_from_energies(energies)
424 log_text = '\n'.join(summary_lines) + '\n'
426 if isinstance(log, str):
427 with paropen(log, 'a') as log_file:
428 log_file.write(log_text)
429 else:
430 log.write(log_text)
432 def get_zero_point_energy(self, freq=None):
433 if freq:
434 raise NotImplementedError()
435 return self.get_vibrations().get_zero_point_energy()
437 def get_mode(self, n):
438 """Get mode number ."""
439 return self.get_vibrations().get_modes(all_atoms=True)[n]
441 def write_mode(self, n=None, kT=units.kB * 300, nimages=30):
442 """Write mode number n to trajectory file. If n is not specified,
443 writes all non-zero modes."""
444 if n is None:
445 for index, energy in enumerate(self.get_energies()):
446 if abs(energy) > 1e-5:
447 self.write_mode(n=index, kT=kT, nimages=nimages)
448 return
450 else:
451 n %= len(self.get_energies())
453 with ase.io.Trajectory('%s.%d.traj' % (self.name, n), 'w') as traj:
454 for image in (self.get_vibrations()
455 .iter_animated_mode(n,
456 temperature=kT, frames=nimages)):
457 traj.write(image)
459 def show_as_force(self, n, scale=0.2, show=True):
460 return self.get_vibrations().show_as_force(n, scale=scale, show=show)
462 def write_jmol(self):
463 """Writes file for viewing of the modes with jmol."""
465 with open(self.name + '.xyz', 'w') as fd:
466 self._write_jmol(fd)
468 def _write_jmol(self, fd):
469 symbols = self.atoms.get_chemical_symbols()
470 freq = self.get_frequencies()
471 for n in range(3 * len(self.indices)):
472 fd.write('%6d\n' % len(self.atoms))
474 if freq[n].imag != 0:
475 c = 'i'
476 freq[n] = freq[n].imag
478 else:
479 freq[n] = freq[n].real
480 c = ' '
482 fd.write('Mode #%d, f = %.1f%s cm^-1'
483 % (n, float(freq[n].real), c))
485 if self.ir:
486 fd.write(', I = %.4f (D/Å)^2 amu^-1.\n' % self.intensities[n])
487 else:
488 fd.write('.\n')
490 mode = self.get_mode(n)
491 for i, pos in enumerate(self.atoms.positions):
492 fd.write('%2s %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n' %
493 (symbols[i], pos[0], pos[1], pos[2],
494 mode[i, 0], mode[i, 1], mode[i, 2]))
496 def fold(self, frequencies, intensities,
497 start=800.0, end=4000.0, npts=None, width=4.0,
498 type='Gaussian', normalize=False):
499 """Fold frequencies and intensities within the given range
500 and folding method (Gaussian/Lorentzian).
501 The energy unit is cm^-1.
502 normalize=True ensures the integral over the peaks to give the
503 intensity.
504 """
506 lctype = type.lower()
507 assert lctype in ['gaussian', 'lorentzian']
508 if not npts:
509 npts = int((end - start) / width * 10 + 1)
510 prefactor = 1
511 if lctype == 'lorentzian':
512 intensities = intensities * width * pi / 2.
513 if normalize:
514 prefactor = 2. / width / pi
515 else:
516 sigma = width / 2. / sqrt(2. * log(2.))
517 if normalize:
518 prefactor = 1. / sigma / sqrt(2 * pi)
520 # Make array with spectrum data
521 spectrum = np.empty(npts)
522 energies = np.linspace(start, end, npts)
523 for i, energy in enumerate(energies):
524 energies[i] = energy
525 if lctype == 'lorentzian':
526 spectrum[i] = (intensities * 0.5 * width / pi /
527 ((frequencies - energy)**2 +
528 0.25 * width**2)).sum()
529 else:
530 spectrum[i] = (intensities *
531 np.exp(-(frequencies - energy)**2 /
532 2. / sigma**2)).sum()
533 return [energies, prefactor * spectrum]
535 def write_dos(self, out='vib-dos.dat', start=800, end=4000,
536 npts=None, width=10,
537 type='Gaussian', method='standard', direction='central'):
538 """Write out the vibrational density of states to file.
540 First column is the wavenumber in cm^-1, the second column the
541 folded vibrational density of states.
542 Start and end points, and width of the Gaussian/Lorentzian
543 should be given in cm^-1."""
544 frequencies = self.get_frequencies(method, direction).real
545 intensities = np.ones(len(frequencies))
546 energies, spectrum = self.fold(frequencies, intensities,
547 start, end, npts, width, type)
549 # Write out spectrum in file.
550 outdata = np.empty([len(energies), 2])
551 outdata.T[0] = energies
552 outdata.T[1] = spectrum
554 with open(out, 'w') as fd:
555 fd.write('# %s folded, width=%g cm^-1\n' % (type.title(), width))
556 fd.write('# [cm^-1] arbitrary\n')
557 for row in outdata:
558 fd.write('%.3f %15.5e\n' %
559 (row[0], row[1]))