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

1# fmt: off

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

5Contacts

6--------

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

8Erik Fransson <erikfr@chalmers.se>

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

12DMol3 environment variables

13----------------------------

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

15to prallelize over

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

20Example

21--------

22>>> from ase.build import bulk

23>>> from ase.calculators import DMol3

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

26>>> calc = DMol3()

27>>> atoms.calc = calc

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

31DMol3 calculator functionality

32-------------------------------

33This calculator does support all the functionality in DMol3.

35Firstly this calculator is limited to only handling either fully

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

38Internal relaxations are not supported by the calculator,

39only support for energy and forces is implemented.

41Reading eigenvalues and kpts are supported.

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

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

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

46format using grd_to_cube().

49DMol3 internal coordinates

50---------------------------

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

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

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

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

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

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

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

59parsed from the .rot file.

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

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

62DMol_atoms * rot_mat = ase_atoms

65DMol3 files

66------------

67The supported DMol3 file formats are:

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

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

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

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

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

74grd outfile for orbitals from DMol3 - cellpar in Angstrom

76"""

78import os

79import re

81import numpy as np

83from ase import Atoms

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

85from ase.io import read

86from ase.io.dmol import write_dmol_car, write_dmol_incoor

87from ase.units import Bohr, Hartree

90class DMol3(FileIOCalculator):

91 """ DMol3 calculator object. """

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

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

95 'symmetry': 'on'}

96 discard_results_on_any_change = True

98 def __init__(self, restart=None,

99 ignore_bad_restart_file=FileIOCalculator._deprecated,

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

101 command=None, **kwargs):

102 """ Construct DMol3 calculator. """

103

104 if command is None:

105 if 'DMOL_COMMAND' in self.cfg:

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

107

108 super().__init__(restart, ignore_bad_restart_file,

109 label, atoms, command=command,

110 **kwargs)

111

112 # tracks if DMol transformed coordinate system

113 self.internal_transformation = False

114

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

116

117 if not np.all(atoms.pbc) and np.any(atoms.pbc):

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

119

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

121 self.internal_transformation = False

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

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

124

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

126

127 if np.all(atoms.pbc):

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

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

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

131

132 self.write_input_file()

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

134

135 def write_input_file(self):

136 """ Writes the input file. """

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

138 self._write_input_file(fd)

139

140 def _write_input_file(self, fd):

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

142

143 # if no key about eigs

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

145

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

147 if isinstance(value, str):

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

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

150 for val in value:

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

152 else:

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

154

155 def read(self, label):

156 FileIOCalculator.read(self, label)

157 geometry = self.label + '.car'

158 output = self.label + '.outmol'

159 force = self.label + '.grad'

160

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

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

163 raise ReadError

164

165 self.atoms = read(geometry)

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

167 self.read_results()

168

169 def read_results(self):

170 finished, message = self.finished_successfully()

171 if not finished:

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

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

174

175 self.find_dmol_transformation()

176 self.read_energy()

177 self.read_forces()

178

179 def finished_successfully(self):

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

181

182 Returns

183 -------

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

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

186

187 """

188 finished = False

189 message = ""

190 for line in self._outmol_lines():

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

192 finished = True

193 if line.startswith('Error'):

194 message += line

195 return finished, message

196

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

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

199 coordinates to ase coordinates.

200

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

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

203

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

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

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

207 function find_transformation(atoms1, atoms2)

208

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

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

211

212 Parameters

213 ----------

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

215 """

216

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

218 dmol_atoms = self.read_atoms_from_outmol()

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

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

221 dmol_atoms.cell) < tol):

222 self.internal_transformation = False

223 else:

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

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

226 raise RuntimeError('Error: transformation matrix does'

227 ' not have determinant 1.0')

228 if err < tol:

229 self.internal_transformation = True

230 self.rotation_matrix = R

231 else:

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

233 ' coordinate transformation')

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

235 try:

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

237 except OSError:

238 self.internal_transformation = False

239 else:

240 self.internal_transformation = True

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

242

243 def read_atoms_from_outmol(self):

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

245 object.

246

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

248 pbc 000.

249

250 Formatting for cell in outmol :

251 translation vector [a0] 1 5.1 0.0 5.1

252 translation vector [a0] 2 5.1 5.1 0.0

253 translation vector [a0] 3 0.0 5.1 5.1

254

255 Formatting for positions in outmol:

256 df ATOMIC COORDINATES (au)

257 df x y z

258 df Si 0.0 0.0 0.0

259 df Si 1.3 3.5 2.2

260 df binding energy -0.2309046Ha

261

262 Returns

263 -------

264 atoms (Atoms object): read atoms object

265 """

266

267 lines = self._outmol_lines()

268 found_cell = False

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

270 symbols = []

271 positions = []

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

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

274

275 for i, line in enumerate(lines):

276 if pattern_translation_vectors.match(line):

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

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

279 found_cell = True

280 if pattern_atomic_coordinates.match(line):

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

282 flds = lines[j].split()

283 symbols.append(flds[1])

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

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

286 atoms.positions *= Bohr

287 atoms.cell *= Bohr

288

289 if found_cell:

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

291 atoms.wrap()

292 else:

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

294 return atoms

295

296 def read_energy(self):

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

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

299 found = False

300 for line in self._outmol_lines():

301 match = energy_regex.match(line)

302 if match:

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

304 found = True

305 if not found:

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

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

308

309 def read_forces(self):

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

