Coverage for /builds/debichem-team/python-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 (

16 CalculatorSetupError,

17 FileIOCalculator,

18 Parameters,

19 ReadError,

20 all_changes,

21 equal,

22)

23from ase.units import Bohr, Hartree

25from .demon_io import parse_xray

27m_e_to_amu = 1822.88839

30class Parameters_deMon(Parameters):

31 """Parameters class for the calculator.

32 Documented in Base_deMon.__init__

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

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

36 input_arguments.

38 """

40 def __init__(

41 self,

42 label='rundir',

43 atoms=None,

44 restart=None,

45 basis_path=None,

46 ignore_bad_restart_file=FileIOCalculator._deprecated,

47 deMon_restart_path='.',

48 title='deMon input file',

49 scftype='RKS',

50 forces=False,

51 dipole=False,

52 xc='VWN',

53 guess='TB',

54 print_out='MOE',

55 basis={},

56 ecps={},

57 mcps={},

58 auxis={},

59 augment={},

60 input_arguments=None):

61 kwargs = locals()

62 kwargs.pop('self')

63 Parameters.__init__(self, **kwargs)

66class Demon(FileIOCalculator):

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

69 implemented_properties = [

70 'energy',

71 'forces',

72 'dipole',

73 'eigenvalues']

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

76 """ASE interface to the deMon code.

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

80 The DEMON_COMMAND environment variable must be set to run the

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

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

84 Parameters:

86 label : str

87 relative path to the run directory

88 atoms : Atoms object

89 the atoms object

90 command : str

91 Command to run deMon. If not present the environment

92 variable DEMON_COMMAND will be used

93 restart : str

94 Relative path to ASE restart directory for parameters and

95 atoms object and results

96 basis_path : str

97 Relative path to the directory containing

98 BASIS, AUXIS, ECPS, MCPS and AUGMENT

99 ignore_bad_restart_file : bool

100 Ignore broken or missing ASE restart files

101 By default, it is an error if the restart

102 file is missing or broken.

103 deMon_restart_path : str

104 Relative path to the deMon restart dir

105 title : str

106 Title in the deMon input file.

107 scftype : str

108 Type of scf

109 forces : bool

110 If True a force calculation will be enforced.

111 dipole : bool

112 If True a dipole calculation will be enforced

113 xc : str

114 xc-functional

115 guess : str

116 guess for initial density and wave functions

117 print_out : str | list

118 Options for the printing in deMon

119 basis : dict

120 Definition of basis sets.

121 ecps : dict

122 Definition of ECPs

123 mcps : dict

124 Definition of MCPs

125 auxis : dict

126 Definition of AUXIS

127 augment : dict

128 Definition of AUGMENT

129 input_arguments : dict

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

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

132 on the same line as the keyword),

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

134 line, the others on following lines.)

135

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

137 the directory ase/test/demon

138

139 """

140

141 parameters = Parameters_deMon(**kwargs)

142

143 # Setup the run command

144 if command is None:

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

146

147 FileIOCalculator.__init__(

148 self,

149 command=command,

150 **parameters)

151

152 def __getitem__(self, key):

153 """Convenience method to retrieve a parameter as

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

155

156 Parameters:

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

158 """

159 return self.parameters[key]

160

161 def set(self, **kwargs):

162 """Set all parameters.

163

164 Parameters:

165 kwargs : Dictionary containing the keywords for deMon

166 """

167 # Put in the default arguments.

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

169

170 if 'parameters' in kwargs:

171 filename = kwargs.pop('parameters')

172 parameters = Parameters.read(filename)

173 parameters.update(kwargs)

174 kwargs = parameters

175

176 changed_parameters = {}

177

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

179 oldvalue = self.parameters.get(key)

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

181 changed_parameters[key] = value

182 self.parameters[key] = value

183

184 return changed_parameters

185

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

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

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

189

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

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

192 todir + '/' + filename)

193 else:

194 raise RuntimeError(

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

196

197 def calculate(self,

198 atoms=None,

199 properties=['energy'],

200 system_changes=all_changes):

201 """Capture the RuntimeError from FileIOCalculator.calculate

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

203

204 See base FileIocalculator for documentation.

205 """

206

207 if atoms is not None:

208 self.atoms = atoms.copy()

209

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

211 command = self.command

212

213 # basis path

214 basis_path = self.parameters['basis_path']

215 if basis_path is None:

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

217

218 if basis_path is None:

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

220 ' or the DEMON_BASIS_PATH' +

221 ' environment variable')

222

223 # link restart file

224 value = self.parameters['guess']

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

226 value2 = self.parameters['deMon_restart_path']

227

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

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

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

231 abspath = op.abspath(value2)

232

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

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

235

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

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

238 else:

