Coverage for /builds/kinetik161/ase/ase/calculators/dmol.py: 28.08%

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

3Contacts

4--------

5Adam Arvidsson <adam.arvidsson@chalmers.se>

6Erik Fransson <erikfr@chalmers.se>

7Anders Hellman <anders.hellman@chalmers.se>

10DMol3 environment variables

11----------------------------

12DMOL_COMMAND should point to the RunDmol script and specify the number of cores

13to prallelize over

15export DMOL_COMMAND="./RunDmol.sh -np 16"

18Example

19--------

20>>> from ase.build import bulk

21>>> from ase.calculators import DMol3

23>>> atoms = bulk('Al','fcc')

24>>> calc = DMol3()

25>>> atoms.calc = calc

26>>> print 'Potential energy %5.5f eV' % atoms.get_potential_energy()

29DMol3 calculator functionality

30-------------------------------

31This calculator does support all the functionality in DMol3.

33Firstly this calculator is limited to only handling either fully

34periodic structures (pbc = [1,1,1]) or non periodic structures (pbc=[0,0,0]).

36Internal relaxations are not supported by the calculator,

37only support for energy and forces is implemented.

39Reading eigenvalues and kpts are supported.

40Be careful with kpts and their directions (see internal coordinates below).

42Outputting the full electron density or specific bands to .grd files can be

43achieved with the plot command. The .grd files can be converted to the cube

44format using grd_to_cube().

47DMol3 internal coordinates

48---------------------------

49DMol3 may change the atomic positions / cell vectors in order to satisfy

50certain criterion ( e.g. molecule symmetry axis along z ). Specifically this

51happens when using Symmetry on/auto. This means the forces read from .grad

52will be in a different coordinates system compared to the atoms object used.

53To solve this the rotation matrix that converts the dmol coordinate system

54to the ase coordinate system is found and applied to the forces.

56For non periodic structures (pbc=[0,0,0]) the rotation matrix can be directly

57parsed from the .rot file.

58For fully periodic structures the rotation matrix is found by reading the

59cell vectors and positions used by dmol and then solving the matrix problem

60DMol_atoms * rot_mat = ase_atoms

63DMol3 files

64------------

65The supported DMol3 file formats are:

67car structure file - Angstrom and cellpar description of cell.

68incoor structure file - Bohr and cellvector describption of cell.

69 Note: incoor file not used if car file present.

70outmol outfile from DMol3 - atomic units (Bohr and Hartree)

71grad outfile for forces from DMol3 - forces in Hartree/Bohr

72grd outfile for orbitals from DMol3 - cellpar in Angstrom

74"""

76import os

77import re

79import numpy as np

81from ase import Atoms

82from ase.calculators.calculator import FileIOCalculator, Parameters, ReadError

83from ase.io import read

84from ase.io.dmol import write_dmol_car, write_dmol_incoor

85from ase.units import Bohr, Hartree

88class DMol3(FileIOCalculator):

89 """ DMol3 calculator object. """

91 implemented_properties = ['energy', 'forces']

92 default_parameters = {'functional': 'pbe',

93 'symmetry': 'on'}

94 discard_results_on_any_change = True

96 def __init__(self, restart=None,

97 ignore_bad_restart_file=FileIOCalculator._deprecated,

98 label='dmol_calc/tmp', atoms=None,

99 command=None, **kwargs):

100 """ Construct DMol3 calculator. """

101

102 if command is None:

103 if 'DMOL_COMMAND' in self.cfg:

104 command = self.cfg['DMOL_COMMAND'] + ' PREFIX > PREFIX.out'

105

106 super().__init__(restart, ignore_bad_restart_file,

107 label, atoms, command=command,

108 **kwargs)

109

110 # tracks if DMol transformed coordinate system

111 self.internal_transformation = False

112

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

114

115 if not (np.all(atoms.pbc) or not np.any(atoms.pbc)):

116 raise RuntimeError('PBC must be all true or all false')

117

118 self.clean() # Remove files from old run

119 self.internal_transformation = False

120 self.ase_positions = atoms.positions.copy()

121 self.ase_cell = atoms.cell.copy()

122

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

124

125 if np.all(atoms.pbc):

126 write_dmol_incoor(self.label + '.incoor', atoms)

127 elif not np.any(atoms.pbc):

128 write_dmol_car(self.label + '.car', atoms)

129

130 self.write_input_file()

131 self.parameters.write(self.label + '.parameters.ase')

132

133 def write_input_file(self):

134 """ Writes the input file. """

135 with open(self.label + '.input', 'w') as fd:

136 self._write_input_file(fd)

137

138 def _write_input_file(self, fd):

139 fd.write('%-32s %s\n' % ('calculate', 'gradient'))

140

141 # if no key about eigs

142 fd.write('%-32s %s\n' % ('print', 'eigval_last_it'))

143

144 for key, value in self.parameters.items():

145 if isinstance(value, str):

146 fd.write('%-32s %s\n' % (key, value))

147 elif isinstance(value, (list, tuple)):

148 for val in value:

149 fd.write('%-32s %s\n' % (key, val))

150 else:

151 fd.write('%-32s %r\n' % (key, value))

152

153 def read(self, label):

154 FileIOCalculator.read(self, label)

155 geometry = self.label + '.car'

156 output = self.label + '.outmol'

157 force = self.label + '.grad'

158

159 for filename in [force, output, geometry]:

160 if not os.path.isfile(filename):

161 raise ReadError

162

163 self.atoms = read(geometry)

164 self.parameters = Parameters.read(self.label + 'parameters.ase')

165 self.read_results()

166

167 def read_results(self):

168 finished, message = self.finished_successfully()

169 if not finished:

170 raise RuntimeError('DMol3 run failed, see outmol file for'

171 ' more info\n\n%s' % message)

172

173 self.find_dmol_transformation()

174 self.read_energy()

175 self.read_forces()

176

177 def finished_successfully(self):

178 """ Reads outmol file and checks if job completed or failed.