311 self.internal_transformation is True. """

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

313 lines = fd.readlines()

314

315 forces = []

316 for i, line in enumerate(lines):

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

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

319 # force = - grad(Epot)

320 forces.append(np.array(

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

322

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

324 if self.internal_transformation:

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

326 self.results['forces'] = forces

327

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

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

330

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

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

333

334 def get_k_point_weights(self):

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

336

337 def get_bz_k_points(self):

338 raise NotImplementedError

339

340 def get_ibz_k_points(self):

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

342

343 def get_spin_polarized(self):

344 return self.read_spin_polarized()

345

346 def get_fermi_level(self):

347 return self.read_fermi()

348

349 def get_energy_contributions(self):

350 return self.read_energy_contributions()

351

352 def get_xc_functional(self):

353 return self.parameters['functional']

354

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

356 """Reads eigenvalues from .outmol file.

357

358 This function splits into two situations:

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

360 2. We have eigenvalues for each k-point

361

362 If calculation is spin_restricted then all eigenvalues

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

364

365 If calculation has no kpts then all eigenvalues

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

367

368 Note DMol does usually NOT print all unoccupied eigenvalues.

369 Meaning number of eigenvalues for different kpts can vary.

370 """

371

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

373 lines = self._outmol_lines()

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

375 for n, line in enumerate(lines):

376

377 # 1. We have no kpts

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

379 spin_key = '+'

380 if self.get_spin_polarized():

381 if spin == 1:

382 spin_key = '-'

383 val_index = -2

384 if mode == 'occupations':

385 val_index = -1

386 values = []

387 m = n + 3

388 while True:

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

390 break

391 flds = lines[m].split()

392 if flds[1] == spin_key:

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

394 m += 1

395 return np.array(values)

396

397 # 2. We have kpts

398 if pattern_kpts.match(line):

399 val_index = 3

400 if self.get_spin_polarized():

401 if spin == 1:

402 val_index = 6

403 if mode == 'occupations':

404 val_index += 1

405 values = []

406 m = n + 2

407 while True:

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

409 break

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

411 m += 1

412 return np.array(values)

413 return None

414

415 def _outmol_lines(self):

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

417 return fd.readlines()

418

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

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

421

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

423 lines = self._outmol_ines()

424

425 values = []

426 for n, line in enumerate(lines):

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

428 if mode == 'ibz_k_points':

429 values.append([float(k_i)

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

431 if mode == 'k_point_weights':

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

433 if values == []:

434 return None

435 return values

436

437 def read_spin_polarized(self):

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

439

440 lines = self._outmol_lines()

441 for n, line in enumerate(lines):

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

443 return False

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

445 return True

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

447

448 def read_fermi(self):

449 """Reads the Fermi level.

450

451 Example line in outmol:

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

453 """

454 lines = self._outmol_lines()

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

456 for line in lines:

457 m = pattern_fermi.match(line)

458 if m:

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

460 return None

461

462 def read_energy_contributions(self):

463 """Reads the different energy contributions."""

464

465 lines = self._outmol_lines()

466 energies = {}

467 for n, line in enumerate(lines):

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

469 m = n + 1

470 while lines[m].strip() != '':

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

472 float(re.findall(

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

474 m += 1

475 return energies

476

477 def clean(self):

478 """ Cleanup after dmol calculation

479

480 Only removes dmol files in self.directory,

481 does not remove the directory itself

482 """

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

484 'incoor', 'kpoints', 'monitor', 'occup', 'outmol',

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

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

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

488

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

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

491

492 for f in files:

493 try:

494 os.remove(f)

495 except OSError:

496 pass

497

498

499# Helper functions

500# ------------------

501

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

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

504

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

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

507

508 Parameters

509 ----------

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

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

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

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

514

515 Returns

516 -------

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

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

519

520 """

521

522 if only_cell:

523 N = 3

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

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

526 else:

527 N = len(atoms1) + 3

528

529 # Setup matrices A and B

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

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

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

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

534 if not only_cell:

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

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

537

538 # Solve least square problem Ax = B

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

540 x = lstsq_fit[0]

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

542

543 # Print comparison between A, B and Ax

544 if verbose:

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

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

547 ax = np.dot(a, x)

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

549 print('(', end='')

550 for a_i in a:

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

552 print(') (', end='')

553 for b_i in b:

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

555 print(') (', end='')

556 for ax_i in ax:

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

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

559

560 return x, error

561

562

563def grd_to_file(atoms, grd_file, new_file):

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

565 cube_file.

566

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

568 same orientation.

569

570 Parameters

571 -----------

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

573 grd_file (str): filename of .grd file

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

575 that supports data argument

576 """

577 from ase.io import write

578

579 atoms_copy = atoms.copy()

580 data, cell, origin = read_grd(grd_file)

581 atoms_copy.cell = cell

582 atoms_copy.positions += origin

583 write(new_file, atoms_copy, data=data)

584

585

586def read_grd(filename):

587 """ Reads .grd file

588

589 Notes

590 -----

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

592 compatible with the cube format

593 Periodic systems is not guaranteed to be oriented correctly

594 """

595 from ase.geometry.cell import cellpar_to_cell

596

597 with open(filename) as fd:

598 lines = fd.readlines()

599

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

601 cell = cellpar_to_cell(cell_data)

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

603 data = np.empty(grid)

604

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

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

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

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

609

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

611 # 1: x , 3: y

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

613 assert fastest_index in [1, 3]

614 if fastest_index == 3:

615 grid[0], grid[1] = grid[1], grid[0]

616

617 dummy_counter = 5

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

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

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

621 if fastest_index == 1:

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

623 elif fastest_index == 3:

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

625 dummy_counter += 1

626

627 return data, cell, origin_xyz

628

629

630if __name__ == '__main__':

631 from ase.build import molecule

632

633 atoms = molecule('H2')

634 calc = DMol3()

635 atoms.calc = calc

636 # ~ 60 sec calculation

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