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 

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

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

12from ase.calculators.genericfileio import (CalculatorTemplate, 

13 GenericFileIOCalculator) 

14 

15 

16class OctopusIOError(IOError): 

17 pass 

18 

19 

20class OctopusProfile: 

21 def __init__(self, argv): 

22 self.argv = argv 

23 

24 def version(self): 

25 from subprocess import check_output 

26 import re 

27 txt = check_output(self.argv + ['--version']).decode('ascii') 

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

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

30 # the first line. 

31 return match.group(1) 

32 

33 def run(self, directory, outputfile): 

34 from subprocess import check_call 

35 with open(directory / outputfile, 'w') as fd: 

36 check_call(self.argv, stdout=fd, cwd=directory) 

37 

38 

39class OctopusTemplate(CalculatorTemplate): 

40 def __init__(self): 

41 super().__init__( 

42 name='octopus', 

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

44 ) 

45 

46 def read_results(self, directory): 

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

48 results = {} 

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

50 results.update(read_static_info(fd)) 

51 

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

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

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

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

56 eigpath = directory / 'static/eigenvalues' 

57 if eigpath.is_file(): 

58 with open(eigpath) as fd: 

59 kpts, eigs, occs = read_eigenvalues_file(fd) 

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

61 # XXX New Octopus probably has symmetry reduction !! 

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

63 ibz_k_points=kpts, 

64 k_point_weights=kpt_weights) 

65 return results 

66 

67 def execute(self, directory, profile): 

68 profile.run(directory, 'octopus.out') 

69 

70 def write_input(self, directory, atoms, parameters, properties): 

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

72 inp = directory / 'inp' 

73 inp.write_text(txt) 

74 

75 

76class Octopus(GenericFileIOCalculator): 

77 """Octopus calculator. 

78 

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

80 

81 def __init__(self, 

82 profile=None, 

83 directory='.', 

84 **kwargs): 

85 """Create Octopus calculator. 

86 

87 Label is always taken as a subdirectory. 

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

89 

90 if profile is None: 

91 profile = OctopusProfile(['octopus']) 

92 

93 super().__init__(profile=profile, template=OctopusTemplate(), 

94 directory=directory, 

95 parameters=kwargs) 

96 

97 @classmethod 

98 def recipe(cls, **kwargs): 

99 from ase import Atoms 

100 system = Atoms() 

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

102 system.calc = calc 

103 try: 

104 system.get_potential_energy() 

105 except OctopusIOError: 

106 pass 

107 else: 

108 raise OctopusIOError('Expected recipe, but found ' 

109 'useful physical output!')