1"""ASE-interface to Octopus. 

2 

3Ask Hjorth Larsen <asklarsen@gmail.com> 

4Carlos de Armas 

5 

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

7""" 

8 

9import os 

10import numpy as np 

11from ase.io.octopus.input import ( 

12 process_special_kwargs, kwargs2atoms, 

13 generate_input, parse_input_file, 

14 normalize_keywords) 

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

16from ase.calculators.calculator import ( 

17 FileIOCalculator, EigenvalOccupationMixin, PropertyNotImplementedError) 

18 

19 

20class OctopusIOError(IOError): 

21 pass 

22 

23 

24class Octopus(FileIOCalculator, EigenvalOccupationMixin): 

25 """Octopus calculator. 

26 

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

28 

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

30 command = 'octopus' 

31 

32 def __init__(self, 

33 restart=None, 

34 label=None, 

35 directory=None, 

36 atoms=None, 

37 command=None, 

38 **kwargs): 

39 """Create Octopus calculator. 

40 

41 Label is always taken as a subdirectory. 

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

43 

44 kwargs.pop('check_keywords', None) # Ignore old keywords 

45 kwargs.pop('troublesome_keywords', None) 

46 

47 if label is not None: 

48 # restart mechanism in Calculator tends to set the label. 

49 #import warnings 

50 #warnings.warn('Please use directory=... instead of label') 

51 directory = label.rstrip('/') 

52 

53 if directory is None: 

54 directory = 'ink-pool' 

55 

56 self.kwargs = {} 

57 

58 FileIOCalculator.__init__(self, restart=restart, 

59 directory=directory, 

60 atoms=atoms, 

61 command=command, **kwargs) 

62 # The above call triggers set() so we can update self.kwargs. 

63 

64 def set(self, **kwargs): 

65 """Set octopus input file parameters.""" 

66 kwargs = normalize_keywords(kwargs) 

67 changes = FileIOCalculator.set(self, **kwargs) 

68 if changes: 

69 self.results.clear() 

70 self.kwargs.update(kwargs) 

71 # XXX should use 'Parameters' but don't know how 

72 

73 def get_xc_functional(self): 

74 """Return the XC-functional identifier. 

75 'LDA', 'PBE', ...""" 

76 return self.kwargs.get('xcfunctional', 'LDA') 

77 

78 def get_bz_k_points(self): 

79 """Return all the k-points in the 1. Brillouin zone. 

80 The coordinates are relative to reciprocal latice vectors.""" 

81 # Have not found nice way of extracting this information 

82 # from Octopus. Thus unimplemented. -askhl 

83 raise NotImplementedError 

84 

85 def get_charges(self, atoms=None): 

86 raise PropertyNotImplementedError 

87 

88 def get_fermi_level(self): 

89 return self.results['efermi'] 

90 

91 def get_potential_energies(self): 

92 raise PropertyNotImplementedError 

93 

94 def get_dipole_moment(self, atoms=None): 

95 if 'dipole' not in self.results: 

96 msg = ('Dipole moment not calculated.\n' 

97 'You may wish to use SCFCalculateDipole=True') 

98 raise OctopusIOError(msg) 

99 return self.results['dipole'] 

100 

101 def get_stresses(self): 

102 raise PropertyNotImplementedError 

103 

104 def get_number_of_spins(self): 

105 """Return the number of spins in the calculation. 

106 Spin-paired calculations: 1, spin-polarized calculation: 2.""" 

107 return 2 if self.get_spin_polarized() else 1 

108 

109 def get_spin_polarized(self): 

110 """Is it a spin-polarized calculation?""" 

111 

112 sc = self.kwargs.get('spincomponents') 

113 if sc is None or sc == 'unpolarized': 

114 return False 

115 elif sc == 'spin_polarized' or sc == 'polarized': 

116 return True 

117 else: 

118 raise NotImplementedError('SpinComponents keyword %s' % sc) 

119 

120 def get_ibz_k_points(self): 

121 """Return k-points in the irreducible part of the Brillouin zone. 

122 The coordinates are relative to reciprocal latice vectors.""" 

123 return self.results['ibz_k_points'] 

124 

125 def get_k_point_weights(self): 

126 return self.results['k_point_weights'] 

127 

128 def get_number_of_bands(self): 

129 return self.results['nbands'] 

130 

131 #def get_magnetic_moments(self, atoms=None): 

132 # if self.results['nspins'] == 1: 

133 # return np.zeros(len(self.atoms)) 

134 # return self.results['magmoms'].copy() 

135 

136 #def get_magnetic_moment(self, atoms=None): 

137 # if self.results['nspins'] == 1: 

138 # return 0.0 

139 # return self.results['magmom'] 

140 

141 def get_occupation_numbers(self, kpt=0, spin=0): 

142 return self.results['occupations'][kpt, spin].copy() 

143 

144 def get_eigenvalues(self, kpt=0, spin=0): 

145 return self.results['eigenvalues'][kpt, spin].copy() 

146 

147 def _getpath(self, path, check=False): 

148 path = os.path.join(self.directory, path) 

149 if check: 

150 if not os.path.exists(path): 

151 raise OctopusIOError('No such file or directory: %s' % path) 

152 return path 

153 

154 def get_atoms(self): 

155 return FileIOCalculator.get_atoms(self) 

156 

157 def read_results(self): 

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

159 with open(self._getpath('static/info', check=True)) as fd: 

160 self.results.update(read_static_info(fd)) 

161 

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

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

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

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

166 try: 

167 eigpath = self._getpath('static/eigenvalues', check=True) 

168 except OctopusIOError: 

169 pass 

170 else: 

171 with open(eigpath) as fd: 

172 kpts, eigs, occs = read_eigenvalues_file(fd) 

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

174 self.results.update(eigenvalues=eigs, occupations=occs, 

175 ibz_k_points=kpts, 

176 k_point_weights=kpt_weights) 

177 

178 def write_input(self, atoms, properties=None, system_changes=None): 

179 FileIOCalculator.write_input(self, atoms, properties=properties, 

180 system_changes=system_changes) 

181 txt = generate_input(atoms, process_special_kwargs(atoms, self.kwargs)) 

182 with open(self._getpath('inp'), 'w') as fd: 

183 fd.write(txt) 

184 

185 def read(self, directory): 

186 # XXX label of restart file may not be the same as actual label! 

187 # This makes things rather tricky. We first set the label to 

188 # that of the restart file and arbitrarily expect the remaining code 

189 # to rectify any consequent inconsistencies. 

190 self.directory = directory 

191 

192 inp_path = self._getpath('inp') 

193 with open(inp_path) as fd: 

194 kwargs = parse_input_file(fd) 

195 

196 self.atoms, kwargs = kwargs2atoms(kwargs) 

197 self.kwargs.update(kwargs) 

198 

199 self.read_results() 

200 

201 @classmethod 

202 def recipe(cls, **kwargs): 

203 from ase import Atoms 

204 system = Atoms() 

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

206 system.calc = calc 

207 try: 

208 system.get_potential_energy() 

209 except OctopusIOError: 

210 pass 

211 else: 

212 raise OctopusIOError('Expected recipe, but found ' 

213 'useful physical output!')