179

180 Returns

181 -------

182 finished (bool): True if job completed, False if something went wrong

183 message (str): If job failed message contains parsed errors, else empty

184

185 """

186 finished = False

187 message = ""

188 for line in self._outmol_lines():

189 if line.rfind('Message: DMol3 job finished successfully') > -1:

190 finished = True

191 if line.startswith('Error'):

192 message += line

193 return finished, message

194

195 def find_dmol_transformation(self, tol=1e-4):

196 """Finds rotation matrix that takes us from DMol internal

197 coordinates to ase coordinates.

198

199 For pbc = [False, False, False] the rotation matrix is parsed from

200 the .rot file, if this file doesnt exist no rotation is needed.

201

202 For pbc = [True, True, True] the Dmol internal cell vectors and

203 positions are parsed and compared to self.ase_cell self.ase_positions.

204 The rotation matrix can then be found by a call to the helper

205 function find_transformation(atoms1, atoms2)

206

207 If a rotation matrix is needed then self.internal_transformation is

208 set to True and the rotation matrix is stored in self.rotation_matrix

209

210 Parameters

211 ----------

212 tol (float): tolerance for check if positions and cell are the same

213 """

214

215 if np.all(self.atoms.pbc): # [True, True, True]

216 dmol_atoms = self.read_atoms_from_outmol()

217 if (np.linalg.norm(self.atoms.positions - dmol_atoms.positions) <

218 tol) and (np.linalg.norm(self.atoms.cell -

219 dmol_atoms.cell) < tol):

220 self.internal_transformation = False

221 else:

222 R, err = find_transformation(dmol_atoms, self.atoms)

223 if abs(np.linalg.det(R) - 1.0) > tol:

224 raise RuntimeError('Error: transformation matrix does'

225 ' not have determinant 1.0')

226 if err < tol:

227 self.internal_transformation = True

228 self.rotation_matrix = R

229 else:

230 raise RuntimeError('Error: Could not find dmol'

231 ' coordinate transformation')

232 elif not np.any(self.atoms.pbc): # [False,False,False]

233 try:

234 data = np.loadtxt(self.label + '.rot')

235 except OSError:

236 self.internal_transformation = False

237 else:

238 self.internal_transformation = True

239 self.rotation_matrix = data[1:].transpose()

240

241 def read_atoms_from_outmol(self):

242 """ Reads atomic positions and cell from outmol file and returns atoms

243 object.

244

245 If no cell vectors are found in outmol the cell is set to np.eye(3) and

246 pbc 000.

247

248 Formatting for cell in outmol :

249 translation vector [a0] 1 5.1 0.0 5.1

250 translation vector [a0] 2 5.1 5.1 0.0

251 translation vector [a0] 3 0.0 5.1 5.1

252

253 Formatting for positions in outmol:

254 df ATOMIC COORDINATES (au)

255 df x y z

