Coverage for /builds/kinetik161/ase/ase/calculators/demon/demon.py: 40.85%

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

3http://www.demon-software.com

5"""

6import os

7import os.path as op

8import shutil

9import subprocess

11import numpy as np

13import ase.data

14import ase.io

15from ase.calculators.calculator import (FileIOCalculator, Parameters,

16 ReadError, all_changes, equal,

17 CalculatorSetupError)

18from ase.units import Bohr, Hartree

20from .demon_io import parse_xray

22m_e_to_amu = 1822.88839

25class Parameters_deMon(Parameters):

26 """Parameters class for the calculator.

27 Documented in Base_deMon.__init__

29 The options here are the most important ones that the user needs to be

30 aware of. Further options accepted by deMon can be set in the dictionary

31 input_arguments.

33 """

35 def __init__(

36 self,

37 label='rundir',

38 atoms=None,

39 restart=None,

40 basis_path=None,

41 ignore_bad_restart_file=FileIOCalculator._deprecated,

42 deMon_restart_path='.',

43 title='deMon input file',

44 scftype='RKS',

45 forces=False,

46 dipole=False,

47 xc='VWN',

48 guess='TB',

49 print_out='MOE',

50 basis={},

51 ecps={},

52 mcps={},

53 auxis={},

54 augment={},

55 input_arguments=None):

56 kwargs = locals()

57 kwargs.pop('self')

58 Parameters.__init__(self, **kwargs)

61class Demon(FileIOCalculator):

62 """Calculator interface to the deMon code. """

64 implemented_properties = [

65 'energy',

66 'forces',

67 'dipole',

68 'eigenvalues']

70 def __init__(self, *, command=None, **kwargs):

71 """ASE interface to the deMon code.

73 The deMon2k code can be obtained from http://www.demon-software.com

75 The DEMON_COMMAND environment variable must be set to run the

76 executable, in bash it would be set along the lines of

77 export DEMON_COMMAND="deMon.4.3.6.std > deMon_ase.out 2>&1"

79 Parameters:

81 label : str

82 relative path to the run directory

83 atoms : Atoms object

84 the atoms object

85 command : str

86 Command to run deMon. If not present the environment

87 variable DEMON_COMMAND will be used

88 restart : str

89 Relative path to ASE restart directory for parameters and

90 atoms object and results

91 basis_path : str

92 Relative path to the directory containing

93 BASIS, AUXIS, ECPS, MCPS and AUGMENT

94 ignore_bad_restart_file : bool

95 Ignore broken or missing ASE restart files

96 By default, it is an error if the restart

97 file is missing or broken.

98 deMon_restart_path : str

99 Relative path to the deMon restart dir

100 title : str

101 Title in the deMon input file.

102 scftype : str

103 Type of scf

104 forces : bool

105 If True a force calculation will be enforced.

106 dipole : bool

107 If True a dipole calculation will be enforced

108 xc : str

109 xc-functional

110 guess : str

111 guess for initial density and wave functions

112 print_out : str | list

113 Options for the printing in deMon

114 basis : dict

115 Definition of basis sets.

116 ecps : dict

117 Definition of ECPs

118 mcps : dict

119 Definition of MCPs

120 auxis : dict

121 Definition of AUXIS

122 augment : dict

123 Definition of AUGMENT

124 input_arguments : dict

125 Explicitly given input arguments. The key is the input keyword

126 and the value is either a str, a list of str (will be written

127 on the same line as the keyword),

128 or a list of lists of str (first list is written on the first

129 line, the others on following lines.)

130

131 For example usage, see the tests h2o.py and h2o_xas_xes.py in

132 the directory ase/test/demon

133

134 """

135

136 parameters = Parameters_deMon(**kwargs)

137

138 # Setup the run command

139 if command is None:

140 command = self.cfg.get('DEMON_COMMAND')

141

142 FileIOCalculator.__init__(

143 self,

144 command=command,

145 **parameters)

146

147 def __getitem__(self, key):

148 """Convenience method to retrieve a parameter as

149 calculator[key] rather than calculator.parameters[key]

150

151 Parameters:

152 key : str, the name of the parameters to get.

153 """

154 return self.parameters[key]

155

156 def set(self, **kwargs):

157 """Set all parameters.

158

159 Parameters:

160 kwargs : Dictionary containing the keywords for deMon