239 raise RuntimeError(

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

241

242 abspath = op.abspath(basis_path)

243

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

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

246

247 if command is None:

248 raise CalculatorSetupError

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

250

251 try:

252 self.read_results()

253 except Exception: # XXX Which kind of exception?

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

255 lines = fd.readlines()

256 debug_lines = 10

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

258 for line in lines[-20:]:

259 print(line.strip())

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

261 raise RuntimeError

262

263 def set_label(self, label):

264 """Set label directory """

265

266 self.label = label

267

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

269 self.directory = self.label

270 if self.directory == '':

271 self.directory = os.curdir

272

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

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

275 See calculator.py for further details.

276

277 Parameters:

278 atoms : The Atoms object to write.

279 properties : The properties which should be calculated.

280 system_changes : List of properties changed since last run.

281

282 """

283 # Call base calculator.

284 FileIOCalculator.write_input(

285 self,

286 atoms=atoms,

287 properties=properties,

288 system_changes=system_changes)

289

290 if system_changes is None and properties is None:

291 return

292

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

294

295 add_print = ''

296

297 # Start writing the file.

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

299

300 # write keyword argument keywords

301 value = self.parameters['title']

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

303

304 fd.write('#\n')

305

306 value = self.parameters['scftype']

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

308

309 value = self.parameters['xc']

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

311

312 value = self.parameters['guess']

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

314

315 # obtain forces through a single BOMD step

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

317 value = self.parameters['forces']

318 if 'forces' in properties or value:

319

320 self._write_argument('DYNAMICS',

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

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

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

324 add_print = add_print + ' ' + 'MD OPT'

325

326 # if dipole is True, enforce dipole calculation.

327 # Otherwise only if asked for

328 value = self.parameters['dipole']

329 if 'dipole' in properties or value:

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

331

332 # print argument, here other options could change this

333 value = self.parameters['print_out']

334 assert isinstance(value, str)

335 value = value + add_print

336

337 if len(value) != 0:

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

339 fd.write('#\n')

340

341 # write general input arguments

342 self._write_input_arguments(fd)

343

344 fd.write('#\n')

345

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

347 basis = self.parameters['basis']

348 if 'all' not in basis:

349 basis['all'] = 'DZVP'

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

351

352 ecps = self.parameters['ecps']

353 if len(ecps) != 0:

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

355

356 mcps = self.parameters['mcps']

357 if len(mcps) != 0:

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

359

360 auxis = self.parameters['auxis']

361 if len(auxis) != 0:

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

363

364 augment = self.parameters['augment']

365 if len(augment) != 0:

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

367

368 # write geometry

369 self._write_atomic_coordinates(fd, atoms)

370

371 # write xyz file for good measure.

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

373

374 def read(self, restart_path):

375 """Read parameters from directory restart_path."""

376

377 self.set_label(restart_path)

378

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

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

381 .format(restart_path))

382

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

384

385 self.read_results()

386

387 def _write_input_arguments(self, fd):

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

389 input_arguments = self.parameters['input_arguments']

390

391 # Early return

392 if input_arguments is None:

393 return

394

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

396 self._write_argument(key, value, fd)

397

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

399 """Write an argument to file.

400 key : a string coresponding to the input keyword

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

402 f : and open file

403 """

404

405 # for only one argument, write on same line

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

407 line = key.upper()

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

409 fd.write(line)

410 fd.write('\n')

411

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

413 # then the rest on new lines

414 else:

415 line = key

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

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

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

419 fd.write(line)

420 fd.write('\n')

421 else:

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

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

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

425 fd.write(line)

426 fd.write('\n')

427 line = ''

428

429 def _write_atomic_coordinates(self, fd, atoms):

430 """Write atomic coordinates.

431

432 Parameters:

433 - f: An open file object.

434 - atoms: An atoms object.

435 """

436

437 fd.write('#\n')

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

439 fd.write('#\n')

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

441

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

443 xyz = atoms.get_positions()[i]

444 chem_symbol = atoms.get_chemical_symbols()[i]

445 chem_symbol += str(i + 1)

446

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

448 # set nuclear charge to 0

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

450 nuc_charge = str(0)

451 else:

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

453

454 mass = atoms.get_masses()[i]

455

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

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

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

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

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

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

462

463 fd.write(line)

464 fd.write('\n')

465

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

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

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

469

470 Parameters:

471 - f: An open file object.

472 - atoms: An atoms object.

473 - basis: A dictionary specifying the basis set

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

475 """

476

477 # basis for all atoms

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

479

480 if 'all' in basis:

481 default_basis = basis['all']

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

483

484 fd.write(line)

485 fd.write('\n')

486

487 # basis for all atomic species

488 chemical_symbols = atoms.get_chemical_symbols()

489 chemical_symbols_set = set(chemical_symbols)

490

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

492 symbol = chemical_symbols_set.pop()

493

494 if symbol in basis:

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

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

497 fd.write(line)

498 fd.write('\n')

499

500 # basis for individual atoms

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

502

503 if i in basis:

504 symbol = str(chemical_symbols[i])

505 symbol += str(i + 1)

506

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

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

509 fd.write(line)

510 fd.write('\n')

511

512 # Analysis routines

513 def read_results(self):

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

515 self.read_energy()

516 self.read_forces(self.atoms)

517 self.read_eigenvalues()

518 self.read_dipole()

519 self.read_xray()

520

521 def read_energy(self):

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

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

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

525

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

527

528 for line in lines:

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

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

531 break

532 else:

533 raise RuntimeError

534

535 def read_forces(self, atoms):

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

537

538 natoms = len(atoms)

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

540

541 if op.isfile(filename):

542 with open(filename) as fd:

543 lines = fd.readlines()

544

545 # find line where the orbitals start

546 flag_found = False

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

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

549 start = i + 4

550 flag_found = True

551 break

552

553 if flag_found:

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

555 for i in range(natoms):

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

557 if len(s) > 0]

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

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

