Coverage for /builds/ase/ase/ase/calculators/demon/demon.py : 10.22%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""This module defines an ASE interface to deMon.
3http://www.demon-software.com
5"""
6import os
7import os.path as op
8import subprocess
9import shutil
11import numpy as np
13from ase.units import Bohr, Hartree
14import ase.data
15from ase.calculators.calculator import FileIOCalculator, ReadError
16from ase.calculators.calculator import Parameters, all_changes
17from ase.calculators.calculator import equal
18import ase.io
19from .demon_io import parse_xray
21m_e_to_amu = 1822.88839
24class Parameters_deMon(Parameters):
25 """Parameters class for the calculator.
26 Documented in Base_deMon.__init__
28 The options here are the most important ones that the user needs to be
29 aware of. Further options accepted by deMon can be set in the dictionary
30 input_arguments.
32 """
34 def __init__(
35 self,
36 label='rundir',
37 atoms=None,
38 command=None,
39 restart=None,
40 basis_path=None,
41 ignore_bad_restart_file=FileIOCalculator._deprecated,
42 deMon_restart_path='.',
43 title='deMon input file',
44 scftype='RKS',
45 forces=False,
46 dipole=False,
47 xc='VWN',
48 guess='TB',
49 print_out='MOE',
50 basis={},
51 ecps={},
52 mcps={},
53 auxis={},
54 augment={},
55 input_arguments=None):
56 kwargs = locals()
57 kwargs.pop('self')
58 Parameters.__init__(self, **kwargs)
61class Demon(FileIOCalculator):
62 """Calculator interface to the deMon code. """
64 implemented_properties = [
65 'energy',
66 'forces',
67 'dipole',
68 'eigenvalues']
70 def __init__(self, **kwargs):
71 """ASE interface to the deMon code.
73 The deMon2k code can be obtained from http://www.demon-software.com
75 The DEMON_COMMAND environment variable must be set to run the
76 executable, in bash it would be set along the lines of
77 export DEMON_COMMAND="deMon.4.3.6.std > deMon_ase.out 2>&1"
79 Parameters:
81 label : str
82 relative path to the run directory
83 atoms : Atoms object
84 the atoms object
85 command : str
86 Command to run deMon. If not present the environment
87 variable DEMON_COMMAND will be used
88 restart : str
89 Relative path to ASE restart directory for parameters and
90 atoms object and results
91 basis_path : str
92 Relative path to the directory containing
93 BASIS, AUXIS, ECPS, MCPS and AUGMENT
94 ignore_bad_restart_file : bool
95 Ignore broken or missing ASE restart files
96 By default, it is an error if the restart
97 file is missing or broken.
98 deMon_restart_path : str
99 Relative path to the deMon restart dir
100 title : str
101 Title in the deMon input file.
102 scftype : str
103 Type of scf
104 forces : bool
105 If True a force calculation will be enforced.
106 dipole : bool
107 If True a dipole calculation will be enforced
108 xc : str
109 xc-functional
110 guess : str
111 guess for initial density and wave functions
112 print_out : str | list
113 Options for the printing in deMon
114 basis : dict
115 Definition of basis sets.
116 ecps : dict
117 Definition of ECPs
118 mcps : dict
119 Definition of MCPs
120 auxis : dict
121 Definition of AUXIS
122 augment : dict
123 Definition of AUGMENT
124 input_arguments : dict
125 Explicitly given input arguments. The key is the input keyword
126 and the value is either a str, a list of str (will be written
127 on the same line as the keyword),
128 or a list of lists of str (first list is written on the first
129 line, the others on following lines.)
131 For example usage, see the tests h2o.py and h2o_xas_xes.py in
132 the directory ase/test/demon
134 """
136 parameters = Parameters_deMon(**kwargs)
138 # Setup the run command
139 command = parameters['command']
140 if command is None:
141 command = os.environ.get('DEMON_COMMAND')
143 if command is None:
144 mess = 'The "DEMON_COMMAND" environment is not defined.'
145 raise ValueError(mess)
146 else:
147 parameters['command'] = command
149 # Call the base class.
150 FileIOCalculator.__init__(
151 self,
152 **parameters)
154 def __getitem__(self, key):
155 """Convenience method to retrieve a parameter as
156 calculator[key] rather than calculator.parameters[key]
158 Parameters:
159 key : str, the name of the parameters to get.
160 """
161 return self.parameters[key]
163 def set(self, **kwargs):
164 """Set all parameters.
166 Parameters:
167 kwargs : Dictionary containing the keywords for deMon
168 """
169 # Put in the default arguments.
170 kwargs = self.default_parameters.__class__(**kwargs)
172 if 'parameters' in kwargs:
173 filename = kwargs.pop('parameters')
174 parameters = Parameters.read(filename)
175 parameters.update(kwargs)
176 kwargs = parameters
178 changed_parameters = {}
180 for key, value in kwargs.items():
181 oldvalue = self.parameters.get(key)
182 if key not in self.parameters or not equal(value, oldvalue):
183 changed_parameters[key] = value
184 self.parameters[key] = value
186 return changed_parameters
188 def link_file(self, fromdir, todir, filename):
190 if op.exists(todir + '/' + filename):
191 os.remove(todir + '/' + filename)
193 if op.exists(fromdir + '/' + filename):
194 os.symlink(fromdir + '/' + filename,
195 todir + '/' + filename)
196 else:
197 raise RuntimeError(
198 "{0} doesn't exist".format(fromdir + '/' + filename))
200 def calculate(self,
201 atoms=None,
202 properties=['energy'],
203 system_changes=all_changes):
204 """Capture the RuntimeError from FileIOCalculator.calculate
205 and add a little debug information from the deMon output.
207 See base FileIocalculator for documentation.
208 """
210 if atoms is not None:
211 self.atoms = atoms.copy()
213 self.write_input(self.atoms, properties, system_changes)
214 if self.command is None:
215 raise RuntimeError('Please set $%s environment variable ' %
216 ('DEMON_COMMAND') +
217 'or supply the command keyword')
218 command = self.command # .replace('PREFIX', self.prefix)
220 # basis path
221 basis_path = self.parameters['basis_path']
222 if basis_path is None:
223 basis_path = os.environ.get('DEMON_BASIS_PATH')
225 if basis_path is None:
226 raise RuntimeError('Please set basis_path keyword,' +
227 ' or the DEMON_BASIS_PATH' +
228 ' environment variable')
230 # link restart file
231 value = self.parameters['guess']
232 if value.upper() == 'RESTART':
233 value2 = self.parameters['deMon_restart_path']
235 if op.exists(self.directory + '/deMon.rst')\
236 or op.islink(self.directory + '/deMon.rst'):
237 os.remove(self.directory + '/deMon.rst')
238 abspath = op.abspath(value2)
240 if op.exists(abspath + '/deMon.mem') \
241 or op.islink(abspath + '/deMon.mem'):
243 shutil.copy(abspath + '/deMon.mem',
244 self.directory + '/deMon.rst')
245 else:
246 raise RuntimeError(
247 "{0} doesn't exist".format(abspath + '/deMon.rst'))
249 abspath = op.abspath(basis_path)
251 for name in ['BASIS', 'AUXIS', 'ECPS', 'MCPS', 'FFDS']:
252 self.link_file(abspath, self.directory, name)
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', 'r') 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 atoms : The Atoms object to write.
284 properties : The properties which should be calculated.
285 system_changes : List of properties changed since last run.
287 """
288 # Call base calculator.
289 FileIOCalculator.write_input(
290 self,
291 atoms=atoms,
292 properties=properties,
293 system_changes=system_changes)
295 if system_changes is None and properties is None:
296 return
298 filename = self.label + '/deMon.inp'
300 add_print = ''
302 # Start writing the file.
303 with open(filename, 'w') as fd:
305 # write keyword argument keywords
306 value = self.parameters['title']
307 self._write_argument('TITLE', value, fd)
309 fd.write('#\n')
311 value = self.parameters['scftype']
312 self._write_argument('SCFTYPE', value, fd)
314 value = self.parameters['xc']
315 self._write_argument('VXCTYPE', value, fd)
317 value = self.parameters['guess']
318 self._write_argument('GUESS', value, fd)
320 # obtain forces through a single BOMD step
321 # only if forces is in properties, or if keyword forces is True
322 value = self.parameters['forces']
323 if 'forces' in properties or value:
325 self._write_argument('DYNAMICS',
326 ['INT=1', 'MAX=0', 'STEP=0'], fd)
327 self._write_argument('TRAJECTORY', 'FORCES', fd)
328 self._write_argument('VELOCITIES', 'ZERO', fd)
329 add_print = add_print + ' ' + 'MD OPT'
331 # if dipole is True, enforce dipole calculation.
332 # Otherwise only if asked for
333 value = self.parameters['dipole']
334 if 'dipole' in properties or value:
335 self._write_argument('DIPOLE', '', fd)
337 # print argument, here other options could change this
338 value = self.parameters['print_out']
339 assert isinstance(value, str)
340 value = value + add_print
342 if not len(value) == 0:
343 self._write_argument('PRINT', value, fd)
344 fd.write('#\n')
346 # write general input arguments
347 self._write_input_arguments(fd)
349 fd.write('#\n')
351 # write basis set, ecps, mcps, auxis, augment
352 basis = self.parameters['basis']
353 if 'all' not in basis:
354 basis['all'] = 'DZVP'
355 self._write_basis(fd, atoms, basis, string='BASIS')
357 ecps = self.parameters['ecps']
358 if not len(ecps) == 0:
359 self._write_basis(fd, atoms, ecps, string='ECPS')
361 mcps = self.parameters['mcps']
362 if not len(mcps) == 0:
363 self._write_basis(fd, atoms, mcps, string='MCPS')
365 auxis = self.parameters['auxis']
366 if not len(auxis) == 0:
367 self._write_basis(fd, atoms, auxis, string='AUXIS')
369 augment = self.parameters['augment']
370 if not len(augment) == 0:
371 self._write_basis(fd, atoms, augment, string='AUGMENT')
373 # write geometry
374 self._write_atomic_coordinates(fd, atoms)
376 # write xyz file for good measure.
377 ase.io.write(self.label + '/deMon_atoms.xyz', self.atoms)
379 def read(self, restart_path):
380 """Read parameters from directory restart_path."""
382 self.set_label(restart_path)
384 if not op.exists(restart_path + '/deMon.inp'):
385 raise ReadError('The restart_path file {0} does not exist'
386 .format(restart_path))
388 self.atoms = self.deMon_inp_to_atoms(restart_path + '/deMon.inp')
390 self.read_results()
392 def _write_input_arguments(self, fd):
393 """Write directly given input-arguments."""
394 input_arguments = self.parameters['input_arguments']
396 # Early return
397 if input_arguments is None:
398 return
400 for key, value in input_arguments.items():
401 self._write_argument(key, value, fd)
403 def _write_argument(self, key, value, fd):
404 """Write an argument to file.
405 key : a string coresponding to the input keyword
406 value : the arguments, can be a string, a number or a list
407 f : and open file
408 """
410 # for only one argument, write on same line
411 if not isinstance(value, (tuple, list)):
412 line = key.upper()
413 line += ' ' + str(value).upper()
414 fd.write(line)
415 fd.write('\n')
417 # for a list, write first argument on the first line,
418 # then the rest on new lines
419 else:
420 line = key
421 if not isinstance(value[0], (tuple, list)):
422 for i in range(len(value)):
423 line += ' ' + str(value[i].upper())
424 fd.write(line)
425 fd.write('\n')
426 else:
427 for i in range(len(value)):
428 for j in range(len(value[i])):
429 line += ' ' + str(value[i][j]).upper()
430 fd.write(line)
431 fd.write('\n')
432 line = ''
434 def _write_atomic_coordinates(self, fd, atoms):
435 """Write atomic coordinates.
437 Parameters:
438 - f: An open file object.
439 - atoms: An atoms object.
440 """
442 fd.write('#\n')
443 fd.write('# Atomic coordinates\n')
444 fd.write('#\n')
445 fd.write('GEOMETRY CARTESIAN ANGSTROM\n')
447 for i in range(len(atoms)):
448 xyz = atoms.get_positions()[i]
449 chem_symbol = atoms.get_chemical_symbols()[i]
450 chem_symbol += str(i + 1)
452 # if tag is set to 1 then we have a ghost atom,
453 # set nuclear charge to 0
454 if atoms.get_tags()[i] == 1:
455 nuc_charge = str(0)
456 else:
457 nuc_charge = str(atoms.get_atomic_numbers()[i])
459 mass = atoms.get_masses()[i]
461 line = '{0:6s}'.format(chem_symbol).rjust(10) + ' '
462 line += '{0:.5f}'.format(xyz[0]).rjust(10) + ' '
463 line += '{0:.5f}'.format(xyz[1]).rjust(10) + ' '
464 line += '{0:.5f}'.format(xyz[2]).rjust(10) + ' '
465 line += '{0:5s}'.format(nuc_charge).rjust(10) + ' '
466 line += '{0:.5f}'.format(mass).rjust(10) + ' '
468 fd.write(line)
469 fd.write('\n')
471 # routine to write basis set inormation, including ecps and auxis
472 def _write_basis(self, fd, atoms, basis={}, string='BASIS'):
473 """Write basis set, ECPs, AUXIS, or AUGMENT basis
475 Parameters:
476 - f: An open file object.
477 - atoms: An atoms object.
478 - basis: A dictionary specifying the basis set
479 - string: 'BASIS', 'ECP','AUXIS' or 'AUGMENT'
480 """
482 # basis for all atoms
483 line = '{0}'.format(string).ljust(10)
485 if 'all' in basis:
486 default_basis = basis['all']
487 line += '({0})'.format(default_basis).rjust(16)
489 fd.write(line)
490 fd.write('\n')
492 # basis for all atomic species
493 chemical_symbols = atoms.get_chemical_symbols()
494 chemical_symbols_set = set(chemical_symbols)
496 for i in range(chemical_symbols_set.__len__()):
497 symbol = chemical_symbols_set.pop()
499 if symbol in basis:
500 line = '{0}'.format(symbol).ljust(10)
501 line += '({0})'.format(basis[symbol]).rjust(16)
502 fd.write(line)
503 fd.write('\n')
505 # basis for individual atoms
506 for i in range(len(atoms)):
508 if i in basis:
509 symbol = str(chemical_symbols[i])
510 symbol += str(i + 1)
512 line = '{0}'.format(symbol).ljust(10)
513 line += '({0})'.format(basis[i]).rjust(16)
514 fd.write(line)
515 fd.write('\n')
517 # Analysis routines
518 def read_results(self):
519 """Read the results from output files."""
520 self.read_energy()
521 self.read_forces(self.atoms)
522 self.read_eigenvalues()
523 self.read_dipole()
524 self.read_xray()
526 def read_energy(self):
527 """Read energy from deMon's text-output file."""
528 with open(self.label + '/deMon.out', 'r') as fd:
529 text = fd.read().upper()
531 lines = iter(text.split('\n'))
533 for line in lines:
534 if line.startswith(' TOTAL ENERGY ='):
535 self.results['energy'] = float(line.split()[-1]) * Hartree
536 break
537 else:
538 raise RuntimeError
540 def read_forces(self, atoms):
541 """Read the forces from the deMon.out file."""
543 natoms = len(atoms)
544 filename = self.label + '/deMon.out'
546 if op.isfile(filename):
547 with open(filename, 'r') as fd:
548 lines = fd.readlines()
550 # find line where the orbitals start
551 flag_found = False
552 for i in range(len(lines)):
553 if lines[i].rfind('GRADIENTS OF TIME STEP 0 IN A.U.') > -1:
554 start = i + 4
555 flag_found = True
556 break
558 if flag_found:
559 self.results['forces'] = np.zeros((natoms, 3), float)
560 for i in range(natoms):
561 line = [s for s in lines[i + start].strip().split(' ')
562 if len(s) > 0]
563 f = -np.array([float(x) for x in line[2:5]])
564 self.results['forces'][i, :] = f * (Hartree / Bohr)
566 def read_eigenvalues(self):
567 """Read eigenvalues from the 'deMon.out' file."""
568 assert os.access(self.label + '/deMon.out', os.F_OK)
570 # Read eigenvalues
571 with open(self.label + '/deMon.out', 'r') as fd:
572 lines = fd.readlines()
574 # try PRINT MOE
575 eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(
576 lines, 'ALPHA MO ENERGIES', 6)
577 eig_beta, occ_beta = self.read_eigenvalues_one_spin(
578 lines, 'BETA MO ENERGIES', 6)
580 # otherwise try PRINT MOS
581 if len(eig_alpha) == 0 and len(eig_beta) == 0:
582 eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(
583 lines, 'ALPHA MO COEFFICIENTS', 5)
584 eig_beta, occ_beta = self.read_eigenvalues_one_spin(
585 lines, 'BETA MO COEFFICIENTS', 5)
587 self.results['eigenvalues'] = np.array([eig_alpha, eig_beta]) * Hartree
588 self.results['occupations'] = np.array([occ_alpha, occ_beta])
590 def read_eigenvalues_one_spin(self, lines, string, neigs_per_line):
591 """Utility method for retreiving eigenvalues after the string "string"
592 with neigs_per_line eigenvlaues written per line
593 """
594 eig = []
595 occ = []
597 skip_line = False
598 more_eigs = False
600 # find line where the orbitals start
601 for i in range(len(lines)):
602 if lines[i].rfind(string) > -1:
603 ii = i
604 more_eigs = True
605 break
607 while more_eigs:
608 # search for two empty lines in a row preceding a line with
609 # numbers
610 for i in range(ii + 1, len(lines)):
611 if len(lines[i].split()) == 0 and \
612 len(lines[i + 1].split()) == 0 and \
613 len(lines[i + 2].split()) > 0:
614 ii = i + 2
615 break
617 # read eigenvalues, occupations
618 line = lines[ii].split()
619 if len(line) < neigs_per_line:
620 # last row
621 more_eigs = False
622 if line[0] != str(len(eig) + 1):
623 more_eigs = False
624 skip_line = True
626 if not skip_line:
627 line = lines[ii + 1].split()
628 for l in line:
629 eig.append(float(l))
630 line = lines[ii + 3].split()
631 for l in line:
632 occ.append(float(l))
633 ii = ii + 3
635 return eig, occ
637 def read_dipole(self):
638 """Read dipole moment."""
639 dipole = np.zeros(3)
640 with open(self.label + '/deMon.out', 'r') as fd:
641 lines = fd.readlines()
643 for i in range(len(lines)):
644 if lines[i].rfind('DIPOLE') > - \
645 1 and lines[i].rfind('XAS') == -1:
646 dipole[0] = float(lines[i + 1].split()[3])
647 dipole[1] = float(lines[i + 2].split()[3])
648 dipole[2] = float(lines[i + 3].split()[3])
650 # debye to e*Ang
651 self.results['dipole'] = dipole * 0.2081943482534
653 break
655 def read_xray(self):
656 """Read deMon.xry if present."""
658 # try to read core IP from, .out file
659 filename = self.label + '/deMon.out'
660 core_IP = None
661 if op.isfile(filename):
662 with open(filename, 'r') as fd:
663 lines = fd.readlines()
665 for i in range(len(lines)):
666 if lines[i].rfind('IONIZATION POTENTIAL') > -1:
667 core_IP = float(lines[i].split()[3])
669 try:
670 mode, ntrans, E_trans, osc_strength, trans_dip = parse_xray(
671 self.label + '/deMon.xry')
672 except ReadError:
673 pass
674 else:
675 xray_results = {'xray_mode': mode,
676 'ntrans': ntrans,
677 'E_trans': E_trans,
678 'osc_strength': osc_strength, # units?
679 'trans_dip': trans_dip, # units?
680 'core_IP': core_IP}
682 self.results['xray'] = xray_results
684 def deMon_inp_to_atoms(self, filename):
685 """Routine to read deMon.inp and convert it to an atoms object."""
687 with open(filename, 'r') as fd:
688 lines = fd.readlines()
690 # find line where geometry starts
691 for i in range(len(lines)):
692 if lines[i].rfind('GEOMETRY') > -1:
693 if lines[i].rfind('ANGSTROM'):
694 coord_units = 'Ang'
695 elif lines.rfind('Bohr'):
696 coord_units = 'Bohr'
697 ii = i
698 break
700 chemical_symbols = []
701 xyz = []
702 atomic_numbers = []
703 masses = []
705 for i in range(ii + 1, len(lines)):
706 try:
707 line = lines[i].split()
709 if len(line) > 0:
710 for symbol in ase.data.chemical_symbols:
711 found = None
712 if line[0].upper().rfind(symbol.upper()) > -1:
713 found = symbol
714 break
716 if found is not None:
717 chemical_symbols.append(found)
718 else:
719 break
721 xyz.append(
722 [float(line[1]), float(line[2]), float(line[3])])
724 if len(line) > 4:
725 atomic_numbers.append(int(line[4]))
727 if len(line) > 5:
728 masses.append(float(line[5]))
730 except Exception: # XXX Which kind of exception?
731 raise RuntimeError
733 if coord_units == 'Bohr':
734 xyz = xyz * Bohr
736 natoms = len(chemical_symbols)
738 # set atoms object
739 atoms = ase.Atoms(symbols=chemical_symbols, positions=xyz)
741 # if atomic numbers were read in, set them
742 if len(atomic_numbers) == natoms:
743 atoms.set_atomic_numbers(atomic_numbers)
745 # if masses were read in, set them
746 if len(masses) == natoms:
747 atoms.set_masses(masses)
749 return atoms