Coverage for /builds/debichem-team/python-ase/ase/calculators/octopus.py: 75.93%

1"""ASE-interface to Octopus. 

2 

3Ask Hjorth Larsen <asklarsen@gmail.com> 

4Carlos de Armas 

5 

6http://tddft.org/programs/octopus/ 

7""" 

8 

9import numpy as np 

10 

11from ase.calculators.genericfileio import ( 

12 BaseProfile, 

13 CalculatorTemplate, 

14 GenericFileIOCalculator, 

15) 

16from ase.io.octopus.input import generate_input, process_special_kwargs 

17from ase.io.octopus.output import read_eigenvalues_file, read_static_info 

18 

19 

20class OctopusIOError(IOError): 

21 pass 

22 

23 

24class OctopusProfile(BaseProfile): 

25 def get_calculator_command(self, inputfile): 

26 return [] 

27 

28 def version(self): 

29 import re 

30 from subprocess import check_output 

31 txt = check_output([*self._split_command, '--version'], 

32 encoding='ascii') 

33 match = re.match(r'octopus\s*(.+)', txt) 

34 # With MPI it prints the line for each rank, but we just match 

35 # the first line. 

36 return match.group(1) 

37 

38 

39class OctopusTemplate(CalculatorTemplate): 

40 _label = 'octopus' 

41 

42 def __init__(self): 

43 super().__init__( 

44 'octopus', 

45 implemented_properties=['energy', 'forces', 'dipole', 'stress'], 

46 ) 

47 self.outputname = f'{self._label}.out' 

48 self.errorname = f'{self._label}.err' 

49 

50 def read_results(self, directory): 

51 """Read octopus output files and extract data.""" 

52 results = {} 

53 with open(directory / 'static/info') as fd: 

54 results.update(read_static_info(fd)) 

55 

56 # If the eigenvalues file exists, we get the eigs/occs from that one. 

57 # This probably means someone ran Octopus in 'unocc' mode to 

58 # get eigenvalues (e.g. for band structures), and the values in 

59 # static/info will be the old (selfconsistent) ones. 

60 eigpath = directory / 'static/eigenvalues' 

61 if eigpath.is_file(): 

62 with open(eigpath) as fd: 

63 kpts, eigs, occs = read_eigenvalues_file(fd) 

64 kpt_weights = np.ones(len(kpts)) # XXX ? Or 1 / len(kpts) ? 

65 # XXX New Octopus probably has symmetry reduction !! 

66 results.update(eigenvalues=eigs, occupations=occs, 

67 ibz_kpoints=kpts, 

68 kpoint_weights=kpt_weights) 

69 return results 

70 

71 def execute(self, directory, profile): 

72 profile.run(directory, None, self.outputname, 

73 errorfile=self.errorname) 

74 

75 def write_input(self, profile, directory, atoms, parameters, properties): 

76 txt = generate_input(atoms, process_special_kwargs(atoms, parameters)) 

77 inp = directory / 'inp' 

78 inp.write_text(txt) 

79 

80 def load_profile(self, cfg, **kwargs): 

81 return OctopusProfile.from_config(cfg, self.name, **kwargs) 

82 

83 

84class Octopus(GenericFileIOCalculator): 

85 """Octopus calculator. 

86 

87 The label is always assumed to be a directory.""" 

88 

89 def __init__(self, profile=None, directory='.', **kwargs): 

90 """Create Octopus calculator. 

91 

92 Label is always taken as a subdirectory. 

93 Restart is taken to be a label.""" 

94 

95 super().__init__(profile=profile, 

96 template=OctopusTemplate(), 

97 directory=directory, 

98 parameters=kwargs) 

99 

100 @classmethod 

101 def recipe(cls, **kwargs): 

102 from ase import Atoms 

103 system = Atoms() 

104 calc = Octopus(CalculationMode='recipe', **kwargs) 

105 system.calc = calc 

106 try: 

107 system.get_potential_energy() 

108 except OctopusIOError: 

109 pass 

110 else: 

111 raise OctopusIOError('Expected recipe, but found ' 

112 'useful physical output!')