1"""This module defines an ASE interface to CRYSTAL14/CRYSTAL17 

2 

3http://www.crystal.unito.it/ 

4 

5Written by: 

6 

7 Daniele Selli, daniele.selli@unimib.it 

8 Gianluca Fazio, g.fazio3@campus.unimib.it 

9 

10The file 'fort.34' contains the input and output geometry 

11and it will be updated during the crystal calculations. 

12The wavefunction is stored in 'fort.20' as binary file. 

13 

14The keywords are given, for instance, as follows: 

15 

16 guess = True, 

17 xc = 'PBE', 

18 kpts = (2,2,2), 

19 otherkeys = [ 'scfdir', 'anderson', ['maxcycles','500'], 

20 ['fmixing','90']], 

21 ... 

22 

23 

24 When used for QM/MM, Crystal calculates coulomb terms 

25 within all point charges. This is wrong and should be corrected by either: 

26 

27 1. Re-calculating the terms and subtracting them 

28 2. Reading in the values from FORCES_CHG.DAT and subtracting 

29 

30 

31 BOTH Options should be available, with 1 as standard, since 2 is 

32 only available in a development version of CRYSTAL 

33 

34""" 

35 

36from ase.units import Hartree, Bohr 

37from ase.io import write 

38import numpy as np 

39import os 

40from ase.calculators.calculator import FileIOCalculator 

41 

42 

43class CRYSTAL(FileIOCalculator): 

44 """ A crystal calculator with ase-FileIOCalculator nomenclature 

45 """ 

46 

47 implemented_properties = ['energy', 'forces', 'stress', 'charges', 

48 'dipole'] 

49 

50 def __init__(self, restart=None, 

51 ignore_bad_restart_file=FileIOCalculator._deprecated, 

52 label='cry', atoms=None, crys_pcc=False, **kwargs): 

53 """Construct a crystal calculator. 

54 

55 """ 

56 # default parameters 

57 self.default_parameters = dict( 

58 xc='HF', 

59 spinpol=False, 

60 oldgrid=False, 

61 neigh=False, 

62 coarsegrid=False, 

63 guess=True, 

64 kpts=None, 

65 isp=1, 

66 basis='custom', 

67 smearing=None, 

68 otherkeys=[]) 

69 

70 self.pcpot = None 

71 self.lines = None 

72 self.atoms = None 

73 self.crys_pcc = crys_pcc # True: Reads Coulomb Correction from file. 

74 self.atoms_input = None 

75 self.outfilename = 'cry.out' 

76 

77 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

78 label, atoms, 

79 **kwargs) 

80 

81 def write_crystal_in(self, filename): 

82 """ Write the input file for the crystal calculation. 

83 Geometry is taken always from the file 'fort.34' 

84 """ 

85 

86 # write BLOCK 1 (only SP with gradients) 

87 with open(filename, 'wt', encoding='latin-1') as outfile: 

88 self._write_crystal_in(outfile) 

89 

90 def _write_crystal_in(self, outfile): 

91 outfile.write('Single point + Gradient crystal calculation \n') 

92 outfile.write('EXTERNAL \n') 

93 outfile.write('NEIGHPRT \n') 

94 outfile.write('0 \n') 

95 

96 if self.pcpot: 

97 outfile.write('POINTCHG \n') 

98 self.pcpot.write_mmcharges('POINTCHG.INP') 

99 

100 # write BLOCK 2 from file (basis sets) 

101 p = self.parameters 

102 if p.basis == 'custom': 

103 outfile.write('END \n') 

104 with open(os.path.join(self.directory, 'basis')) as basisfile: 

105 basis_ = basisfile.readlines() 

106 for line in basis_: 

107 outfile.write(line) 

108 outfile.write('99 0 \n') 

109 outfile.write('END \n') 

110 else: 

111 outfile.write('BASISSET \n') 

112 outfile.write(p.basis.upper() + '\n') 

113 

114 # write BLOCK 3 according to parameters set as input 

115 # ----- write hamiltonian 

116 

117 if self.atoms.get_initial_magnetic_moments().any(): 

118 p.spinpol = True 

119 

120 if p.xc == 'HF': 

121 if p.spinpol: 

122 outfile.write('UHF \n') 

123 else: 

124 outfile.write('RHF \n') 

125 elif p.xc == 'MP2': 

126 outfile.write('MP2 \n') 

127 outfile.write('ENDMP2 \n') 

128 else: 

129 outfile.write('DFT \n') 

130 # Standalone keywords and LDA are given by a single string. 

131 if isinstance(p.xc, str): 

