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

1# fmt: off

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

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

7Written by:

9 Daniele Selli, daniele.selli@unimib.it

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

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

13and it will be updated during the crystal calculations.

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

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

18 guess = True,

19 xc = 'PBE',

20 kpts = (2,2,2),

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

22 ['fmixing','90']],

23 ...

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

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

29 1. Re-calculating the terms and subtracting them

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

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

34 only available in a development version of CRYSTAL

36"""

38import os

40import numpy as np

42from ase.calculators.calculator import FileIOCalculator

43from ase.io import write

44from ase.units import Bohr, Hartree

47class CRYSTAL(FileIOCalculator):

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

49 """

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

52 'dipole']

54 def __init__(self, restart=None,

55 ignore_bad_restart_file=FileIOCalculator._deprecated,

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

57 """Construct a crystal calculator.

59 """

60 # default parameters

61 self.default_parameters = dict(

62 xc='HF',

63 spinpol=False,

64 oldgrid=False,

65 neigh=False,

66 coarsegrid=False,

67 guess=True,

68 kpts=None,

69 isp=1,

70 basis='custom',

71 smearing=None,

72 otherkeys=[])

74 self.pcpot = None

75 self.lines = None

76 self.atoms = None

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

78 self.atoms_input = None

79 self.outfilename = 'cry.out'

81 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,

82 label, atoms,

83 **kwargs)

85 def write_crystal_in(self, filename):

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

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

88 """

90 # write BLOCK 1 (only SP with gradients)

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

92 self._write_crystal_in(outfile)

94 def _write_crystal_in(self, outfile):

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

96 outfile.write('EXTERNAL \n')

97 outfile.write('NEIGHPRT \n')

98 outfile.write('0 \n')

100 if self.pcpot:

101 outfile.write('POINTCHG \n')

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

103

104 # write BLOCK 2 from file (basis sets)

105 p = self.parameters

106 if p.basis == 'custom':

107 outfile.write('END \n')

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

109 basis_ = basisfile.readlines()

110 for line in basis_:

111 outfile.write(line)

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

113 outfile.write('END \n')

114 else:

115 outfile.write('BASISSET \n')

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

117

118 # write BLOCK 3 according to parameters set as input

119 # ----- write hamiltonian

120

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

122 p.spinpol = True

123

124 if p.xc == 'HF':

125 if p.spinpol:

126 outfile.write('UHF \n')

127 else:

128 outfile.write('RHF \n')

129 elif p.xc == 'MP2':

130 outfile.write('MP2 \n')

131 outfile.write('ENDMP2 \n')

132 else:

133 outfile.write('DFT \n')

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

135 if isinstance(p.xc, str):

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

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

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

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

140 else:

141 x, c = p.xc

142 outfile.write('EXCHANGE \n')

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

144 outfile.write('CORRELAT \n')

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

146 if p.spinpol:

147 outfile.write('SPIN \n')

148 if p.oldgrid:

149 outfile.write('OLDGRID \n')

150 if p.coarsegrid:

151 outfile.write('RADIAL\n')

152 outfile.write('1\n')

153 outfile.write('4.0\n')

154 outfile.write('20\n')

155 outfile.write('ANGULAR\n')

156 outfile.write('5\n')

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

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

159 outfile.write('END \n')

160 # When guess=True, wf is read.

161 if p.guess:

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

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

164 outfile.write('GUESSP \n')

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

166 outfile.write('GUESSP \n')

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

168

169 # smearing

170 if p.smearing is not None:

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

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

173 else:

174 outfile.write('SMEAR \n')

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

176

177 # ----- write other CRYSTAL keywords

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

179

180 for keyword in p.otherkeys:

181 if isinstance(keyword, str):

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

183 else:

184 for key in keyword:

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

186

187 ispbc = self.atoms.get_pbc()

188 self.kpts = p.kpts

189

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

191 if any(ispbc):

192 if self.kpts is None:

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

194 else:

195 self.kpts = None

196

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

198 # Pack net and k-point density definition are

199 # not allowed.

200 if self.kpts is not None:

201 if isinstance(self.kpts, float):

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

203 if isinstance(self.kpts, list):

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

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

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

207 outfile.write('SHRINK \n')

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

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

210 if ispbc[2]:

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

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

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

214 elif ispbc[1]:

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

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

217 + ' 1 \n')

218 elif ispbc[0]:

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

220 + ' 1 1 \n')

221

222 # GRADCAL command performs a single

223 # point and prints out the forces

224 # also on the charges

225 outfile.write('GRADCAL \n')

226 outfile.write('END \n')

227

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

229 FileIOCalculator.write_input(

230 self, atoms, properties, system_changes)

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

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

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

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

235 self.atoms_input = atoms

236 self.atoms = None

237

238 def read_results(self):

239 """ all results are read from OUTPUT file

240 It will be destroyed after it is read to avoid

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

242

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

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

245 self.lines = myfile.readlines()

246