256 df Si 0.0 0.0 0.0

257 df Si 1.3 3.5 2.2

258 df binding energy -0.2309046Ha

259

260 Returns

261 -------

262 atoms (Atoms object): read atoms object

263 """

264

265 lines = self._outmol_lines()

266 found_cell = False

267 cell = np.zeros((3, 3))

268 symbols = []

269 positions = []

270 pattern_translation_vectors = re.compile(r'\s+translation\s+vector')

271 pattern_atomic_coordinates = re.compile(r'df\s+ATOMIC\s+COORDINATES')

272

273 for i, line in enumerate(lines):

274 if pattern_translation_vectors.match(line):

275 cell[int(line.split()[3]) - 1, :] = \

276 np.array([float(x) for x in line.split()[-3:]])

277 found_cell = True

278 if pattern_atomic_coordinates.match(line):

279 for ind, j in enumerate(range(i + 2, i + 2 + len(self.atoms))):

280 flds = lines[j].split()

281 symbols.append(flds[1])

282 positions.append(flds[2:5])

283 atoms = Atoms(symbols=symbols, positions=positions, cell=cell)

284 atoms.positions *= Bohr

285 atoms.cell *= Bohr

286

287 if found_cell:

288 atoms.pbc = [True, True, True]

289 atoms.wrap()

290 else:

291 atoms.pbc = [False, False, False]

292 return atoms

293

294 def read_energy(self):

295 """ Find and return last occurrence of Ef in outmole file. """

296 energy_regex = re.compile(r'^Ef\s+(\S+)Ha')

297 found = False

298 for line in self._outmol_lines():

299 match = energy_regex.match(line)

300 if match:

301 energy = float(match.group(1))

302 found = True

303 if not found:

304 raise RuntimeError('Could not read energy from outmol')

305 self.results['energy'] = energy * Hartree

306

307 def read_forces(self):

308 """ Read forces from .grad file. Applies self.rotation_matrix if

309 self.internal_transformation is True. """

310 with open(self.label + '.grad') as fd:

311 lines = fd.readlines()

312

313 forces = []

314 for i, line in enumerate(lines):

315 if line.startswith('$gradients'):

316 for j in range(i + 1, i + 1 + len(self.atoms)):

317 # force = - grad(Epot)

318 forces.append(np.array(

319 [-float(x) for x in lines[j].split()[1:4]]))

320

321 forces = np.array(forces) * Hartree / Bohr

322 if self.internal_transformation:

323 forces = np.dot(forces, self.rotation_matrix)

324 self.results['forces'] = forces

325

326 def get_eigenvalues(self, kpt=0, spin=0):

327 return self.read_eigenvalues(kpt, spin, 'eigenvalues')

328

329 def get_occupations(self, kpt=0, spin=0):

330 return self.read_eigenvalues(kpt, spin, 'occupations')

331

332 def get_k_point_weights(self):

333 return self.read_kpts(mode='k_point_weights')

334

335 def get_bz_k_points(self):

336 raise NotImplementedError

337

338 def get_ibz_k_points(self):

339 return self.read_kpts(mode='ibz_k_points')

340

341 def get_spin_polarized(self):

342 return self.read_spin_polarized()

343

344 def get_fermi_level(self):

345 return self.read_fermi()

346

347 def get_energy_contributions(self):

348 return self.read_energy_contributions()

349

350 def get_xc_functional(self):

351 return self.parameters['functional']

352

353 def read_eigenvalues(self, kpt=0, spin=0, mode='eigenvalues'):

354 """Reads eigenvalues from .outmol file.

355

356 This function splits into two situations:

357 1. We have no kpts just the raw eigenvalues ( Gamma point )

358 2. We have eigenvalues for each k-point

359

360 If calculation is spin_restricted then all eigenvalues

361 will be returned no matter what spin parameter is set to.

362

363 If calculation has no kpts then all eigenvalues

364 will be returned no matter what kpts parameter is set to.

365

366 Note DMol does usually NOT print all unoccupied eigenvalues.

367 Meaning number of eigenvalues for different kpts can vary.

