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

1# fmt: off

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

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

7"""

8import os

9import os.path as op

10import shutil

11import subprocess

13import numpy as np

15import ase.data

16import ase.io

17from ase.calculators.calculator import (

18 CalculatorSetupError,

19 FileIOCalculator,

20 Parameters,

21 ReadError,

22 all_changes,

23 equal,

24)

25from ase.units import Bohr, Hartree

27from .demon_io import parse_xray

29m_e_to_amu = 1822.88839

32class Parameters_deMon(Parameters):

33 """Parameters class for the calculator.

34 Documented in Base_deMon.__init__

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

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

38 input_arguments.

40 """

42 def __init__(

43 self,

44 label='rundir',

45 atoms=None,

46 restart=None,

47 basis_path=None,

48 ignore_bad_restart_file=FileIOCalculator._deprecated,

49 deMon_restart_path='.',

50 title='deMon input file',

51 scftype='RKS',

52 forces=False,

53 dipole=False,

54 xc='VWN',

55 guess='TB',

56 print_out='MOE',

57 basis={},

58 ecps={},

59 mcps={},

60 auxis={},

61 augment={},

62 input_arguments=None):

63 kwargs = locals()

64 kwargs.pop('self')

65 Parameters.__init__(self, **kwargs)

68class Demon(FileIOCalculator):

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

71 implemented_properties = [

72 'energy',

73 'forces',

74 'dipole',

75 'eigenvalues']

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

78 """ASE interface to the deMon code.

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

82 The DEMON_COMMAND environment variable must be set to run the

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

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

86 Parameters:

88 label : str

89 relative path to the run directory

90 atoms : Atoms object

91 the atoms object

92 command : str

93 Command to run deMon. If not present the environment

94 variable DEMON_COMMAND will be used

95 restart : str

96 Relative path to ASE restart directory for parameters and

97 atoms object and results

98 basis_path : str

99 Relative path to the directory containing

100 BASIS, AUXIS, ECPS, MCPS and AUGMENT

101 ignore_bad_restart_file : bool

102 Ignore broken or missing ASE restart files

103 By default, it is an error if the restart

104 file is missing or broken.

105 deMon_restart_path : str

106 Relative path to the deMon restart dir

107 title : str

108 Title in the deMon input file.

109 scftype : str

110 Type of scf

111 forces : bool

112 If True a force calculation will be enforced.

113 dipole : bool

114 If True a dipole calculation will be enforced

115 xc : str

116 xc-functional

117 guess : str

118 guess for initial density and wave functions

119 print_out : str | list

120 Options for the printing in deMon

121 basis : dict

122 Definition of basis sets.

123 ecps : dict

124 Definition of ECPs

125 mcps : dict

126 Definition of MCPs

127 auxis : dict

128 Definition of AUXIS

129 augment : dict

130 Definition of AUGMENT

131 input_arguments : dict

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

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

134 on the same line as the keyword),

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

136 line, the others on following lines.)

137

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

139 the directory ase/test/demon

140

141 """

142

143 parameters = Parameters_deMon(**kwargs)

144

145 # Setup the run command

146 if command is None:

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

148

149 FileIOCalculator.__init__(

150 self,

151 command=command,

152 **parameters)

153

154 def __getitem__(self, key):

155 """Convenience method to retrieve a parameter as

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

157

158 Parameters:

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

160 """

161 return self.parameters[key]

162

163 def set(self, **kwargs):

164 """Set all parameters.

165

166 Parameters:

167 kwargs : Dictionary containing the keywords for deMon

168 """

169 # Put in the default arguments.

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

171

172 if 'parameters' in kwargs:

173 filename = kwargs.pop('parameters')

174 parameters = Parameters.read(filename)

175 parameters.update(kwargs)

176 kwargs = parameters

177

178 changed_parameters = {}

179

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

181 oldvalue = self.parameters.get(key)

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

183 changed_parameters[key] = value

184 self.parameters[key] = value

185

186 return changed_parameters

187

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

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

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

191

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

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

194 todir + '/' + filename)

195 else:

196 raise RuntimeError(

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

198

199 def calculate(self,

200 atoms=None,

201 properties=['energy'],

202 system_changes=all_changes):

203 """Capture the RuntimeError from FileIOCalculator.calculate

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

205

206 See base FileIocalculator for documentation.