247 self.atoms = self.atoms_input

248 # Energy line index

249 estring1 = 'SCF ENDED'

250 estring2 = 'TOTAL ENERGY + DISP'

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

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

253 index_energy = iline

254 pos_en = 8

255 break

256 else:

257 raise RuntimeError('Problem in reading energy')

258 # Check if there is dispersion corrected

259 # energy value.

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

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

262 index_energy = iline

263 pos_en = 5

264

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

266 e_coul = 0

267 if self.pcpot:

268 if self.crys_pcc:

269 self.pcpot.read_pc_corrections()

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

271 self.pcpot.crys_pcc = True

272 else:

273 self.pcpot.manual_pc_correct()

274 e_coul, _f_coul = self.pcpot.coulomb_corrections

275

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

277 energy -= e_coul # e_coul already in eV.

278

279 self.results['energy'] = energy

280 # Force line indexes

281 fstring = 'CARTESIAN FORCES'

282 gradients = []

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

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

285 index_force_begin = iline + 2

286 break

287 else:

288 raise RuntimeError('Problem in reading forces')

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

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

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

292 if len(word) == 5:

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

294 elif len(word) == 4:

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

296 else:

297 raise RuntimeError('Problem in reading forces')

298

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

300

301 self.results['forces'] = forces

302

303 # stress stuff begins

304 sstring = 'STRESS TENSOR, IN'

305 have_stress = False

306 stress = []

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

308 if sstring in line:

309 have_stress = True

310 start = iline + 4

311 end = start + 3

312 for i in range(start, end):

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

314 stress.append(cell)

315 if have_stress:

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

317 self.results['stress'] = stress

318

319 # stress stuff ends

320

321 # Get partial charges on atoms.

322 # In case we cannot find charges

323 # they are set to None

324 qm_charges = []

325

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

327 # ----- string search, which corresponds

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

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

330 if 'TOTAL ATOMIC CHARGE' in line:

331 chargestart = n + 1

332 lines1 = self.lines[chargestart:(chargestart

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

334 atomnum = self.atoms.get_atomic_numbers()

335 words = []

336 for line in lines1:

337 for el in line.split():

338 words.append(float(el))

339 i = 0

340 for atn in atomnum:

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

342 i = i + 1

343 charges = np.array(qm_charges)

344 self.results['charges'] = charges

345

346 # Read dipole moment.

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

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

349 if 'DIPOLE MOMENT ALONG' in line:

350 dipolestart = n + 2

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

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

353 break

354 # debye to e*Ang

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

356

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

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

359 """

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

361 return self.pcpot

362

363

364class PointChargePotential:

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

366 """Point-charge potential for CRYSTAL.

367 """

368 self.mmcharges = mmcharges

369 self.directory = directory

370 self.mmpositions = None

371 self.mmforces = None

372 self.coulomb_corrections = None

373 self.crys_pcc = False

374

375 def set_positions(self, mmpositions):

376 self.mmpositions = mmpositions

377

378 def set_charges(self, mmcharges):

379 self.mmcharges = mmcharges

380

381 def write_mmcharges(self, filename):

382 """ mok all

383 write external charges as monopoles for CRYSTAL.

384

385 """

386 if self.mmcharges is None:

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

388 return

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

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

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

392 [x, y, z] = pos

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

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

395

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

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

398 if get_forces:

399 return self.read_forces_on_pointcharges()

400 else:

401 return np.zeros_like(self.mmpositions)

402

403 def read_forces_on_pointcharges(self):

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

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

406 lines = infile.readlines()

407

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

409 # read in force and energy Coulomb corrections

410 if self.crys_pcc:

411 self.read_pc_corrections()

412 else:

413 self.manual_pc_correct()

414 _e_coul, f_coul = self.coulomb_corrections

415

416 external_forces = []

417 for n, line in enumerate(lines):

418 if ('RESULTANT FORCE' in line):

419 chargeend = n - 1

420 break

421 else:

422 raise RuntimeError(

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

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

425 for line in lines1:

426 external_forces.append(

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

428

429 f = np.array(external_forces) - f_coul

430 f *= (Hartree / Bohr)

431

432 return f

433

434 def read_pc_corrections(self):

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

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

437 to be subtracted again.

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

439

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

441 'FORCES_CHG.DAT')) as infile:

442 lines = infile.readlines()

443

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

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

446

447 e *= Hartree

448

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

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

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

452

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

454 for i, l in enumerate(f_lines):

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

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

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

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

459

460 self.coulomb_corrections = (e, pc_forces)

461

462 def manual_pc_correct(self):

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

464

465 R = self.mmpositions / Bohr

466 charges = self.mmcharges

467

468 forces = np.zeros_like(R)

469 energy = 0.0

470

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

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

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

474 d = d2**0.5

475

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

477

478 energy += np.sum(e_c)

479

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

481

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

483 forces[m + 1:] += F

484

485 energy *= Hartree

486

487 self.coulomb_corrections = (energy, forces)