From 64541f25ca5dc1264b8ede947a1db61f4694a305 Mon Sep 17 00:00:00 2001 From: Edward Linscott Date: Thu, 30 Nov 2023 13:30:15 +0100 Subject: [PATCH] Bugfixes to MT corrections and allowing for calculations without empty states --- src/koopmans/calculators/_projwfc.py | 3 ++- src/koopmans/pseudopotentials.py | 7 ++++--- src/koopmans/utils/_io.py | 2 +- src/koopmans/workflows/_dft.py | 18 ++++++++-------- src/koopmans/workflows/_koopmans_dscf.py | 22 +++++++++++--------- src/koopmans/workflows/_unfold_and_interp.py | 22 +++++++++++++++----- src/koopmans/workflows/_workflow.py | 22 ++++++++++---------- 7 files changed, 56 insertions(+), 40 deletions(-) diff --git a/src/koopmans/calculators/_projwfc.py b/src/koopmans/calculators/_projwfc.py index 1338cdbd1..81181bf07 100644 --- a/src/koopmans/calculators/_projwfc.py +++ b/src/koopmans/calculators/_projwfc.py @@ -75,7 +75,8 @@ def generate_dos(self) -> GridDOSCollection: # The filename does not encode the principal quantum number n. In order to recover this number, we compare # the reported angular momentum quantum number l against the list of expected orbitals, and infer n # assuming only that the file corresponding to nl will come before (n+1)l - orbitals = copy.copy(self._expected_orbitals[atom.symbol]) + label = atom.symbol + str(atom.tag) if atom.tag > 0 else atom.symbol + orbitals = copy.copy(self._expected_orbitals[label]) for filename in filenames: # Infer l from the filename subshell = filename.name[-2] diff --git a/src/koopmans/pseudopotentials.py b/src/koopmans/pseudopotentials.py index 45b1c2e9d..b56046190 100644 --- a/src/koopmans/pseudopotentials.py +++ b/src/koopmans/pseudopotentials.py @@ -160,14 +160,15 @@ def expected_subshells(atoms: Atoms, pseudopotentials: Dict[str, str], expected_orbitals = {} for atom in atoms: - if atom.symbol in expected_orbitals: + label = atom.symbol + str(atom.tag) if atom.tag > 0 else atom.symbol + if label in expected_orbitals: continue - pseudo_file = pseudopotentials[atom.symbol] + pseudo_file = pseudopotentials[label] z_core = atom.number - valence_from_pseudo(pseudo_file, pseudo_dir) if z_core in z_core_to_first_orbital: first_orbital = z_core_to_first_orbital[z_core] else: raise ValueError(f'Failed to identify the subshells of the valence of {pseudo_file}') all_orbitals = list(z_core_to_first_orbital.values()) + ['5d', '6p', '6d'] - expected_orbitals[atom.symbol] = sorted(all_orbitals[all_orbitals.index(first_orbital):]) + expected_orbitals[label] = sorted(all_orbitals[all_orbitals.index(first_orbital):]) return expected_orbitals diff --git a/src/koopmans/utils/_io.py b/src/koopmans/utils/_io.py index 79c88b485..f3ecac6b1 100644 --- a/src/koopmans/utils/_io.py +++ b/src/koopmans/utils/_io.py @@ -52,7 +52,7 @@ def construct_cell_parameters_block(atoms: Atoms) -> Dict[str, Any]: params = dict(**cell_to_parameters(atoms.cell)) else: params = {'vectors': [list(row) for row in atoms.cell[:]], 'units': 'angstrom'} - params['periodic'] = atoms.pbc + params['periodic'] = list([bool(x) for x in atoms.pbc]) return params diff --git a/src/koopmans/workflows/_dft.py b/src/koopmans/workflows/_dft.py index d342284a1..ccd5908a3 100644 --- a/src/koopmans/workflows/_dft.py +++ b/src/koopmans/workflows/_dft.py @@ -36,15 +36,14 @@ def _run(self): if self.parameters.from_scratch: utils.system_call(f'rm -r {calc.parameters.outdir} 2>/dev/null', False) - calc.prefix = 'dft' - calc.directory = '.' - calc.parameters.ndr = 50 - calc.parameters.ndw = 51 - calc.parameters.restart_mode = 'from_scratch' - calc.parameters.do_orbdep = False - calc.parameters.fixed_state = False - calc.parameters.do_outerloop = True - calc.parameters.which_compensation = 'tcc' + # Setting defaults + defaults = {'prefix': 'dft', 'directory': '.', 'ndr': 50, + 'ndw': 51, 'restart_mode': 'from_scratch', 'do_orbdep': False, + 'fixed_state': False, 'do_outerloop': True} + + for key, val in defaults.items(): + if calc.parameters[key] is None: + calc.parameters[key] = val if calc.parameters.maxiter is None: calc.parameters.maxiter = 300 @@ -53,6 +52,7 @@ def _run(self): if calc.parameters.empty_states_maxstep is None: calc.parameters.empty_states_maxstep = 300 + # Running the calculator self.run_calculator(calc, enforce_ss=self.parameters.fix_spin_contamination) return calc diff --git a/src/koopmans/workflows/_koopmans_dscf.py b/src/koopmans/workflows/_koopmans_dscf.py index 213504e08..f5de038ad 100644 --- a/src/koopmans/workflows/_koopmans_dscf.py +++ b/src/koopmans/workflows/_koopmans_dscf.py @@ -39,7 +39,7 @@ def __init__(self, *args, redo_smooth_dft: Optional[bool] = None, # If periodic, convert the kcp calculation into a Γ-only supercell calculation kcp_params = self.calculator_parameters['kcp'] - if all(self.atoms.pbc): + if any(self.atoms.pbc): spins: List[Optional[str]] if self.parameters.spin_polarized: spins = ['up', 'down'] @@ -268,7 +268,7 @@ def _run(self) -> None: self.perform_final_calculations() # Postprocessing - if all(self.atoms.pbc): + if any(self.atoms.pbc): if self.parameters.calculate_bands in [None, True] and self.projections and self.kpoints.path is not None: # Calculate interpolated band structure and DOS with UI from koopmans import workflows @@ -394,23 +394,25 @@ def perform_initialization(self) -> None: self.run_calculator(calc, enforce_ss=False) # Check the consistency between the PW and CP band gaps - pw_calc = [c for c in self.calculations if isinstance( - c, calculators.PWCalculator) and c.parameters.calculation == 'nscf'][-1] - pw_gap = pw_calc.results['lumo_energy'] - pw_calc.results['homo_energy'] - cp_gap = calc.results['lumo_energy'] - calc.results['homo_energy'] - if abs(pw_gap - cp_gap) > 2e-2 * pw_gap: - raise ValueError(f'PW and CP band gaps are not consistent: {pw_gap} {cp_gap}') + if calc.results['lumo_energy']: + pw_calc = [c for c in self.calculations if isinstance( + c, calculators.PWCalculator) and c.parameters.calculation == 'nscf'][-1] + pw_gap = pw_calc.results['lumo_energy'] - pw_calc.results['homo_energy'] + cp_gap = calc.results['lumo_energy'] - calc.results['homo_energy'] + if abs(pw_gap - cp_gap) > 2e-2 * pw_gap: + raise ValueError(f'PW and CP band gaps are not consistent: {pw_gap} {cp_gap}') # The CP restarting from Wannier functions must be already converged Eini = calc.results['convergence']['filled'][0]['Etot'] Efin = calc.results['energy'] if abs(Efin - Eini) > 1e-6 * abs(Efin): - raise ValueError(f'Too much difference between the initial and final CP energies: {Eini} {Efin}') + utils.warn(f'There is a large difference ({Efin - Eini:.5f} Ha) between the initial and final CP energies ({Eini:.3f} and {Efin:.3f} Ha). This should not happen.') # Add to the outdir of dft_init a link to the files containing the Wannier functions dst = Path(f'{calc.parameters.outdir}/{calc.parameters.prefix}_{calc.parameters.ndw}.save/K00001/') for file in ['evc_occupied1.dat', 'evc_occupied2.dat', 'evc0_empty1.dat', 'evc0_empty2.dat']: - utils.symlink(f'{restart_dir}/{file}', dst, force=True) + if (restart_dir / file).exists(): + utils.symlink(f'{restart_dir}/{file}', dst, force=True) elif self.parameters.functional in ['ki', 'pkipz']: calc = self.new_kcp_calculator('dft_init') diff --git a/src/koopmans/workflows/_unfold_and_interp.py b/src/koopmans/workflows/_unfold_and_interp.py index 77745e0cc..5c9108c9f 100644 --- a/src/koopmans/workflows/_unfold_and_interp.py +++ b/src/koopmans/workflows/_unfold_and_interp.py @@ -67,6 +67,7 @@ def _run(self) -> None: else: spins = [None] + energies_to_merge = [] for spin, band_filling in zip(spins, self.bands.filling): # Extract the centers and spreads corresponding to this particular spin centers = np.array([center for c, p in zip(w90_calcs, self.projections) @@ -74,7 +75,17 @@ def _run(self) -> None: spreads = np.array([spread for c, p in zip(w90_calcs, self.projections) for spread in c.results['spreads'] if p.spin == spin]) - for filled, filling in zip([True, False], ['occ', 'emp']): + energies_to_merge.append([]) + + filleds = [True] + if not all(band_filling): + filleds.append(False) + + for filled in filleds: + if filled: + filling = 'occ' + else: + filling = 'emp' calc_presets = filling if spin: calc_presets += '_' + spin @@ -97,17 +108,18 @@ def _run(self) -> None: # Run the calculator self.run_calculator(calc, enforce_ss=False) + # Store the bandstructure energies that we will merge below + energies_to_merge[-1].append(self.calculations[-1].results['band structure'].energies) + # Merge the two calculations to print out the DOS and bands calc = self.new_ui_calculator('merge') # Merge the bands if self.parameters.spin_polarized: - energies = [[c.results['band structure'].energies for c in subset] - for subset in [self.calculations[-4:-2], self.calculations[-2:]]] - reference = np.max([e[0] for e in energies]) + reference = np.max([e[0] for e in energies_to_merge]) energies_np = np.concatenate([np.concatenate(e, axis=2) for e in energies], axis=0) else: - energies = [c.results['band structure'].energies for c in self.calculations[-2:]] + energies = energies_to_merge[0] reference = np.max(energies[0]) energies_np = np.concatenate(energies, axis=2) calc.results['band structure'] = BandStructure(self.kpoints.path, energies_np, reference=reference) diff --git a/src/koopmans/workflows/_workflow.py b/src/koopmans/workflows/_workflow.py index 8591714f0..65a4de092 100644 --- a/src/koopmans/workflows/_workflow.py +++ b/src/koopmans/workflows/_workflow.py @@ -117,7 +117,6 @@ def __init__(self, atoms: Atoms, autogenerate_settings: bool = True, version: Optional[str] = None, **kwargs: Dict[str, Any]): - # Parsing parameters self.parameters = settings.WorkflowSettingsDict(**parameters) for key, value in kwargs.items(): @@ -229,12 +228,13 @@ def __init__(self, atoms: Atoms, # Before saving the calculator_parameters, automatically generate some keywords and perform some sanity checks if self.parameters.task != 'ui' and autogenerate_settings: + # Automatically calculate nelec/nelup/neldw/etc using information contained in the pseudopotential files # and the kcp settings nelec = nelec_from_pseudos(self.atoms, self.pseudopotentials, self.parameters.pseudo_directory) - tot_charge = calculator_parameters['kcp'].get('tot_charge', 0) + tot_charge = kwargs.get('tot_charge', calculator_parameters['kcp'].get('tot_charge', 0)) nelec -= tot_charge - tot_mag = calculator_parameters['kcp'].get('tot_magnetization', nelec % 2) + tot_mag = kwargs.get('tot_magnetization', calculator_parameters['kcp'].get('tot_magnetization', nelec % 2)) nelup = int(nelec / 2 + tot_mag / 2) neldw = int(nelec / 2 - tot_mag / 2) @@ -435,11 +435,11 @@ def _run_sanity_checks(self): self.parameters.calculate_alpha = False # Checking periodic image correction schemes - if not self.parameters.calculate_alpha: - # If we are not calculating alpha, we do not consider charged systems and therefore we don't need image - # corrections, so we skip the following checks - pass - elif all(self.atoms.pbc): + # if not self.parameters.calculate_alpha: + # # If we are not calculating alpha, we do not consider charged systems and therefore we don't need image + # # corrections, so we skip the following checks + # pass + if all(self.atoms.pbc): if self.parameters.method == 'dfpt': # For DPFT, we use gb_correction if self.parameters.gb_correction is None: @@ -680,11 +680,11 @@ def new_calculator(self, # Add Martyna-Tuckerman settings for kcp calculators if calc.parameters.is_valid('which_compensation'): if self.parameters.mt_correction: - if all(self.atoms.pbc): + if not any(self.atoms.pbc): calc.parameters.which_compensation = 'tcc' - elif self.atoms.pbc == [False, False, True]: + elif all(self.atoms.pbc == [False, False, True]): calc.parameters.which_compensation = 'tcc1d' - elif self.atoms.pbc == [True, True, False]: + elif all(self.atoms.pbc == [True, True, False]): calc.parameters.which_compensation = 'tcc2d' else: raise ValueError('Martyna-Tuckerman correction not implemented for this geometry; for 1D and 2D '