207 """

208

209 if atoms is not None:

210 self.atoms = atoms.copy()

211

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

213 command = self.command

214

215 # basis path

216 basis_path = self.parameters['basis_path']

217 if basis_path is None:

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

219

220 if basis_path is None:

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

222 ' or the DEMON_BASIS_PATH' +

223 ' environment variable')

224

225 # link restart file

226 value = self.parameters['guess']

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

228 value2 = self.parameters['deMon_restart_path']

229

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

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

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

233 abspath = op.abspath(value2)

234

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

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

237

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

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

240 else:

241 raise RuntimeError(

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

243

244 abspath = op.abspath(basis_path)

245

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

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

248

249 if command is None:

250 raise CalculatorSetupError

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

252

253 try:

254 self.read_results()

255 except Exception: # XXX Which kind of exception?

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

257 lines = fd.readlines()

258 debug_lines = 10

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

260 for line in lines[-20:]:

261 print(line.strip())

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

263 raise RuntimeError

264

265 def set_label(self, label):

266 """Set label directory """

267

268 self.label = label

269

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

271 self.directory = self.label

272 if self.directory == '':

273 self.directory = os.curdir

274

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

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

277 See calculator.py for further details.

278

279 Parameters:

280 atoms : The Atoms object to write.

281 properties : The properties which should be calculated.

282 system_changes : List of properties changed since last run.

283

284 """

285 # Call base calculator.

286 FileIOCalculator.write_input(

287 self,

288 atoms=atoms,

289 properties=properties,

290 system_changes=system_changes)

291

292 if system_changes is None and properties is None:

293 return

294

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

296

297 add_print = ''

298

299 # Start writing the file.

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

301

302 # write keyword argument keywords

303 value = self.parameters['title']

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

305

306 fd.write('#\n')

307

308 value = self.parameters['scftype']

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

310

311 value = self.parameters['xc']

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

313

314 value = self.parameters['guess']

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

316

317 # obtain forces through a single BOMD step

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

319 value = self.parameters['forces']

320 if 'forces' in properties or value:

321

322 self._write_argument('DYNAMICS',

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

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

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

326 add_print = add_print + ' ' + 'MD OPT'

327

328 # if dipole is True, enforce dipole calculation.

329 # Otherwise only if asked for

330 value = self.parameters['dipole']

331 if 'dipole' in properties or value:

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

333

334 # print argument, here other options could change this

335 value = self.parameters['print_out']

336 assert isinstance(value, str)

337 value = value + add_print

338

339 if len(value) != 0:

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

341 fd.write('#\n')

342

343 # write general input arguments

344 self._write_input_arguments(fd)

345

346 fd.write('#\n')

347

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

349 basis = self.parameters['basis']

350 if 'all' not in basis:

351 basis['all'] = 'DZVP'

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

353

354 ecps = self.parameters['ecps']

355 if len(ecps) != 0:

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

357

358 mcps = self.parameters['mcps']

359 if len(mcps) != 0:

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

361

362 auxis = self.parameters['auxis']

363 if len(auxis) != 0:

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

365

366 augment = self.parameters['augment']

367 if len(augment) != 0:

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

369

370 # write geometry

371 self._write_atomic_coordinates(fd, atoms)

372

373 # write xyz file for good measure.

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

375

376 def read(self, restart_path):

377 """Read parameters from directory restart_path."""

378

379 self.set_label(restart_path)

380

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

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

383 .format(restart_path))

384

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

386

387 self.read_results()

388

389 def _write_input_arguments(self, fd):

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

391 input_arguments = self.parameters['input_arguments']

392

393 # Early return

394 if input_arguments is None:

395 return

396

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

398 self._write_argument(key, value, fd)

399

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

401 """Write an argument to file.

402 key : a string coresponding to the input keyword

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

404 f : and open file

405 """

406

407 # for only one argument, write on same line

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

409 line = key.upper()

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

411 fd.write(line)

412 fd.write('\n')

413

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

415 # then the rest on new lines

416 else:

417 line = key

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

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

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

421 fd.write(line)

422 fd.write('\n')

423 else:

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

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

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

427 fd.write(line)

428 fd.write('\n')

429 line = ''

430

431 def _write_atomic_coordinates(self, fd, atoms):

432 """Write atomic coordinates.

433

434 Parameters:

435 - f: An open file object.

436 - atoms: An atoms object.

437 """

438

439 fd.write('#\n')

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

441 fd.write('#\n')

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

443

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

445 xyz = atoms.get_positions()[i]

446 chem_symbol = atoms.get_chemical_symbols()[i]

447 chem_symbol += str(i + 1)

448

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

450 # set nuclear charge to 0

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

452 nuc_charge = str(0)

453 else:

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

455

456 mass = atoms.get_masses()[i]

457

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

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

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

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

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

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

464

465 fd.write(line)

466 fd.write('\n')

467

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

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

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

471

472 Parameters:

473 - f: An open file object.

474 - atoms: An atoms object.

475 - basis: A dictionary specifying the basis set

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

