Coverage for /builds/kinetik161/ase/ase/calculators/demonnano.py: 44.37%

1# flake8: noqa 

2"""This module defines an ASE interface to deMon-nano. 

3 

4Link to the open-source DFTB code deMon-nano: 

5http://demon-nano.ups-tlse.fr/ 

6 

7export ASE_DEMONNANO_COMMAND="/path/to/bin/deMon.username.x" 

8export DEMONNANO_BASIS_PATH="/path/to/basis/" 

9 

10The file 'deMon.inp' contains the input geometry and parameters 

11The file 'deMon.out' contains the results 

12 

13""" 

14import os 

15import os.path as op 

16# import subprocess 

17import pathlib as pl 

18 

19import numpy as np 

20 

21import ase.data 

22import ase.io 

23from ase.calculators.calculator import FileIOCalculator, Parameters, ReadError 

24from ase.units import Bohr, Hartree 

25 

26 

27class DemonNanoParameters(Parameters): 

28 """Parameters class for the calculator. 

29 

30 The options here are the most important ones that the user needs to be 

31 aware of. Further options accepted by deMon can be set in the dictionary 

32 input_arguments. 

33 

34 """ 

35 

36 def __init__( 

37 self, 

38 label='.', 

39 atoms=None, 

40 command=None, 

41 basis_path=None, 

42 restart_path='.', 

43 print_out='ASE', 

44 title='deMonNano input file', 

45 forces=False, 

46 input_arguments=None): 

47 kwargs = locals() 

48 kwargs.pop('self') 

49 Parameters.__init__(self, **kwargs) 

50 

51 

52class DemonNano(FileIOCalculator): 

53 """Calculator interface to the deMon-nano code. """ 

54 

55 implemented_properties = ['energy', 'forces'] 

56 

57 def __init__(self, **kwargs): 

58 """ASE interface to the deMon-nano code. 

59 

60 The deMon-nano code can be obtained from http://demon-nano.ups-tlse.fr/ 

61 

62 The ASE_DEMONNANO_COMMAND environment variable must be set to run the executable, in bash it would be set along the lines of 

63 export ASE_DEMONNANO_COMMAND="pathway-to-deMon-binary/deMon.username.x" 

64 

65 Parameters: 

66 

67 label : str 

68 relative path to the run directory 

69 atoms : Atoms object 

70 the atoms object 

71 command : str 

72 Command to run deMon. If not present, the environment variable ASE_DEMONNANO_COMMAND is used 

73 basis_path : str 

74 Relative path to the directory containing DFTB-SCC or DFTB-0 parameters 

75 If not present, the environment variable DEMONNANO_BASIS_PATH is used 

76 restart_path : str 

77 Relative path to the deMon restart dir 

78 title : str 

79 Title in the deMon input file. 

80 forces : bool 

81 If True a force calculation is enforced 

82 print_out : str | list 

83 Options for the printing in deMon 

84 input_arguments : dict 

85 Explicitly given input arguments. The key is the input keyword 

86 and the value is either a str, a list of str (will be written on the same line as the keyword), 

87 or a list of lists of str (first list is written on the first line, the others on following lines.) 

88 """ 

89 

90 parameters = DemonNanoParameters(**kwargs) 

91 

92 # basis path 

93 basis_path = parameters['basis_path'] 

94 if basis_path is None: 

95 basis_path = self.cfg.get('DEMONNANO_BASIS_PATH') 

96 

97 if basis_path is None: 

98 mess = 'The "DEMONNANO_BASIS_PATH" environment is not defined.' 

99 raise ValueError(mess) 

100 else: 

101 parameters['basis_path'] = basis_path 

102 

103 # Call the base class. 

104 FileIOCalculator.__init__( 

105 self, 

106 **parameters) 

107 

108 def __getitem__(self, key): 

109 """Convenience method to retrieve a parameter as 

110 calculator[key] rather than calculator.parameters[key] 

111 

112 Parameters: 

113 key : str, the name of the parameters to get. 

114 """ 

