1"""This module defines an ASE interface to GULP. 

2 

3Written by: 

4 

5Andy Cuko <andi.cuko@upmc.fr> 

6Antoni Macia <tonimacia@gmail.com> 

7 

8EXPORT ASE_GULP_COMMAND="/path/to/gulp < PREFIX.gin > PREFIX.got" 

9 

10Keywords 

11Options 

12 

13""" 

14import os 

15import re 

16import numpy as np 

17from ase.units import eV, Ang 

18from ase.calculators.calculator import FileIOCalculator, ReadError 

19 

20 

21class GULPOptimizer: 

22 def __init__(self, atoms, calc): 

23 self.atoms = atoms 

24 self.calc = calc 

25 

26 def todict(self): 

27 return {'type': 'optimization', 

28 'optimizer': 'GULPOptimizer'} 

29 

30 def run(self, fmax=None, steps=None, **gulp_kwargs): 

31 if fmax is not None: 

32 gulp_kwargs['gmax'] = fmax 

33 if steps is not None: 

34 gulp_kwargs['maxcyc'] = steps 

35 

36 self.calc.set(**gulp_kwargs) 

37 self.atoms.calc = self.calc 

38 self.atoms.get_potential_energy() 

39 self.atoms.cell = self.calc.get_atoms().cell 

40 self.atoms.positions[:] = self.calc.get_atoms().positions 

41 

42 

43class GULP(FileIOCalculator): 

44 implemented_properties = ['energy', 'free_energy', 'forces', 'stress'] 

45 command = 'gulp < PREFIX.gin > PREFIX.got' 

46 discard_results_on_any_change = True 

47 default_parameters = dict( 

48 keywords='conp gradients', 

49 options=[], 

50 shel=[], 

51 library="ffsioh.lib", 

52 conditions=None) 

53 

54 def get_optimizer(self, atoms): 

55 gulp_keywords = self.parameters.keywords.split() 

56 if 'opti' not in gulp_keywords: 

57 raise ValueError('Can only create optimizer from GULP calculator ' 

58 'with "opti" keyword. Current keywords: {}' 

59 .format(gulp_keywords)) 

60 

61 opt = GULPOptimizer(atoms, self) 

62 return opt 

63 

64 def __init__(self, restart=None, 

65 ignore_bad_restart_file=FileIOCalculator._deprecated, 

66 label='gulp', atoms=None, optimized=None, 

67 Gnorm=1000.0, steps=1000, conditions=None, **kwargs): 

68 """Construct GULP-calculator object.""" 

69 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

70 label, atoms, **kwargs) 

71 self.optimized = optimized 

72 self.Gnorm = Gnorm 

73 self.steps = steps 

74 self.conditions = conditions 

75 self.library_check() 

76 self.atom_types = [] 

77 

78 # GULP prints the fractional coordinates before the Final 

79 # lattice vectors so they need to be stored and then atoms 

80 # positions need to be set after we get the Final lattice 

81 # vectors 

82 self.fractional_coordinates = None 

83 

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

85 FileIOCalculator.write_input(self, atoms, properties, system_changes) 

86 p = self.parameters 

87 

88 # Build string to hold .gin input file: 

89 s = p.keywords 

90 s += '\ntitle\nASE calculation\nend\n\n' 

91 

92 if all(self.atoms.pbc): 

93 cell_params = self.atoms.cell.cellpar() 

94 # Formating is necessary since Gulp max-line-length restriction 

95 s += 'cell\n{0:9.6f} {1:9.6f} {2:9.6f} ' \ 

96 '{3:8.5f} {4:8.5f} {5:8.5f}\n'.format(*cell_params) 

97 s += 'frac\n' 

98 coords = self.atoms.get_scaled_positions() 

99 else: 

100 s += 'cart\n' 

101 coords = self.atoms.get_positions() 

102 

103 if self.conditions is not None: 

104 c = self.conditions 

105 labels = c.get_atoms_labels() 

106 self.atom_types = c.get_atom_types() 

107 else: 

108 labels = self.atoms.get_chemical_symbols() 

109 

110 for xyz, symbol in zip(coords, labels): 

111 s += ' {0:2} core' \ 

112 ' {1:10.7f} {2:10.7f} {3:10.7f}\n' .format(symbol, *xyz) 

113 if symbol in p.shel: 

114 s += ' {0:2} shel' \ 

115 ' {1:10.7f} {2:10.7f} {3:10.7f}\n' .format(symbol, *xyz) 

116 

117 if p.library: 

118 s += '\nlibrary {0}\n'.format(p.library) 

119 

120 if p.options: 

121 for t in p.options: 

122 s += '%s\n' % t 

123 with open(self.prefix + '.gin', 'w') as fd: 

124 fd.write(s) 

125 

126 def read_results(self): 