560

561 def read_eigenvalues(self):

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

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

564

565 # Read eigenvalues

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

567 lines = fd.readlines()

568

569 # try PRINT MOE

570 eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(

571 lines, 'ALPHA MO ENERGIES', 6)

572 eig_beta, occ_beta = self.read_eigenvalues_one_spin(

573 lines, 'BETA MO ENERGIES', 6)

574

575 # otherwise try PRINT MOS

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

577 eig_alpha, occ_alpha = self.read_eigenvalues_one_spin(

578 lines, 'ALPHA MO COEFFICIENTS', 5)

579 eig_beta, occ_beta = self.read_eigenvalues_one_spin(

580 lines, 'BETA MO COEFFICIENTS', 5)

581

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

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

584

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

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

587 with neigs_per_line eigenvlaues written per line

588 """

589 eig = []

590 occ = []

591

592 skip_line = False

593 more_eigs = False

594

595 # find line where the orbitals start

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

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

598 ii = i

599 more_eigs = True

600 break

601

602 while more_eigs:

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

604 # numbers

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

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

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

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

609 ii = i + 2

610 break

611

612 # read eigenvalues, occupations

613 line = lines[ii].split()

614 if len(line) < neigs_per_line:

615 # last row

616 more_eigs = False

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

618 more_eigs = False

619 skip_line = True

620

621 if not skip_line:

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

623 for l in line:

624 eig.append(float(l))

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

626 for l in line:

627 occ.append(float(l))

628 ii = ii + 3

629

630 return eig, occ

631

632 def read_dipole(self):

633 """Read dipole moment."""

634 dipole = np.zeros(3)

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

636 lines = fd.readlines()

637

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

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

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

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

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

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

644

645 # debye to e*Ang

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

647

648 break

649

650 def read_xray(self):

651 """Read deMon.xry if present."""

652

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

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

655 core_IP = None

656 if op.isfile(filename):

657 with open(filename) as fd:

658 lines = fd.readlines()

659

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

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

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

663

664 try:

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

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

667 except ReadError:

668 pass

669 else:

670 xray_results = {'xray_mode': mode,

671 'ntrans': ntrans,

672 'E_trans': E_trans,

673 'osc_strength': osc_strength, # units?

674 'trans_dip': trans_dip, # units?

675 'core_IP': core_IP}

676

677 self.results['xray'] = xray_results

678

679 def deMon_inp_to_atoms(self, filename):

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

681

682 with open(filename) as fd:

683 lines = fd.readlines()

684

685 # find line where geometry starts

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

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

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

689 coord_units = 'Ang'

690 elif lines.rfind('Bohr'):

691 coord_units = 'Bohr'

692 ii = i

693 break

694

695 chemical_symbols = []

696 xyz = []

697 atomic_numbers = []

698 masses = []

699

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

701 try:

702 line = lines[i].split()

703

704 if len(line) > 0:

705 for symbol in ase.data.chemical_symbols:

706 found = None

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

708 found = symbol

709 break

710

711 if found is not None:

712 chemical_symbols.append(found)

713 else:

714 break

715

716 xyz.append(

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

718

719 if len(line) > 4:

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

721

722 if len(line) > 5:

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

724

725 except Exception: # XXX Which kind of exception?

726 raise RuntimeError

727

728 if coord_units == 'Bohr':

729 xyz *= Bohr

730

731 natoms = len(chemical_symbols)

732

733 # set atoms object

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

735

736 # if atomic numbers were read in, set them

737 if len(atomic_numbers) == natoms:

738 atoms.set_atomic_numbers(atomic_numbers)

739

740 # if masses were read in, set them

741 if len(masses) == natoms:

742 atoms.set_masses(masses)

743

744 return atoms