Coverage for /builds/hweiske/ase/ase/calculators/nwchem.py: 72.73%
44 statements
« prev ^ index » next coverage.py v7.2.7, created at 2024-04-22 11:22 +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 self.calc = None
175 def input_filename(self):
176 return f'{self.prefix}.nwi'
178 def output_filename(self):
179 return f'{self.prefix}.nwo'
181 def write_input(self, atoms, properties=None, system_changes=None):
182 FileIOCalculator.write_input(self, atoms, properties, system_changes)
184 # Prepare perm and scratch directories
185 perm = os.path.abspath(self.parameters.get('perm', self.label))
186 scratch = os.path.abspath(self.parameters.get('scratch', self.label))
187 os.makedirs(perm, exist_ok=True)
188 os.makedirs(scratch, exist_ok=True)
190 io.write(Path(self.directory) / self.input_filename(),
191 atoms, properties=properties,
192 label=self.label, **self.parameters)
194 def read_results(self):
195 output = io.read(self.output_filename())
196 self.calc = output.calc
197 self.results = output.calc.results
199 def band_structure(self):
200 self.calculate()
201 perm = self.parameters.get('perm', self.label)
202 if self.calc.get_spin_polarized():
203 alpha = np.loadtxt(os.path.join(perm, self.label + '.alpha_band'))
204 beta = np.loadtxt(os.path.join(perm, self.label + '.beta_band'))
205 energies = np.array([alpha[:, 1:], beta[:, 1:]]) * Hartree
206 else:
207 data = np.loadtxt(os.path.join(perm,
208 self.label + '.restricted_band'))
209 energies = data[np.newaxis, :, 1:] * Hartree
210 eref = self.calc.get_fermi_level()
211 if eref is None:
212 eref = 0.
213 return BandStructure(self.parameters.bandpath, energies, eref)