127 FileIOCalculator.read(self, self.label) 

128 if not os.path.isfile(self.label + '.got'): 

129 raise ReadError 

130 

131 with open(self.label + '.got') as fd: 

132 lines = fd.readlines() 

133 

134 cycles = -1 

135 self.optimized = None 

136 for i, line in enumerate(lines): 

137 m = re.match(r'\s*Total lattice energy\s*=\s*(\S+)\s*eV', line) 

138 if m: 

139 energy = float(m.group(1)) 

140 self.results['energy'] = energy 

141 self.results['free_energy'] = energy 

142 

143 elif line.find('Optimisation achieved') != -1: 

144 self.optimized = True 

145 

146 elif line.find('Final Gnorm') != -1: 

147 self.Gnorm = float(line.split()[-1]) 

148 

149 elif line.find('Cycle:') != -1: 

150 cycles += 1 

151 

152 elif line.find('Final Cartesian derivatives') != -1: 

153 s = i + 5 

154 forces = [] 

155 while True: 

156 s = s + 1 

157 if lines[s].find("------------") != -1: 

158 break 

159 if lines[s].find(" s ") != -1: 

160 continue 

161 g = lines[s].split()[3:6] 

162 G = [-float(x) * eV / Ang for x in g] 

163 forces.append(G) 

164 forces = np.array(forces) 

165 self.results['forces'] = forces 

166 

167 elif line.find('Final internal derivatives') != -1: 

168 s = i + 5 

169 forces = [] 

170 while True: 

171 s = s + 1 

172 if lines[s].find("------------") != -1: 

173 break 

174 g = lines[s].split()[3:6] 

175 

176 # Uncomment the section below to separate the numbers when 

177 # there is no space between them, in the case of long 

178 # numbers. This prevents the code to break if numbers are 

179 # too big. 

180 

181 '''for t in range(3-len(g)): 

182 g.append(' ') 

183 for j in range(2): 

184 min_index=[i+1 for i,e in enumerate(g[j][1:]) 

185 if e == '-'] 

186 if j==0 and len(min_index) != 0: 

187 if len(min_index)==1: 

188 g[2]=g[1] 

189 g[1]=g[0][min_index[0]:] 

190 g[0]=g[0][:min_index[0]] 

191 else: 

192 g[2]=g[0][min_index[1]:] 

193 g[1]=g[0][min_index[0]:min_index[1]] 

194 g[0]=g[0][:min_index[0]] 

195 break 

196 if j==1 and len(min_index) != 0: 

197 g[2]=g[1][min_index[0]:] 

198 g[1]=g[1][:min_index[0]]''' 

199 

200 G = [-float(x) * eV / Ang for x in g] 

201 forces.append(G) 

202 forces = np.array(forces) 

203 self.results['forces'] = forces 

204 

205 elif line.find('Final cartesian coordinates of atoms') != -1: 

206 s = i + 5 

207 positions = [] 

208 while True: 

209 s = s + 1 

210 if lines[s].find("------------") != -1: 

211 break 

212 if lines[s].find(" s ") != -1: 

213 continue 

214 xyz = lines[s].split()[3:6] 

215 XYZ = [float(x) * Ang for x in xyz] 

216 positions.append(XYZ) 

217 positions = np.array(positions) 

218 self.atoms.set_positions(positions) 

219 

220 elif line.find('Final stress tensor components') != -1: 

221 res = [0., 0., 0., 0., 0., 0.] 

222 for j in range(3): 

223 var = lines[i + j + 3].split()[1] 

224 res[j] = float(var) 

225 var = lines[i + j + 3].split()[3] 

226 res[j + 3] = float(var) 

227 stress = np.array(res) 

228 self.results['stress'] = stress 

229 

230 elif line.find('Final Cartesian lattice vectors') != -1: 

231 lattice_vectors = np.zeros((3, 3)) 

232 s = i + 2 

233 for j in range(s, s + 3): 

234 temp = lines[j].split() 

235 for k in range(3): 

236 lattice_vectors[j - s][k] = float(temp[k]) 

237 self.atoms.set_cell(lattice_vectors) 

238 if self.fractional_coordinates is not None: 

239 self.fractional_coordinates = np.array( 

240 self.fractional_coordinates) 

241 self.atoms.set_scaled_positions( 

242 self.fractional_coordinates) 

243 

244 elif line.find('Final fractional coordinates of atoms') != -1: 

245 s = i + 5 

246 scaled_positions = [] 

247 while True: 

248 s = s + 1 

249 if lines[s].find("------------") != -1: 

250 break 

251 if lines[s].find(" s ") != -1: 

252 continue 

253 xyz = lines[s].split()[3:6] 

254 XYZ = [float(x) for x in xyz] 

255 scaled_positions.append(XYZ) 

256 self.fractional_coordinates = scaled_positions 