132 xc = {'LDA': 'EXCHANGE\nLDA\nCORRELAT\nVWN', 

133 'PBE': 'PBEXC'}.get(p.xc, p.xc) 

134 outfile.write(xc.upper() + '\n') 

135 # Custom xc functional are given by a tuple of string 

136 else: 

137 x, c = p.xc 

138 outfile.write('EXCHANGE \n') 

139 outfile.write(x + ' \n') 

140 outfile.write('CORRELAT \n') 

141 outfile.write(c + ' \n') 

142 if p.spinpol: 

143 outfile.write('SPIN \n') 

144 if p.oldgrid: 

145 outfile.write('OLDGRID \n') 

146 if p.coarsegrid: 

147 outfile.write('RADIAL\n') 

148 outfile.write('1\n') 

149 outfile.write('4.0\n') 

150 outfile.write('20\n') 

151 outfile.write('ANGULAR\n') 

152 outfile.write('5\n') 

153 outfile.write('0.1667 0.5 0.9 3.05 9999.0\n') 

154 outfile.write('2 6 8 13 8\n') 

155 outfile.write('END \n') 

156 # When guess=True, wf is read. 

157 if p.guess: 

158 # wf will be always there after 2nd step. 

159 if os.path.isfile('fort.20'): 

160 outfile.write('GUESSP \n') 

161 elif os.path.isfile('fort.9'): 

162 outfile.write('GUESSP \n') 

163 os.system('cp fort.9 fort.20') 

164 

165 # smearing 

166 if p.smearing is not None: 

167 if p.smearing[0] != 'Fermi-Dirac': 

168 raise ValueError('Only Fermi-Dirac smearing is allowed.') 

169 else: 

170 outfile.write('SMEAR \n') 

171 outfile.write(str(p.smearing[1] / Hartree) + ' \n') 

172 

173 # ----- write other CRYSTAL keywords 

174 # ----- in the list otherkey = ['ANDERSON', ...] . 

175 

176 for keyword in p.otherkeys: 

177 if isinstance(keyword, str): 

178 outfile.write(keyword.upper() + '\n') 

179 else: 

180 for key in keyword: 

181 outfile.write(key.upper() + '\n') 

182 

183 ispbc = self.atoms.get_pbc() 

184 self.kpts = p.kpts 

185 

186 # if it is periodic, gamma is the default. 

187 if any(ispbc): 

188 if self.kpts is None: 

189 self.kpts = (1, 1, 1) 

190 else: 

191 self.kpts = None 

192 

193 # explicit lists of K-points, shifted Monkhorst- 

194 # Pack net and k-point density definition are 

195 # not allowed. 

196 if self.kpts is not None: 

197 if isinstance(self.kpts, float): 

198 raise ValueError('K-point density definition not allowed.') 

199 if isinstance(self.kpts, list): 

200 raise ValueError('Explicit K-points definition not allowed.') 

201 if isinstance(self.kpts[-1], str): 

202 raise ValueError('Shifted Monkhorst-Pack not allowed.') 

203 outfile.write('SHRINK \n') 

204 # isp is by default 1, 2 is suggested for metals. 

205 outfile.write('0 ' + str(p.isp * max(self.kpts)) + ' \n') 

206 if ispbc[2]: 

207 outfile.write(str(self.kpts[0]) 

208 + ' ' + str(self.kpts[1]) 

209 + ' ' + str(self.kpts[2]) + ' \n') 

210 elif ispbc[1]: 

211 outfile.write(str(self.kpts[0]) 

212 + ' ' + str(self.kpts[1]) 

213 + ' 1 \n') 

214 elif ispbc[0]: 

215 outfile.write(str(self.kpts[0]) 

216 + ' 1 1 \n') 

217 

218 # GRADCAL command performs a single 

219 # point and prints out the forces 

220 # also on the charges 

221 outfile.write('GRADCAL \n') 

222 outfile.write('END \n') 

223 

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

225 FileIOCalculator.write_input( 

226 self, atoms, properties, system_changes) 

227 self.write_crystal_in(os.path.join(self.directory, 'INPUT')) 

228 write(os.path.join(self.directory, 'fort.34'), atoms) 

229 # self.atoms is none until results are read out, 

230 # then it is set to the ones at writing input 

231 self.atoms_input = atoms 

232 self.atoms = None 

233 

234 def read_results(self): 

235 """ all results are read from OUTPUT file 

236 It will be destroyed after it is read to avoid 

237 reading it once again after some runtime error """ 

238 

239 with open(os.path.join(self.directory, 'OUTPUT'), 

240 'rt', 

241 encoding='latin-1') as myfile: 

242 self.lines = myfile.readlines() 

