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