161 """

162 # Put in the default arguments.

163 kwargs = self.default_parameters.__class__(**kwargs)

164

165 if 'parameters' in kwargs:

166 filename = kwargs.pop('parameters')

167 parameters = Parameters.read(filename)

168 parameters.update(kwargs)

169 kwargs = parameters

170

171 changed_parameters = {}

172

173 for key, value in kwargs.items():

174 oldvalue = self.parameters.get(key)

175 if key not in self.parameters or not equal(value, oldvalue):

176 changed_parameters[key] = value

177 self.parameters[key] = value

178

179 return changed_parameters

180

181 def link_file(self, fromdir, todir, filename):

182 if op.exists(todir + '/' + filename):

183 os.remove(todir + '/' + filename)

184

185 if op.exists(fromdir + '/' + filename):

186 os.symlink(fromdir + '/' + filename,

187 todir + '/' + filename)

188 else:

189 raise RuntimeError(

190 "{} doesn't exist".format(fromdir + '/' + filename))

191

192 def calculate(self,

193 atoms=None,

194 properties=['energy'],

195 system_changes=all_changes):

196 """Capture the RuntimeError from FileIOCalculator.calculate

197 and add a little debug information from the deMon output.

198

199 See base FileIocalculator for documentation.

200 """

201

202 if atoms is not None:

203 self.atoms = atoms.copy()

204

205 self.write_input(self.atoms, properties, system_changes)

206 command = self.command

207

208 # basis path

209 basis_path = self.parameters['basis_path']

210 if basis_path is None:

211 basis_path = self.cfg.get('DEMON_BASIS_PATH')

212

213 if basis_path is None:

214 raise RuntimeError('Please set basis_path keyword,' +

215 ' or the DEMON_BASIS_PATH' +

216 ' environment variable')

217

218 # link restart file

219 value = self.parameters['guess']

220 if value.upper() == 'RESTART':

221 value2 = self.parameters['deMon_restart_path']

222

223 if op.exists(self.directory + '/deMon.rst')\

224 or op.islink(self.directory + '/deMon.rst'):

225 os.remove(self.directory + '/deMon.rst')

226 abspath = op.abspath(value2)

227

228 if op.exists(abspath + '/deMon.mem') \

229 or op.islink(abspath + '/deMon.mem'):

230

231 shutil.copy(abspath + '/deMon.mem',

232 self.directory + '/deMon.rst')

233 else:

234 raise RuntimeError(

235 "{} doesn't exist".format(abspath + '/deMon.rst'))

236

237 abspath = op.abspath(basis_path)

238

239 for name in ['BASIS', 'AUXIS', 'ECPS', 'MCPS', 'FFDS']:

240 self.link_file(abspath, self.directory, name)

241

242 if command is None:

243 raise CalculatorSetupError

244 subprocess.check_call(command, shell=True, cwd=self.directory)

245

246 try:

247 self.read_results()

248 except Exception: # XXX Which kind of exception?

249 with open(self.directory + '/deMon.out') as fd:

250 lines = fd.readlines()

251 debug_lines = 10

252 print('##### %d last lines of the deMon.out' % debug_lines)

253 for line in lines[-20:]:

254 print(line.strip())

255 print('##### end of deMon.out')

256 raise RuntimeError

257

258 def set_label(self, label):

259 """Set label directory """

260

261 self.label = label

262

263 # in our case self.directory = self.label

264 self.directory = self.label

265 if self.directory == '':

266 self.directory = os.curdir

267

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

269 """Write input (in)-file.

270 See calculator.py for further details.

271

272 Parameters:

273 atoms : The Atoms object to write.

274 properties : The properties which should be calculated.

275 system_changes : List of properties changed since last run.

276