243 

244 self.atoms = self.atoms_input 

245 # Energy line index 

246 estring1 = 'SCF ENDED' 

247 estring2 = 'TOTAL ENERGY + DISP' 

248 for iline, line in enumerate(self.lines): 

249 if line.find(estring1) >= 0: 

250 index_energy = iline 

251 pos_en = 8 

252 break 

253 else: 

254 raise RuntimeError('Problem in reading energy') 

255 # Check if there is dispersion corrected 

256 # energy value. 

257 for iline, line in enumerate(self.lines): 

258 if line.find(estring2) >= 0: 

259 index_energy = iline 

260 pos_en = 5 

261 

262 # If there's a point charge potential (QM/MM), read corrections 

263 e_coul = 0 

264 if self.pcpot: 

265 if self.crys_pcc: 

266 self.pcpot.read_pc_corrections() 

267 # also pass on to pcpot that it should read in from file 

268 self.pcpot.crys_pcc = True 

269 else: 

270 self.pcpot.manual_pc_correct() 

271 e_coul, f_coul = self.pcpot.coulomb_corrections 

272 

273 energy = float(self.lines[index_energy].split()[pos_en]) * Hartree 

274 energy -= e_coul # e_coul already in eV. 

275 

276 self.results['energy'] = energy 

277 # Force line indexes 

278 fstring = 'CARTESIAN FORCES' 

279 gradients = [] 

280 for iline, line in enumerate(self.lines): 

281 if line.find(fstring) >= 0: 

282 index_force_begin = iline + 2 

283 break 

284 else: 

285 raise RuntimeError('Problem in reading forces') 

286 for j in range(index_force_begin, index_force_begin + len(self.atoms)): 

287 word = self.lines[j].split() 

288 # If GHOST atoms give problems, have a close look at this 

289 if len(word) == 5: 

290 gradients.append([float(word[k + 2]) for k in range(0, 3)]) 

291 elif len(word) == 4: 

292 gradients.append([float(word[k + 1]) for k in range(0, 3)]) 

293 else: 

294 raise RuntimeError('Problem in reading forces') 

295 

296 forces = np.array(gradients) * Hartree / Bohr 

297 

298 self.results['forces'] = forces 

299 

300 # stress stuff begins 

301 sstring = 'STRESS TENSOR, IN' 

302 have_stress = False 

303 stress = [] 

304 for iline, line in enumerate(self.lines): 

305 if sstring in line: 

306 have_stress = True 

307 start = iline + 4 

308 end = start + 3 

309 for i in range(start, end): 

310 cell = [float(x) for x in self.lines[i].split()] 

311 stress.append(cell) 

312 if have_stress: 

313 stress = -np.array(stress) * Hartree / Bohr**3 

314 self.results['stress'] = stress 

315 

316 # stress stuff ends 

317 

318 # Get partial charges on atoms. 

319 # In case we cannot find charges 

320 # they are set to None 

321 qm_charges = [] 

322 

323 # ----- this for cycle finds the last entry of the 

324 # ----- string search, which corresponds 

325 # ----- to the charges at the end of the SCF. 

326 for n, line in enumerate(self.lines): 

327 if 'TOTAL ATOMIC CHARGE' in line: 

328 chargestart = n + 1 

