1# -*- coding: utf-8 -*- 

2# 

3# Copyright (c) 2017, the cclib development team 

4# 

5# This file is part of cclib (http://cclib.github.io) and is distributed under 

6# the terms of the BSD 3-Clause License. 

7 

8"""A writer for MOLDEN format files.""" 

9 

10import os.path 

11import math 

12import decimal 

13import numpy 

14 

15from cclib.parser import utils 

16from cclib.io import filewriter 

17 

18 

19def round_molden(num, p=6): 

20 """Molden style number rounding in [Atoms] section.""" 

21 # Digit at pth position after dot. 

22 p_digit = math.floor(abs(num) * 10 ** p) % 10 

23 # If the 6th digit after dot is greater than 5, but is not 7, 

24 # round the number upto 6th place. 

25 # Else truncate at 6th digit after dot. 

26 if p_digit > 5 and p_digit != 7: 

27 return round(num, p) 

28 if num >= 0: 

29 return math.floor(num * 10 ** p) / 10 ** p 

30 else: 

31 return math.ceil(num * 10 ** p) / 10 ** p 

32 

33 

34class MOLDEN(filewriter.Writer): 

35 """A writer for MOLDEN files.""" 

36 

37 required_attrs = ('atomcoords', 'atomnos', 'natom') 

38 

39 def _title(self, path): 

40 """Return filename without extension to be used as title.""" 

41 title = os.path.basename(os.path.splitext(path)[0]) 

42 return title 

43 

44 def _coords_from_ccdata(self, index): 

45 """Create [Atoms] section using geometry at the given index.""" 

46 elements = [self.pt.element[Z] for Z in self.ccdata.atomnos] 

47 if self.ghost is not None: 

48 elements = [self.ghost if e is None else e for e in elements] 

49 elif None in elements: 

50 raise ValueError('It seems that there is at least one ghost atom ' + 

51 'in these elements. Please use the ghost flag to'+ 

52 ' specify a label for the ghost atoms.') 

53 atomcoords = self.ccdata.atomcoords[index] 

54 atomnos = self.ccdata.atomnos 

55 nos = range(self.ccdata.natom) 

56 

57 # element_name number atomic_number x y z 

58 atom_template = '{:2s} {:5d} {:2d} {:12.6f} {:12.6f} {:12.6f}' 

59 lines = [] 

60 for element, no, atomno, coord in zip(elements, nos, atomnos, 

61 atomcoords): 

62 x, y, z = map(round_molden, coord) 

63 lines.append(atom_template.format(element, no + 1, atomno, 

64 x, y, z)) 

65 

66 return lines 

67 

68 def _gto_from_ccdata(self): 

69 """Create [GTO] section using gbasis. 

70 

71 atom_sequence_number1 0 

72 shell_label number_of_primitives 1.00 

73 exponent_primitive_1 contraction_coeff_1 (contraction_coeff_1) 

74 ... 

75 empty line 

76 atom_sequence__number2 0 

77 """ 

78 

79 gbasis = self.ccdata.gbasis 

80 label_template = '{:s} {:5d} 1.00' 

81 basis_template = '{:15.9e} {:15.9e}' 

82 lines = [] 

83 

84 for no, basis in enumerate(gbasis): 

85 lines.append('{:3d} 0'.format(no + 1)) 

86 for prims in basis: 

87 lines.append(label_template.format(prims[0].lower(), 

88 len(prims[1]))) 

89 for prim in prims[1]: 

90 lines.append(basis_template.format(prim[0], prim[1])) 

91 lines.append('') 

92 lines.append('') 

93 return lines 

94 

95 def _scfconv_from_ccdata(self): 

96 """Create [SCFCONV] section using gbasis. 

97 

98 scf-first 1 THROUGH 12 

99 -672.634394 

100 ... 

101 -673.590571 

102 -673.590571 

103 """ 

104 

105 lines = ["scf-first 1 THROUGH %d" % len(self.ccdata.scfenergies)] 

106 