368 """

369

370 assert mode in ['eigenvalues', 'occupations']

371 lines = self._outmol_lines()

372 pattern_kpts = re.compile(r'Eigenvalues for kvector\s+%d' % (kpt + 1))

373 for n, line in enumerate(lines):

374

375 # 1. We have no kpts

376 if line.split() == ['state', 'eigenvalue', 'occupation']:

377 spin_key = '+'

378 if self.get_spin_polarized():

379 if spin == 1:

380 spin_key = '-'

381 val_index = -2

382 if mode == 'occupations':

383 val_index = -1

384 values = []

385 m = n + 3

386 while True:

387 if lines[m].strip() == '':

388 break

389 flds = lines[m].split()

390 if flds[1] == spin_key:

391 values.append(float(flds[val_index]))

392 m += 1

393 return np.array(values)

394

395 # 2. We have kpts

396 if pattern_kpts.match(line):

397 val_index = 3

398 if self.get_spin_polarized():

399 if spin == 1:

400 val_index = 6

401 if mode == 'occupations':

402 val_index += 1

403 values = []

404 m = n + 2

405 while True:

406 if lines[m].strip() == '':

407 break

408 values.append(float(lines[m].split()[val_index]))

409 m += 1

410 return np.array(values)

411 return None

412

413 def _outmol_lines(self):

414 with open(self.label + '.outmol') as fd:

415 return fd.readlines()

416

417 def read_kpts(self, mode='ibz_k_points'):

418 """ Returns list of kpts coordinates or kpts weights. """

419

420 assert mode in ['ibz_k_points', 'k_point_weights']

421 lines = self._outmol_ines()

422

423 values = []

424 for n, line in enumerate(lines):

425 if line.startswith('Eigenvalues for kvector'):

426 if mode == 'ibz_k_points':

427 values.append([float(k_i)

428 for k_i in lines[n].split()[4:7]])

429 if mode == 'k_point_weights':

430 values.append(float(lines[n].split()[8]))

431 if values == []:

432 return None

433 return values

434

435 def read_spin_polarized(self):

436 """Reads, from outmol file, if calculation is spin polarized."""

437

438 lines = self._outmol_lines()

439 for n, line in enumerate(lines):

440 if line.rfind('Calculation is Spin_restricted') > -1:

441 return False

442 if line.rfind('Calculation is Spin_unrestricted') > -1:

443 return True

444 raise OSError('Could not read spin restriction from outmol')

445

446 def read_fermi(self):

447 """Reads the Fermi level.

448

449 Example line in outmol:

450 Fermi Energy: -0.225556 Ha -6.138 eV xyz text

451 """

452 lines = self._outmol_lines()

453 pattern_fermi = re.compile(r'Fermi Energy:\s+(\S+)\s+Ha')

454 for line in lines:

455 m = pattern_fermi.match(line)

456 if m:

457 return float(m.group(1)) * Hartree

458 return None

459

460 def read_energy_contributions(self):

461 """Reads the different energy contributions."""

462

463 lines = self._outmol_lines()

464 energies = {}

465 for n, line in enumerate(lines):

466 if line.startswith('Energy components'):

467 m = n + 1

468 while not lines[m].strip() == '':

469 energies[lines[m].split('=')[0].strip()] = \

470 float(re.findall(

471 r"[-+]?\d*\.\d+|\d+", lines[m])[0]) * Hartree

472 m += 1

473 return energies

474

475 def clean(self):

476 """ Cleanup after dmol calculation

477

478 Only removes dmol files in self.directory,

479 does not remove the directory itself

480 """

481 file_extensions = ['basis', 'car', 'err', 'grad', 'input', 'inatm',

482 'incoor', 'kpoints', 'monitor', 'occup', 'outmol',

483 'outatom', 'rot', 'sdf', 'sym', 'tpotl', 'tpdensk',

484 'torder', 'out', 'parameters.ase']

485 files_to_clean = ['DMol3.log', 'stdouterr.txt', 'mpd.hosts']

486

487 files = [os.path.join(self.directory, f) for f in files_to_clean]

488 files += [''.join((self.label, '.', ext)) for ext in file_extensions]

489

490 for f in files:

491 try:

492 os.remove(f)

493 except OSError:

494 pass

495

496

497# Helper functions

498# ------------------

499

500def find_transformation(atoms1, atoms2, verbose=False, only_cell=False):

501 """ Solves Ax = B where A and B are cell and positions from atoms objects.

502

503 Uses numpys least square solver to solve the problem Ax = B where A and

504 B are cell vectors and positions for atoms1 and atoms2 respectively.

505

506 Parameters

507 ----------

508 atoms1 (Atoms object): First atoms object (A)

509 atoms2 (Atoms object): Second atoms object (B)

510 verbose (bool): If True prints for each i A[i], B[i], Ax[i]

