• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

pyiron / pyiron_atomistics / 10502979427

22 Aug 2024 06:26AM UTC coverage: 70.935% (+0.005%) from 70.93%
10502979427

Pull #1419

github

web-flow
Merge 8e363d3eb into 29fb911e0
Pull Request #1419: Always overwrite eigenvalues

1 of 1 new or added line in 1 file covered. (100.0%)

2 existing lines in 1 file now uncovered.

10680 of 15056 relevant lines covered (70.94%)

0.71 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

80.13
/pyiron_atomistics/sphinx/base.py
1
# coding: utf-8
2
# Copyright (c) Max-Planck-Institut für Eisenforschung GmbH -Computational Materials Design (CM) Department
3
# Distributed under the terms of "New BSD License", see the LICENSE file.
4

5
from __future__ import division, print_function
1✔
6

7
import os
1✔
8
import posixpath
1✔
9
import stat
1✔
10
import subprocess
1✔
11
import tarfile
1✔
12
import warnings
1✔
13
from shutil import move as movefile
1✔
14
from subprocess import PIPE
1✔
15
from tempfile import TemporaryDirectory
1✔
16

17
import numpy as np
1✔
18
import scipy.constants
1✔
19
import spglib
1✔
20
from pyiron_base import DataContainer, job_status_successful_lst, state
1✔
21
from pyiron_snippets.deprecate import deprecate
1✔
22

23
from pyiron_atomistics.dft.job.generic import GenericDFTJob
1✔
24
from pyiron_atomistics.dft.waves.electronic import ElectronicStructure
1✔
25
from pyiron_atomistics.sphinx.input_writer import (
1✔
26
    Group,
27
    copy_potentials,
28
    get_structure_group,
29
    write_spin_constraints,
30
)
31
from pyiron_atomistics.sphinx.output_parser import (
1✔
32
    SphinxLogParser,
33
    collect_energy_dat,
34
    collect_energy_struct,
35
    collect_eps_dat,
36
    collect_relaxed_hist,
37
    collect_residue_dat,
38
    collect_spins_dat,
39
)
40
from pyiron_atomistics.sphinx.potential import SphinxJTHPotentialFile
1✔
41
from pyiron_atomistics.sphinx.structure import read_atoms
1✔
42
from pyiron_atomistics.sphinx.util import sxversions
1✔
43
from pyiron_atomistics.sphinx.volumetric_data import SphinxVolumetricData
1✔
44
from pyiron_atomistics.vasp.potential import (
1✔
45
    VaspPotentialFile,
46
    VaspPotentialSetter,
47
    strip_xc_from_potential_name,
48
)
49

50
__author__ = "Osamu Waseda, Jan Janssen"
1✔
51
__copyright__ = (
1✔
52
    "Copyright 2021, Max-Planck-Institut für Eisenforschung GmbH - "
53
    "Computational Materials Design (CM) Department"
54
)
55
__version__ = "1.0"
1✔
56
__maintainer__ = "Jan Janssen"
1✔
57
__email__ = "janssen@mpie.de"
1✔
58
__status__ = "development"
1✔
59
__date__ = "Sep 1, 2017"
1✔
60

61
BOHR_TO_ANGSTROM = (
1✔
62
    scipy.constants.physical_constants["Bohr radius"][0] / scipy.constants.angstrom
63
)
64
HARTREE_TO_EV = scipy.constants.physical_constants["Hartree energy in eV"][0]
1✔
65
RYDBERG_TO_EV = HARTREE_TO_EV / 2
1✔
66
HARTREE_OVER_BOHR_TO_EV_OVER_ANGSTROM = HARTREE_TO_EV / BOHR_TO_ANGSTROM
1✔
67

68

69
class SphinxBase(GenericDFTJob):
1✔
70
    """
71
    Class to setup and run SPHInX simulations.
72

73
    Inherits pyiron_atomistics.atomistics.job.generic.GenericJob. The functions in
74
    these modules are written such that the function names and attributes
75
    are very pyiron-generic (get_structure(), molecular_dynamics(),
76
    version) but internally handle SPHInX specific input and output.
77

78
    Alternatively, because SPHInX inputs are built on a group-based
79
    format, users have the option to set specific groups and parameters
80
    directly, e.g.
81

82
    ```python
83
    # Modify/add a new parameter via
84
    job.input.sphinx.PAWHamiltonian.nEmptyStates = 15
85
    job.input.sphinx.PAWHamiltonian.dipoleCorrection = True
86
    # or
87
    job.input.sphinx.PAWHamiltonian.set("nEmptyStates", 15)
88
    job.input.sphinx.PAWHamiltonian.set("dipoleCorrection", True)
89
    # Modify/add a sub-group via
90
    job.input.sphinx.initialGuess.rho.charged = {"charge": 2, "z": 25}
91
    # or
92
    job.input.sphinx.initialGuess.rho.set("charged", {"charge": 2, "z": 25})
93
    ```
94

95
    Args:
96
        project: Project object (defines path where job will be
97
                 created and stored)
98
        job_name (str): name of the job (must be unique within
99
                        this project path)
100
    """
101

102
    """Version of the data format in hdf5"""
1✔
103
    __hdf_version__ = "0.1.0"
1✔
104

105
    def __init__(self, project, job_name):
1✔
106
        super(SphinxBase, self).__init__(project, job_name)
1✔
107

108
        # keeps both the generic parameters as well as the sphinx specific
109
        # input groups
110
        self.input = Group(table_name="parameters", lazy=True)
1✔
111
        self.load_default_input()
1✔
112
        self.output = Output(job=self)
1✔
113
        self._potential = VaspPotentialSetter([])
1✔
114
        if self.check_vasp_potentials():
1✔
115
            self.input["VaspPot"] = True  # use VASP potentials if available
×
116
        self._generic_input["restart_for_band_structure"] = False
1✔
117
        self._generic_input["path_name"] = None
1✔
118
        self._generic_input["n_path"] = None
1✔
119
        self._generic_input["fix_spin_constraint"] = False
1✔
120

121
    def update_sphinx(self):
1✔
122
        if self.output.old_version:
×
123
            _update_datacontainer(self)
×
124

125
    def __getitem__(self, item):
1✔
126
        if not isinstance(item, str):
1✔
127
            return super().__getitem__(item)
×
128
        result = None
1✔
129
        if item[-1] == "/":
1✔
130
            item = item[:-1]
×
131
        for tag in item.split("/"):
1✔
132
            try:  # horrible workaround to be removed when hdf output becomes consistent
1✔
133
                if result is None:
1✔
134
                    result = super().__getitem__(tag)
1✔
135
                else:
136
                    result = result[tag]
1✔
137
                if hasattr(result, "list_nodes") and "TYPE" in result.list_nodes():
1✔
138
                    result = result.to_object()
1✔
139
            except (ValueError, KeyError):
×
140
                return None
×
141
        return result
1✔
142

143
    @property
1✔
144
    def structure(self):
1✔
145
        """
146

147
        Returns:
148

149
        """
150
        return GenericDFTJob.structure.fget(self)
×
151

152
    @structure.setter
1✔
153
    def structure(self, structure):
1✔
154
        """
155

156
        Args:
157
            structure:
158

159
        Returns:
160

161
        """
162
        GenericDFTJob.structure.fset(self, structure)
×
163
        if structure is not None:
×
164
            self._potential = VaspPotentialSetter(
×
165
                element_lst=structure.get_species_symbols().tolist()
166
            )
167

168
    @property
1✔
169
    def id_pyi_to_spx(self):
1✔
170
        if self.structure is None:
1✔
171
            raise ValueError("Structure not set")
×
172
        # Translate the chemical symbols into indices
173
        indices = np.unique(self.structure.get_chemical_symbols(), return_inverse=True)[
1✔
174
            1
175
        ]
176
        # Add small ramp to ensure order uniqueness
177
        indices = indices + np.arange(len(indices)) / len(indices)
1✔
178
        return np.argsort(indices)
1✔
179

180
    @property
1✔
181
    def id_spx_to_pyi(self):
1✔
182
        if self.structure is None:
1✔
183
            raise ValueError("Structure not set")
×
184
        return np.argsort(self.id_pyi_to_spx)
1✔
185

186
    @property
1✔
187
    def plane_wave_cutoff(self):
1✔
188
        if "eCut" in self.input.sphinx.basis.keys():
1✔
189
            return self.input.sphinx.basis["eCut"] * RYDBERG_TO_EV
1✔
190
        else:
191
            return self.input["EnCut"]
×
192

193
    @property
1✔
194
    def fix_spin_constraint(self):
1✔
195
        return self._generic_input["fix_spin_constraint"]
1✔
196

197
    @fix_spin_constraint.setter
1✔
198
    def fix_spin_constraint(self, boolean):
1✔
199
        if not isinstance(boolean, bool):
1✔
200
            raise ValueError("fix_spin_constraint has to be a boolean")
1✔
201
        self._generic_input["fix_spin_constraint"] = boolean
1✔
202
        self.structure.set_array(
1✔
203
            "spin_constraint", np.array(len(self.structure) * [boolean])
204
        )
205

206
    @plane_wave_cutoff.setter
1✔
207
    def plane_wave_cutoff(self, val):
1✔
208
        """
209
        Function to setup the energy cut-off for the SPHInX job in eV.
210

211
        Args:
212
            val (int): energy cut-off in eV
213
        """
214
        if val <= 0:
1✔
215
            raise ValueError("Cutoff radius value not valid")
1✔
216
        elif val < 50:
1✔
217
            warnings.warn(
1✔
218
                "The given cutoff is either very small (probably "
219
                + "too small) or was accidentally given in Ry. "
220
                + "Please make sure it is in eV (1eV = 13.606 Ry)."
221
            )
222
        self.input["EnCut"] = val
1✔
223
        self.input.sphinx.basis.eCut = self.input["EnCut"] / RYDBERG_TO_EV
1✔
224

225
    @property
1✔
226
    def exchange_correlation_functional(self):
1✔
227
        return self.input["Xcorr"]
1✔
228

229
    @exchange_correlation_functional.setter
1✔
230
    def exchange_correlation_functional(self, val):
1✔
231
        """
232
        Args:
233
            val:
234

235
        Returns:
236
        """
237
        if val.upper() in ["PBE", "LDA"]:
1✔
238
            self.input["Xcorr"] = val.upper()
1✔
239
        else:
240
            warnings.warn(
1✔
241
                "Exchange correlation function not recognized (\
242
                    recommended: PBE or LDA)",
243
                SyntaxWarning,
244
            )
245
            self.input["Xcorr"] = val
1✔
246
        if "xc" in self.input.sphinx.PAWHamiltonian.keys():
1✔
247
            self.input.sphinx.PAWHamiltonian.xc = self.input["Xcorr"]
1✔
248

249
    @property
1✔
250
    def potential_view(self):
1✔
251
        if self.structure is None:
1✔
252
            raise ValueError("Can't list potentials unless a structure is set")
×
253
        else:
254
            if self.input["VaspPot"]:
1✔
255
                potentials = VaspPotentialFile(xc=self.input["Xcorr"])
×
256
            else:
257
                potentials = SphinxJTHPotentialFile(xc=self.input["Xcorr"])
1✔
258
            df = potentials.find(self.structure.get_species_symbols().tolist())
1✔
259
            if len(df) > 0:
1✔
260
                df["Name"] = [
1✔
261
                    strip_xc_from_potential_name(n) for n in df["Name"].values
262
                ]
263
            return df
1✔
264

265
    @property
1✔
266
    def potential_list(self):