107 for scfenergy in self.ccdata.scfenergies: 

108 lines.append('{:15.6f}'.format(scfenergy)) 

109 

110 return lines 

111 

112 def _rearrange_mocoeffs(self, mocoeffs): 

113 """Rearrange cartesian F functions in mocoeffs. 

114 

115 Molden's order: 

116 xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz 

117 cclib's order: 

118 XXX, YYY, ZZZ, XXY, XXZ, YYX, YYZ, ZZX, ZZY, XYZ 

119 cclib's order can be converted by: 

120 moving YYX two indexes ahead, and 

121 moving YYZ two indexes back. 

122 """ 

123 aonames = self.ccdata.aonames 

124 mocoeffs = mocoeffs.tolist() 

125 

126 pos_yyx = [key for key, val in enumerate(aonames) if '_YYX' in val] 

127 pos_yyz = [key for key, val in enumerate(aonames) if '_YYZ' in val] 

128 

129 if pos_yyx: 

130 for pos in pos_yyx: 

131 mocoeffs.insert(pos-2, mocoeffs.pop(pos)) 

132 if pos_yyz: 

133 for pos in pos_yyz: 

134 mocoeffs.insert(pos+2, mocoeffs.pop(pos)) 

135 

136 return mocoeffs 

137 

138 

139 def _mo_from_ccdata(self): 

140 """Create [MO] section. 

141 

142 Sym= symmetry_label_1 

143 Ene= mo_energy_1 

144 Spin= (Alpha|Beta) 

145 Occup= mo_occupation_number_1 

146 ao_number_1 mo_coefficient_1 

147 ... 

148 ao_number_n mo_coefficient_n 

149 ... 

150 """ 

151 

152 moenergies = self.ccdata.moenergies 

153 mocoeffs = self.ccdata.mocoeffs 

154 homos = self.ccdata.homos 

155 mult = self.ccdata.mult 

156 

157 has_syms = False 

158 lines = [] 

159 

160 # Sym attribute is optional in [MO] section. 

161 if hasattr(self.ccdata, 'mosyms'): 

162 has_syms = True 

163 syms = self.ccdata.mosyms 

164 else: 

165 syms = numpy.full_like(moenergies, 'A', dtype=str) 

166 unres = len(moenergies) > 1 

167 openshell = len(homos) > 1 

168 

169 spin = 'Alpha' 

170 for i in range(len(moenergies)): 

171 for j in range(len(moenergies[i])): 

172 lines.append(' Sym= %s' % syms[i][j]) 

173 moenergy = utils.convertor(moenergies[i][j], 'eV', 'hartree') 

174 lines.append(' Ene= {:10.4f}'.format(moenergy)) 

175 lines.append(' Spin= %s' % spin) 

176 if unres and openshell: 

177 if j <= homos[i]: 

178 lines.append(' Occup= {:10.6f}'.format(1.0)) 

179 else: 

180 lines.append(' Occup= {:10.6f}'.format(0.0)) 

181 elif not unres and openshell: 

182 occ = numpy.sum(j <= homos) 

183 if j <= homos[i]: 

184 lines.append(' Occup= {:10.6f}'.format(occ)) 

185 else: 

186 lines.append(' Occup= {:10.6f}'.format(0.0)) 

187 else: 

188 if j <= homos[i]: 

189 lines.append(' Occup= {:10.6f}'.format(2.0)) 

190 else: 

191 lines.append(' Occup= {:10.6f}'.format(0.0)) 

192 # Rearrange mocoeffs according to Molden's lexicographical order. 

193 mocoeffs[i][j] = self._rearrange_mocoeffs(mocoeffs[i][j]) 

194 for k, mocoeff in enumerate(mocoeffs[i][j]): 

195 lines.append('{:4d} {:10.6f}'.format(k + 1, mocoeff)) 

196 

197 spin = 'Beta' 

198 

199 return lines 

200 

201 def generate_repr(self): 

202 """Generate the MOLDEN representation of the logfile data.""" 