511 only_cell (bool): If True only cell in used, otherwise cell and positions.

512

513 Returns

514 -------

515 x (np.array((3,3))): Least square solution to Ax = B

516 error (float): The error calculated as np.linalg.norm(Ax-b)

517

518 """

519

520 if only_cell:

521 N = 3

522 elif len(atoms1) != len(atoms2):

523 raise RuntimeError('Atoms object must be of same length')

524 else:

525 N = len(atoms1) + 3

526

527 # Setup matrices A and B

528 A = np.zeros((N, 3))

529 B = np.zeros((N, 3))

530 A[0:3, :] = atoms1.cell

531 B[0:3, :] = atoms2.cell

532 if not only_cell:

533 A[3:, :] = atoms1.positions

534 B[3:, :] = atoms2.positions

535

536 # Solve least square problem Ax = B

537 lstsq_fit = np.linalg.lstsq(A, B, rcond=-1)

538 x = lstsq_fit[0]

539 error = np.linalg.norm(np.dot(A, x) - B)

540

541 # Print comparison between A, B and Ax

542 if verbose:

543 print('%17s %33s %35s %24s' % ('A', 'B', 'Ax', '|Ax-b|'))

544 for a, b in zip(A, B):

545 ax = np.dot(a, x)

546 loss = np.linalg.norm(ax - b)

547 print('(', end='')

548 for a_i in a:

549 print('%8.5f' % a_i, end='')

550 print(') (', end='')

551 for b_i in b:

552 print('%8.5f ' % b_i, end='')

553 print(') (', end='')

554 for ax_i in ax:

555 print('%8.5f ' % ax_i, end='')

556 print(') %8.5f' % loss)

557

558 return x, error

559

560

561def grd_to_file(atoms, grd_file, new_file):

562 """ Reads grd_file and converts data to cube format and writes to

563 cube_file.

564

565 Note: content of grd_file and atoms object are assumed to match with the

566 same orientation.

567

568 Parameters

569 -----------

570 atoms (Atoms object): atoms object grd_file data is for

571 grd_file (str): filename of .grd file

572 new_file (str): filename to write grd-data to, must be ASE format

573 that supports data argument

574 """

575 from ase.io import write

576

577 atoms_copy = atoms.copy()

578 data, cell, origin = read_grd(grd_file)

579 atoms_copy.cell = cell

580 atoms_copy.positions += origin

581 write(new_file, atoms_copy, data=data)

582

583

584def read_grd(filename):

585 """ Reads .grd file

586

587 Notes

588 -----

589 origin_xyz is offset with half a grid point in all directions to be

590 compatible with the cube format

591 Periodic systems is not guaranteed to be oriented correctly

592 """

593 from ase.geometry.cell import cellpar_to_cell

594

595 with open(filename) as fd:

596 lines = fd.readlines()

597

598 cell_data = np.array([float(fld) for fld in lines[2].split()])

599 cell = cellpar_to_cell(cell_data)

600 grid = [int(fld) + 1 for fld in lines[3].split()]

601 data = np.empty(grid)

602

603 origin_data = [int(fld) for fld in lines[4].split()[1:]]

604 origin_xyz = cell[0] * (-float(origin_data[0]) - 0.5) / (grid[0] - 1) + \

605 cell[1] * (-float(origin_data[2]) - 0.5) / (grid[1] - 1) + \

606 cell[2] * (-float(origin_data[4]) - 0.5) / (grid[2] - 1)

607

608 # Fastest index describes which index ( x or y ) varies fastest

609 # 1: x , 3: y

610 fastest_index = int(lines[4].split()[0])

611 assert fastest_index in [1, 3]

612 if fastest_index == 3:

613 grid[0], grid[1] = grid[1], grid[0]

614

615 dummy_counter = 5

616 for i in range(grid[2]):

617 for j in range(grid[1]):

618 for k in range(grid[0]): # Fastest index

619 if fastest_index == 1:

620 data[k, j, i] = float(lines[dummy_counter])

621 elif fastest_index == 3:

622 data[j, k, i] = float(lines[dummy_counter])

623 dummy_counter += 1

624

625 return data, cell, origin_xyz

626

627

628if __name__ == '__main__':

629 from ase.build import molecule

630

631 atoms = molecule('H2')

632 calc = DMol3()

633 atoms.calc = calc

634 # ~ 60 sec calculation

635 print('Potential energy %5.5f eV' % atoms.get_potential_energy())