Coverage for /builds/kinetik161/ase/ase/calculators/crystal.py: 13.60%

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

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

5Written by:

7 Daniele Selli, daniele.selli@unimib.it

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

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.

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

16 guess = True,

17 xc = 'PBE',

18 kpts = (2,2,2),

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

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

21 ...

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

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

27 1. Re-calculating the terms and subtracting them

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

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

32 only available in a development version of CRYSTAL

34"""

36import os

38import numpy as np

40from ase.calculators.calculator import FileIOCalculator

41from ase.io import write

42from ase.units import Bohr, Hartree

45class CRYSTAL(FileIOCalculator):

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

47 """

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

50 'dipole']

52 def __init__(self, restart=None,

53 ignore_bad_restart_file=FileIOCalculator._deprecated,

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

55 """Construct a crystal calculator.

57 """

58 # default parameters

59 self.default_parameters = dict(

60 xc='HF',

61 spinpol=False,

62 oldgrid=False,

63 neigh=False,

64 coarsegrid=False,

65 guess=True,

66 kpts=None,

67 isp=1,

68 basis='custom',

69 smearing=None,

70 otherkeys=[])

72 self.pcpot = None

73 self.lines = None

74 self.atoms = None

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

76 self.atoms_input = None

77 self.outfilename = 'cry.out'

79 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,

80 label, atoms,

81 **kwargs)

83 def write_crystal_in(self, filename):

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

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

86 """

88 # write BLOCK 1 (only SP with gradients)

89 with open(filename, 'w', encoding='latin-1') as outfile:

90 self._write_crystal_in(outfile)

92 def _write_crystal_in(self, outfile):

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

94 outfile.write('EXTERNAL \n')

95 outfile.write('NEIGHPRT \n')

96 outfile.write('0 \n')

98 if self.pcpot:

99 outfile.write('POINTCHG \n')

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

101

102 # write BLOCK 2 from file (basis sets)

103 p = self.parameters

104 if p.basis == 'custom':

105 outfile.write('END \n')

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

107 basis_ = basisfile.readlines()

108 for line in basis_:

109 outfile.write(line)

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

111 outfile.write('END \n')

112 else:

113 outfile.write('BASISSET \n')

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

115

116 # write BLOCK 3 according to parameters set as input

117 # ----- write hamiltonian

118

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

120 p.spinpol = True

121

122 if p.xc == 'HF':

123 if p.spinpol:

124 outfile.write('UHF \n')

125 else:

126 outfile.write('RHF \n')

127 elif p.xc == 'MP2':

128 outfile.write('MP2 \n')

129 outfile.write('ENDMP2 \n')

130 else:

131 outfile.write('DFT \n')

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

133 if isinstance(p.xc, str):

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

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

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

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

138 else:

139 x, c = p.xc

140 outfile.write('EXCHANGE \n')

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

142 outfile.write('CORRELAT \n')

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

144 if p.spinpol:

145 outfile.write('SPIN \n')

146 if p.oldgrid:

147 outfile.write('OLDGRID \n')

148 if p.coarsegrid:

149 outfile.write('RADIAL\n')

150 outfile.write('1\n')

151 outfile.write('4.0\n')

152 outfile.write('20\n')

153 outfile.write('ANGULAR\n')

154 outfile.write('5\n')

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

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

157 outfile.write('END \n')

158 # When guess=True, wf is read.

159 if p.guess:

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

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

162 outfile.write('GUESSP \n')

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

164 outfile.write('GUESSP \n')

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

166

167 # smearing

168 if p.smearing is not None:

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

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

171 else:

172 outfile.write('SMEAR \n')

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

174

175 # ----- write other CRYSTAL keywords

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

177

178 for keyword in p.otherkeys:

179 if isinstance(keyword, str):

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

181 else:

182 for key in keyword:

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

184

185 ispbc = self.atoms.get_pbc()

186 self.kpts = p.kpts

187

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

189 if any(ispbc):

190 if self.kpts is None:

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

192 else:

193 self.kpts = None

194

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

196 # Pack net and k-point density definition are

197 # not allowed.

198 if self.kpts is not None:

199 if isinstance(self.kpts, float):

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

201 if isinstance(self.kpts, list):

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

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

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

205 outfile.write('SHRINK \n')

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

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

208 if ispbc[2]:

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

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

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

212 elif ispbc[1]:

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

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

215 + ' 1 \n')

216 elif ispbc[0]:

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

218 + ' 1 1 \n')

219

220 # GRADCAL command performs a single

221 # point and prints out the forces

222 # also on the charges

223 outfile.write('GRADCAL \n')

224 outfile.write('END \n')

225

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

227 FileIOCalculator.write_input(

228 self, atoms, properties, system_changes)

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

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

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

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

233 self.atoms_input = atoms

234 self.atoms = None

235

236 def read_results(self):

237 """ all results are read from OUTPUT file

238 It will be destroyed after it is read to avoid

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

240

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

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

243 self.lines = myfile.readlines()