115 return self.parameters[key] 

116 

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

118 """Write input (in)-file. 

119 See calculator.py for further details. 

120 

121 Parameters: 

122 atoms : The Atoms object to write. 

123 properties : The properties which should be calculated. 

124 system_changes : List of properties changed since last run. 

125 

126 """ 

127 # Call base calculator. 

128 FileIOCalculator.write_input( 

129 self, 

130 atoms=atoms, 

131 properties=properties, 

132 system_changes=system_changes) 

133 

134 if system_changes is None and properties is None: 

135 return 

136 

137 filename = self.label + '/deMon.inp' 

138 

139 # Start writing the file. 

140 with open(filename, 'w') as fd: 

141 # write keyword argument keywords 

142 value = self.parameters['title'] 

143 self._write_argument('TITLE', value, fd) 

144 fd.write('\n') 

145 

146 # obtain forces through a single BOMD step 

147 # only if forces is in properties, or if keyword forces is True 

148 value = self.parameters['forces'] 

149 if 'forces' in properties or value: 

150 self._write_argument('MDYNAMICS', 'ZERO', fd) 

151 self._write_argument('MDSTEP', 'MAX=1', fd) 

152 # default timestep is 0.25 fs if not enough - uncomment the line below 

153 # self._write_argument('TIMESTEP', '0.1', fd) 

154 

155 # print argument, here other options could change this 

156 value = self.parameters['print_out'] 

157 assert (isinstance(value, str)) 

158 

159 if not len(value) == 0: 

160 self._write_argument('PRINT', value, fd) 

161 fd.write('\n') 

162 

163 # write general input arguments 

164 self._write_input_arguments(fd) 

165 

166 if 'BASISPATH' not in self.parameters['input_arguments']: 

167 value = self.parameters['basis_path'] 

168 fd.write(value) 

169 fd.write('\n') 

170 

171 # write geometry 

172 self._write_atomic_coordinates(fd, atoms) 

173 

174 # write xyz file for good measure. 

175 ase.io.write(self.label + '/deMon_atoms.xyz', self.atoms) 

176 

177 def read(self, restart_path): 

178 """Read parameters from directory restart_path.""" 

179 

180 self.set_label(restart_path) 

181 rpath = pl.Path(restart_path) 

182 

183 if not (rpath / 'deMon.inp').exists(): 

184 raise ReadError('The restart_path file {} does not exist' 

185 .format(rpath)) 

186 

187 self.atoms = self.deMon_inp_to_atoms(rpath / 'deMon.inp') 

188 

189 self.read_results() 

190 

191 def _write_input_arguments(self, fd): 

192 """Write directly given input-arguments.""" 

193 input_arguments = self.parameters['input_arguments'] 

194 

195 # Early return 

196 if input_arguments is None: 

197 return 

198 

199 for key, value in input_arguments.items(): 

200 self._write_argument(key, value, fd) 

201 

202 def _write_argument(self, key, value, fd): 

203 """Write an argument to file. 

204 key : a string coresponding to the input keyword 

205 value : the arguments, can be a string, a number or a list 

206 fd : and open file 

207 """ 

208 if key == 'BASISPATH': 

209 # Write a basis path to file. 

210 # Has to be in lowercase for deMon-nano to work 

211 line = value.lower() 

212 fd.write(line) 

213 fd.write('\n') 

214 elif not isinstance(value, (tuple, list)): 

215 # for only one argument, write on same line 

216 line = key.upper() 

217 line += ' ' + str(value).upper() 

218 fd.write(line) 

219 fd.write('\n') 

220 

221 # for a list, write first argument on the first line, 

222 # then the rest on new lines 

223 else: 

224 line = key 

225 if not isinstance(value[0], (tuple, list)): 

226 for i in range(len(value)): 

227 line += ' ' + str(value[i].upper()) 

228 fd.write(line) 

229 fd.write('\n') 

230 else: 

231 for i in range(len(value)): 

232 for j in range(len(value[i])): 

