Coverage for /builds/debichem-team/python-ase/ase/calculators/dftd3.py: 86.91%
275 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1import os
2import subprocess
3from pathlib import Path
4from warnings import warn
6import numpy as np
8from ase.calculators.calculator import (
9 BaseCalculator,
10 Calculator,
11 FileIOCalculator,
12)
13from ase.io import write
14from ase.io.vasp import write_vasp
15from ase.parallel import world
16from ase.units import Bohr, Hartree
19def dftd3_defaults():
20 default_parameters = {'xc': None, # PBE if no custom damping parameters
21 'grad': True, # calculate forces/stress
22 'abc': False, # ATM 3-body contribution
23 'cutoff': 95 * Bohr, # Cutoff for 2-body calcs
24 'cnthr': 40 * Bohr, # Cutoff for 3-body and CN calcs
25 'old': False, # use old DFT-D2 method instead
26 'damping': 'zero', # Default to zero-damping
27 'tz': False, # 'triple zeta' alt. parameters
28 's6': None, # damping parameters start here
29 'sr6': None,
30 's8': None,
31 'sr8': None,
32 'alpha6': None,
33 'a1': None,
34 'a2': None,
35 'beta': None}
36 return default_parameters
39class DFTD3(BaseCalculator):
40 """Grimme DFT-D3 calculator"""
42 def __init__(self,
43 label='ase_dftd3', # Label for dftd3 output files
44 command=None, # Command for running dftd3
45 dft=None, # DFT calculator
46 comm=world,
47 **kwargs):
49 # Convert from 'func' keyword to 'xc'. Internally, we only store
50 # 'xc', but 'func' is also allowed since it is consistent with the
51 # CLI dftd3 interface.
52 func = kwargs.pop('func', None)
53 if func is not None:
54 if kwargs.get('xc') is not None:
55 raise RuntimeError('Both "func" and "xc" were provided! '
56 'Please provide at most one of these '
57 'two keywords. The preferred keyword '
58 'is "xc"; "func" is allowed for '
59 'consistency with the CLI dftd3 '
60 'interface.')
61 kwargs['xc'] = func
63 # If the user did not supply an XC functional, but did attach a
64 # DFT calculator that has XC set, then we will use that. Note that
65 # DFTD3's spelling convention is different from most, so in general
66 # you will have to explicitly set XC for both the DFT calculator and
67 # for DFTD3 (and DFTD3's will likely be spelled differently...)
68 if dft is not None and kwargs.get('xc') is None:
69 dft_xc = dft.parameters.get('xc')
70 if dft_xc is not None:
71 kwargs['xc'] = dft_xc
73 dftd3 = PureDFTD3(label=label, command=command, comm=comm, **kwargs)
75 # dftd3 only implements energy, forces, and stresses (for periodic
76 # systems). But, if a DFT calculator is attached, and that calculator
77 # implements more properties, we expose those properties.
78 # dftd3 contributions for those properties will be zero.
79 if dft is None:
80 self.implemented_properties = list(dftd3.dftd3_properties)
81 else:
82 self.implemented_properties = list(dft.implemented_properties)
84 # Should our arguments be "parameters" (passed to superclass)
85 # or are they not really "parameters"?
86 #
87 # That's not really well defined. Let's not do anything then.
88 super().__init__()
90 self.dftd3 = dftd3
91 self.dft = dft
93 def todict(self):
94 return {}
96 def calculate(self, atoms, properties, system_changes):
97 common_props = set(self.dftd3.dftd3_properties) & set(properties)
98 dftd3_results = self._get_properties(atoms, common_props, self.dftd3)
100 if self.dft is None:
101 results = dftd3_results
102 else:
103 dft_results = self._get_properties(atoms, properties, self.dft)
104 results = dict(dft_results)
105 for name in set(results) & set(dftd3_results):
106 assert np.shape(results[name]) == np.shape(dftd3_results[name])
107 results[name] += dftd3_results[name]
109 # Although DFTD3 may have calculated quantities not provided
110 # by the calculator (e.g. stress), it would be wrong to
111 # return those! Return only what corresponds to the DFT calc.
112 assert set(results) == set(dft_results)
113 self.results = results
115 def _get_properties(self, atoms, properties, calc):
116 # We want any and all properties that the calculator
117 # normally produces. So we intend to rob the calc.results
118 # dictionary instead of only getting the requested properties.
120 import copy
121 for name in properties:
122 calc.get_property(name, atoms)
123 assert name in calc.results
125 # XXX maybe use get_properties() when that makes sense.
126 results = copy.deepcopy(calc.results)
127 assert set(properties) <= set(results)
128 return results
131class PureDFTD3(FileIOCalculator):
132 """DFTD3 calculator without corresponding DFT contribution.
134 This class is an implementation detail."""
136 name = 'puredftd3'
138 dftd3_properties = {'energy', 'free_energy', 'forces', 'stress'}
139 implemented_properties = list(dftd3_properties)
140 default_parameters = dftd3_defaults()
141 damping_methods = {'zero', 'bj', 'zerom', 'bjm'}
142 _legacy_default_command = 'dftd3'
144 def __init__(self,
145 *,
146 label='ase_dftd3', # Label for dftd3 output files
147 command=None, # Command for running dftd3
148 comm=world,
149 **kwargs):
151 # FileIOCalculator would default to self.name to get the envvar
152 # which determines the command.
153 # We'll have to overrule that if we want to keep scripts working:
154 command = command or self.cfg.get('ASE_DFTD3_COMMAND')
156 super().__init__(label=label,
157 command=command,
158 **kwargs)
160 # TARP: This is done because the calculator does not call
161 # FileIOCalculator.calculate, but Calculator.calculate and does not
162 # use the profile defined in FileIOCalculator.__init__
163 self.comm = comm
165 def set(self, **kwargs):
166 changed_parameters = {}
168 # Check for unknown arguments. Don't raise an error, just let the
169 # user know that we don't understand what they're asking for.
170 unknown_kwargs = set(kwargs) - set(self.default_parameters)
171 if unknown_kwargs:
172 warn('WARNING: Ignoring the following unknown keywords: {}'
173 ''.format(', '.join(unknown_kwargs)))
175 changed_parameters.update(FileIOCalculator.set(self, **kwargs))
177 # Ensure damping method is valid (zero, bj, zerom, bjm).
178 damping = self.parameters['damping']
179 if damping is not None:
180 damping = damping.lower()
181 if damping not in self.damping_methods:
182 raise ValueError(f'Unknown damping method {damping}!')
184 # d2 only is valid with 'zero' damping
185 elif self.parameters['old'] and damping != 'zero':
186 raise ValueError('Only zero-damping can be used with the D2 '
187 'dispersion correction method!')
189 # If cnthr (cutoff for three-body and CN calculations) is greater
190 # than cutoff (cutoff for two-body calculations), then set the former
191 # equal to the latter, since that doesn't make any sense.
192 if self.parameters['cnthr'] > self.parameters['cutoff']:
193 warn('WARNING: CN cutoff value of {cnthr} is larger than '
194 'regular cutoff value of {cutoff}! Reducing CN cutoff '
195 'to {cutoff}.'
196 ''.format(cnthr=self.parameters['cnthr'],
197 cutoff=self.parameters['cutoff']))
198 self.parameters['cnthr'] = self.parameters['cutoff']
200 # If you only care about the energy, gradient calculations (forces,
201 # stresses) can be bypassed. This will greatly speed up calculations
202 # in dense 3D-periodic systems with three-body corrections. But, we
203 # can no longer say that we implement forces and stresses.
204 # if not self.parameters['grad']:
205 # for val in ['forces', 'stress']:
206 # if val in self.implemented_properties:
207 # self.implemented_properties.remove(val)
209 # Check to see if we're using custom damping parameters.
210 zero_damppars = {'s6', 'sr6', 's8', 'sr8', 'alpha6'}
211 bj_damppars = {'s6', 'a1', 's8', 'a2', 'alpha6'}
212 zerom_damppars = {'s6', 'sr6', 's8', 'beta', 'alpha6'}
213 all_damppars = zero_damppars | bj_damppars | zerom_damppars
215 self.custom_damp = False
217 damppars = set(kwargs) & all_damppars
218 if damppars:
219 self.custom_damp = True
220 if damping == 'zero':
221 valid_damppars = zero_damppars
222 elif damping in ['bj', 'bjm']:
223 valid_damppars = bj_damppars
224 elif damping == 'zerom':
225 valid_damppars = zerom_damppars
227 # If some but not all damping parameters are provided for the
228 # selected damping method, raise an error. We don't have "default"
229 # values for damping parameters, since those are stored in the
230 # dftd3 executable & depend on XC functional.
231 missing_damppars = valid_damppars - damppars
232 if missing_damppars and missing_damppars != valid_damppars:
233 raise ValueError('An incomplete set of custom damping '
234 'parameters for the {} damping method was '
235 'provided! Expected: {}; got: {}'
236 ''.format(damping,
237 ', '.join(valid_damppars),
238 ', '.join(damppars)))
240 # If a user provides damping parameters that are not used in the
241 # selected damping method, let them know that we're ignoring them.
242 # If the user accidentally provided the *wrong* set of parameters,
243 # (e.g., the BJ parameters when they are using zero damping), then
244 # the previous check will raise an error, so we don't need to
245 # worry about that here.
246 if damppars - valid_damppars:
247 warn('WARNING: The following damping parameters are not '
248 'valid for the {} damping method and will be ignored: {}'
249 ''.format(damping,
250 ', '.join(damppars)))
252 # The default XC functional is PBE, but this is only set if the user
253 # did not provide their own value for xc or any custom damping
254 # parameters.
255 if self.parameters['xc'] and self.custom_damp:
256 warn('WARNING: Custom damping parameters will be used '
257 'instead of those parameterized for {}!'
258 ''.format(self.parameters['xc']))
260 if changed_parameters:
261 self.results.clear()
262 return changed_parameters
264 def calculate(self, atoms, properties, system_changes):
265 # We don't call FileIOCalculator.calculate here, because that method
266 # calls subprocess.call(..., shell=True), which we don't want to do.
267 # So, we reproduce some content from that method here.
268 Calculator.calculate(self, atoms, properties, system_changes)
270 # If a parameter file exists in the working directory, delete it
271 # first. If we need that file, we'll recreate it later.
272 localparfile = os.path.join(self.directory, '.dftd3par.local')
273 if self.comm.rank == 0 and os.path.isfile(localparfile):
274 os.remove(localparfile)
276 # Write XYZ or POSCAR file and .dftd3par.local file if we are using
277 # custom damping parameters.
278 self.write_input(self.atoms, properties, system_changes)
279 # command = self._generate_command()
281 inputs = DFTD3Inputs(command=self.command, prefix=self.label,
282 atoms=self.atoms, parameters=self.parameters)
283 command = inputs.get_argv(custom_damp=self.custom_damp)
285 # Finally, call dftd3 and parse results.
286 # DFTD3 does not run in parallel
287 # so we only need it to run on 1 core
288 errorcode = 0
289 if self.comm.rank == 0:
290 with open(self.label + '.out', 'w') as fd:
291 errorcode = subprocess.call(command,
292 cwd=self.directory, stdout=fd)
294 errorcode = self.comm.sum_scalar(errorcode)
296 if errorcode:
297 raise RuntimeError('%s returned an error: %d' %
298 (self.name, errorcode))
300 self.read_results()
302 def write_input(self, atoms, properties=None, system_changes=None):
303 FileIOCalculator.write_input(self, atoms, properties=properties,
304 system_changes=system_changes)
305 # dftd3 can either do fully 3D periodic or non-periodic calculations.
306 # It cannot do calculations that are only periodic in 1 or 2
307 # dimensions. If the atoms object is periodic in only 1 or 2
308 # dimensions, then treat it as a fully 3D periodic system, but warn
309 # the user.
311 if self.custom_damp:
312 damppars = _get_damppars(self.parameters)
313 else:
314 damppars = None
316 pbc = any(atoms.pbc)
317 if pbc and not all(atoms.pbc):
318 warn('WARNING! dftd3 can only calculate the dispersion energy '
319 'of non-periodic or 3D-periodic systems. We will treat '
320 'this system as 3D-periodic!')
322 if self.comm.rank == 0:
323 self._actually_write_input(
324 directory=Path(self.directory), atoms=atoms,
325 properties=properties, prefix=self.label,
326 damppars=damppars, pbc=pbc)
328 def _actually_write_input(self, directory, prefix, atoms, properties,
329 damppars, pbc):
330 if pbc:
331 fname = directory / f'{prefix}.POSCAR'
332 # We sort the atoms so that the atomtypes list becomes as
333 # short as possible. The dftd3 program can only handle 10
334 # atomtypes
335 write_vasp(fname, atoms, sort=True)
336 else:
337 fname = directory / f'{prefix}.xyz'
338 write(fname, atoms, format='xyz', parallel=False)
340 # Generate custom damping parameters file. This is kind of ugly, but
341 # I don't know of a better way of doing this.
342 if damppars is not None:
343 damp_fname = directory / '.dftd3par.local'
344 with open(damp_fname, 'w') as fd:
345 fd.write(' '.join(damppars))
347 def _outname(self):
348 return Path(self.directory) / f'{self.label}.out'
350 def _read_and_broadcast_results(self):
351 from ase.parallel import broadcast
352 if self.comm.rank == 0:
353 output = DFTD3Output(directory=self.directory,
354 stdout_path=self._outname())
355 dct = output.read(atoms=self.atoms,
356 read_forces=bool(self.parameters['grad']))
357 else:
358 dct = None
360 dct = broadcast(dct, root=0, comm=self.comm)
361 return dct
363 def read_results(self):
364 results = self._read_and_broadcast_results()
365 self.results = results
368class DFTD3Inputs:
369 dftd3_flags = {'grad', 'pbc', 'abc', 'old', 'tz'}
371 def __init__(self, command, prefix, atoms, parameters):
372 self.command = command
373 self.prefix = prefix
374 self.atoms = atoms
375 self.parameters = parameters
377 @property
378 def pbc(self):
379 return any(self.atoms.pbc)
381 @property
382 def inputformat(self):
383 if self.pbc:
384 return 'POSCAR'
385 else:
386 return 'xyz'
388 def get_argv(self, custom_damp):
389 argv = self.command.split()
391 argv.append(f'{self.prefix}.{self.inputformat}')
393 if not custom_damp:
394 xc = self.parameters.get('xc')
395 if xc is None:
396 xc = 'pbe'
397 argv += ['-func', xc.lower()]
399 for arg in self.dftd3_flags:
400 if self.parameters.get(arg):
401 argv.append('-' + arg)
403 if self.pbc:
404 argv.append('-pbc')
406 argv += ['-cnthr', str(self.parameters['cnthr'] / Bohr)]
407 argv += ['-cutoff', str(self.parameters['cutoff'] / Bohr)]
409 if not self.parameters['old']:
410 argv.append('-' + self.parameters['damping'])
412 return argv
415class DFTD3Output:
416 def __init__(self, directory, stdout_path):
417 self.directory = Path(directory)
418 self.stdout_path = Path(stdout_path)
420 def read(self, *, atoms, read_forces):
421 results = {}
423 energy = self.read_energy()
424 results['energy'] = energy
425 results['free_energy'] = energy
427 if read_forces:
428 results['forces'] = self.read_forces(atoms)
430 if any(atoms.pbc):
431 results['stress'] = self.read_stress(atoms.cell)
433 return results
435 def read_forces(self, atoms):
436 forcename = self.directory / 'dftd3_gradient'
437 with open(forcename) as fd:
438 forces = self.parse_forces(fd)
440 assert len(forces) == len(atoms)
442 forces *= -Hartree / Bohr
443 # XXXX ordering!
444 if any(atoms.pbc):
445 # This seems to be due to vasp file sorting.
446 # If that sorting rule changes, we will get garbled
447 # forces!
448 ind = np.argsort(atoms.symbols)
449 forces[ind] = forces.copy()
450 return forces
452 def read_stress(self, cell):
453 volume = cell.volume
454 assert volume > 0
456 stress = self.read_cellgradient()
457 stress *= Hartree / Bohr / volume
458 stress = stress.T @ cell
459 return stress.flat[[0, 4, 8, 5, 2, 1]]
461 def read_cellgradient(self):
462 with (self.directory / 'dftd3_cellgradient').open() as fd:
463 return self.parse_cellgradient(fd)
465 def read_energy(self) -> float:
466 with self.stdout_path.open() as fd:
467 return self.parse_energy(fd, self.stdout_path)
469 def parse_energy(self, fd, outname):
470 for line in fd:
471 if line.startswith(' program stopped'):
472 if 'functional name unknown' in line:
473 message = ('Unknown DFTD3 functional name. '
474 'Please check the dftd3.f source file '
475 'for the list of known functionals '
476 'and their spelling.')
477 else:
478 message = ('dftd3 failed! Please check the {} '
479 'output file and report any errors '
480 'to the ASE developers.'
481 ''.format(outname))
482 raise RuntimeError(message)
484 if line.startswith(' Edisp'):
485 # line looks something like this:
486 #
487 # Edisp /kcal,au,ev: xxx xxx xxx
488 #
489 parts = line.split()
490 assert parts[1][0] == '/'
491 index = 2 + parts[1][1:-1].split(',').index('au')
492 e_dftd3 = float(parts[index]) * Hartree
493 return e_dftd3
495 raise RuntimeError('Could not parse energy from dftd3 '
496 'output, see file {}'.format(outname))
498 def parse_forces(self, fd):
499 forces = []
500 for i, line in enumerate(fd):
501 forces.append(line.split())
502 return np.array(forces, dtype=float)
504 def parse_cellgradient(self, fd):
505 stress = np.zeros((3, 3))
506 for i, line in enumerate(fd):
507 for j, x in enumerate(line.split()):
508 stress[i, j] = float(x)
509 # Check if all stress elements are present?
510 # Check if file is longer?
511 return stress
514def _get_damppars(par):
515 damping = par['damping']
517 damppars = []
519 # s6 is always first
520 damppars.append(str(float(par['s6'])))
522 # sr6 is the second value for zero{,m} damping, a1 for bj{,m}
523 if damping in ['zero', 'zerom']:
524 damppars.append(str(float(par['sr6'])))
525 elif damping in ['bj', 'bjm']:
526 damppars.append(str(float(par['a1'])))
528 # s8 is always third
529 damppars.append(str(float(par['s8'])))
531 # sr8 is fourth for zero, a2 for bj{,m}, beta for zerom
532 if damping == 'zero':
533 damppars.append(str(float(par['sr8'])))
534 elif damping in ['bj', 'bjm']:
535 damppars.append(str(float(par['a2'])))
536 elif damping == 'zerom':
537 damppars.append(str(float(par['beta'])))
538 # alpha6 is always fifth
539 damppars.append(str(int(par['alpha6'])))
541 # last is the version number
542 if par['old']:
543 damppars.append('2')
544 elif damping == 'zero':
545 damppars.append('3')
546 elif damping == 'bj':
547 damppars.append('4')
548 elif damping == 'zerom':
549 damppars.append('5')
550 elif damping == 'bjm':
551 damppars.append('6')
552 return damppars