Coverage for ase / calculators / demon / demon.py: 40.85%
355 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3"""This module defines an ASE interface to deMon.
5http://www.demon-software.com
7"""
8import os
9import os.path as op
10import shutil
11import subprocess
13import numpy as np
15import ase.data
16import ase.io
17from ase.calculators.calculator import (
18 CalculatorSetupError,
19 FileIOCalculator,
20 Parameters,
21 ReadError,
22 all_changes,
23 equal,
24)
25from ase.units import Bohr, Hartree
27from .demon_io import parse_xray
29m_e_to_amu = 1822.88839
32class Parameters_deMon(Parameters):
33 """Parameters class for the calculator.
34 Documented in Base_deMon.__init__
36 The options here are the most important ones that the user needs to be
37 aware of. Further options accepted by deMon can be set in the dictionary
38 input_arguments.
40 """
42 def __init__(
43 self,
44 label='rundir',
45 atoms=None,
46 restart=None,
47 basis_path=None,
48 ignore_bad_restart_file=FileIOCalculator._deprecated,
49 deMon_restart_path='.',
50 title='deMon input file',
51 scftype='RKS',
52 forces=False,
53 dipole=False,
54 xc='VWN',
55 guess='TB',
56 print_out='MOE',
57 basis={},
58 ecps={},
59 mcps={},
60 auxis={},
61 augment={},
62 input_arguments=None):
63 kwargs = locals()
64 kwargs.pop('self')
65 Parameters.__init__(self, **kwargs)
68class Demon(FileIOCalculator):
69 """Calculator interface to the deMon code. """
71 implemented_properties = [
72 'energy',
73 'forces',
74 'dipole',
75 'eigenvalues']
77 def __init__(self, *, command=None, **kwargs):
78 """ASE interface to the deMon code.
80 The deMon2k code can be obtained from http://www.demon-software.com
82 The DEMON_COMMAND environment variable must be set to run the
83 executable, in bash it would be set along the lines of
84 export DEMON_COMMAND="deMon.4.3.6.std > deMon_ase.out 2>&1"
86 Parameters
87 ----------
89 label : str
90 relative path to the run directory
91 atoms : Atoms object
92 the atoms object
93 command : str
94 Command to run deMon. If not present the environment
95 variable DEMON_COMMAND will be used
96 restart : str
97 Relative path to ASE restart directory for parameters and
98 atoms object and results
99 basis_path : str
100 Relative path to the directory containing
101 BASIS, AUXIS, ECPS, MCPS and AUGMENT
102 ignore_bad_restart_file : bool
103 Ignore broken or missing ASE restart files
104 By default, it is an error if the restart
105 file is missing or broken.
106 deMon_restart_path : str
107 Relative path to the deMon restart dir
108 title : str
109 Title in the deMon input file.
110 scftype : str
111 Type of scf
112 forces : bool
113 If True a force calculation will be enforced.
114 dipole : bool
115 If True a dipole calculation will be enforced
116 xc : str
117 xc-functional
118 guess : str
119 guess for initial density and wave functions
120 print_out : str | list
121 Options for the printing in deMon
122 basis : dict
123 Definition of basis sets.
124 ecps : dict
125 Definition of ECPs
126 mcps : dict
127 Definition of MCPs
128 auxis : dict
129 Definition of AUXIS
130 augment : dict
131 Definition of AUGMENT
132 input_arguments : dict
133 Explicitly given input arguments. The key is the input keyword
134 and the value is either a str, a list of str (will be written
135 on the same line as the keyword),
136 or a list of lists of str (first list is written on the first
137 line, the others on following lines.)
139 For example usage, see the tests h2o.py and h2o_xas_xes.py in
140 the directory ase/test/demon
142 """
144 parameters = Parameters_deMon(**kwargs)
146 # Setup the run command
147 if command is None:
148 command = self.cfg.get('DEMON_COMMAND')
150 FileIOCalculator.__init__(
151 self,
152 command=command,
153 **parameters)
155 def __getitem__(self, key):
156 """Convenience method to retrieve a parameter as
157 calculator[key] rather than calculator.parameters[key]
159 Parameters
160 ----------
161 key : str, the name of the parameters to get.
162 """
163 return self.parameters[key]
165 def set(self, **kwargs):
166 """Set all parameters.
168 Parameters
169 ----------
170 kwargs : Dictionary containing the keywords for deMon
171 """
172 # Put in the default arguments.
173 kwargs = self.default_parameters.__class__(**kwargs)
175 if 'parameters' in kwargs:
176 filename = kwargs.pop('parameters')
177 parameters = Parameters.read(filename)
178 parameters.update(kwargs)
179 kwargs = parameters
181 changed_parameters = {}
183 for key, value in kwargs.items():
184 oldvalue = self.parameters.get(key)
185 if key not in self.parameters or not equal(value, oldvalue):
186 changed_parameters[key] = value
187 self.parameters[key] = value
189 return changed_parameters
191 def link_file(self, fromdir, todir, filename):
192 if op.exists(todir + '/' + filename):
193 os.remove(todir + '/' + filename)
195 if op.exists(fromdir + '/' + filename):
196 os.symlink(fromdir + '/' + filename,
197 todir + '/' + filename)
198 else:
199 raise RuntimeError(
200 "{} doesn't exist".format(fromdir + '/' + filename))
202 def calculate(self,
203 atoms=None,
204 properties=['energy'],
205 system_changes=all_changes):
206 """Capture the RuntimeError from FileIOCalculator.calculate
207 and add a little debug information from the deMon output.
209 See base FileIocalculator for documentation.
210 """
212 if atoms is not None:
213 self.atoms = atoms.copy()
215 self.write_input(self.atoms, properties, system_changes)
216 command = self.command
218 # basis path
219 basis_path = self.parameters['basis_path']
220 if basis_path is None:
221 basis_path = self.cfg.get('DEMON_BASIS_PATH')
223 if basis_path is None:
224 raise RuntimeError('Please set basis_path keyword,' +
225 ' or the DEMON_BASIS_PATH' +
226 ' environment variable')
228 # link restart file
229 value = self.parameters['guess']
230 if value.upper() == 'RESTART':
231 value2 = self.parameters['deMon_restart_path']
233 if op.exists(self.directory + '/deMon.rst')\
234 or op.islink(self.directory + '/deMon.rst'):
235 os.remove(self.directory + '/deMon.rst')
236 abspath = op.abspath(value2)
238 if op.exists(abspath + '/deMon.mem') \
239 or op.islink(abspath + '/deMon.mem'):
241 shutil.copy(abspath + '/deMon.mem',
242 self.directory + '/deMon.rst')
243 else:
244 raise RuntimeError(
245 "{} doesn't exist".format(abspath + '/deMon.rst'))
247 abspath = op.abspath(basis_path)
249 for name in ['BASIS', 'AUXIS', 'ECPS', 'MCPS', 'FFDS']:
250 self.link_file(abspath, self.directory, name)
252 if command is None:
253 raise CalculatorSetupError
254 subprocess.check_call(command, shell=True, cwd=self.directory)
256 try:
257 self.read_results()
258 except Exception: # XXX Which kind of exception?
259 with open(self.directory + '/deMon.out') as fd:
260 lines = fd.readlines()
261 debug_lines = 10
262 print('##### %d last lines of the deMon.out' % debug_lines)
263 for line in lines[-20:]:
264 print(line.strip())
265 print('##### end of deMon.out')
266 raise RuntimeError
268 def set_label(self, label):
269 """Set label directory """
271 self.label = label
273 # in our case self.directory = self.label
274 self.directory = self.label
275 if self.directory == '':
276 self.directory = os.curdir
278 def write_input(self, atoms, properties=None, system_changes=None):
279 """Write input (in)-file.
280 See calculator.py for further details.
282 Parameters
283 ----------
284 atoms : The Atoms object to write.
285 properties : The properties which should be calculated.
286 system_changes : List of properties changed since last run.
288 """
289 # Call base calculator.
290 FileIOCalculator.write_input(
291 self,
292 atoms=atoms,
293 properties=properties,
294 system_changes=system_changes)
296 if system_changes is None and properties is None:
297 return
299 filename = f'{self.directory}/deMon.inp'
301 add_print = ''
303 # Start writing the file.
304 with open(filename, 'w') as fd:
306 # write keyword argument keywords
307 value = self.parameters['title']
308 self._write_argument('TITLE', value, fd)
310 fd.write('#\n')
312 value = self.parameters['scftype']
313 self._write_argument('SCFTYPE', value, fd)
315 value = self.parameters['xc']
316 self._write_argument('VXCTYPE', value, fd)
318 value = self.parameters['guess']
319 self._write_argument('GUESS', value, fd)
321 # obtain forces through a single BOMD step
322 # only if forces is in properties, or if keyword forces is True
323 value = self.parameters['forces']
324 if 'forces' in properties or value:
326 self._write_argument('DYNAMICS',
327 ['INT=1', 'MAX=0', 'STEP=0'], fd)
328 self._write_argument('TRAJECTORY', 'FORCES', fd)
329 self._write_argument('VELOCITIES', 'ZERO', fd)
330 add_print = add_print + ' ' + 'MD OPT'
332 # if dipole is True, enforce dipole calculation.
333 # Otherwise only if asked for
334 value = self.parameters['dipole']
335 if 'dipole' in properties or value:
336 self._write_argument('DIPOLE', '', fd)
338 # print argument, here other options could change this
339 value = self.parameters['print_out']
340 assert isinstance(value, str)
341 value = value + add_print
343 if len(value) != 0:
344 self._write_argument('PRINT', value, fd)
345 fd.write('#\n')
347 # write general input arguments
348 self._write_input_arguments(fd)
350 fd.write('#\n')
352 # write basis set, ecps, mcps, auxis, augment
353 basis = self.parameters['basis']
354 if 'all' not in basis:
355 basis['all'] = 'DZVP'
356 self._write_basis(fd, atoms, basis, string='BASIS')
358 ecps = self.parameters['ecps']
359 if len(ecps) != 0:
360 self._write_basis(fd, atoms, ecps, string='ECPS')
362 mcps = self.parameters['mcps']
363 if len(mcps) != 0:
364 self._write_basis(fd, atoms, mcps, string='MCPS')
366 auxis = self.parameters['auxis']
367 if len(auxis) != 0:
368 self._write_basis(fd, atoms, auxis, string='AUXIS')
370 augment = self.parameters['augment']
371 if len(augment) != 0:
372 self._write_basis(fd, atoms, augment, string='AUGMENT')
374 # write geometry
375 self._write_atomic_coordinates(fd, atoms)
377 # write xyz file for good measure.
378 ase.io.write(f'{self.directory}/deMon_atoms.xyz', self.atoms)
380 def read(self, restart_path):
381 """Read parameters from directory restart_path."""
383 self.set_label(restart_path)
385 if not op.exists(restart_path + '/deMon.inp'):
386 raise ReadError('The restart_path file {} does not exist'
387 .format(restart_path))
389 self.atoms = self.deMon_inp_to_atoms(restart_path + '/deMon.inp')
391 self.read_results()
393 def _write_input_arguments(self, fd):
394 """Write directly given input-arguments."""
395 input_arguments = self.parameters['input_arguments']
397 # Early return
398 if input_arguments is None:
399 return
401 for key, value in input_arguments.items():
402 self._write_argument(key, value, fd)
404 def _write_argument(self, key, value, fd):
405 """Write an argument to file.
406 key : a string coresponding to the input keyword
407 value : the arguments, can be a string, a number or a list
408 f : and open file
409 """
411 # for only one argument, write on same line
412 if not isinstance(value, (tuple, list)):
413 line = key.upper()
414 line += ' ' + str(value).upper()
415 fd.write(line)
416 fd.write('\n')
418 # for a list, write first argument on the first line,
419 # then the rest on new lines
420 else:
421 line = key
422 if not isinstance(value[0], (tuple, list)):
423 for i in range(len(value)):
424 line += ' ' + str(value[i].upper())
425 fd.write(line)
426 fd.write('\n')
427 else:
428 for i in range(len(value)):
429 for j in range(len(value[i])):
430 line += ' ' + str(value[i][j]).upper()
431 fd.write(line)
432 fd.write('\n')
433 line = ''
435 def _write_atomic_coordinates(self, fd, atoms):
436 """Write atomic coordinates.
438 Parameters
439 ----------
440 - f: An open file object.
441 - atoms: An atoms object.
442 """
444 fd.write('#\n')
445 fd.write('# Atomic coordinates\n')
446 fd.write('#\n')
447 fd.write('GEOMETRY CARTESIAN ANGSTROM\n')
449 for i in range(len(atoms)):
450 xyz = atoms.get_positions()[i]
451 chem_symbol = atoms.get_chemical_symbols()[i]
452 chem_symbol += str(i + 1)
454 # if tag is set to 1 then we have a ghost atom,
455 # set nuclear charge to 0
456 if atoms.get_tags()[i] == 1:
457 nuc_charge = str(0)
458 else:
459 nuc_charge = str(atoms.get_atomic_numbers()[i])
461 mass = atoms.get_masses()[i]
463 line = f'{chem_symbol:6s}'.rjust(10) + ' '
464 line += f'{xyz[0]:.5f}'.rjust(10) + ' '
465 line += f'{xyz[1]:.5f}'.rjust(10) + ' '
466 line += f'{xyz[2]:.5f}'.rjust(10) + ' '
467 line += f'{nuc_charge:5s}'.rjust(10) + ' '
468 line += f'{mass:.5f}'.rjust(10) + ' '
470 fd.write(line)
471 fd.write('\n')
473 # routine to write basis set inormation, including ecps and auxis
474 def _write_basis(self, fd, atoms, basis={}, string='BASIS'):
475 """Write basis set, ECPs, AUXIS, or AUGMENT basis
477 Parameters
478 ----------
479 - f: An open file object.
480 - atoms: An atoms object.
481 - basis: A dictionary specifying the basis set
482 - string: 'BASIS', 'ECP','AUXIS' or 'AUGMENT'
483 """
485 # basis for all atoms
486 line = f'{string}'.ljust(10)
488 if 'all' in basis:
489 default_basis = basis['all']
490 line += f'({default_basis})'.rjust(16)
492 fd.write(line)
493 fd.write('\n')
495 # basis for all atomic species
496 chemical_symbols = atoms.get_chemical_symbols()
497 chemical_symbols_set = set(chemical_symbols)
499 for _ in range(chemical_symbols_set.__len__()):
500 symbol = chemical_symbols_set.pop()
502 if symbol in basis:
503 line = f'{symbol}'.ljust(10)
504 line += f'({basis[symbol]})'.rjust(16)
505 fd.write(line)
506 fd.write('\n')
508 # basis for individual atoms
509 for i in range(len(atoms)):
511 if i in basis:
512 symbol = str(chemical_symbols[i])
513 symbol += str(i + 1)
515 line = f'{symbol}'.ljust(10)
516 line += f'({basis[i]})'.rjust(16)
517 fd.write(line)
518 fd.write('\n')
520 # Analysis routines
521 def read_results(self):
522 """Read the results from output files."""
523 self.read_energy()
524 self.read_forces(self.atoms)
525 self.read_eigenvalues()
526 self.read_dipole()
527 self.read_xray()
529 def read_energy(self):
530 """Read energy from deMon's text-output file."""
531 with open(self.label + '/deMon.out') as fd:
532 text = fd.read().upper()
534 lines = iter(text.split('\n'))
536 for line in lines:
537 if line.startswith(' TOTAL ENERGY ='):
538 self.results['energy'] = float(line.split()[-1]) * Hartree
539 break
540 else:
541 raise RuntimeError
543 def read_forces(self, atoms):
544 """Read the forces from the deMon.out file."""
546 natoms = len(atoms)
547 filename = self.label + '/deMon.out'
549 if op.isfile(filename):
550 with open(filename) as fd:
551 lines = fd.readlines()
553 # find line where the orbitals start
554 flag_found = False
555 for i in range(len(lines)):
556 if lines[i].rfind('GRADIENTS OF TIME STEP 0 IN A.U.') > -1:
557 start = i + 4
558 flag_found = True
559 break
561 if flag_found:
562 self.results['forces'] = np.zeros((natoms, 3), float)
563 for i in range(natoms):
564 line = [s for s in lines[i + start].strip().split(' ')
565 if len(s) > 0]
566 f = -np.array([float(x) for x in line[2:5]])
567 self.results['forces'][i, :] = f * (Hartree / Bohr)
569 def read_eigenvalues(self):
570 """Read eigenvalues from the 'deMon.out' file."""
571 assert os.access(self.label + '/deMon.out', os.F_OK)
573 # Read eigenvalues
574 with open(self.label + '/deMon.out') as fd:
575 lines = fd.readlines()
577 # try PRINT MOE
578 eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(
579 lines, 'ALPHA MO ENERGIES', 6)
580 eig_beta, occ_beta = self.read_eigenvalues_one_spin(
581 lines, 'BETA MO ENERGIES', 6)
583 # otherwise try PRINT MOS
584 if len(eig_alpha) == 0 and len(eig_beta) == 0:
585 eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(
586 lines, 'ALPHA MO COEFFICIENTS', 5)
587 eig_beta, occ_beta = self.read_eigenvalues_one_spin(
588 lines, 'BETA MO COEFFICIENTS', 5)
590 self.results['eigenvalues'] = np.array([eig_alpha, eig_beta]) * Hartree
591 self.results['occupations'] = np.array([occ_alpha, occ_beta])
593 def read_eigenvalues_one_spin(self, lines, string, neigs_per_line):
594 """Utility method for retreiving eigenvalues after the string "string"
595 with neigs_per_line eigenvlaues written per line
596 """
597 eig = []
598 occ = []
600 skip_line = False
601 more_eigs = False
603 # find line where the orbitals start
604 for i in range(len(lines)):
605 if lines[i].rfind(string) > -1:
606 ii = i
607 more_eigs = True
608 break
610 while more_eigs:
611 # search for two empty lines in a row preceding a line with
612 # numbers
613 for i in range(ii + 1, len(lines)):
614 if len(lines[i].split()) == 0 and \
615 len(lines[i + 1].split()) == 0 and \
616 len(lines[i + 2].split()) > 0:
617 ii = i + 2
618 break
620 # read eigenvalues, occupations
621 line = lines[ii].split()
622 if len(line) < neigs_per_line:
623 # last row
624 more_eigs = False
625 if line[0] != str(len(eig) + 1):
626 more_eigs = False
627 skip_line = True
629 if not skip_line:
630 line = lines[ii + 1].split()
631 for l in line:
632 eig.append(float(l))
633 line = lines[ii + 3].split()
634 for l in line:
635 occ.append(float(l))
636 ii = ii + 3
638 return eig, occ
640 def read_dipole(self):
641 """Read dipole moment."""
642 dipole = np.zeros(3)
643 with open(self.label + '/deMon.out') as fd:
644 lines = fd.readlines()
646 for i in range(len(lines)):
647 if lines[i].rfind('DIPOLE') > - \
648 1 and lines[i].rfind('XAS') == -1:
649 dipole[0] = float(lines[i + 1].split()[3])
650 dipole[1] = float(lines[i + 2].split()[3])
651 dipole[2] = float(lines[i + 3].split()[3])
653 # debye to e*Ang
654 self.results['dipole'] = dipole * 0.2081943482534
656 break
658 def read_xray(self):
659 """Read deMon.xry if present."""
661 # try to read core IP from, .out file
662 filename = self.label + '/deMon.out'
663 core_IP = None
664 if op.isfile(filename):
665 with open(filename) as fd:
666 lines = fd.readlines()
668 for i in range(len(lines)):
669 if lines[i].rfind('IONIZATION POTENTIAL') > -1:
670 core_IP = float(lines[i].split()[3])
672 try:
673 mode, ntrans, E_trans, osc_strength, trans_dip = parse_xray(
674 self.label + '/deMon.xry')
675 except ReadError:
676 pass
677 else:
678 xray_results = {'xray_mode': mode,
679 'ntrans': ntrans,
680 'E_trans': E_trans,
681 'osc_strength': osc_strength, # units?
682 'trans_dip': trans_dip, # units?
683 'core_IP': core_IP}
685 self.results['xray'] = xray_results
687 def deMon_inp_to_atoms(self, filename):
688 """Routine to read deMon.inp and convert it to an atoms object."""
690 with open(filename) as fd:
691 lines = fd.readlines()
693 # find line where geometry starts
694 for i in range(len(lines)):
695 if lines[i].rfind('GEOMETRY') > -1:
696 if lines[i].rfind('ANGSTROM'):
697 coord_units = 'Ang'
698 elif lines.rfind('Bohr'):
699 coord_units = 'Bohr'
700 ii = i
701 break
703 chemical_symbols = []
704 xyz = []
705 atomic_numbers = []
706 masses = []
708 for i in range(ii + 1, len(lines)):
709 try:
710 line = lines[i].split()
712 if len(line) > 0:
713 for symbol in ase.data.chemical_symbols:
714 found = None
715 if line[0].upper().rfind(symbol.upper()) > -1:
716 found = symbol
717 break
719 if found is not None:
720 chemical_symbols.append(found)
721 else:
722 break
724 xyz.append(
725 [float(line[1]), float(line[2]), float(line[3])])
727 if len(line) > 4:
728 atomic_numbers.append(int(line[4]))
730 if len(line) > 5:
731 masses.append(float(line[5]))
733 except Exception: # XXX Which kind of exception?
734 raise RuntimeError
736 if coord_units == 'Bohr':
737 xyz *= Bohr
739 natoms = len(chemical_symbols)
741 # set atoms object
742 atoms = ase.Atoms(symbols=chemical_symbols, positions=xyz)
744 # if atomic numbers were read in, set them
745 if len(atomic_numbers) == natoms:
746 atoms.set_atomic_numbers(atomic_numbers)
748 # if masses were read in, set them
749 if len(masses) == natoms:
750 atoms.set_masses(masses)
752 return atoms