244

245 self.atoms = self.atoms_input

246 # Energy line index

247 estring1 = 'SCF ENDED'

248 estring2 = 'TOTAL ENERGY + DISP'

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

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

251 index_energy = iline

252 pos_en = 8

253 break

254 else:

255 raise RuntimeError('Problem in reading energy')

256 # Check if there is dispersion corrected

257 # energy value.

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

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

260 index_energy = iline

261 pos_en = 5

262

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

264 e_coul = 0

265 if self.pcpot:

266 if self.crys_pcc:

267 self.pcpot.read_pc_corrections()

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

269 self.pcpot.crys_pcc = True

270 else:

271 self.pcpot.manual_pc_correct()

272 e_coul, f_coul = self.pcpot.coulomb_corrections

273

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

275 energy -= e_coul # e_coul already in eV.

276

277 self.results['energy'] = energy

278 # Force line indexes

279 fstring = 'CARTESIAN FORCES'

280 gradients = []

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

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

283 index_force_begin = iline + 2

284 break

285 else:

286 raise RuntimeError('Problem in reading forces')

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

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

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

290 if len(word) == 5:

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

292 elif len(word) == 4:

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

294 else:

295 raise RuntimeError('Problem in reading forces')

296

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

298

299 self.results['forces'] = forces

300

301 # stress stuff begins

302 sstring = 'STRESS TENSOR, IN'

303 have_stress = False

304 stress = []

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

306 if sstring in line:

307 have_stress = True

308 start = iline + 4

309 end = start + 3

310 for i in range(start, end):

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

312 stress.append(cell)

313 if have_stress:

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

315 self.results['stress'] = stress

316

317 # stress stuff ends

318

319 # Get partial charges on atoms.

320 # In case we cannot find charges

321 # they are set to None

322 qm_charges = []

323

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

325 # ----- string search, which corresponds

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

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

328 if 'TOTAL ATOMIC CHARGE' in line:

329 chargestart = n + 1

330 lines1 = self.lines[chargestart:(chargestart

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

332 atomnum = self.atoms.get_atomic_numbers()

333 words = []

334 for line in lines1:

335 for el in line.split():

336 words.append(float(el))

337 i = 0

338 for atn in atomnum:

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

340 i = i + 1

341 charges = np.array(qm_charges)

342 self.results['charges'] = charges

343

344 # Read dipole moment.

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

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

347 if 'DIPOLE MOMENT ALONG' in line:

348 dipolestart = n + 2

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

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

351 break

352 # debye to e*Ang

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

354

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

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

357 """

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

359 return self.pcpot

360

361

362class PointChargePotential:

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

364 """Point-charge potential for CRYSTAL.

365 """

366 self.mmcharges = mmcharges

367 self.directory = directory

368 self.mmpositions = None

369 self.mmforces = None

370 self.coulomb_corrections = None

371 self.crys_pcc = False

372

373 def set_positions(self, mmpositions):

374 self.mmpositions = mmpositions

375

376 def set_charges(self, mmcharges):

377 self.mmcharges = mmcharges

378

379 def write_mmcharges(self, filename):

380 """ mok all

381 write external charges as monopoles for CRYSTAL.

382

383 """

384 if self.mmcharges is None:

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

386 return

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

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

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

390 [x, y, z] = pos

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

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

393

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

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

396 if get_forces:

397 return self.read_forces_on_pointcharges()

398 else:

399 return np.zeros_like(self.mmpositions)

400

401 def read_forces_on_pointcharges(self):

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

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

404 lines = infile.readlines()

405

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

407 # read in force and energy Coulomb corrections

408 if self.crys_pcc:

409 self.read_pc_corrections()

410 else:

411 self.manual_pc_correct()

412 e_coul, f_coul = self.coulomb_corrections

413

414 external_forces = []

415 for n, line in enumerate(lines):

416 if ('RESULTANT FORCE' in line):

417 chargeend = n - 1

418 break

419 else:

420 raise RuntimeError(

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

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

423 for line in lines1:

424 external_forces.append(

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

426

427 f = np.array(external_forces) - f_coul

428 f *= (Hartree / Bohr)

429

430 return f

431

432 def read_pc_corrections(self):

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

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

435 to be subtracted again.

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

437

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

439 'FORCES_CHG.DAT')) as infile:

440 lines = infile.readlines()

441

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

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

444

445 e *= Hartree

446

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

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

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

450

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

452 for i, l in enumerate(f_lines):

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

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

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

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

457

458 self.coulomb_corrections = (e, pc_forces)

459

460 def manual_pc_correct(self):

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

462

463 R = self.mmpositions / Bohr

464 charges = self.mmcharges

465

466 forces = np.zeros_like(R)

467 energy = 0.0

468

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

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

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

472 d = d2**0.5

473

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

475

476 energy += np.sum(e_c)

477

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

479

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

481 forces[m + 1:] += F

482

483 energy *= Hartree

484

485 self.coulomb_corrections = (energy, forces)