1✔
267
        return list(self.potential_view["Name"].values)
1✔
268

269
    @property
1✔
270
    def potential(self):
1✔
271
        return self._potential
1✔
272

273
    def get_kpoints(self):
1✔
274
        return self.input.KpointFolding
×
275

276
    def get_version_float(self):
1✔
277
        version_str = self.executable.version.split("_")[0]
1✔
278
        version_float = float(version_str.split(".")[0])
1✔
279
        if len(version_str.split(".")) > 1:
1✔
280
            version_float += float("0." + "".join(version_str.split(".")[1:]))
1✔
281
        return version_float
1✔
282

283
    def set_input_to_read_only(self):
1✔
284
        """
285
        This function enforces read-only mode for the input classes,
286
        but it has to be implemented in the individual classes.
287
        """
288
        super(SphinxBase, self).set_input_to_read_only()
×
289
        self.input.read_only = True
×
290

291
    def get_scf_group(
1✔
292
        self, maxSteps=None, keepRhoFixed=False, dEnergy=None, algorithm="blockCCG"
293
    ):
294
        """
295
        SCF group setting for SPHInX
296
        for all args refer to calc_static or calc_minimize
297
        """
298

299
        scf_group = Group()
1✔
300
        if algorithm.upper() == "CCG":
1✔
301
            algorithm = "CCG"
×
302
        elif algorithm.upper() != "BLOCKCCG":
1✔
303
            warnings.warn(
1✔
304
                "Algorithm not recognized -> setting to blockCCG. \
305
                    Alternatively, choose algorithm=CCG",
306
                SyntaxWarning,
307
            )
308
            algorithm = "blockCCG"
1✔
309

310
        if keepRhoFixed:
1✔
311
            scf_group["keepRhoFixed"] = True
1✔
312
        else:
313
            scf_group["rhoMixing"] = str(self.input["rhoMixing"])
1✔
314
            scf_group["spinMixing"] = str(self.input["spinMixing"])
1✔
315
            if "nPulaySteps" in self.input:
1✔
316
                scf_group["nPulaySteps"] = str(self.input["nPulaySteps"])
1✔
317
        if dEnergy is None:
1✔
318
            scf_group["dEnergy"] = self.input["Ediff"] / HARTREE_TO_EV
1✔
319
        else:
320
            scf_group["dEnergy"] = str(dEnergy)
1✔
321
        if maxSteps is None:
1✔
322
            scf_group["maxSteps"] = str(self.input["Estep"])
1✔
323
        else:
324
            scf_group["maxSteps"] = str(maxSteps)
1✔
325
        if "preconditioner" in self.input and self.input["preconditioner"] != "KERKER":
1✔
326
            scf_group.create_group("preconditioner")["type"] = self.input[
1✔
327
                "preconditioner"
328
            ]
329
        else:
330
            scf_group.create_group("preconditioner")["type"] = "KERKER"
1✔
331
            scf_group.preconditioner["scaling"] = self.input["rhoResidualScaling"]
1✔
332
            scf_group.preconditioner["spinScaling"] = self.input["spinResidualScaling"]
1✔
333
        scf_group.create_group(algorithm)
1✔
334
        if "maxStepsCCG" in self.input:
1✔
335
            scf_group[algorithm]["maxStepsCCG"] = self.input["maxStepsCCG"]
1✔
336
        if "blockSize" in self.input and algorithm == "blockCCG":
1✔
337
            scf_group[algorithm]["blockSize"] = self.input["blockSize"]
1✔
338
        if "nSloppy" in self.input and algorithm == "blockCCG":
1✔
339
            scf_group[algorithm]["nSloppy"] = self.input["nSloppy"]
1✔
340
        if self.input["WriteWaves"] is False:
1✔
341
            scf_group["noWavesStorage"] = True
1✔
342
        return scf_group
1✔
343

344
    def get_structure_group(self, keep_angstrom=False):
1✔
345
        """
346
        create a SPHInX Group object based on self.structure
347

348
        Args:
349
            keep_angstrom (bool): Store distances in Angstroms or Bohr
350
        """
351
        return get_structure_group(
1✔
352
            positions=self.structure.positions,
353
            cell=self.structure.cell,
354
            elements=self.structure.get_chemical_symbols(),
355
            movable=self.structure.arrays.get("selective_dynamics", None),
356
            labels=self.structure.get_initial_magnetic_moments(),
357
            use_symmetry=self.fix_symmetry,
358
            keep_angstrom=keep_angstrom,
359
        )
360

361
    def load_default_input(self):
1✔
362
        """
363
        Set defaults for generic parameters and create SPHInX input groups.
364
        """
365

366
        sph = self.input.create_group("sphinx")
1✔
367
        sph.create_group("pawPot")
1✔
368
        sph.create_group("structure")
1✔
369
        sph.create_group("basis")
1✔
370
        sph.create_group("PAWHamiltonian")
1✔
371
        sph.create_group("initialGuess")
1✔
372
        sph.create_group("main")
1✔
373

374
        self.input.EnCut = 340
1✔
375
        self.input.KpointCoords = [0.5, 0.5, 0.5]
1✔
376
        self.input.KpointFolding = [4, 4, 4]
1✔
377
        self.input.EmptyStates = "auto"
1✔
378
        self.input.MethfesselPaxton = 1
1✔
379
        self.input.Sigma = 0.2
1✔
380
        self.input.Xcorr = "PBE"
1✔
381
        self.input.VaspPot = False
1✔
382
        self.input.Estep = 100
1✔
383
        self.input.Ediff = 1.0e-4
1✔
384
        self.input.WriteWaves = True
1✔
385
        self.input.KJxc = False
1✔
386
        self.input.SaveMemory = True
1✔
387
        self.input.rhoMixing = 1.0
1✔
388
        self.input.spinMixing = 1.0
1✔
389
        self.input.rhoResidualScaling = 1.0
1✔
390
        self.input.spinResidualScaling = 1.0
1✔
391
        self.input.CheckOverlap = True
1✔
392
        self.input.THREADS = 1
1✔
393
        self.input.use_on_the_fly_cg_optimization = True
1✔
394

395
    def load_structure_group(self, keep_angstrom=False):
1✔
396
        """
397
        Build + load the structure group based on self.structure
398

399
        Args:
400
            keep_angstrom (bool): Store distances in Angstroms or Bohr
401
        """
402
        self.input.sphinx.structure = self.get_structure_group(
1✔
403
            keep_angstrom=keep_angstrom
404
        )
405

406
    def load_species_group(self, check_overlap=True, potformat="VASP"):
1✔
407
        """
408
        Build the species Group object based on self.structure
409

410
        Args:
411
            check_overlap (bool): Whether to check overlap
412
                (see set_check_overlap)
413
            potformat (str): type of pseudopotentials that will be
414
                read. Options are JTH or VASP.
415
        """
416

417
        self.input.sphinx.pawPot = Group({"species": []})
1✔
418
        for species_obj in self.structure.get_species_objects():
1✔
419
            if species_obj.Parent is not None:
1✔
420
                elem = species_obj.Parent
×
421
            else:
422
                elem = species_obj.Abbreviation
1✔
423
            if potformat == "JTH":
1✔
424
                self.input.sphinx.pawPot["species"].append(
1✔
425
                    Group(
426
                        {
427
                            "name": '"' + elem + '"',
428
                            "potType": '"AtomPAW"',
429
                            "element": '"' + elem + '"',
430
                            "potential": f'"{elem}_GGA.atomicdata"',
431
                        }
432
                    )
433
                )
434
            elif potformat == "VASP":
×
435
                self.input.sphinx.pawPot["species"].append(
×
436
                    Group(
437
                        {
438
                            "name": '"' + elem + '"',
439
                            "potType": '"VASP"',
440
                            "element": '"' + elem + '"',
441
                            "potential": '"' + elem + "_POTCAR" + '"',
442
                        }
443
                    )
444
                )
445
            else:
446
                raise ValueError("Potential must be JTH or VASP")
×
447
        if not check_overlap:
1✔
448
            self.input.sphinx.pawPot["species"][-1]["checkOverlap"] = "false"
×
449
        if self.input["KJxc"]:
1✔
450
            self.input.sphinx.pawPot["kjxc"] = True
×
451

452
    def load_main_group(self):
1✔
453
        """
454
        Load the main Group.
455

456
        The group is populated based on the type of calculation and settings in
457
        the self.input.
458
        """
459

460
        if (
1✔
461
            len(self.restart_file_list) != 0
462
            and not self._generic_input["restart_for_band_structure"]
463
        ):
464
            self.input.sphinx.main.get("scfDiag", create=True).append(
1✔
465
                self.get_scf_group(maxSteps=10, keepRhoFixed=True, dEnergy=1.0e-4)
466
            )
467
        if "Istep" in self.input:
1✔
468
            optimizer = "linQN"
1✔
469
            if self.input.use_on_the_fly_cg_optimization:
1✔
470
                optimizer = "ricQN"
1✔
471
            self.input.sphinx.main[optimizer] = Group(table_name="input")
1✔
472
            self.input.sphinx.main[optimizer]["maxSteps"] = str(self.input["Istep"])
1✔
473
            self.input.sphinx.main[optimizer]["maxStepLength"] = str(
1✔
474
                0.1 / BOHR_TO_ANGSTROM
475
            )
476
            if "dE" in self.input:
1✔
477
                self.input.sphinx.main[optimizer]["dEnergy"] = str(
1✔
478
                    self.input["dE"] / HARTREE_TO_EV
479
                )
480
            if "dF" in self.input:
1✔
481
                self.input.sphinx.main[optimizer]["dF"] = str(
×
482
                    self.input["dF"] / HARTREE_OVER_BOHR_TO_EV_OVER_ANGSTROM
483
                )
484
            bgroup = self.input.sphinx.main[optimizer].create_group("bornOppenheimer")
1✔
485
            bgroup["scfDiag"] = self.get_scf_group()
1✔
486
        else:
487
            scf = self.input.sphinx.main.get("scfDiag", create=True)
1✔
488
            if self._generic_input["restart_for_band_structure"]:
1✔
489
                scf.append(self.get_scf_group(keepRhoFixed=True))
×
490
            else:
491
                scf.append(self.get_scf_group())
1✔
492
            if self.executable.version is not None:
1✔
493
                if self.get_version_float() > 2.5:
1✔
494
                    efgroup = self.input.sphinx.main.create_group("evalForces")
1✔
495
                    efgroup["file"] = '"relaxHist.sx"'
1✔
496
            else:
497
                warnings.warn("executable version could not be identified")
×
498

499
    def load_basis_group(self):
1✔
500
        """
501
        Load the basis Group.
502

503
        The group is populated using setdefault to avoid
504
        overwriting values that were previously (intentionally)
505
        modified.
506
        """
507
        self.input.sphinx.basis.setdefault("eCut", self.input["EnCut"] / RYDBERG_TO_EV)
1✔
508
        self.input.sphinx.basis.get("kPoint", create=True)
1✔
509
        if "KpointCoords" in self.input:
1✔
510
            self.input.sphinx.basis.kPoint.setdefault(
1✔
511
                "coords", np.array(self.input["KpointCoords"])
512
            )
513
        self.input.sphinx.basis.kPoint.setdefault("weight", 1)
1✔
514
        self.input.sphinx.basis.kPoint.setdefault("relative", True)
1✔
515
        if "KpointFolding" in self.input:
1✔
516
            self.input.sphinx.basis.setdefault(
1✔
517
                "folding", np.array(self.input["KpointFolding"])
518
            )