277 """

278 # Call base calculator.

279 FileIOCalculator.write_input(

280 self,

281 atoms=atoms,

282 properties=properties,

283 system_changes=system_changes)

284

285 if system_changes is None and properties is None:

286 return

287

288 filename = f'{self.directory}/deMon.inp'

289

290 add_print = ''

291

292 # Start writing the file.

293 with open(filename, 'w') as fd:

294

295 # write keyword argument keywords

296 value = self.parameters['title']

297 self._write_argument('TITLE', value, fd)

298

299 fd.write('#\n')

300

301 value = self.parameters['scftype']

302 self._write_argument('SCFTYPE', value, fd)

303

304 value = self.parameters['xc']

305 self._write_argument('VXCTYPE', value, fd)

306

307 value = self.parameters['guess']

308 self._write_argument('GUESS', value, fd)

309

310 # obtain forces through a single BOMD step

311 # only if forces is in properties, or if keyword forces is True

312 value = self.parameters['forces']

313 if 'forces' in properties or value:

314

315 self._write_argument('DYNAMICS',

316 ['INT=1', 'MAX=0', 'STEP=0'], fd)

317 self._write_argument('TRAJECTORY', 'FORCES', fd)

318 self._write_argument('VELOCITIES', 'ZERO', fd)

319 add_print = add_print + ' ' + 'MD OPT'

320

321 # if dipole is True, enforce dipole calculation.

322 # Otherwise only if asked for

323 value = self.parameters['dipole']

324 if 'dipole' in properties or value:

325 self._write_argument('DIPOLE', '', fd)

326

327 # print argument, here other options could change this

328 value = self.parameters['print_out']

329 assert isinstance(value, str)

330 value = value + add_print

331

332 if not len(value) == 0:

333 self._write_argument('PRINT', value, fd)

334 fd.write('#\n')

335

336 # write general input arguments

337 self._write_input_arguments(fd)

338

339 fd.write('#\n')

340

341 # write basis set, ecps, mcps, auxis, augment

342 basis = self.parameters['basis']

343 if 'all' not in basis:

344 basis['all'] = 'DZVP'

345 self._write_basis(fd, atoms, basis, string='BASIS')

346

347 ecps = self.parameters['ecps']

348 if not len(ecps) == 0:

349 self._write_basis(fd, atoms, ecps, string='ECPS')

350

351 mcps = self.parameters['mcps']

352 if not len(mcps) == 0:

353 self._write_basis(fd, atoms, mcps, string='MCPS')

354

355 auxis = self.parameters['auxis']

356 if not len(auxis) == 0:

357 self._write_basis(fd, atoms, auxis, string='AUXIS')

358

359 augment = self.parameters['augment']

360 if not len(augment) == 0:

361 self._write_basis(fd, atoms, augment, string='AUGMENT')

362

363 # write geometry

364 self._write_atomic_coordinates(fd, atoms)

365

366 # write xyz file for good measure.

367 ase.io.write(f'{self.directory}/deMon_atoms.xyz', self.atoms)

368

369 def read(self, restart_path):

370 """Read parameters from directory restart_path."""

371

372 self.set_label(restart_path)

373

374 if not op.exists(restart_path + '/deMon.inp'):

375 raise ReadError('The restart_path file {} does not exist'

376 .format(restart_path))

377

378 self.atoms = self.deMon_inp_to_atoms(restart_path + '/deMon.inp')

379

380 self.read_results()

381

382 def _write_input_arguments(self, fd):

383 """Write directly given input-arguments."""

384 input_arguments = self.parameters['input_arguments']

385

386 # Early return

387 if input_arguments is None:

388 return

389

390 for key, value in input_arguments.items():

391 self._write_argument(key, value, fd)

392

393 def _write_argument(self, key, value, fd):

394 """Write an argument to file.

395 key : a string coresponding to the input keyword

396 value : the arguments, can be a string, a number or a list

397 f : and open file

398 """

399

400 # for only one argument, write on same line

401 if not isinstance(value, (tuple, list)):

402 line = key.upper()

403 line += ' ' + str(value).upper()

404 fd.write(line)

405 fd.write('\n')

406

407 # for a list, write first argument on the first line,

408 # then the rest on new lines

409 else:

410 line = key

411 if not isinstance(value[0], (tuple, list)):

412 for i in range(len(value)):

413 line += ' ' + str(value[i].upper())

414 fd.write(line)

415 fd.write('\n')

416 else:

417 for i in range(len(value)):

418 for j in range(len(value[i])):

419 line += ' ' + str(value[i][j]).upper()

420 fd.write(line)

421 fd.write('\n')

422 line = ''

423

424 def _write_atomic_coordinates(self, fd, atoms):

425 """Write atomic coordinates.

426

427 Parameters:

428 - f: An open file object.

429 - atoms: An atoms object.

430 """

431

432 fd.write('#\n')

433 fd.write('# Atomic coordinates\n')

434 fd.write('#\n')

435 fd.write('GEOMETRY CARTESIAN ANGSTROM\n')

436

437 for i in range(len(atoms)):

438 xyz = atoms.get_positions()[i]

439 chem_symbol = atoms.get_chemical_symbols()[i]

440 chem_symbol += str(i + 1)

441

442 # if tag is set to 1 then we have a ghost atom,

443 # set nuclear charge to 0

444 if atoms.get_tags()[i] == 1:

445 nuc_charge = str(0)

446 else:

447 nuc_charge = str(atoms.get_atomic_numbers()[i])

448

449 mass = atoms.get_masses()[i]

450

451 line = f'{chem_symbol:6s}'.rjust(10) + ' '

452 line += f'{xyz[0]:.5f}'.rjust(10) + ' '

453 line += f'{xyz[1]:.5f}'.rjust(10) + ' '

454 line += f'{xyz[2]:.5f}'.rjust(10) + ' '

455 line += f'{nuc_charge:5s}'.rjust(10) + ' '

456 line += f'{mass:.5f}'.rjust(10) + ' '

457

458 fd.write(line)

459 fd.write('\n')

460

461 # routine to write basis set inormation, including ecps and auxis

462 def _write_basis(self, fd, atoms, basis={}, string='BASIS'):

463 """Write basis set, ECPs, AUXIS, or AUGMENT basis