233 line += ' ' + str(value[i][j]).upper() 

234 fd.write(line) 

235 fd.write('\n') 

236 line = '' 

237 

238 def _write_atomic_coordinates(self, fd, atoms): 

239 """Write atomic coordinates. 

240 Parameters: 

241 - fd: An open file object. 

242 - atoms: An atoms object. 

243 """ 

244 # fd.write('#\n') 

245 # fd.write('# Atomic coordinates\n') 

246 # fd.write('#\n') 

247 fd.write('GEOMETRY CARTESIAN ANGSTROM\n') 

248 

249 for sym, pos in zip(atoms.symbols, atoms.positions): 

250 fd.write('{:9s} {:10.5f} {:10.5f} {:10.5f}\n'.format(sym, *pos)) 

251 

252 fd.write('\n') 

253 

254# Analysis routines 

255 def read_results(self): 

256 """Read the results from output files.""" 

257 self.read_energy() 

258 self.read_forces(self.atoms) 

259 # self.read_eigenvalues() 

260 

261 def read_energy(self): 

262 """Read energy from deMon.ase output file.""" 

263 

264 epath = pl.Path(self.label) 

265 

266 if not (epath / 'deMon.ase').exists(): 

267 raise ReadError('The deMonNano output file for ASE {} does not exist' 

268 .format(epath)) 

269 

270 filename = self.label + '/deMon.ase' 

271 

272 if op.isfile(filename): 

273 with open(filename) as fd: 

274 lines = fd.readlines() 

275 

276 for i in range(len(lines)): 

277 if lines[i].startswith(' DFTB total energy [Hartree]'): 

278 self.results['energy'] = float(lines[i + 1]) * Hartree 

279 break 

280 

281 def read_forces(self, atoms): 

282 """Read forces from the deMon.ase file.""" 

283 

284 natoms = len(atoms) 

285 epath = pl.Path(self.label) 

286 

287 if not (epath / 'deMon.ase').exists(): 

288 raise ReadError('The deMonNano output file for ASE {} does not exist' 

289 .format(epath)) 

290 

291 filename = self.label + '/deMon.ase' 

292 

293 with open(filename) as fd: 

294 lines = fd.readlines() 

295 

296 # find line where the forces start 

297 flag_found = False 

298 for i in range(len(lines)): 

299 if 'DFTB gradients at 0 time step in a.u.' in lines[i]: 

300 start = i + 1 

301 flag_found = True 

302 break 

303 

304 if flag_found: 

305 self.results['forces'] = np.zeros((natoms, 3), float) 

306 for i in range(natoms): 

307 line = [s for s in lines[i + start].strip().split(' ') 

308 if len(s) > 0] 

309 f = -np.array([float(x) for x in line[1:4]]) 

310 # output forces in a.u. 

311 # self.results['forces'][i, :] = f 

312 # output forces with real dimension 

313 self.results['forces'][i, :] = f * (Hartree / Bohr) 

314 

315 def deMon_inp_to_atoms(self, filename): 

316 """Routine to read deMon.inp and convert it to an atoms object.""" 

317 

318 read_flag = False 

319 chem_symbols = [] 

320 xyz = [] 

321 

322 with open(filename) as fd: 

323 for line in fd: 

324 if 'GEOMETRY' in line: 

325 read_flag = True 

326 if 'ANGSTROM' in line: 

327 coord_units = 'Ang' 

328 elif 'BOHR' in line: 

329 coord_units = 'Bohr' 

330 

331 if read_flag: 

332 tokens = line.split() 

333 symbol = tokens[0] 

334 xyz_loc = np.array(tokens[1:4]).astype(float) 

335 if read_flag and tokens: 

336 chem_symbols.append(symbol) 

337 xyz.append(xyz_loc) 

338 

339 if coord_units == 'Bohr': 

340 xyz *= Bohr 

341 

342 # set atoms object 

343 atoms = ase.Atoms(symbols=chem_symbols, positions=xyz) 

344 

345 return atoms