519
        self.input.sphinx.basis.setdefault("saveMemory", self.input["SaveMemory"])
1✔
520

521
    def load_hamilton_group(self):
1✔
522
        """
523
        Load the PAWHamiltonian Group.
524

525
        The group is populated using setdefault to avoid
526
        overwriting values that were previously (intentionally)
527
        modified.
528
        """
529
        self.input.sphinx.PAWHamiltonian.setdefault(
1✔
530
            "nEmptyStates", self.input["EmptyStates"]
531
        )
532
        self.input.sphinx.PAWHamiltonian.setdefault("ekt", self.input["Sigma"])
1✔
533
        for k in ["MethfesselPaxton", "FermiDirac"]:
1✔
534
            if k in self.input.list_nodes():
1✔
535
                self.input.sphinx.PAWHamiltonian.setdefault(k, self.input[k])
1✔
536
                break
1✔
537
        self.input.sphinx.PAWHamiltonian.setdefault("xc", self.input["Xcorr"])
1✔
538
        self.input.sphinx.PAWHamiltonian["spinPolarized"] = self._spin_enabled
1✔
539

540
    def load_guess_group(self, update_spins=True):
1✔
541
        """
542
        Load the initialGuess Group.
543

544
        The group is populated using setdefault to avoid
545
        overwriting values that were previously (intentionally)
546
        modified.
547

548
        Args:
549
            update_spins (bool): whether or not to reload the
550
                atomicSpin groups based on the latest structure.
551
                Defaults to True.
552
        """
553

554
        charge_density_file = None
1✔
555
        for ff in self.restart_file_list:
1✔
556
            if "rho.sxb" in ff.split("/")[-1]:
1✔
557
                charge_density_file = ff
1✔
558
        wave_function_file = None
1✔
559
        for ff in self.restart_file_list:
1✔
560
            if "waves.sxb" in ff.split("/")[-1]:
1✔
561
                wave_function_file = ff
1✔
562

563
        # introduce short alias for initialGuess group
564
        guess = self.input.sphinx.initialGuess
1✔
565

566
        guess.setdefault("waves", Group())
1✔
567
        guess.waves.setdefault("pawBasis", True)
1✔
568
        if wave_function_file is None:
1✔
569
            guess.waves.setdefault("lcao", Group())
1✔
570
        else:
571
            guess.waves.setdefault("file", '"' + wave_function_file + '"')
1✔
572
            # TODO: only for hybrid functionals
573
            guess.setdefault("exchange", Group())
1✔
574
            guess.exchange.setdefault("file", '"' + wave_function_file + '"')
1✔
575

576
        guess.setdefault("rho", Group())
1✔
577
        if charge_density_file is None:
1✔
578
            if wave_function_file is None:
1✔
579
                guess.rho.setdefault("atomicOrbitals", True)
1✔
580
                if self._spin_enabled:
1✔
581
                    init_spins = self.structure.get_initial_magnetic_moments()
1✔
582
                    # --- validate that initial spin moments are scalar
583
                    for spin in init_spins:
1✔
584
                        if isinstance(spin, list) or isinstance(spin, np.ndarray):
1✔
585
                            raise ValueError("SPHInX only supports collinear spins.")
×
586
                    guess.rho.get("atomicSpin", create=True)
1✔
587
                    if update_spins:
1✔
588
                        guess.rho.atomicSpin.clear()
1✔
589
                    # --- create initial spins if needed
590
                    if len(guess.rho.atomicSpin) == 0:
1✔
591
                        # set initial spin via label for each unique value of spin
592
                        # dict.from_keys (...).keys () deduplicates
593
                        for spin in dict.fromkeys(init_spins).keys():
1✔
594
                            guess.rho["atomicSpin"].append(
1✔
595
                                Group(
596
                                    {
597
                                        "label": '"spin_' + str(spin) + '"',
598
                                        "spin": str(spin),
599
                                    }
600
                                )
601
                            )
602
            else:
603
                guess.rho.setdefault("fromWaves", True)
1✔
604
        else:
605
            guess.rho.setdefault("file", '"' + charge_density_file + '"')
1✔
606

607
        if "noWavesStorage" not in guess:
1✔
608
            guess["noWavesStorage"] = not self.input["WriteWaves"]
1✔
609

610
    def calc_static(self, electronic_steps=100):
1✔
611
        """
612
        Setup the hamiltonian to perform a static SCF run.
613

614
        Loads defaults for all SPHInX input groups, including a static
615
        main Group.
616

617
        Args:
618
            electronic_steps (int): max # of electronic steps
619
        """
620
        if electronic_steps is not None:
1✔
621
            self.input["Estep"] = electronic_steps
1✔
622
        for arg in ["Istep", "dF", "dE"]:
1✔
623
            if arg in self.input:
1✔
624
                del self.input[arg]
1✔
625
        super().calc_static(electronic_steps=electronic_steps)
1✔
626
        self.load_default_groups()
1✔
627

628
    def calc_minimize(
1✔
629
        self,
630
        electronic_steps=60,
631
        ionic_steps=None,
632
        max_iter=None,
633
        pressure=None,
634
        ionic_energy=None,
635
        ionic_forces=None,
636
        ionic_energy_tolerance=None,
637
        ionic_force_tolerance=None,
638
        volume_only=False,
639
    ):
640
        """
641
        Setup the hamiltonian to perform ionic relaxations.
642

643
        The convergence goal can be set using either the
644
        ionic_energy_tolerance as a limit for fluctuations in energy or the
645
        ionic_force_tolerance.
646

647
        Loads defaults for all SPHInX input groups, including a
648
        ricQN-based main Group.
649

650
        .. warning::
651
            Sphinx does not support volume minimizations!  Calling this method with `pressure` or `volume_only` results
652
            in an error.
653

654
        Args:
655
            pressure:
656
            max_iter:
657
            electronic_steps (int): maximum number of electronic steps
658
                                    per electronic convergence
659
            ionic_steps (int): maximum number of ionic steps
660
            ionic_energy (float): convergence goal in terms of
661
                                  energy (depreciated use ionic_energy_tolerance instead)
662
            ionic_energy_tolerance (float): convergence goal in terms of
663
                                  energy (optional)
664
            ionic_forces (float): convergence goal in terms of
665
                                  forces (depreciated use ionic_force_tolerance instead)
666
            ionic_force_tolerance (float): convergence goal in terms of
667
                                  forces (optional)
668
            volume_only (bool):
669
        """
670
        if pressure is not None or volume_only:
1✔
671
            raise NotImplementedError(
×
672
                "pressure minimization is not implemented in SPHInX"
673
            )
674
        if electronic_steps is not None:
1✔
675
            self.input["Estep"] = electronic_steps
1✔
676
        if ionic_steps is not None:
1✔
677
            self.input["Istep"] = ionic_steps
1✔
678
        elif "Istep" not in self.input:
1✔
679
            self.input["Istep"] = 100
1✔
680
        if ionic_force_tolerance is not None:
1✔
681
            if ionic_force_tolerance < 0:
×
682
                raise ValueError("ionic_force_tolerance must be a positive integer")
×
683
            self.input["dF"] = float(ionic_force_tolerance)
×
684
        if ionic_energy_tolerance is not None:
1✔
685
            if ionic_energy_tolerance < 0:
1✔
686
                raise ValueError("ionic_force_tolerance must be a positive integer")
×
687
            self.input["dE"] = float(ionic_energy_tolerance)
1✔
688
        super(SphinxBase, self).calc_minimize(
1✔
689
            electronic_steps=electronic_steps,
690
            ionic_steps=ionic_steps,
691
            max_iter=max_iter,
692
            pressure=pressure,
693
            ionic_energy_tolerance=ionic_energy_tolerance,
694
            ionic_force_tolerance=ionic_force_tolerance,
695
            volume_only=volume_only,
696
        )
697
        self.load_default_groups()
1✔
698

699
    def calc_md(
1✔
700
        self, temperature=None, n_ionic_steps=1000, n_print=1, time_step=1.0, **kwargs
701
    ):
702
        raise NotImplementedError("calc_md() not implemented in SPHInX.")
×
703

704
    def restart_for_band_structure_calculations(self, job_name=None):
1✔
705
        """
706
        Restart a new job created from an existing calculation
707
        by reading the charge density for band structures.
708

709
        Args:
710
            job_name (str/None): Job name
711

712
        Returns:
713
            pyiron_atomistics.sphinx.sphinx.sphinx: new job instance
714
        """
715
        return self.restart_from_charge_density(
×
716
            job_name=job_name, band_structure_calc=True
717
        )
718

719
    def restart_from_charge_density(
1✔
720
        self, job_name=None, job_type="Sphinx", band_structure_calc=False
721
    ):
722
        """
723
        Restart a new job created from an existing calculation
724
        by reading the charge density.
725

726
        Args:
727
            job_name (str/None): Job name
728
            job_type (str/None): Job type. If not specified a SPHInX job type is assumed (actually
729
                this is all that's currently supported)
730
            band_structure_calc (bool): has to be True for band structure calculations.
731

732
        Returns:
733
            pyiron_atomistics.sphinx.sphinx.sphinx: new job instance
734
        """
735
        ham_new = self.restart(
×
736
            job_name=job_name,
737
            job_type=job_type,
738
            from_wave_functions=False,
739
            from_charge_density=True,
740
        )
741
        if band_structure_calc:
×
742
            ham_new._generic_input["restart_for_band_structure"] = True
×
743
            # --- clean up minimization related settings
744
            for setting in ["Istep", "dF", "dE"]:
×
745
                if setting in ham_new.input:
×
746
                    del ham_new.input[setting]
×
747
            # remove optimization-related stuff from GenericDFTJob
748
            super(SphinxBase, self).calc_static()
×
749
            # --- recreate main group
750
            del ham_new.input.sphinx["main"]
×
751
            ham_new.input.sphinx.create_group("main")
×
752
            ham_new.load_main_group()
×
753
        return ham_new
×
754

755
    def restart_from_wave_functions(
1✔
756
        self,
757
        job_name=None,
758
        job_type="Sphinx",
759
    ):
760
        """
761
        Restart a new job created from an existing calculation
762
        by reading the wave functions.
763

764
        Args:
765
            job_name (str): Job name
766
            job_type (str): Job type. If not specified a SPHInX job type is assumed (actually
767
                this is all that's currently supported.)
768

769
        Returns:
770
            pyiron_atomistics.sphinx.sphinx.sphinx: new job instance
771
        """
772
        return self.restart(
×
773
            job_name=job_name,
774
            job_type=job_type,
775
            from_wave_functions=True,
776
            from_charge_density=False,
777
        )
778

779
    def restart(
1✔
780
        self,
781
        job_name=None,
782
        job_type=None,
783
        from_charge_density=True,
784
        from_wave_functions=True,
785
    ):
786
        if not self.status.finished and not self.is_compressed():
×
787
            # self.decompress()
788
            with warnings.catch_warnings(record=True) as w:
×
789
                try:
×
790
                    self.collect_output()
×
791
                except AssertionError as orig_error:
×
792
                    if from_charge_density or from_wave_functions:
×
793
                        raise AssertionError(
×
794
                            orig_error.message
795
                            + "\nCowardly refusing to use density or wavefunctions for restart.\n"
796
                            + "Solution: set from_charge_density and from_wave_functions to False."
797
                        )
798
                if len(w) > 0:
×
799
                    self.status.not_converged = True
×
800
        new_job = super(SphinxBase, self).restart(job_name=job_name, job_type=job_type)
×
801

802
        new_job.input = self.input.copy()
×
803

804
        recreate_guess = False
