Coverage for /builds/debichem-team/python-ase/ase/calculators/nwchem.py: 71.74%
46 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1"""This module defines an ASE interface to NWchem
3https://nwchemgit.github.io
4"""
5import os
6from pathlib import Path
8import numpy as np
10from ase import io
11from ase.calculators.calculator import FileIOCalculator
12from ase.spectrum.band_structure import BandStructure
13from ase.units import Hartree
16class NWChem(FileIOCalculator):
17 implemented_properties = ['energy', 'free_energy',
18 'forces', 'stress', 'dipole']
19 _legacy_default_command = 'nwchem PREFIX.nwi > PREFIX.nwo'
20 accepts_bandpath_keyword = True
21 discard_results_on_any_change = True
23 fileio_rules = FileIOCalculator.ruleset(
24 extend_argv=['{prefix}.nwi'],
25 stdout_name='{prefix}.nwo')
27 def __init__(self, restart=None,
28 ignore_bad_restart_file=FileIOCalculator._deprecated,
29 label='nwchem', atoms=None, command=None, **kwargs):
30 """
31 NWChem keywords are specified using (potentially nested)
32 dictionaries. Consider the following input file block::
34 dft
35 odft
36 mult 2
37 convergence energy 1e-9 density 1e-7 gradient 5e-6
38 end
40 This can be generated by the NWChem calculator by using the
41 following settings:
42 >>> from ase.calculators.nwchem import NWChem
43 >>> calc = NWChem(dft={'odft': None,
44 ... 'mult': 2,
45 ... 'convergence': {'energy': 1e-9,
46 ... 'density': 1e-7,
47 ... 'gradient': 5e-6,
48 ... },
49 ... },
50 ... )
52 In addition, the calculator supports several special keywords:
54 theory: str
55 Which NWChem module should be used to calculate the
56 energies and forces. Supported values are ``'dft'``,
57 ``'scf'``, ``'mp2'``, ``'ccsd'``, ``'tce'``, ``'tddft'``,
58 ``'pspw'``, ``'band'``, and ``'paw'``. If not provided, the
59 calculator will attempt to guess which theory to use based
60 on the keywords provided by the user.
61 xc: str
62 The exchange-correlation functional to use. Only relevant
63 for DFT calculations.
64 task: str
65 What type of calculation is to be performed, e.g.
66 ``'energy'``, ``'gradient'``, ``'optimize'``, etc. When
67 using ``'SocketIOCalculator'``, ``task`` should be set
68 to ``'optimize'``. In most other circumstances, ``task``
69 should not be set manually.
70 basis: str or dict
71 Which basis set to use for gaussian-type orbital
72 calculations. Set to a string to use the same basis for all
73 elements. To use a different basis for different elements,
74 provide a dict of the form:
76 >>> calc = NWChem(...,
77 ... basis={'O': '3-21G',
78 ... 'Si': '6-31g'})
80 basispar: str
81 Additional keywords to put in the NWChem ``basis`` block,
82 e.g. ``'rel'`` for relativistic bases.
83 symmetry: int or str
84 The point group (for gaussian-type orbital calculations) or
85 space group (for plane-wave calculations) of the system.
86 Supports both group names (e.g. ``'c2v'``, ``'Fm3m'``) and
87 numbers (e.g. ``225``).
88 autosym: bool
89 Whether NWChem should automatically determine the symmetry
90 of the structure (defaults to ``False``).
91 center: bool
92 Whether NWChem should automatically center the structure
93 (defaults to ``False``). Enable at your own risk.
94 autoz: bool
95 Whether NWChem should automatically construct a Z-matrix
96 for your molecular system (defaults to ``False``).
97 geompar: str
98 Additional keywords to put in the NWChem `geometry` block,
99 e.g. ``'nucleus finite'`` for gaussian-shaped nuclear
100 charges. Do not set ``'autosym'``, ``'center'``, or
101 ``'autoz'`` in this way; instead, use the appropriate
102 keyword described above for these settings.
103 set: dict
104 Used to manually create or modify entries in the NWChem
105 rtdb. For example, the following settings enable
106 pseudopotential filtering for plane-wave calculations::
108 set nwpw:kbpp_ray .true.
109 set nwpw:kbpp_filter .true.
111 These settings are generated by the NWChem calculator by
112 passing the arguments:
114 >>> calc = NWChem(...,
115 >>> set={'nwpw:kbpp_ray': True,
116 >>> 'nwpw:kbpp_filter': True})
118 kpts: (int, int, int), or dict
119 Indicates which k-point mesh to use. Supported syntax is
120 similar to that of GPAW. Implies ``theory='band'``.
121 bandpath: BandPath object
122 The band path to use for a band structure calculation.
123 Implies ``theory='band'``.
124 pretasks: list of dict
125 Tasks used to produce a better initial guess
126 for the wavefunction.
127 These task typically use a cheaper level of theory
128 or smaller basis set (but not both).
129 The output energy and forces should remain unchanged
130 regardless of the number of tasks or their parameters,
131 but the runtime may be significantly improved.
133 For example, a MP2 calculation preceded by guesses at the
134 DFT and HF levels would be
136 >>> calc = NWChem(theory='mp2', basis='aug-cc-pvdz',
137 >>> pretasks=[
138 >>> {'dft': {'xc': 'hfexch'},
139 >>> 'set': {'lindep:n_dep': 0}},
140 >>> {'theory': 'scf', 'set': {'lindep:n_dep': 0}},
141 >>> ])
143 Each dictionary could contain any of the other parameters,
144 except those which pertain to global configurations
145 (e.g., geometry details, scratch dir).
147 The default basis set is that of the final step in the calculation,
148 or that of the previous step that which defines a basis set.
149 For example, all steps in the example will use aug-cc-pvdz
150 because the last step is the only one which defines a basis.
152 Steps which change basis set must use the same theory.
153 The following specification would perform SCF using the 3-21G
154 basis set first, then B3LYP//3-21g, and then B3LYP//6-31G(2df,p)
156 >>> calc = NWChem(theory='dft', xc='b3lyp', basis='6-31g(2df,p)',
157 >>> pretasks=[
158 >>> {'theory': 'scf', 'basis': '3-21g',
159 >>> 'set': {'lindep:n_dep': 0}},
160 >>> {'dft': {'xc': 'b3lyp'}},
161 >>> ])
163 The :code:`'set': {'lindep:n_dep': 0}` option is highly suggested
164 as a way to avoid errors relating to symmetry changes between tasks.
166 The calculator will configure appropriate options for saving
167 and loading intermediate wavefunctions, and
168 place an "ignore" task directive between each step so that
169 convergence errors in intermediate steps do not halt execution.
170 """
171 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
172 label, atoms, command, **kwargs)
173 if self.prefix is None:
174 self.prefix = 'nwchem'
175 self.calc = None
177 def input_filename(self):
178 return f'{self.prefix}.nwi'
180 def output_filename(self):
181 return f'{self.prefix}.nwo'
183 def write_input(self, atoms, properties=None, system_changes=None):
184 FileIOCalculator.write_input(self, atoms, properties, system_changes)
186 # Prepare perm and scratch directories
187 perm = os.path.abspath(self.parameters.get('perm', self.label))
188 scratch = os.path.abspath(self.parameters.get('scratch', self.label))
189 os.makedirs(perm, exist_ok=True)
190 os.makedirs(scratch, exist_ok=True)
192 io.write(Path(self.directory) / self.input_filename(),
193 atoms, properties=properties,
194 label=self.label, **self.parameters)
196 def read_results(self):
197 output = io.read(Path(self.directory) / self.output_filename())
198 self.calc = output.calc
199 self.results = output.calc.results
201 def band_structure(self):
202 self.calculate()
203 perm = self.parameters.get('perm', self.label)
204 if self.calc.get_spin_polarized():
205 alpha = np.loadtxt(os.path.join(perm, self.label + '.alpha_band'))
206 beta = np.loadtxt(os.path.join(perm, self.label + '.beta_band'))
207 energies = np.array([alpha[:, 1:], beta[:, 1:]]) * Hartree
208 else:
209 data = np.loadtxt(os.path.join(perm,
210 self.label + '.restricted_band'))
211 energies = data[np.newaxis, :, 1:] * Hartree
212 eref = self.calc.get_fermi_level()
213 if eref is None:
214 eref = 0.
215 return BandStructure(self.parameters.bandpath, energies, eref)