464

465 Parameters:

466 - f: An open file object.

467 - atoms: An atoms object.

468 - basis: A dictionary specifying the basis set

469 - string: 'BASIS', 'ECP','AUXIS' or 'AUGMENT'

470 """

471

472 # basis for all atoms

473 line = f'{string}'.ljust(10)

474

475 if 'all' in basis:

476 default_basis = basis['all']

477 line += f'({default_basis})'.rjust(16)

478

479 fd.write(line)

480 fd.write('\n')

481

482 # basis for all atomic species

483 chemical_symbols = atoms.get_chemical_symbols()

484 chemical_symbols_set = set(chemical_symbols)

485

486 for i in range(chemical_symbols_set.__len__()):

487 symbol = chemical_symbols_set.pop()

488

489 if symbol in basis:

490 line = f'{symbol}'.ljust(10)

491 line += f'({basis[symbol]})'.rjust(16)

492 fd.write(line)

493 fd.write('\n')

494

495 # basis for individual atoms

496 for i in range(len(atoms)):

497

498 if i in basis:

499 symbol = str(chemical_symbols[i])

500 symbol += str(i + 1)

501

502 line = f'{symbol}'.ljust(10)

503 line += f'({basis[i]})'.rjust(16)

504 fd.write(line)

505 fd.write('\n')

506

507 # Analysis routines

508 def read_results(self):

509 """Read the results from output files."""

510 self.read_energy()

511 self.read_forces(self.atoms)

512 self.read_eigenvalues()

513 self.read_dipole()

514 self.read_xray()

515

516 def read_energy(self):

517 """Read energy from deMon's text-output file."""

518 with open(self.label + '/deMon.out') as fd:

519 text = fd.read().upper()

520

521 lines = iter(text.split('\n'))

522

523 for line in lines:

524 if line.startswith(' TOTAL ENERGY ='):

525 self.results['energy'] = float(line.split()[-1]) * Hartree

526 break

527 else:

528 raise RuntimeError

529

530 def read_forces(self, atoms):

531 """Read the forces from the deMon.out file."""

532

533 natoms = len(atoms)

534 filename = self.label + '/deMon.out'

535

536 if op.isfile(filename):

537 with open(filename) as fd:

538 lines = fd.readlines()

539

540 # find line where the orbitals start

541 flag_found = False

542 for i in range(len(lines)):

543 if lines[i].rfind('GRADIENTS OF TIME STEP 0 IN A.U.') > -1:

544 start = i + 4

545 flag_found = True

546 break

547

548 if flag_found:

549 self.results['forces'] = np.zeros((natoms, 3), float)

550 for i in range(natoms):

551 line = [s for s in lines[i + start].strip().split(' ')

552 if len(s) > 0]

553 f = -np.array([float(x) for x in line[2:5]])

554 self.results['forces'][i, :] = f * (Hartree / Bohr)

555

556 def read_eigenvalues(self):

557 """Read eigenvalues from the 'deMon.out' file."""

558 assert os.access(self.label + '/deMon.out', os.F_OK)

559

560 # Read eigenvalues

561 with open(self.label + '/deMon.out') as fd:

562 lines = fd.readlines()

563

564 # try PRINT MOE

565 eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(

566 lines, 'ALPHA MO ENERGIES', 6)

567 eig_beta, occ_beta = self.read_eigenvalues_one_spin(

568 lines, 'BETA MO ENERGIES', 6)

569

570 # otherwise try PRINT MOS

571 if len(eig_alpha) == 0 and len(eig_beta) == 0:

572 eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(

573 lines, 'ALPHA MO COEFFICIENTS', 5)

574 eig_beta, occ_beta = self.read_eigenvalues_one_spin(

575 lines, 'BETA MO COEFFICIENTS', 5)

576

577 self.results['eigenvalues'] = np.array([eig_alpha, eig_beta]) * Hartree

578 self.results['occupations'] = np.array([occ_alpha, occ_beta])

579

580 def read_eigenvalues_one_spin(self, lines, string, neigs_per_line):