×
805
        if from_charge_density and os.path.isfile(
×
806
            posixpath.join(self.working_directory, "rho.sxb")
807
        ):
808
            new_job.restart_file_list.append(
×
809
                posixpath.join(self.working_directory, "rho.sxb")
810
            )
811
            del new_job.input.sphinx.initialGuess["rho"]
×
812
            recreate_guess = True
×
813

814
        elif from_charge_density:
×
815
            self._logger.warning(
×
816
                msg=f"A charge density from job: {self.job_name} "
817
                + "is not generated and therefore it can't be read."
818
            )
819
        if from_wave_functions and os.path.isfile(
×
820
            posixpath.join(self.working_directory, "waves.sxb")
821
        ):
822
            new_job.restart_file_list.append(
×
823
                posixpath.join(self.working_directory, "waves.sxb")
824
            )
825
            try:
×
826
                del new_job.input.sphinx.initialGuess["rho"]
×
827
            except KeyError:
×
828
                pass
×
829
            del new_job.input.sphinx.initialGuess["waves"]
×
830
            recreate_guess = True
×
831

832
        elif from_wave_functions:
×
833
            self._logger.warning(
×
834
                msg="No wavefunction file (waves.sxb) was found for "
835
                + f"job {self.job_name} in {self.working_directory}."
836
            )
837
        if recreate_guess:
×
838
            new_job.load_guess_group()
×
839
        return new_job
×
840

841
    def relocate_hdf5(self, h5_path=None):
1✔
842
        self.input._force_load()
×
843
        super().relocate_hdf5(h5_path=h5_path)
×
844

845
    def to_hdf(self, hdf=None, group_name=None):
1✔
846
        """
847
        Stores the instance attributes into the hdf5 file
848

849
        Args:
850
            hdf (str): Path to the hdf5 file
851
            group_name (str): Name of the group which contains the object
852
        """
853
        super(SphinxBase, self).to_hdf(hdf=hdf, group_name=group_name)
1✔
854
        self._structure_to_hdf()
1✔
855
        with self._hdf5.open("input") as hdf:
1✔
856
            self.input.to_hdf(hdf)
1✔
857
        self.output.to_hdf(self._hdf5)
1✔
858

859
    def from_hdf(self, hdf=None, group_name=None):
1✔
860
        """
861
        Recreates instance from the hdf5 file
862

863
        Args:
864
            hdf (str): Path to the hdf5 file
865
            group_name (str): Name of the group which contains the object
866
        """
867
        if "HDF_VERSION" not in self._hdf5.keys():
×
868
            from pyiron_base import GenericParameters
×
869

870
            super(SphinxBase, self).from_hdf(hdf=hdf, group_name=group_name)
×
871
            self._structure_from_hdf()
×
872
            gp = GenericParameters(table_name="input")
×
873
            gp.from_hdf(self._hdf5)
×
874
            for k in gp.keys():
×
875
                self.input[k] = gp[k]
×
876
        elif self._hdf5["HDF_VERSION"] == "0.1.0":
×
877
            super(SphinxBase, self).from_hdf(hdf=hdf, group_name=group_name)
×
878
            self._structure_from_hdf()
×
879
            with self._hdf5.open("input") as hdf:
×
880
                self.input.from_hdf(hdf, group_name="parameters")
×
881
        self.output.from_hdf(self._hdf5)
×
882

883
    def from_directory(self, directory, file_name="structure.sx"):
1✔
884
        try:
×
885
            if not self.status.finished:
×
886
                file_path = posixpath.join(directory, file_name)
×
887
                if os.path.isfile(file_path):
×
888
                    self.structure = read_atoms(file_path)
×
889
                else:
890
                    raise ValueError(
×
891
                        f"File {file_path} not found. "
892
                        "Please double check the directory and file name."
893
                    )
894

895
                self.output.collect(directory=directory)
×
896
                self.to_hdf(self._hdf5)
×
897
            else:
898
                self.output.from_hdf(self._hdf5)
×
899
            self.status.finished = True
×
900
        except Exception as err:
×
901
            print(err)
×
902
            self.status.aborted = True
×
903

904
    def set_check_overlap(self, check_overlap=True):
1✔
905
        """
906
        Args:
907
            check_overlap (bool): Whether to check overlap
908

909
        Comments:
910
            Certain PAW-pseudo-potentials have an intrinsic pathology:
911
            their PAW overlap operator is not generally positive definite
912
            (i.e., the PAW-corrected norm of a wavefunction could become
913
            negative). SPHInX usually refuses to use such problematic
914
            potentials. This behavior can be overridden by setting
915
            check_overlap to False.
916
        """
917
        if not isinstance(check_overlap, bool):
1✔
918
            raise TypeError("check_overlap has to be a boolean")
1✔
919

920
        if self.get_version_float() < 2.51 and not check_overlap:
×
921
            warnings.warn(
×
922
                "SPHInX executable version has to be 2.5.1 or above "
923
                + "in order for the overlap to be considered. "
924
                + "Change it via job.executable.version"
925
            )
926
        self.input.CheckOverlap = check_overlap
×
927

928
    def set_mixing_parameters(
1✔
929
        self,
930
        method=None,
931
        n_pulay_steps=None,
932
        density_mixing_parameter=None,
933
        spin_mixing_parameter=None,
934
        density_residual_scaling=None,
935
        spin_residual_scaling=None,
936
    ):
937
        """
938
        Further information can be found on the website:
939
        https://sxrepo.mpie.de
940
        """
941
        method_list = ["PULAY", "KERKER", "LINEAR"]
1✔
942
        if method is not None and method.upper() not in method_list:
1✔
943
            raise ValueError("Mixing method has to be PULAY or KERKER")
1✔
944
        if method is not None:
1✔
945
            self.input["mixingMethod"] = method.upper().replace("KERKER", "LINEAR")
1✔
946
        if n_pulay_steps is not None:
1✔
947
            self.input["nPulaySteps"] = int(n_pulay_steps)
1✔
948
        if density_mixing_parameter is not None:
1✔
949
            if density_mixing_parameter > 1.0 or density_mixing_parameter < 0:
1✔
950
                raise ValueError(
1✔
951
                    "density_mixing_parameter has to be between 0 and 1 "
952
                    + "(default value is 1)"
953
                )
954
            self.input["rhoMixing"] = density_mixing_parameter
1✔
955
        if spin_mixing_parameter is not None:
1✔
956
            if spin_mixing_parameter > 1.0 or spin_mixing_parameter < 0:
1✔
957
                raise ValueError(
1✔
958
                    "spin_mixing_parameter has to be between 0 and 1 "
959
                    + "(default value is 1)"
960
                )
961
            self.input["spinMixing"] = spin_mixing_parameter
1✔
962
        if density_residual_scaling is not None:
1✔
963
            if density_residual_scaling <= 0:
1✔
964
                raise ValueError("density_residual_scaling must be a positive value")
1✔
965
            self.input["rhoResidualScaling"] = density_residual_scaling
1✔
966
        if spin_residual_scaling is not None:
1✔
967
            if spin_residual_scaling <= 0:
1✔
968
                raise ValueError("spin_residual_scaling   must be a positive value")
1✔
969
            self.input["spinResidualScaling"] = spin_residual_scaling
1✔
970

971
    set_mixing_parameters.__doc__ = (
1✔
972
        GenericDFTJob.set_mixing_parameters.__doc__ + set_mixing_parameters.__doc__
973
    )
974

975
    def set_occupancy_smearing(self, smearing=None, width=None, order=1):
1✔
976
        """
977
        Set how the finite temperature smearing is applied in
978
        determining partial occupancies
979

980
        Args:
981
            smearing (str): Type of smearing, 'FermiDirac' or 'MethfesselPaxton'
982
            width (float): Smearing width (eV) (default: 0.2)
983
            order (int): Smearing order
984
        """
985
        if smearing is not None:
1✔
986
            if not isinstance(smearing, str):
1✔
987
                raise ValueError("Smearing must be a string")
1✔
988
            if smearing.lower().startswith("meth"):
1✔
989
                self.input.MethfesselPaxton = order
1✔
990
                if "FermiDirac" in self.input.list_nodes():
1✔
991
                    del self.input["FermiDirac"]
1✔
992
            elif smearing.lower().startswith("fermi"):
1✔
993
                self.input.FermiDirac = order
1✔
994
                if "MethfesselPaxton" in self.input.list_nodes():
1✔
995
                    del self.input["MethfesselPaxton"]
1✔
996
        if width is not None and width < 0:
1✔
997
            raise ValueError("Smearing value must be a float >= 0")
1✔
998
        if width is not None:
1✔
999
            self.input["Sigma"] = width
1✔
1000

1001
    @deprecate(
1✔
1002
        ionic_forces="Use ionic_force_tolerance",
1003
        ionic_energy="use ionic_energy_tolerance",
1004
    )
1005
    def set_convergence_precision(
1✔
1006
        self,
1007
        ionic_energy_tolerance=None,
1008
        ionic_force_tolerance=None,
1009
        ionic_energy=None,
1010
        electronic_energy=None,
1011
        ionic_forces=None,
1012
    ):
1013
        """
1014
        Sets the electronic and ionic convergence precision.
1015

1016
        For ionic convergence either the energy or the force
1017
        precision is required.
1018

1019
        Args:
1020
            ionic_energy (float): Ionic energy convergence precision
1021
                                  (depreciated use ionic_energy_tolerance instead)
1022
            ionic_energy_tolerance (float): Ionic energy convergence precision
1023
            electronic_energy (float): Electronic energy convergence
1024
                                       precision
1025
            ionic_forces (float): Ionic force convergence precision
1026
                                  (depreciated use ionic_force_tolerance instead)
1027
            ionic_force_tolerance (float): Ionic force convergence precision
1028
        """
1029
        if ionic_forces is not None:
1✔
1030
            ionic_force_tolerance = ionic_forces
×
1031
        if ionic_energy is not None:
1✔
1032
            ionic_energy_tolerance = ionic_energy
×
1033
        assert (
1✔
1034
            ionic_energy_tolerance is None or ionic_energy_tolerance > 0
1035
        ), "ionic_energy_tolerance must be a positive float"
1036
        assert (
1✔
1037
            ionic_force_tolerance is None or ionic_force_tolerance > 0
1038
        ), "ionic_force_tolerance must be a positive float"
1039
        assert (
1✔
1040
            electronic_energy is None or electronic_energy > 0
1041
        ), "electronic_energy must be a positive float"
1042
        if ionic_energy_tolerance is not None or ionic_force_tolerance is not None:
1✔
1043
            # self.input["dE"] = ionic_energy_tolerance
1044
            # self.input["dF"] = ionic_force_tolerance
1045
            print("Setting calc_minimize")
1✔
1046
            self.calc_minimize(
1✔
1047
                ionic_energy_tolerance=ionic_energy_tolerance,
1048
                ionic_force_tolerance=ionic_force_tolerance,
1049
            )
1050
        if electronic_energy is not None:
1✔
1051
            self.input["Ediff"] = electronic_energy
1✔
1052

1053
    def set_empty_states(self, n_empty_states=None):
1✔
1054
        """
1055
        Function to set the number of empty states.
1056

1057
        Args:
1058
            n_empty_states (int/None): Number of empty states.
1059
            If None, sets it to 'auto'.
1060

1061
        Comments:
1062
            If this number is too low, the algorithm will not be
1063
            able to able to swap wave functions near the chemical
1064
            potential. If the number is too high, computation time
1065
            will be wasted for the higher energy states and
1066
            potentially lead to a memory problem.
1067

1068
            In contrast to VASP, this function sets only the number
1069
            of empty states and not the number of total states.
1070

1071
            The default value is 0.5*NIONS+3 for non-magnetic systems
1072
            and 1.5*NIONS+3 for magnetic systems
1073
        """