477 """

478

479 # basis for all atoms

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

481

482 if 'all' in basis:

483 default_basis = basis['all']

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

485

486 fd.write(line)

487 fd.write('\n')

488

489 # basis for all atomic species

490 chemical_symbols = atoms.get_chemical_symbols()

491 chemical_symbols_set = set(chemical_symbols)

492

493 for _ in range(chemical_symbols_set.__len__()):

494 symbol = chemical_symbols_set.pop()

495

496 if symbol in basis:

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

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

499 fd.write(line)

500 fd.write('\n')

501

502 # basis for individual atoms

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

504

505 if i in basis:

506 symbol = str(chemical_symbols[i])

507 symbol += str(i + 1)

508

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

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

511 fd.write(line)

512 fd.write('\n')

513

514 # Analysis routines

515 def read_results(self):

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

517 self.read_energy()

518 self.read_forces(self.atoms)

519 self.read_eigenvalues()

520 self.read_dipole()

521 self.read_xray()

522

523 def read_energy(self):

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

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

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

527

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

529

530 for line in lines:

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

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

533 break

534 else:

535 raise RuntimeError

536

537 def read_forces(self, atoms):

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

539

540 natoms = len(atoms)

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

542

543 if op.isfile(filename):

544 with open(filename) as fd:

545 lines = fd.readlines()

546

547 # find line where the orbitals start

548 flag_found = False

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

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

551 start = i + 4

552 flag_found = True

553 break

554

555 if flag_found:

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

557 for i in range(natoms):

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

559 if len(s) > 0]

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

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

562

563 def read_eigenvalues(self):

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

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

566

567 # Read eigenvalues

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

569 lines = fd.readlines()

570

571 # try PRINT MOE

572 eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(

573 lines, 'ALPHA MO ENERGIES', 6)

574 eig_beta, occ_beta = self.read_eigenvalues_one_spin(

575 lines, 'BETA MO ENERGIES', 6)

576

577 # otherwise try PRINT MOS

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

579 eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(

580 lines, 'ALPHA MO COEFFICIENTS', 5)

581 eig_beta, occ_beta = self.read_eigenvalues_one_spin(

582 lines, 'BETA MO COEFFICIENTS', 5)

583

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

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

586

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

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

589 with neigs_per_line eigenvlaues written per line

590 """

591 eig = []

592 occ = []

593

594 skip_line = False

595 more_eigs = False

596

597 # find line where the orbitals start

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

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

600 ii = i

601 more_eigs = True

602 break

603

604 while more_eigs:

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

606 # numbers

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

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

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

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

611 ii = i + 2

612 break

613

614 # read eigenvalues, occupations

615 line = lines[ii].split()

616 if len(line) < neigs_per_line:

617 # last row

618 more_eigs = False

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

620 more_eigs = False

621 skip_line = True

622

623 if not skip_line:

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

625 for l in line:

626 eig.append(float(l))

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

628 for l in line:

629 occ.append(float(l))

630 ii = ii + 3

631

632 return eig, occ

633

634 def read_dipole(self):

635 """Read dipole moment."""

636 dipole = np.zeros(3)

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

638 lines = fd.readlines()

639

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

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

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

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

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

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

646

647 # debye to e*Ang

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

649

650 break

651

652 def read_xray(self):

653 """Read deMon.xry if present."""

654

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

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

657 core_IP = None

658 if op.isfile(filename):

659 with open(filename) as fd:

660 lines = fd.readlines()

661

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

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

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

665

666 try:

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

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

669 except ReadError:

670 pass

671 else:

672 xray_results = {'xray_mode': mode,

673 'ntrans': ntrans,

674 'E_trans': E_trans,

675 'osc_strength': osc_strength, # units?

676 'trans_dip': trans_dip, # units?

677 'core_IP': core_IP}

678

679 self.results['xray'] = xray_results

680

681 def deMon_inp_to_atoms(self, filename):

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

683

684 with open(filename) as fd:

685 lines = fd.readlines()

686

687 # find line where geometry starts

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

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

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

691 coord_units = 'Ang'

692 elif lines.rfind('Bohr'):

693 coord_units = 'Bohr'

694 ii = i

695 break

696

697 chemical_symbols = []

698 xyz = []

699 atomic_numbers = []

700 masses = []

701

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

703 try:

704 line = lines[i].split()

705

706 if len(line) > 0:

707 for symbol in ase.data.chemical_symbols:

708 found = None

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

710 found = symbol

711 break

712

713 if found is not None:

714 chemical_symbols.append(found)

715 else:

716 break

717

718 xyz.append(

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

720

721 if len(line) > 4:

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

723

724 if len(line) > 5:

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

726

727 except Exception: # XXX Which kind of exception?

728 raise RuntimeError

729

730 if coord_units == 'Bohr':

731 xyz *= Bohr

732

733 natoms = len(chemical_symbols)

734

735 # set atoms object

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

737

738 # if atomic numbers were read in, set them

739 if len(atomic_numbers) == natoms:

740 atoms.set_atomic_numbers(atomic_numbers)

741

742 # if masses were read in, set them

743 if len(masses) == natoms:

744 atoms.set_masses(masses)

745

746 return atoms