257 

258 self.steps = cycles 

259 

260 def get_opt_state(self): 

261 return self.optimized 

262 

263 def get_opt_steps(self): 

264 return self.steps 

265 

266 def get_Gnorm(self): 

267 return self.Gnorm 

268 

269 def library_check(self): 

270 if self.parameters['library'] is not None: 

271 if 'GULP_LIB' not in os.environ: 

272 raise RuntimeError("Be sure to have set correctly $GULP_LIB " 

273 "or to have the force field library.") 

274 

275 

276class Conditions: 

277 """Atomic labels for the GULP calculator. 

278 

279 This class manages an array similar to 

280 atoms.get_chemical_symbols() via get_atoms_labels() method, but 

281 with atomic labels in stead of atomic symbols. This is useful 

282 when you need to use calculators like GULP or lammps that use 

283 force fields. Some force fields can have different atom type for 

284 the same element. In this class you can create a set_rule() 

285 function that assigns labels according to structural criteria.""" 

286 

287 def __init__(self, atoms): 

288 self.atoms = atoms 

289 self.atoms_symbols = atoms.get_chemical_symbols() 

290 self.atoms_labels = atoms.get_chemical_symbols() 

291 self.atom_types = [] 

292 

293 def min_distance_rule(self, sym1, sym2, 

294 ifcloselabel1=None, ifcloselabel2=None, 

295 elselabel1=None, max_distance=3.0): 

296 """Find pairs of atoms to label based on proximity. 

297 

298 This is for, e.g., the ffsioh or catlow force field, where we 

299 would like to identify those O atoms that are close to H 

300 atoms. For each H atoms, we must specially label one O atom. 

301 

302 This function is a rule that allows to define atom labels (like O1, 

303 O2, O_H etc..) starting from element symbols of an Atoms 

304 object that a force field can use and according to distance 

305 parameters. 

306 

307 Example: 

308 atoms = read('some_xyz_format.xyz') 

309 a = Conditions(atoms) 

310 a.set_min_distance_rule('O', 'H', ifcloselabel1='O2', 

311 ifcloselabel2='H', elselabel1='O1') 

312 new_atoms_labels = a.get_atom_labels() 

313 

314 In the example oxygens O are going to be labeled as O2 if they 

315 are close to a hydrogen atom othewise are labeled O1. 

316 

317 """ 

318 

319 if ifcloselabel1 is None: 

320 ifcloselabel1 = sym1 

321 if ifcloselabel2 is None: 

322 ifcloselabel2 = sym2 

323 if elselabel1 is None: 

324 elselabel1 = sym1 

325 

326 # self.atom_types is a list of element types used instead of element 

327 # symbols in orger to track the changes made. Take care of this because 

328 # is very important.. gulp_read function that parse the output 

329 # has to know which atom_type it has to associate with which 

330 # atom_symbol 

331 # 

332 # Example: [['O','O1','O2'],['H', 'H_C', 'H_O']] 

333 # this because Atoms oject accept only atoms symbols 

334 self.atom_types.append([sym1, ifcloselabel1, elselabel1]) 

335 self.atom_types.append([sym2, ifcloselabel2]) 

336 

337 dist_mat = self.atoms.get_all_distances() 

338 index_assigned_sym1 = [] 

339 index_assigned_sym2 = [] 

340 

341 for i in range(len(self.atoms_symbols)): 

342 if self.atoms_symbols[i] == sym2: 

343 dist_12 = 1000 

344 index_assigned_sym2.append(i) 

345 for t in range(len(self.atoms_symbols)): 

346 if (self.atoms_symbols[t] == sym1 

347 and dist_mat[i, t] < dist_12 

348 and t not in index_assigned_sym1): 

349 dist_12 = dist_mat[i, t] 

350 closest_sym1_index = t 

351 index_assigned_sym1.append(closest_sym1_index) 

352 

353 for i1, i2 in zip(index_assigned_sym1, index_assigned_sym2): 

354 if dist_mat[i1, i2] > max_distance: 

355 raise ValueError('Cannot unambiguously apply minimum-distance ' 

356 'rule because pairings are not obvious. ' 

357 'If you wish to ignore this, then increase ' 

358 'max_distance.') 

359 

360 for s in range(len(self.atoms_symbols)): 

361 if s in index_assigned_sym1: 

362 self.atoms_labels[s] = ifcloselabel1 

363 elif s not in index_assigned_sym1 and self.atoms_symbols[s] == sym1: 

364 self.atoms_labels[s] = elselabel1 

365 elif s in index_assigned_sym2: 

366 self.atoms_labels[s] = ifcloselabel2 

367 

368 def get_atom_types(self): 

369 return self.atom_types 

370 

371 def get_atoms_labels(self): 

372 labels = np.array(self.atoms_labels) 

373 return labels