1074
        if n_empty_states is None:
1✔
1075
            # will be converted later; see load_default_groups
1076
            self.input["EmptyStates"] = "auto"
1✔
1077
        else:
1078
            if n_empty_states < 0:
1✔
1079
                raise ValueError("Number of empty states must be greater than 0")
1✔
1080
            self.input["EmptyStates"] = n_empty_states
1✔
1081
        self.input.sphinx.PAWHamiltonian.nEmptyStates = self.input["EmptyStates"]
1✔
1082

1083
    def _set_kpoints(
1✔
1084
        self,
1085
        mesh=None,
1086
        scheme="MP",
1087
        center_shift=None,
1088
        symmetry_reduction=True,
1089
        manual_kpoints=None,
1090
        weights=None,
1091
        reciprocal=True,
1092
        n_path=None,
1093
        path_name=None,
1094
    ):
1095
        """
1096
        Function to setup the k-points for the SPHInX job
1097

1098
        Args:
1099
            reciprocal (bool): Tells if the supplied values are in
1100
                               reciprocal (direct) or cartesian coordinates
1101
                               (in reciprocal space) (not implemented)
1102
            weights (list): Manually supplied weights to each k-point in
1103
                            case of the manual mode (not implemented)
1104
            manual_kpoints (list): Manual list of k-points (not implemented)
1105
            symmetry_reduction (bool): Tells if the symmetry reduction is
1106
                                       to be applied to the k-points
1107
            scheme (str): Type of k-point generation scheme ('MP' or 'Line')
1108
            mesh (list): Size of the mesh (in the MP scheme)
1109
            center_shift (list): Shifts the center of the mesh from the
1110
                                 gamma point by the given vector
1111
            n_path (int): Number of points per trace part for line mode
1112
            path_name (str): Name of high symmetry path used for band
1113
                             structure calculations.
1114
        """
1115
        if not isinstance(symmetry_reduction, bool):
1✔
1116
            raise ValueError("symmetry_reduction has to be a boolean")
1✔
1117
        if manual_kpoints is not None:
1✔
1118
            raise ValueError(
×
1119
                "manual_kpoints is not yet implemented in pyiron for SPHInX"
1120
            )
1121
        if weights is not None:
1✔
1122
            raise ValueError(
×
1123
                "manual weights are not yet implmented in Pyiron for " + "SPHInX"
1124
            )
1125

1126
        if scheme == "MP":
1✔
1127
            # Remove kPoints and set kPoint
1128
            if "kPoints" in self.input.sphinx.basis:
1✔
1129
                del self.input.sphinx.basis.kPoints
1✔
1130
            self.input.sphinx.basis.get("kPoint", create=True)
1✔
1131
            if mesh is not None:
1✔
1132
                self.input["KpointFolding"] = list(mesh)
1✔
1133
                self.input.sphinx.basis["folding"] = np.array(
1✔
1134
                    self.input["KpointFolding"]
1135
                )
1136
            if center_shift is not None:
1✔
1137
                self.input["KpointCoords"] = list(center_shift)
1✔
1138
                self.input.sphinx.basis["kPoint"]["coords"] = np.array(
1✔
1139
                    self.input["KpointCoords"]
1140
                )
1141
                self.input.sphinx.basis.kPoint["weight"] = 1
1✔
1142
                self.input.sphinx.basis.kPoint["relative"] = True
1✔
1143

1144
        elif scheme == "Line":
1✔
1145
            # Remove Kpoint and set Kpoints
1146

1147
            if "kPoint" in self.input.sphinx.basis:
1✔
1148
                del self.input.sphinx.basis["kPoint"]
1✔
1149
                del self.input["KpointFolding"]
1✔
1150
                del self.input["KpointCoords"]
1✔
1151
                if "folding" in self.input.sphinx.basis:
1✔
1152
                    del self.input.sphinx.basis["folding"]
1✔
1153
            if n_path is None and self._generic_input["n_path"] is None:
1✔
1154
                raise ValueError("'n_path' has to be defined")
1✔
1155
            if n_path is None:
1✔
1156
                n_path = self._generic_input["n_path"]
×
1157
            else:
1158
                self._generic_input["n_path"] = n_path
1✔
1159

1160
            if self.structure.get_high_symmetry_points() is None:
1✔
1161
                raise ValueError("no 'high_symmetry_points' defined for 'structure'.")
×
1162

1163
            if path_name is None and self._generic_input["path_name"] is None:
1✔
1164
                raise ValueError("'path_name' has to be defined")
1✔
1165
            if path_name is None:
1✔
1166
                path_name = self._generic_input["path_name"]
×
1167
            else:
1168
                self._generic_input["path_name"] = path_name
1✔
1169

1170
            try:
1✔
1171
                path = self.structure.get_high_symmetry_path()[path_name]
1✔
1172
            except KeyError:
1✔
1173
                raise AssertionError("'{}' is not a valid path!".format(path_name))
1✔
1174

1175
            def make_point(point, n_path):
1✔
1176
                return Group(
1✔
1177
                    {
1178
                        "coords": np.array(
1179
                            self.structure.get_high_symmetry_points()[point]
1180
                        ),
1181
                        "nPoints": n_path,
1182
                        "label": '"{}"'.format(point.replace("'", "p")),
1183
                    }
1184
                )
1185

1186
            kpoints = Group({"relative": True})
1✔
1187
            kpoints["from"] = make_point(path[0][0], None)
1✔
1188
            # from nodes are not supposed to have a nPoints attribute
1189
            del kpoints["from/nPoints"]
1✔
1190

1191
            kpoints.create_group("to").append(make_point(path[0][1], n_path))
1✔
1192

1193
            for segment in path[1:]:
1✔
1194
                # if the last node on the so far is not the same as the first
1195
                # node of this path segment, then we need to insert another
1196
                # node into the path to alert sphinx that we want a cut in our
1197
                # band structure (n_path = 0)
1198
                if '"{}"'.format(segment[0]) != kpoints.to[-1].label:
1✔
1199
                    kpoints["to"].append(make_point(segment[0], 0))
1✔
1200

1201
                kpoints["to"].append(make_point(segment[1], n_path))
1✔
1202

1203
            self.input.sphinx.basis["kPoints"] = kpoints
1✔
1204
        else:
1205
            raise ValueError(
1✔
1206
                "only Monkhorst-Pack mesh and Line mode\
1207
                are currently implemented in Pyiron for SPHInX"
1208
            )
1209

1210
    def load_default_groups(self):
1✔
1211
        """
1212
        Populates input groups with the default values.
1213

1214
        Nearly every default simply points to a variable stored in
1215
        self.input.
1216

1217
        Does not load job.input.structure or job.input.pawPot.
1218
        These groups should usually be modified via job.structure,
1219
        in which case they will be set at the last minute when
1220
        the job is run. These groups can be synced to job.structure
1221
        at any time using job.load_structure_group() and
1222
        job.load_species_group().
1223
        """
1224

1225
        if self.structure is None:
1✔
1226
            raise AssertionError(
1✔
1227
                f"{self.job_name} has not been assigned "
1228
                + "a structure. Please load one first (e.g. "
1229
                + f"{self.job_name}.structure = ...)"
1230
            )
1231

1232
        if self.input["EmptyStates"] == "auto":
1✔
1233
            if self._spin_enabled:
1✔
1234
                self.input["EmptyStates"] = int(1.5 * len(self.structure) + 3)
1✔
1235
            else:
1236
                self.input["EmptyStates"] = int(len(self.structure) + 3)
1✔
1237

1238
        if not self.input.sphinx.basis.read_only:
1✔
1239
            self.load_basis_group()
1✔
1240
        if not self.input.sphinx.structure.read_only:
1✔
1241
            self.load_structure_group()
1✔
1242
        if self.input["VaspPot"]:
1✔
1243
            potformat = "VASP"
×
1244
        else:
1245
            potformat = "JTH"
1✔
1246
        if not self.input.sphinx.pawPot.read_only:
1✔
1247
            self.load_species_group(
1✔
1248
                check_overlap=self.input.CheckOverlap, potformat=potformat
1249
            )
1250
        if not self.input.sphinx.initialGuess.read_only:
1✔
1251
            self.load_guess_group()
1✔
1252
        if not self.input.sphinx.PAWHamiltonian.read_only:
1✔
1253
            self.load_hamilton_group()
1✔
1254
        if not self.input.sphinx.main.read_only:
1✔
1255
            self.load_main_group()
1✔
1256

1257
    def list_potentials(self):
1✔
1258
        """
1259
        Lists all the possible POTCAR files for the elements in the structure depending on the XC functional
1260

1261
        Returns:
1262
           list: a list of available potentials
1263
        """
1264
        return self.potential_list
1✔
1265

1266
    def _get_potential_path(
1✔
1267
        self,
1268
        potformat="JTH",
1269
        xc=None,
1270
        cwd=None,
1271
        pot_path_dict=None,
1272
        modified_elements=None,
1273
    ):
1274
        """
1275
        Copy potential files
1276

1277
        Args:
1278
            potformat (str):
1279
            xc (str/None):
1280
            cwd (str/None):
1281
            pot_path_dict (dict):
1282
            modified_elements (dict):
1283
        """
1284

1285
        if pot_path_dict is None:
1✔
1286
            pot_path_dict = {}
1✔
1287

1288
        if potformat == "JTH":
1✔
1289
            potentials = SphinxJTHPotentialFile(xc=xc)
1✔
1290
            pot_path_dict.setdefault("PBE", "jth-gga-pbe")
1✔
1291
        elif potformat == "VASP":
×
1292
            potentials = VaspPotentialFile(xc=xc)
×
1293
            pot_path_dict.setdefault("PBE", "paw-gga-pbe")
×
1294
            pot_path_dict.setdefault("LDA", "paw-lda")
×
1295
        else:
1296
            raise ValueError("Only JTH and VASP potentials are supported!")
×
1297

1298
        ori_paths, des_paths = [], []
1✔
1299
        for species_obj in self.structure.get_species_objects():
1✔
1300
            if species_obj.Parent is not None:
1✔
1301
                elem = species_obj.Parent
×
1302
            else:
1303
                elem = species_obj.Abbreviation
1✔
1304

1305
            if "pseudo_potcar_file" in species_obj.tags.keys():
1✔
1306
                new_element = species_obj.tags["pseudo_potcar_file"]
×
1307
                potentials.add_new_element(parent_element=elem, new_element=new_element)
×
1308
                potential_path = potentials.find_potential_file(
×
1309
                    path=potentials.find_default(new_element)["Filename"].values[0][0]
1310
                )
1311
                assert os.path.isfile(
×
1312
                    potential_path
1313
                ), "such a file does not exist in the pp directory"
1314
            elif elem in modified_elements.keys():
1✔
1315
                new_element = modified_elements[elem]
×
1316
                if os.path.isabs(new_element):
×
1317
                    potential_path = new_element
×
1318
                else:
1319
                    potentials.add_new_element(
×
1320
                        parent_element=elem, new_element=new_element
1321
                    )
1322
                    potential_path = potentials.find_potential_file(
×
1323
                        path=potentials.find_default(new_element)["Filename"].values[0][
1324
                            0
1325
                        ]
1326
                    )
1327
            else:
1328
                ori_paths.append(
1✔
1329
                    potentials.find_potential_file(
1330
                        path=potentials.find_default(elem)["Filename"].values[0][0]
1331
                    )
1332
                )
1333
            if potformat == "JTH":
1✔
1334
                des_paths.append(posixpath.join(cwd, elem + "_GGA.atomicdata"))
1✔
1335
            else:
1336
                des_paths.append(posixpath.join(cwd, elem + "_POTCAR"))
×
1337
        return {"origins": ori_paths, "destinations": des_paths}
1✔
1338

1339
    def write_input(self):
1✔
1340
        """
1341
        Generate all the required input files for the SPHInX job.
1342

1343
        Creates:
1344
        structure.sx: structure associated w/ job
1345
        all pseudopotential files
1346
        spins.in (if necessary): constrained spin moments
1347
        input.sx: main input file with all sub-groups
1348

1349
        Automatically called by job.run()
1350
        """
1351
        super().write_input()
1✔
1352

1353
        # If the structure group was not modified directly by the
1354
        # user, via job.input.structure (which is likely True),
1355
        # load it based on job.structure.
1356
        structure_sync = str(self.input.sphinx.structure) == str(
1✔
1357
            self.get_structure_group()
1358
        )
1359
        if not structure_sync and not self.input.sphinx.structure.read_only:
1✔
1360
            self.load_structure_group()
×
1361

1362
        # copy potential files to working directory
1363
        if self.input["VaspPot"]:
1✔
1364
            potformat = "VASP"
×
1365
        else:
1366
            potformat = "JTH"
1✔
1367
        # If the species group was not modified directly by the user,
1368
        # via job.input.pawPot (which is likely True),
1369
        # load it based on job.structure.
1370
        if not structure_sync and not self.input.sphinx.pawPot.read_only:
1✔
1371
            self.load_species_group(
×
1372
                check_overlap=self.input.CheckOverlap, potformat=potformat
1373
            )
1374

1375
        modified_elements = {
1✔
1376
            key: value
1377
            for key, value in self._potential.to_dict().items()
1378
            if value is not None
1379
        }
1380

1381
        copy_potentials(
1✔
1382
            **self._get_potential_path(
1383
                potformat=potformat,
1384
                xc=self.input["Xcorr"],
1385
                cwd=self.working_directory,
1386
                modified_elements=modified_elements,
1387
            )
1388
        )
1389

1390
        # Write spin constraints, if set via _generic_input.
1391
        all_groups = [
1✔
1392
            self.input.sphinx.pawPot,
1393
            self.input.sphinx.structure,
1394
            self.input.sphinx.basis,
1395
            self.input.sphinx.PAWHamiltonian,
1396
            self.input.sphinx.initialGuess,
1397
            self.input.sphinx.main,
1398
        ]
1399

1400
        if self._generic_input["fix_spin_constraint"] and self.structure.has(
1✔
1401
            "initial_magmoms"
1402
        ):
1403
            self.input.sphinx.spinConstraint = Group()
1✔
1404
            all_groups.append(self.input.sphinx.spinConstraint)
1✔
1405
            write_spin_constraints(
1✔
1406
                cwd=self.working_directory,
1407
                magmoms=self.structure.get_initial_magnetic_moments()[
1408
                    self.id_pyi_to_spx
1409
                ],
1410
                constraints=self.structure.spin_constraint[self.id_pyi_to_spx],
1411
            )
1412
            self.input.sphinx.spinConstraint.setdefault("file", '"spins.in"')
1✔
1413

1414
        # In case the entire group was
1415
        # set/overwritten as a normal dict.
1416
        for group in all_groups:
1✔
1417
            group = Group(group)
1✔
1418

1419
        # write input.sx
1420
        file_name = posixpath.join(self.working_directory, "input.sx")
1✔
1421
        with open(file_name, "w") as f:
1✔
1422
            f.write(f"//{self.job_name}\n")
1✔
1423
            f.write("//SPHInX input file generated by pyiron\n\n")
1✔
1424
            f.write("format paw;\n")
1✔
1425
            f.write("include <parameters.sx>;\n\n")
1✔
1426
            f.write(self.input.sphinx.to_sphinx(indent=0))
1✔
1427

1428
    @property
1✔
1429
    def _spin_enabled(self):
1✔
1430
        return self.structure.has("initial_magmoms")
1✔
1431

1432
    def get_charge_density(self):
1✔
1433
        """
1434
        Gets the charge density from the hdf5 file. This value is normalized by the volume
1435

1436
        Returns:
1437
            pyiron_atomistics.atomistics.volumetric.generic.VolumetricData
1438
        """
1439
        if self.status not in job_status_successful_lst:
×
1440
            return
×
1441
        else:
1442
            with self.project_hdf5.open("output") as ho:
×
1443
                cd_obj = SphinxVolumetricData()
×
1444
                cd_obj.from_hdf(ho, "charge_density")
×
1445
            cd_obj.atoms = self.get_structure(-1)
×
1446
            return cd_obj
×
1447

1448
    def get_electrostatic_potential(self):
1✔
1449
        """
1450
        Gets the electrostatic potential from the hdf5 file.
1451

1452
        Returns:
1453
            pyiron_atomistics.atomistics.volumetric.generic.VolumetricData
1454
        """
1455
        if self.status not in job_status_successful_lst:
×
1456
            return
×
1457
        else:
1458
            with self.project_hdf5.open("output") as ho:
×
1459
                es_obj = SphinxVolumetricData()
×
1460
                es_obj.from_hdf(ho, "electrostatic_potential")
×
1461
            es_obj.atoms = self.get_structure(-1)
×
1462
            return es_obj
×
1463

1464
    def collect_output(self, force_update=False, compress_files=True):
1✔
1465
        """
1466
        Collects the outputs and stores them to the hdf file
1467
        """
1468
        if self.is_compressed():
1✔
1469
            warnings.warn("Job already compressed - output not collected")
×
1470
            return
×
1471
        self.output.collect(directory=self.working_directory)
1✔
1472
        self.output.to_hdf(self._hdf5, force_update=force_update)
1✔
1473
        if compress_files:
1✔
1474
            self.compress()
1✔
1475

1476
    def convergence_check(self):
1✔
1477
        """
1478
        Checks for electronic and ionic convergence according to the user specified tolerance
1479

1480
        Returns:
1481

1482
            bool: True if converged
1483

1484
        """
1485
        # Checks if sufficient empty states are present
1486
        if not self.nbands_convergence_check():
×
1487
            return False
×
1488
        return self.output.generic.dft.scf_convergence[-1]
×
1489

1490
    def collect_logfiles(self):
1✔
1491
        """
1492
        Collect errors and warnings.
1493
        """
1494
        self.collect_errors()
×
1495
        self.collect_warnings()
×
1496

1497
    def collect_warnings(self):
1✔
1498
        """
1499
        Collects warnings from the SPHInX run
1500
        """
1501
        self._logger.info(
×
1502
            "collect_warnings() is not yet \
1503
            implemented for SPHInX"
1504
        )
1505

1506
    def collect_errors(self):
1✔
1507
        """
1508
        Collects errors from the SPHInX run
1509
        """
1510
        self._logger.info("collect_errors() is not yet implemented for SPHInX")
×
1511

1512
    def get_n_ir_reciprocal_points(
1✔
1513
        self, is_time_reversal=True, symprec=1e-5, ignore_magmoms=False
1514
    ):
1515
        lattice = self.structure.cell
1✔
1516
        positions = self.structure.get_scaled_positions()
1✔
1517
        numbers = self.structure.get_atomic_numbers()
1✔
1518
        if ignore_magmoms:
1✔
1519
            magmoms = np.zeros(len(positions))
×
1520
        else:
1521
            magmoms = self.structure.get_initial_magnetic_moments()
1✔
1522
        mag_num = np.array(list(zip(magmoms, numbers)))
1✔
1523
        satz = np.unique(mag_num, axis=0)
1✔
1524
        numbers = []
1✔
1525
        for nn in np.all(satz == mag_num[:, np.newaxis], axis=-1):
1✔
1526
            numbers.append(np.arange(len(satz))[nn][0])
1✔
1527
        mapping, _ = spglib.get_ir_reciprocal_mesh(
1✔
1528
            mesh=[int(self.input["KpointFolding"][k]) for k in range(3)],
1529
            cell=(lattice, positions, numbers),
1530
            is_shift=np.dot(self.structure.cell, np.array(self.input["KpointCoords"])),
1531
            is_time_reversal=is_time_reversal,
1532
            symprec=symprec,
1533
        )
1534
        return len(np.unique(mapping))
1✔
1535

1536
    def check_setup(self):
1✔
1537
        with warnings.catch_warnings(record=True) as w:
1✔
1538
            # Check for parameters that were not modified but
1539
            # possibly should have (encut, kpoints, smearing, etc.),
1540
            # or were set to nonsensical values.
1541

1542
            if (
1✔
1543
                not (
1544
                    isinstance(self.input.sphinx.basis["eCut"], int)
1545
                    or isinstance(self.input.sphinx.basis["eCut"], float)
1546
                )
1547
                or round(self.input.sphinx.basis["eCut"] * RYDBERG_TO_EV, 0) == 340
1548
            ):
1549
                warnings.warn(
1✔
1550
                    "Energy cut-off value wrong or not modified from default "
1551
                    + "340 eV; change it via job.set_encut()"
1552
                )
1553
            if "kPoint" in self.input.sphinx.basis:
1✔
1554
                if not (
1✔
1555
                    isinstance(self.input.sphinx.basis["kPoint"]["coords"], np.ndarray)
1556
                    or len(self.input.sphinx.basis["kPoint"]["coords"]) != 3
1557
                ):
1558
                    warnings.warn("K point coordinates seem to be inappropriate")
×
1559
            if (
1✔
1560
                not (
1561
                    isinstance(self.input.sphinx.PAWHamiltonian["ekt"], int)
1562
                    or isinstance(self.input.sphinx.PAWHamiltonian["ekt"], float)
1563
                )
1564
                or round(self.input.sphinx.PAWHamiltonian["ekt"], 1) == 0.2
1565
            ):
1566
                warnings.warn(
1✔
1567
                    "Fermi smearing value wrong or not modified from default "
1568
                    + "0.2 eV; change it via job.set_occupancy_smearing()"
1569
                )
1570
            if not (
1✔
1571
                isinstance(self.input.sphinx.basis["folding"], np.ndarray)
1572
                or len(self.input.sphinx.basis["folding"]) != 3
1573
            ) or self.input.sphinx.basis["folding"].tolist() == [4, 4, 4]:
1574
                warnings.warn(
1✔
1575
                    "K point folding wrong or not modified from default "
1576
                    + "[4,4,4]; change it via job.set_kpoints()"
1577
                )
1578
            if self.get_n_ir_reciprocal_points() < self.server.cores:
1✔
1579
                warnings.warn(
1✔
1580
                    "Number of cores exceed number of irreducible "
1581
                    + "reciprocal points: "
1582
                    + str(self.get_n_ir_reciprocal_points())
1583
                )
1584
            if self.input["EmptyStates"] == "auto":
1✔
1585
                if self._spin_enabled:
1✔
1586
                    warnings.warn(
×
1587
                        "Number of empty states was not specified. Default: "
1588
                        + "3+NIONS*1.5 for magnetic systems. "
1589
                    )
1590
                else:
1591
                    warnings.warn(
1✔
1592
                        "Number of empty states was not specified. Default: "
1593
                        + "3+NIONS for non-magnetic systems"
1594
                    )
1595

1596
            return len(w) == 0
1✔
1597

1598
    def validate_ready_to_run(self):