581 """Utility method for retreiving eigenvalues after the string "string"

582 with neigs_per_line eigenvlaues written per line

583 """

584 eig = []

585 occ = []

586

587 skip_line = False

588 more_eigs = False

589

590 # find line where the orbitals start

591 for i in range(len(lines)):

592 if lines[i].rfind(string) > -1:

593 ii = i

594 more_eigs = True

595 break

596

597 while more_eigs:

598 # search for two empty lines in a row preceding a line with

599 # numbers

600 for i in range(ii + 1, len(lines)):

601 if len(lines[i].split()) == 0 and \

602 len(lines[i + 1].split()) == 0 and \

603 len(lines[i + 2].split()) > 0:

604 ii = i + 2

605 break

606

607 # read eigenvalues, occupations

608 line = lines[ii].split()

609 if len(line) < neigs_per_line:

610 # last row

611 more_eigs = False

612 if line[0] != str(len(eig) + 1):

613 more_eigs = False

614 skip_line = True

615

616 if not skip_line:

617 line = lines[ii + 1].split()

618 for l in line:

619 eig.append(float(l))

620 line = lines[ii + 3].split()

621 for l in line:

622 occ.append(float(l))

623 ii = ii + 3

624

625 return eig, occ

626

627 def read_dipole(self):

628 """Read dipole moment."""

629 dipole = np.zeros(3)

630 with open(self.label + '/deMon.out') as fd:

631 lines = fd.readlines()

632

633 for i in range(len(lines)):

634 if lines[i].rfind('DIPOLE') > - \

635 1 and lines[i].rfind('XAS') == -1:

636 dipole[0] = float(lines[i + 1].split()[3])

637 dipole[1] = float(lines[i + 2].split()[3])

638 dipole[2] = float(lines[i + 3].split()[3])

639

640 # debye to e*Ang

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

642

643 break

644

645 def read_xray(self):

646 """Read deMon.xry if present."""

647

648 # try to read core IP from, .out file

649 filename = self.label + '/deMon.out'

650 core_IP = None

651 if op.isfile(filename):

652 with open(filename) as fd:

653 lines = fd.readlines()

654

655 for i in range(len(lines)):

656 if lines[i].rfind('IONIZATION POTENTIAL') > -1:

657 core_IP = float(lines[i].split()[3])

658

659 try:

660 mode, ntrans, E_trans, osc_strength, trans_dip = parse_xray(

661 self.label + '/deMon.xry')

662 except ReadError:

663 pass

664 else:

665 xray_results = {'xray_mode': mode,

666 'ntrans': ntrans,

667 'E_trans': E_trans,

668 'osc_strength': osc_strength, # units?

669 'trans_dip': trans_dip, # units?

670 'core_IP': core_IP}

671

672 self.results['xray'] = xray_results

673

674 def deMon_inp_to_atoms(self, filename):

675 """Routine to read deMon.inp and convert it to an atoms object."""

676

677 with open(filename) as fd:

678 lines = fd.readlines()

679

680 # find line where geometry starts

681 for i in range(len(lines)):

682 if lines[i].rfind('GEOMETRY') > -1:

683 if lines[i].rfind('ANGSTROM'):

684 coord_units = 'Ang'

685 elif lines.rfind('Bohr'):

686 coord_units = 'Bohr'

687 ii = i

688 break

689

690 chemical_symbols = []

691 xyz = []

692 atomic_numbers = []

693 masses = []

694

695 for i in range(ii + 1, len(lines)):

696 try:

697 line = lines[i].split()

698

699 if len(line) > 0:

700 for symbol in ase.data.chemical_symbols:

701 found = None

702 if line[0].upper().rfind(symbol.upper()) > -1:

703 found = symbol

704 break

705

706 if found is not None:

707 chemical_symbols.append(found)

708 else:

709 break

710

711 xyz.append(

712 [float(line[1]), float(line[2]), float(line[3])])

713

714 if len(line) > 4:

715 atomic_numbers.append(int(line[4]))

716

717 if len(line) > 5:

718 masses.append(float(line[5]))

719

720 except Exception: # XXX Which kind of exception?

721 raise RuntimeError

722

723 if coord_units == 'Bohr':

724 xyz *= Bohr

725

726 natoms = len(chemical_symbols)

727

728 # set atoms object

729 atoms = ase.Atoms(symbols=chemical_symbols, positions=xyz)

730

731 # if atomic numbers were read in, set them

732 if len(atomic_numbers) == natoms:

733 atoms.set_atomic_numbers(atomic_numbers)

734

735 # if masses were read in, set them

736 if len(masses) == natoms:

737 atoms.set_masses(masses)

738

739 return atoms