Coverage for strongcoca/response/nwchem.py: 99%
178 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-26 18:44 +0000
1import logging
2import os
3from collections import OrderedDict
4from typing import Any, Dict, IO, List, Union
6import numpy as np
7from ase import Atom, Atoms
8from ase.calculators.singlepoint import (SinglePointDFTCalculator,
9 SinglePointKPoint)
11from ..units import eV_to_au, au_to_eV, au_to_eVA
12from ..types import Path_str, Vector
13from .casida import CasidaResponse
14from .utilities import Broadening, NoArtificialBroadening
17logger = logging.getLogger(__name__)
20def read_nwchem_casida(filename: Path_str,
21 broadening: Broadening = NoArtificialBroadening(),
22 name: str = 'NWChem Casida Response') -> CasidaResponse:
23 """Returns a freshly initialized :class:`~strongcoca.response.CasidaResponse` object
24 using an NWChem output file to set up the components of the Casida representation.
26 Parameters
27 ----------
28 filename
29 Name of NWChem output file; can be compressed.
30 broadening
31 Artificial broadening used for the continuous response;
32 defaults to no broadening.
33 name
34 Name of subsystem.
35 """
36 logger.debug(f'Entering read_nwchem_casida: filename= {filename}')
37 try:
38 atoms = _read_nwchem_logfile(filename)
39 except Exception as e:
40 raise Exception(f'Failed to parse NWChem output file {filename}\n{str(e)}')
41 if 'excited_states' not in atoms.info or len(atoms.info['excited_states']) is False:
42 raise ValueError(f'NWChem output file {filename} does'
43 ' not contain excited states information')
45 # include only singlet excitations
46 excited_states = atoms.info['excited_states']
47 excited_states = [e for e in excited_states if e['type'] == 'singlet']
48 n = len(excited_states)
50 D_n = np.zeros(n)
51 K_nn = np.zeros((n, n))
52 mu_nv = np.zeros((n, 3))
54 for i, state in enumerate(excited_states):
55 D_n[i] = state['energy'] * eV_to_au
56 mu_nv[i, :] = state['dipole_moment']
58 response = CasidaResponse(D_n, K_nn, mu_nv, broadening=broadening,
59 name=name, units='au')
60 response._atoms = atoms
61 response._pbc = atoms.pbc
62 return response
65def _read_nwchem_logfile(filename: Path_str) -> Atoms:
66 """This function parses the log file from a NWChem calculation and stores the
67 data as an :class`ASE Atoms <~ase.Atoms>` object. Each element of the list
68 corresponds to one ionic step. Additional information pertaining e.g., to
69 eigenvalues or excitations are stored in the 'info' dictionary of the ASE
70 atoms class.
72 Parameters
73 ----------
74 filename
75 Input file name.
77 Returns
78 -------
79 Atomic configuration with associated information.
80 """
81 import gzip
82 # GzipFile doesn't appear as IO type. If it does at some point,
83 # it should suffice to do just
84 # fd: IO
85 fd: Union[IO, gzip.GzipFile]
87 filename = os.path.expanduser(filename)
88 if filename.endswith('.gz'):
89 fd = gzip.open(filename)
90 elif filename.endswith('.bz2'):
91 import bz2
92 fd = bz2.open(filename, 'rb')
93 else:
94 fd = open(filename, 'rb')
95 data = [row.decode() if isinstance(row, bytes) else row for row in fd]
97 atoms = Atoms(pbc=False)
98 energy: float = np.nan
99 forces: List[Vector] = []
100 eKS_energies: List[float] = []
101 eKS_occupations: List[float] = []
102 parameters: Dict[str, Any] = OrderedDict()
104 for n, line in enumerate(data):
106 if 'echo of input deck' in line:
107 input_deck = []
108 k = 1
109 while '=====================' not in data[n+k]:
110 line = data[n+k].strip()
111 k += 1
112 if len(line) == 0:
113 continue
114 input_deck.append(line)
115 parameters['input_deck'] = input_deck
117 # Computational parameters
118 elif 'XC Information' in line:
119 parameters['xc'] = data[n+2].strip()
121 elif 'Charge :' in line:
122 flds = line.split()
123 assert len(flds) == 3, \
124 f'NWChem parser failed to read charge from line {n+1}\n --> {line}'
125 parameters['charge'] = float(flds[2])
127 elif 'Spin multiplicity:' in line:
128 flds = line.split()
129 assert len(flds) == 3, \
130 f'NWChem parser failed to read spin multiplicity from line {n+1}\n --> {line}'
131 parameters['spin_multiplicity'] = int(flds[2])
133 # Atomic configuration
134 elif 'Output coordinates in angstroms' in line:
135 atoms = Atoms(pbc=False)
136 k = 4
137 while len(data[n+k].strip()) > 0:
138 flds = data[n+k].split()
139 assert len(flds) == 6, \
140 'NWChem parser failed to read atomic number and' \
141 f' coordinates from line {n+k+1}\n --> {data[n+k]}'
142 atom_type = flds[1]
143 position = np.array([float(flds[3]), float(flds[4]), float(flds[5])])
144 atoms.append(Atom(atom_type, position))
145 k += 1
147 # Total energy and initialization of Atoms object
148 elif 'Total DFT energy' in line:
149 flds = line.split()
150 assert len(flds) == 5, \
151 f'NWChem parser failed to read energy from line {n+1}\n --> {line}'
152 energy = float(flds[4]) * au_to_eV
154 # Forces
155 elif 'DFT ENERGY GRADIENTS' in line:
156 k = 4
157 forces = []
158 while len(data[n+k].strip()) > 0:
159 flds = data[n+k].split()
160 assert len(flds) == 8, \
161 f'NWChem parser failed to read force vector from line {n+k}\n --> {data[n+k]}'
162 force = np.array([float(flds[5]), float(flds[6]), float(flds[7])])
163 force *= au_to_eVA
164 forces.append(force)
165 k += 1
167 # Kohn-Sham eigen energies
168 elif 'DFT Final Molecular Orbital Analysis' in line:
169 eKS_energies = []
170 eKS_occupations = []
172 elif 'Vector' in line and 'Occ' in line:
173 s = line
174 s = s.replace('Vector', ' ')
175 s = s.replace('Occ=', ' ')
176 s = s.replace('E=', ' ')
177 s = s.replace('D+', 'E+').replace('D-', 'E-')
178 flds = s.split()
179 assert len(flds) >= 3, \
180 f'NWChem parser failed to read KS eigenstate from line {n+1}\n --> {line}'
181 eKS_occupations.append(float(flds[1]))
182 eKS_energies.append(float(flds[2]) * au_to_eV)
184 # Excitations
185 elif 'NWChem TDDFT Module' in line:
186 maintag = 'excited_states'
187 atoms.info[maintag] = []
189 elif 'Root' in line:
190 flds = line.replace('=', ' ').split()
191 assert len(flds) > 5, \
192 f'NWChem parser failed to read Casida root from line {n+1}\n --> {line}'
194 maintag = 'excited_states'
195 if maintag not in atoms.info:
196 continue
198 transition: Dict[str, Any] = {}
199 atoms.info[maintag].append(transition)
201 # basic information (energy, type, symmetry, ...)
202 transition['type'] = flds[2]
203 transition['symmetry'] = flds[3]
204 transition['energy'] = float(flds[4]) * au_to_eV
206 transition['dipole_moment'] = np.zeros(3)
207 transition['quadrupole_moment'] = np.zeros((3, 3))
209 if 'Spin forbidden' in data[n+2]:
210 continue
212 # transition moment, dipole
213 lineno = n + 2
214 flds = data[lineno].split()
215 assert len(flds) == 8, \
216 'NWChem parser failed to read transition dipole moment' \
217 f' from line {lineno+1}\n --> {data[lineno]}'
218 transition['dipole_moment'][0] = float(flds[3])
219 transition['dipole_moment'][1] = float(flds[5])
220 transition['dipole_moment'][2] = float(flds[7])
222 # transition moment, quadrupole
223 lineno = n + 3
224 flds = data[lineno].split()
225 assert len(flds) == 8, \
226 'NWChem parser failed to read transition quadrupole moment' \
227 f' from line {lineno+1}\n --> {data[lineno]} (1)'
228 transition['quadrupole_moment'][0, 0] = float(flds[3])
229 transition['quadrupole_moment'][0, 1] = float(flds[5])
230 transition['quadrupole_moment'][0, 2] = float(flds[7])
231 lineno = n + 4
232 flds = data[lineno].split()
233 assert len(flds) == 8, \
234 'NWChem parser failed to read transition quadrupole moment' \
235 f' from line {lineno+1}\n --> {data[lineno]} (2)'
236 transition['quadrupole_moment'][1, 1] = float(flds[3])
237 transition['quadrupole_moment'][1, 2] = float(flds[5])
238 transition['quadrupole_moment'][2, 2] = float(flds[7])
239 transition['quadrupole_moment'][1, 0] = transition['quadrupole_moment'][0, 1]
240 transition['quadrupole_moment'][2, 0] = transition['quadrupole_moment'][0, 2]
241 transition['quadrupole_moment'][2, 1] = transition['quadrupole_moment'][1, 2]
243 # Parsing Hessian from 'frequencies' run
244 elif 'P.Frequency' in line:
245 maintag = 'normal_modes'
246 atoms.info[maintag] = {}
247 offset = 0
248 if len(atoms.info[maintag].keys()) > 0: 248 ↛ 249line 248 didn't jump to line 249 because the condition on line 248 was never true
249 offset = sorted(atoms.info[maintag].keys())[-1]
251 flds = line.split()
252 for m, fld in enumerate(flds[1:]):
253 mode = offset + m + 1
254 atoms.info[maintag][mode] = {}
255 atoms.info[maintag][mode]['frequency'] = float(fld)
256 atoms.info[maintag][mode]['normal_mode'] = np.zeros((len(atoms), 3))
258 for k in range(3 * len(atoms)):
259 lineno = n + k + 2
260 line = data[lineno].strip()
261 flds = line.split()
262 assert len(flds) > 2, \
263 f'NWChem parser failed to read normal mode from line {lineno+1}\n' \
264 f' --> {line} (1)'
265 iat = int(k / 3)
266 assert iat <= len(atoms), \
267 f'NWChem parser failed to read normal mode from line {lineno+1}\n' \
268 f' --> {line} (2)'
269 alph = k % 3
270 for m, fld in enumerate(flds[1:]):
271 mode = offset + m + 1
272 assert mode <= 3 * len(atoms), \
273 f'NWChem parser failed to read normal mode from line {lineno+1}\n' \
274 f' --> {line} (3)'
275 atoms.info[maintag][mode]['normal_mode'][iat, alph] = float(fld)
277 # Compile Atoms object
278 if energy is np.nan or len(atoms) == 0:
279 raise ValueError(f'Failed to extract energy and configuration from input file {filename}')
280 atoms.calc = SinglePointDFTCalculator(atoms, energy=energy, forces=forces)
281 atoms.calc.name = 'nwchem'
282 atoms.calc.parameters = parameters
283 if len(eKS_energies) > 0 and len(eKS_occupations) > 0:
284 atoms.calc.kpts = [SinglePointKPoint(1, 0, 0, eKS_energies, eKS_occupations)]
286 return atoms