1✔
1599
        """
1600
        Checks whether parameters are set appropriately. It does not
1601
        mean the simulation won't run if it returns False.
1602
        """
1603

1604
        all_groups = {
1✔
1605
            "job.input.pawPot": self.input.sphinx.pawPot,
1606
            "job.input.structure": self.input.sphinx.structure,
1607
            "job.input.basis": self.input.sphinx.basis,
1608
            "job.input.PAWHamiltonian": self.input.sphinx.PAWHamiltonian,
1609
            "job.input.initialGuess": self.input.sphinx.initialGuess,
1610
            "job.input.main": self.input.sphinx.main,
1611
        }
1612

1613
        if np.any([len(all_groups[group]) == 0 for group in all_groups]):
1✔
1614
            self.load_default_groups()
1✔
1615

1616
        if self.structure is None:
1✔
1617
            raise AssertionError(
1✔
1618
                "Structure not set; set it via job.structure = "
1619
                + "Project().create_structure()"
1620
            )
1621
        if self.input["THREADS"] > self.server.cores:
1✔
1622
            raise AssertionError(
1✔
1623
                "Number of cores cannot be smaller than the number "
1624
                + "of OpenMP threads"
1625
            )
1626
        with warnings.catch_warnings(record=True) as w:
1✔
1627
            # Warn about discrepancies between values in
1628
            # self.input and individual groups, in case
1629
            # a user modified them directly
1630
            if round(self.input["EnCut"], 0) != round(
1✔
1631
                self.input.sphinx.basis.eCut * RYDBERG_TO_EV, 0
1632
            ):
1633
                warnings.warn(
1✔
1634
                    "job.input.basis.eCut was modified directly. "
1635
                    "It is recommended to set it via job.set_encut()"
1636
                )
1637

1638
            if round(self.input["Sigma"], 1) != round(
1✔
1639
                self.input.sphinx.PAWHamiltonian.ekt, 1
1640
            ):
1641
                warnings.warn(
1✔
1642
                    "job.input.PAWHamiltonian.ekt was modified directly. "
1643
                    "It is recommended to set it via "
1644
                    "job.set_occupancy_smearing()"
1645
                )
1646

1647
            if self.input["Xcorr"] != self.input.sphinx.PAWHamiltonian.xc:
1✔
1648
                warnings.warn(
1✔
1649
                    "job.input.PAWHamiltonian.xc was modified directly. "
1650
                    "It is recommended to set it via "
1651
                    "job.exchange_correlation_functional()"
1652
                )
1653

1654
            if (
1✔
1655
                self.input["EmptyStates"]
1656
                != self.input.sphinx.PAWHamiltonian.nEmptyStates
1657
            ):
1658
                warnings.warn(
1✔
1659
                    "job.input.PAWHamiltonian.nEmptyStates was modified "
1660
                    "directly. It is recommended to set it via "
1661
                    "job.set_empty_states()"
1662
                )
1663

1664
            if (
1✔
1665
                "KpointCoords" in self.input
1666
                and np.array(self.input.KpointCoords).tolist()
1667
                != np.array(self.input.sphinx.basis.kPoint.coords).tolist()
1668
            ) or (
1669
                "KpointFolding" in self.input
1670
                and np.array(self.input.KpointFolding).tolist()
1671
                != np.array(self.input.sphinx.basis.folding).tolist()
1672
            ):
1673
                warnings.warn(
1✔
1674
                    "job.input.basis.kPoint was modified directly. "
1675
                    "It is recommended to set all k-point settings via "
1676
                    "job.set_kpoints()"
1677
                )
1678

1679
            structure_sync = str(self.input.sphinx.structure) == str(
1✔
1680
                self.get_structure_group()
1681
            )
1682
            if not structure_sync and not self.input.sphinx.structure.read_only:
1✔
1683
                warnings.warn(
1✔
1684
                    "job.input.structure != job.structure. "
1685
                    "The current job.structure will overwrite "
1686
                    "any changes you may might have made to "
1687
                    "job.input.structure in the meantime. "
1688
                    "To disable this overwrite, "
1689
                    "set job.input.structure.read_only = True. "
1690
                    "To disable this warning, call "
1691
                    "job.load_structure_group() after making changes "
1692
                    "to job.structure."
1693
                )
1694

1695
            if len(w) > 0:
1✔
1696
                print("WARNING:")
1✔
1697
                for ww in w:
1✔
1698
                    print(ww.message)
1✔
1699
                return False
1✔
1700
            else:
1701
                return True
1✔
1702

1703
    def compress(self, files_to_compress=None):
1✔
1704
        """
1705
        Compress the output files of a job object.
1706

1707
        Args:
1708
            files_to_compress (list): A list of files to compress (optional)
1709
        """
1710
        # delete empty files
1711
        if files_to_compress is None:
1✔
1712
            files_to_compress = [
1✔
1713
                f
1714
                for f in self.files.list()
1715
                if (
1716
                    f not in ["rho.sxb", "waves.sxb"]
1717
                    and not stat.S_ISFIFO(
1718
                        os.stat(os.path.join(self.working_directory, f)).st_mode
1719
                    )
1720
                )
1721
            ]
1722
        for f in self.files.list():
1✔
1723
            filename = os.path.join(self.working_directory, f)
1✔
1724
            if (
1✔
1725
                f not in files_to_compress
1726
                and os.path.exists(filename)
1727
                and os.stat(filename).st_size == 0
1728
            ):
1729
                os.remove(filename)
×
1730
        super(SphinxBase, self).compress(files_to_compress=files_to_compress)
1✔
1731

1732
    @staticmethod
1✔
1733
    def check_vasp_potentials():
1✔
1734
        return any(
1✔
1735
            [
1736
                os.path.exists(
1737
                    os.path.join(p, "vasp", "potentials", "potpaw", "Fe", "POTCAR")
1738
                )
1739
                for p in state.settings.resource_paths
1740
            ]
1741
        )
1742

1743
    def run_addon(
1✔
1744
        self,
1745
        addon,
1746
        args=None,
1747
        from_tar=None,
1748
        silent=False,
1749
        log=True,
1750
        version=None,
1751
        debug=False,
1752
    ):
1753
        """Run a SPHInX addon
1754

1755
        addon          - name of addon (str)
1756
        args           - arguments (str or list)
1757
        from_tar       - if job is compressed, extract these files (list)
1758
        silent         - do not print output for successful runs?
1759
        log            - produce log file?
1760
        version        - which sphinx version to load (str or None)
1761
        debug          - return subprocess.CompletedProcess ?
1762

1763
        """
1764
        if self.is_compressed() and from_tar is None:
1✔
1765
            raise FileNotFoundError(
1✔
1766
                "Cannot run add-on on compressed job without 'from_tar' parameter.\n"
1767
                + "   Solution 1: Run .decompress () first.\n"
1768
                + "   Solution 2: specify from_tar list to run in temporary directory\n"
1769
                + "   Solution 3: run with from_tar=[] if no files from tar are needed\n"
1770
            )
1771

1772
        # prepare argument list
1773
        if args is None:
1✔
1774
            args = ""
1✔
1775
        elif isinstance(args, list):
1✔
1776
            args = " ".join(args)
1✔
1777
        if log:
1✔
1778
            args += " --log"
1✔
1779

1780
        cmd = addon + " " + args
1✔
1781

1782
        # --- handle versions
1783
        sxv = sxversions()
1✔
1784
        if (
1✔
1785
            version is None
1786
            and self.executable.version is not None
1787
            and self.executable.version in sxv.keys()
1788
        ):
1789
            version = self.executable.version
×
1790
            if not silent:
×
1791
                print("Taking version '" + version + "' from job._executable version")
×
1792
        if isinstance(version, str):
1✔
1793
            if version in sxv.keys():
1✔
1794
                cmd = sxv[version] + " && " + cmd
1✔
1795
            elif version != "":
1✔
1796
                raise KeyError(
1✔
1797
                    "version '"
1798
                    + version
1799
                    + "' not found. Available versions are: '"
1800
                    + "', '".join(sxv.keys())
1801
                    + "'."
1802
                )
1803
            # version="" overrides job.executable_version
1804
        elif version is not None:
1✔
1805
            raise TypeError("version must be str or None")
1✔
1806

1807
        if isinstance(from_tar, str):
1✔
1808
            from_tar = [from_tar]
1✔
1809
        if self.is_compressed() and isinstance(from_tar, list):
1✔
1810
            # run addon in temporary directory
1811
            with TemporaryDirectory() as tempd:
1✔
1812
                if not silent:
1✔
1813
                    print("Running {} in temporary directory {}".format(addon, tempd))
1✔
1814

1815
                # --- extract files from list
1816
                # note: tf should be obtained from JobCore to ensure encapsulation
1817
                tarfilename = os.path.join(
1✔
1818
                    self.working_directory, self.job_name + ".tar.bz2"
1819
                )
1820
                with tarfile.open(tarfilename, "r:bz2") as tf:
1✔
1821
                    for file in from_tar:
1✔
1822
                        try:
1✔
1823
                            tf.extract(file, path=tempd)
1✔
1824
                        except:
1✔
1825
                            print("Cannot extract " + file + " from " + tarfilename)
1✔
1826

1827
                # --- link other files
1828
                linkfiles = []
1✔
1829
                for file in self.files.list():
1✔
1830
                    linkfile = os.path.join(tempd, file)
1✔
1831
                    if not os.path.isfile(linkfile):
1✔
1832
                        os.symlink(
1✔
1833
                            os.path.join(self.working_directory, file),
1834
                            linkfile,
1835
                            target_is_directory=True,
1836
                        )
1837
                        linkfiles.append(linkfile)
1✔
1838
                # now run
1839
                out = subprocess.run(
1✔
1840
                    cmd, cwd=tempd, shell=True, stdout=PIPE, stderr=PIPE, text=True
1841
                )
1842

1843
                # now clean tempdir
1844
                for file in from_tar:
1✔
1845
                    try:
1✔
1846
                        os.remove(os.path.join(tempd, file))
1✔
1847
                    except FileNotFoundError:
1✔
1848
                        pass
1✔
1849
                for linkfile in linkfiles:
1✔
1850
                    if os.path.islink:
1✔
1851
                        os.remove(linkfile)
1✔
1852

1853
                # move output to working directory for successful runs
1854
                if out.returncode == 0:
1✔
1855
                    for file in os.listdir(tempd):
1✔
1856
                        movefile(os.path.join(tempd, file), self.working_directory)
1✔
1857
                        if not silent:
1✔
1858
                            print("Copying " + file + " to " + self.working_directory)
×
1859
                else:
1860
                    print(addon + " crashed - potential output files are not kept.")
1✔
1861

1862
        else:
1863
            out = subprocess.run(
1✔
1864
                cmd,
1865
                cwd=self.working_directory,
1866
                shell=True,
1867
                stdout=PIPE,
1868
                stderr=PIPE,
1869
                text=True,
1870
            )
1871
            if out.returncode != 0:
1✔
1872
                print(addon + " crashed.")
×
1873

1874
        # print output
1875
        if not silent or out.returncode != 0:
1✔
1876
            if out.returncode != 0:
1✔
1877
                print(addon + " output:\n\n")
1✔
1878
            print(out.stdout)
1✔
1879
            if out.returncode != 0:
1✔
1880
                print(out.stderr)
1✔
1881
        return out if debug else None
1✔
1882

1883

1884
class Output:
1✔
1885
    """
1886
    Handles the output from a SPHInX simulation.
1887
    """
1888

1889
    def __init__(self, job):