329 lines1 = self.lines[chargestart:(chargestart 

330 + (len(self.atoms) - 1) // 6 + 1)] 

331 atomnum = self.atoms.get_atomic_numbers() 

332 words = [] 

333 for line in lines1: 

334 for el in line.split(): 

335 words.append(float(el)) 

336 i = 0 

337 for atn in atomnum: 

338 qm_charges.append(-words[i] + atn) 

339 i = i + 1 

340 charges = np.array(qm_charges) 

341 self.results['charges'] = charges 

342 

343 # Read dipole moment. 

344 dipole = np.zeros([1, 3]) 

345 for n, line in enumerate(self.lines): 

346 if 'DIPOLE MOMENT ALONG' in line: 

347 dipolestart = n + 2 

348 dipole = np.array([float(f) for f in 

349 self.lines[dipolestart].split()[2:5]]) 

350 break 

351 # debye to e*Ang 

352 self.results['dipole'] = dipole * 0.2081943482534 

353 

354 def embed(self, mmcharges=None, directory='./'): 

355 """Embed atoms in point-charges (mmcharges) 

356 """ 

357 self.pcpot = PointChargePotential(mmcharges, self.directory) 

358 return self.pcpot 

359 

360 

361class PointChargePotential: 

362 def __init__(self, mmcharges, directory='./'): 

363 """Point-charge potential for CRYSTAL. 

364 """ 

365 self.mmcharges = mmcharges 

366 self.directory = directory 

367 self.mmpositions = None 

368 self.mmforces = None 

369 self.coulomb_corrections = None 

370 self.crys_pcc = False 

371 

372 def set_positions(self, mmpositions): 

373 self.mmpositions = mmpositions 

374 

375 def set_charges(self, mmcharges): 

376 self.mmcharges = mmcharges 

377 

378 def write_mmcharges(self, filename): 

379 """ mok all 

380 write external charges as monopoles for CRYSTAL. 

381 

382 """ 

383 if self.mmcharges is None: 

384 print("CRYSTAL: Warning: not writing external charges ") 

385 return 

386 with open(os.path.join(self.directory, filename), 'w') as charge_file: 

387 charge_file.write(str(len(self.mmcharges)) + ' \n') 

388 for [pos, charge] in zip(self.mmpositions, self.mmcharges): 

389 [x, y, z] = pos 

390 charge_file.write('%12.6f %12.6f %12.6f %12.6f \n' 

391 % (x, y, z, charge)) 

392 

393 def get_forces(self, calc, get_forces=True): 

394 """ returns forces on point charges if the flag get_forces=True """ 

395 if get_forces: 

396 return self.read_forces_on_pointcharges() 

397 else: 

398 return np.zeros_like(self.mmpositions) 

399 

400 def read_forces_on_pointcharges(self): 

401 """Read Forces from CRYSTAL output file (OUTPUT).""" 

402 with open(os.path.join(self.directory, 'OUTPUT'), 'r') as infile: 

403 lines = infile.readlines() 

404 

405 print('PCPOT crys_pcc: ' + str(self.crys_pcc)) 

406 # read in force and energy Coulomb corrections 

407 if self.crys_pcc: 

408 self.read_pc_corrections() 

409 else: 

410 self.manual_pc_correct() 

411 e_coul, f_coul = self.coulomb_corrections 

412 

413 external_forces = [] 

414 for n, line in enumerate(lines): 

415 if ('RESULTANT FORCE' in line): 

416 chargeend = n - 1 

417 break 

418 else: 

419 raise RuntimeError( 

420 'Problem in reading forces on MM external-charges') 

421 lines1 = lines[(chargeend - len(self.mmcharges)):chargeend] 

422 for line in lines1: 

423 external_forces.append( 

424 [float(i) for i in line.split()[2:]]) 

425 

426 f = np.array(external_forces) - f_coul 

427 f *= (Hartree / Bohr) 

428 

429 return f 

430 

431 def read_pc_corrections(self): 

432 ''' Crystal calculates Coulomb forces and energies between all 

433 point charges, and adds that to the QM subsystem. That needs 

434 to be subtracted again. 

435 This will be standard in future CRYSTAL versions .''' 

436 

437 with open(os.path.join(self.directory, 

438 'FORCES_CHG.DAT'), 'r') as infile: 

439 lines = infile.readlines() 

440 

441 e = [float(x.split()[-1]) 

442 for x in lines if 'SELF-INTERACTION ENERGY(AU)' in x][0] 

443 

444 e *= Hartree 

445 

446 f_lines = [s for s in lines if '199' in s] 

447 assert len(f_lines) == len(self.mmcharges), \ 

448 'Mismatch in number of point charges from FORCES_CHG.dat' 

449 

450 pc_forces = np.zeros((len(self.mmcharges), 3)) 

451 for i, l in enumerate(f_lines): 

452 first = l.split(str(i + 1) + ' 199 ') 

453 assert len(first) == 2, 'Problem reading FORCES_CHG.dat' 

454 f = first[-1].split() 

455 pc_forces[i] = [float(x) for x in f] 

456 

457 self.coulomb_corrections = (e, pc_forces) 

458 

459 def manual_pc_correct(self): 

460 ''' For current versions of CRYSTAL14/17, manual Coulomb correction ''' 

461 

462 R = self.mmpositions / Bohr 

463 charges = self.mmcharges 

464 

465 forces = np.zeros_like(R) 

466 energy = 0.0 

467 

468 for m in range(len(charges)): 

469 D = R[m + 1:] - R[m] 

470 d2 = (D**2).sum(1) 

471 d = d2**0.5 

472 

473 e_c = charges[m + 1:] * charges[m] / d 

474 

475 energy += np.sum(e_c) 

476 

477 F = (e_c / d2)[:, None] * D 

478 

479 forces[m] -= F.sum(0) 

480 forces[m + 1:] += F 

481 

482 energy *= Hartree 

483 

484 self.coulomb_corrections = (energy, forces)