203 

204 molden_lines = ['[Molden Format]'] 

205 

206 # Title of file. 

207 if self.jobfilename is not None: 

208 molden_lines.append('[Title]') 

209 molden_lines.append(self._title(self.jobfilename)) 

210 

211 # Coordinates for the Electron Density/Molecular orbitals. 

212 # [Atoms] (Angs|AU) 

213 unit = "Angs" 

214 molden_lines.append('[Atoms] %s' % unit) 

215 # Last set of coordinates for geometry optimization runs. 

216 index = -1 

217 molden_lines.extend(self._coords_from_ccdata(index)) 

218 

219 # Either both [GTO] and [MO] should be present or none of them. 

220 if hasattr(self.ccdata, 'gbasis') and hasattr(self.ccdata, 'mocoeffs')\ 

221 and hasattr(self.ccdata, 'moenergies'): 

222 

223 molden_lines.append('[GTO]') 

224 molden_lines.extend(self._gto_from_ccdata()) 

225 

226 molden_lines.append('[MO]') 

227 molden_lines.extend(self._mo_from_ccdata()) 

228 

229 # Omitting until issue #390 is resolved. 

230 # https://github.com/cclib/cclib/issues/390 

231 # if hasattr(self.ccdata, 'scfenergies'): 

232 # if len(self.ccdata.scfenergies) > 1: 

233 # molden_lines.append('[SCFCONV]') 

234 # molden_lines.extend(self._scfconv_from_ccdata()) 

235 

236 # molden_lines.append('') 

237 

238 return '\n'.join(molden_lines) 

239 

240 

241class MoldenReformatter: 

242 """Reformat Molden output files.""" 

243 

244 def __init__(self, filestring): 

245 self.filestring = filestring 

246 

247 def scinotation(self, num): 

248 """Convert Molden style number formatting to scientific notation. 

249 0.9910616900D+02 --> 9.910617e+01 

250 """ 

251 num = num.replace('D', 'e') 

252 return str('%.9e' % decimal.Decimal(num)) 

253 

254 def reformat(self): 

255 """Reformat Molden output file to: 

256 - use scientific notation, 

257 - split sp molecular orbitals to s and p, and 

258 - replace multiple spaces with single.""" 

259 filelines = iter(self.filestring.split("\n")) 

260 lines = [] 

261 

262 for line in filelines: 

263 line = line.replace('\n', '') 

264 # Replace multiple spaces with single spaces. 

265 line = ' '.join(line.split()) 

266 

267 # Check for [Title] section. 

268 if '[title]' in line.lower(): 

269 # skip the title 

270 line = next(filelines) 

271 line = next(filelines) 

272 

273 # Exclude SCFCONV section until issue #390 is resolved. 

274 # https://github.com/cclib/cclib/issues/390 

275 if '[scfconv]' in line.lower(): 

276 break 

277 

278 # Although Molden format specifies Sym in [MO] section, 

279 # the Molden program does not print it. 

280 if 'sym' in line.lower(): 

281 continue 

282 

283 # Convert D notation to scientific notation. 

284 if 'D' in line: 

285 vals = line.split() 

286 vals = [self.scinotation(i) for i in vals] 

287 lines.append(' '.join(vals)) 

288 

289 # Convert sp to s and p orbitals. 

290 elif 'sp' in line: 

291 n_prim = int(line.split()[1]) 

292 new_s = ['s ' + str(n_prim) + ' 1.00'] 

293 new_p = ['p ' + str(n_prim) + ' 1.00'] 

294 while n_prim > 0: 

295 n_prim -= 1 

296 line = next(filelines).split() 

297 new_s.append(self.scinotation(line[0]) + ' ' 

298 + self.scinotation(line[1])) 

299 new_p.append(self.scinotation(line[0]) + ' ' 

300 + self.scinotation(line[2])) 

301 lines.extend(new_s) 

302 lines.extend(new_p) 

303 else: 

304 lines.append(line) 

305 

306 return '\n'.join(lines)