1✔
1890
        self._job = job
1✔
1891
        self.generic = DataContainer(table_name="output/generic")
1✔
1892
        self.charge_density = SphinxVolumetricData()
1✔
1893
        self.electrostatic_potential = SphinxVolumetricData()
1✔
1894
        self.generic.create_group("dft")
1✔
1895
        self.old_version = False
1✔
1896

1897
    def collect_spins_dat(self, file_name="spins.dat", cwd=None):
1✔
1898
        """
1899

1900
        Args:
1901
            file_name:
1902
            cwd:
1903

1904
        Returns:
1905

1906
        """
1907
        try:
1✔
1908
            results = collect_spins_dat(
1✔
1909
                file_name=file_name, cwd=cwd, index_permutation=self._job.id_spx_to_pyi
1910
            )
1911
        except FileNotFoundError:
1✔
1912
            return
1✔
1913
        for k, v in results.items():
1✔
1914
            self.generic.dft[k] = v
1✔
1915

1916
    def collect_energy_dat(self, file_name="energy.dat", cwd=None):
1✔
1917
        """
1918

1919
        Args:
1920
            file_name:
1921
            cwd:
1922

1923
        Returns:
1924

1925
        """
1926
        try:
1✔
1927
            results = collect_energy_dat(file_name=file_name, cwd=cwd)
1✔
1928
        except FileNotFoundError:
1✔
1929
            return
1✔
1930
        for k, v in results.items():
1✔
1931
            self.generic.dft[k] = v
1✔
1932

1933
    def collect_residue_dat(self, file_name="residue.dat", cwd=None):
1✔
1934
        """
1935

1936
        Args:
1937
            file_name:
1938
            cwd:
1939

1940
        Returns:
1941

1942
        """
1943
        try:
1✔
1944
            results = collect_residue_dat(file_name=file_name, cwd=cwd)
1✔
1945
        except FileNotFoundError:
1✔
1946
            return
1✔
1947
        for k, v in results.items():
1✔
1948
            self.generic.dft[k] = v
1✔
1949

1950
    def collect_eps_dat(self, file_name=None, cwd=None):
1✔
1951
        """
1952

1953
        Args:
1954
            file_name:
1955
            cwd:
1956

1957
        Returns:
1958

1959
        """
1960
        try:
1✔
1961
            results = collect_eps_dat(
1✔
1962
                file_name=file_name, cwd=cwd, spins=self._job._spin_enabled
1963
            )
1964
            for k, v in results.items():
1✔
1965
                self.generic.dft[k] = v
1✔
1966
        except FileNotFoundError:
1✔
1967
            return
1✔
1968

1969
    def collect_energy_struct(self, file_name="energy-structOpt.dat", cwd=None):
1✔
1970
        """
1971

1972
        Args:
1973
            file_name:
1974
            cwd:
1975

1976
        Returns:
1977

1978
        """
1979
        try:
1✔
1980
            results = collect_energy_struct(file_name=file_name, cwd=cwd)
1✔
1981
        except FileNotFoundError:
1✔
1982
            return
1✔
1983
        for k, v in results.items():
1✔
1984
            self.generic.dft[k] = v
1✔
1985

1986
    def collect_sphinx_log(self, file_name="sphinx.log", cwd=None):
1✔
1987
        """
1988

1989
        Args:
1990
            file_name:
1991
            cwd:
1992

1993
        Returns:
1994

1995
        """
1996

1997
        try:
1✔
1998
            results = SphinxLogParser(
1✔
1999
                file_name=file_name, cwd=cwd, index_permutation=self._job.id_spx_to_pyi
2000
            ).results
2001
        except FileNotFoundError:
1✔
2002
            return None
1✔
2003
        if len(results) == 0:
1✔
2004
            self._job.status.aborted = True
1✔
2005
            return None
1✔
2006
        if not results["generic"].pop("job_finished"):
1✔
2007
            self._job.status.aborted = True
×
2008
        for key, value in results["generic"].items():
1✔
2009
            if key not in self.generic:
1✔
2010
                self.generic[key] = value
1✔
2011
        for key, value in results["dft"].items():
1✔
2012
            if key not in self.generic.dft:
1✔
2013
                self.generic.dft[key] = value
1✔
2014

2015
    def collect_relaxed_hist(self, file_name="relaxHist.sx", cwd=None):
1✔
2016
        """
2017

2018
        Args:
2019
            file_name:
2020
            cwd:
2021

2022
        Returns:
2023

2024
        """
2025
        try:
1✔
2026
            results = collect_relaxed_hist(
1✔
2027
                file_name=file_name, cwd=cwd, index_permutation=self._job.id_spx_to_pyi
2028
            )
2029
        except FileNotFoundError:
1✔
2030
            return
1✔
2031
        for k, v in results.items():
1✔
2032
            self.generic[k] = v
1✔
2033

2034
    def collect_charge_density(self, file_name, cwd):
1✔
2035
        if (
1✔
2036
            file_name in os.listdir(cwd)
2037
            and os.stat(posixpath.join(cwd, file_name)).st_size != 0
2038
        ):
2039
            self.charge_density.from_file(
1✔
2040
                filename=posixpath.join(cwd, file_name), normalize=True
2041
            )
2042

2043
    def collect_electrostatic_potential(self, file_name, cwd):
1✔
2044
        if (
1✔
2045
            file_name in os.listdir(cwd)
2046
            and os.stat(posixpath.join(cwd, file_name)).st_size != 0
2047
        ):
2048
            self.electrostatic_potential.from_file(
1✔
2049
                filename=posixpath.join(cwd, file_name), normalize=False
2050
            )
2051

2052
    def _get_electronic_structure_object(self):
1✔
2053
        es = ElectronicStructure()
1✔
2054
        # This line won't be necessary in the newer version of sphinx parser
2055
        if len(self.generic.dft.bands_eigen_values) == 0:
1✔
2056
            return es
×
2057
        eig_mat = self.generic.dft.bands_eigen_values[-1]
1✔
2058
        occ_mat = self.generic.dft.bands_occ[-1]
1✔
2059
        if len(eig_mat.shape) == 3:
1✔
2060
            es.eigenvalue_matrix = eig_mat
1✔
2061
            es.occupancy_matrix = occ_mat
1✔
2062
        else:
UNCOV
2063
            es.eigenvalue_matrix = np.array([eig_mat])
×
UNCOV
2064
            es.occupancy_matrix = np.array([occ_mat])
×
2065
        es.efermi = self.generic.dft.bands_e_fermi[-1]
1✔
2066
        es.n_spins = len(es.occupancy_matrix)
1✔
2067
        es.kpoint_list = self.generic.dft.kpoints_cartesian
1✔
2068
        es.kpoint_weights = self.generic.dft.bands_k_weights
1✔
2069
        es.generate_from_matrices()
1✔
2070
        return es
1✔
2071

2072
    def collect(self, directory=None):
1✔
2073
        """
2074
        The collect function, collects all the output from a SPHInX simulation.
2075

2076
        Args:
2077
            directory (str): the directory to collect the output from.
2078
        """
2079
        if directory is None:
1✔
2080
            directory = self._job.working_directory
×
2081
        self.collect_energy_struct(file_name="energy-structOpt.dat", cwd=directory)
1✔
2082
        self.collect_sphinx_log(file_name="sphinx.log", cwd=directory)
1✔
2083
        self.collect_energy_dat(file_name="energy.dat", cwd=directory)
1✔
2084
        self.collect_residue_dat(file_name="residue.dat", cwd=directory)
1✔
2085
        self.collect_eps_dat(file_name=None, cwd=directory)
1✔
2086
        self.collect_spins_dat(file_name="spins.dat", cwd=directory)
1✔
2087
        self.collect_relaxed_hist(file_name="relaxHist.sx", cwd=directory)
1✔
2088
        self.collect_electrostatic_potential(file_name="vElStat-eV.sxb", cwd=directory)
1✔
2089
        self.collect_charge_density(file_name="rho.sxb", cwd=directory)
1✔
2090

2091
    def to_hdf(self, hdf, force_update=False):
1✔
2092
        """
2093
        Store output in an HDF5 file
2094

2095
        Args:
2096
            hdf: HDF5 group
2097
            force_update(bool):
2098
        """
2099

2100
        def get_last(arr):
1✔
2101
            return np.array([vv[-1] for vv in arr])
1✔
2102

2103
        for k in self.generic.dft.list_nodes():
1✔
2104
            if "scf" in k and k != "scf_convergence":
1✔
2105
                self.generic.dft[k.replace("scf_", "")] = get_last(self.generic.dft[k])
1✔
2106
        if "energy_free" in self.generic.dft.list_nodes():
1✔
2107
            self.generic.energy_tot = self.generic.dft.energy_free
1✔
2108
            self.generic.energy_pot = self.generic.dft.energy_free
1✔
2109
        if "positions" not in self.generic.list_nodes():
1✔
2110
            self.generic.positions = np.array([self._job.structure.positions])
1✔
2111
        if (
1✔
2112
            "cells" not in self.generic.list_nodes()
2113
            and "cell" in self.generic.list_nodes()
2114
        ):
2115
            self.generic.cells = self.generic.cell
1✔
2116
        self.generic.to_hdf(hdf=hdf)
1✔
2117

2118
        with hdf.open("output") as hdf5_output:
1✔
2119
            if self.electrostatic_potential.total_data is not None:
1✔
2120
                self.electrostatic_potential.to_hdf(
1✔
2121
                    hdf5_output, group_name="electrostatic_potential"
2122
                )
2123
            if self.charge_density.total_data is not None:
1✔
2124
                self.charge_density.to_hdf(hdf5_output, group_name="charge_density")
1✔
2125
            if "bands_occ" in self.generic.dft:
1✔
2126
                try:
1✔
2127
                    es = self._get_electronic_structure_object()
1✔
2128
                    if len(es.kpoint_list) > 0:
1✔
2129
                        es.to_hdf(hdf5_output)
1✔
2130
                except IndexError:
1✔
2131
                    warnings.warn("Electronic structure parsing failed")
1✔
2132
            with hdf5_output.open("electronic_structure") as hdf5_es:
1✔
2133
                if "dos" not in hdf5_es.list_groups():
1✔
2134
                    hdf5_es.create_group("dos")
1✔
2135
                with hdf5_es.open("dos") as hdf5_dos:
1✔
2136
                    warning_message = " is not stored in SPHInX; use job.get_density_of_states instead"
1✔
2137
                    for k in ["energies", "int_densities", "tot_densities"]:
1✔
2138
                        hdf5_dos[k] = k + warning_message
1✔
2139

2140
    def from_hdf(self, hdf):
1✔
2141
        """
2142
        Load output from an HDF5 file
2143
        """
2144
        try:
×
2145
            self.generic.from_hdf(hdf=hdf)
×
2146
        except ValueError:
×
2147
            warnings.warn(
×
2148
                "You are using an old version of SPHInX output - update via job.update_sphinx()"
2149
            )
2150
            self.old_version = True
×
2151
            pass
×
2152

2153

2154
def _update_datacontainer(job):
1✔
2155
    job.output.generic.create_group("dft")
×
2156
    for node in job["output/generic/dft"].list_nodes():
×
2157
        job.output.generic.dft[node] = job["output/generic/dft"][node]
×
2158
    for node in job["output/generic"].list_nodes():
×
2159
        job.output.generic[node] = job["output/generic"][node]
×
2160
    job["output/generic"].remove_group()
×
2161
    job.output.generic.to_hdf(hdf=job.project_hdf5)
×
2162
    job.output.old_version = False
×
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc