diff --git a/LICENSE b/LICENSE new file mode 100644 index 00000000..5ca21d97 --- /dev/null +++ b/LICENSE @@ -0,0 +1,12 @@ +pseudo_dojo is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation, either version 2.1 of the License, or +(at your option) any later version. + +pseudo_dojo is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License along with pseudo_dojo. +If not, see . diff --git a/PROJECT.rst b/PROJECT.rst new file mode 100644 index 00000000..343bf3ad --- /dev/null +++ b/PROJECT.rst @@ -0,0 +1,197 @@ +Global long-term objectives: +---------------------------- + + #. To have, on the Web, sets of validated pseudopotentials, for the whole periodic table, + for different exchange-correlation functionals, with different possibilities of + semi-core electrons (e.g. for GW), different cut-off radii (e.g. high pressure application), + with an optimal cut-off energy. + + #. To have a Web portal for the generation/validation of new pseudopotentials. + +Intermediate objectives (more and more difficult steps): +-------------------------------------------------------- + + #. Robust generation of one pseudopotential, be given the atomic number + (for a given exchange-correlation functional, definition of semi-core electrons, + definition of cut-off radii). The cut-off energy does not matter at this stage. + + #. Computation of selected physical properties for selected systems, associated with one + given pseudopotential (automatic computation of the cut-off energy, computation of the + total energy, the interatomic distance, the lattice parameter of the elemental solid + and one oxide, also dimer). Results presented on the Web. + + #. Validation of pseudopotentials with respect to a reference. + + #. Semi-automatic improvement of the generation of the pseudopotential: + more accurate physical properties, better energy cut-off. + + #. Automatic generation of one set of pseudopotentials, and associated automatic procedure + of calculation, validation. + + #. The transfer to the Web of the different results of the objectives. + +Components (database): +---------------------- + + #. A set of repositories (one for each atomic number). Repository is placed in the + directory where Z is the atomic number (three digits) followed by the symbol, + with a dash in between, e.g. 001-H, or 092-U. + Each of these repository is an "ATOM repository", and will contain subdirectories of two + kinds, see sections below. + + #. Structuration inside one ATOM repository: + + - A ReferenceData subdirectory, that is not tight to one pseudopotential. + + - For each pseudopotential: + /[PAW|NC|...]/num_valence_electrons/xc_type/ID , which might called a "Pseudo-atom Box" or patbox. + + The num_valence_electrons might be 4e, or 22e ... + The xc_type might be GGA-PBE or LDA-PW91 or the libcx ID: ex. LIBXC-X001-C012. + ID will be a digit, e.g. 1, or 2, etc ... + These IDs will not have any predefined meaning. Some of the pseudo-atom boxes might be good + for a specific purpose (e.g. GW or High-pressure), but this will be determined from the + database of results for each pseudo-atom box, by a script, at the demand of one user + (or for populating a Web page). + + #. Content of the ReferenceData subdirectory of the "ATOM repository": + + - A (xml, json?) file with the atomic configuration for each possible + num_valence_electrons, and other data needed for pseudo-atom generators that are not + specific to a pseudo-atom generator. Standard name: atomic_config.xml. + + - Possibly, some CIF files for a elemental solid (or more than one), and for oxide(s) + or hydride(s), or potassium-based compounds. + + - A set of master data file (xml) containing the description of the different test systems + for the specific atom. Each test system belongs to a test system class: + + - atom + + - dimer + + - elemental + + - oxide + + - hydride or potassium-based compound. + + Within each class it is labelled with an ID (number starting from 1). + This description contains insulator/metal and magnetism information, and either + the name of the cif file to be used, or the reference length for the dimer. + Standard name: .description__ID.xml + + #. Content of one pseudo-atom Box: + + - Subdirectories: + + - + + - atom_X, dimer_X + + - elemental_X, oxide_X + + - possibly hydride_X or K_X . + + Where might be atompaw, or ape, or fhi98pp, or ... + And where X is the ID defined in D3c. + + - A (xml) summary file containing metadata concerning this pseudo-atom box, obtained by + running the applications in the different subdirectories, and also describing the + validation criteria (this implies a set of runs). + Standard name: .summary.[PAW|NC|...].#valence_electrons..ID.xml + + #. Content of the subdirectory of the pseudo-atom Box: + + - Optionally, the specific input data needed for the generator (PAW or NC), to complement + the content of the atomic_config file. Typically cut-off radii. + Standard name atomic_data_.xml, e.g. atomic_data_atompaw.xml + + - A pseudo-atom data generator input file (PAW or NC) - might have been automatically generated + from atomc_config and the file in D5a. Standard name .in + + - A pseudo-atom data file (PAW or NC) - has been automatically generated (output of the atomic generator). + Standard name .pseudoatom.[PAW|NC|...].#valence_electrons..ID. + + - + This is the pseudopotential file, or the PAW atomic data file. + The might be .fhi or .pawps , or other postfix. + + #. Content of the _X (where class is atom, dimer, elemental, oxide, hydride ...) + D6a Subdirectories abinit_runY and elk_runY, where Y is an integer starting from 1. + + #. Content of the abinit_runY directory + D7a This is a working directory for one abinit run. It contains an ABINIT input file + usually automatically generated from D3c, specialized for the pseudo-atom box and the system. + D7b For Y=1 : determination of a basic k point grid, using kptrlen and prtkpt. + Can be used by elk, see elk_1. + D7c For Y=2 : computation of total energy as a function of ecut, for basic k point grid, + and, for metals, using the tsmear determined by elk_1. + + #. Content of the elk_runY subdirectory: + + - This is a working directory for one elk run. It contains an ELK input file + usually automatically generated from D3c, specialized for the pseudo-atom box + and the system, and using the k point grid determined by D7b . + + - For Y=1, determination of the tsmear. + +Components (software). +---------------------- + +They should be placed inside ABINIT package psps/script, for testing/coherency purposes +across the different directories. + + #. A "pseudo-atom box" creator (init_patbox.py), to be called inside the psps/ directory. + (propose options for the path described in D2, then create the path, + and the directories of the pseudo-atom box, and also bzr add the dirs) + + #. A cif2cml translator, to go from D3b to D3c. + + #. A script to initialize the file .in mentioned in D5b from D3a + atomic_config.xml and D5a atomic_data_.xml + + #. A pseudopotential generator, e.g. ATOMPAW (already placed inside the ABINIT package) + + #. A driver of abinit: generation of abinit input files, running of abinit, gathering of the data in D4b. + ACTUALLY NEED A LIST OF TASKS / VALIDATION CRITERIA / to be defined. + + #. A driver of elk, and a binary for elk. + + #. A validator. + + #. More scripts to be added ... + +Miscellaneous +------------- + + #. Reference oxygen PAW data files for different XC functionals, reference hydrogen PAW data files + for different XC functionals. Placed in the abinit/psps/RefPseudoAtoms subdirectory of + the ABINIT package. And to be copied in the patbox at init time. + +Strategy +-------- + + - Work component by component, by placing these components under version management and + automatic testing, with appropriate hardware. + + - Define the files and their format (including metadata) in an iterative way, with possibilities + to regenerate them in an automatic way + + - Gradual understanding of the CPU constraints, memory constraint, and human time needed. + + - Adjust the objectives to stay realist. + +TO BE KEPT IN MIND FOR FURTHER SPECIFICATION +-------------------------------------------- + + - set up of a bot (on the machine nazgul): be given the ABINIT branch, and the pseudopotentials => + computation of the physical characteristics of this pseudopotential + + - set up of the corresponding "on-demand" mechanism + + - set up a new waterfall: the list of files that will be provided will be quite different + from the usual bots + + - set up of a new Web window to visualize the files (to be discussed). + diff --git a/README.md b/README.md deleted file mode 100644 index 1671b213..00000000 --- a/README.md +++ /dev/null @@ -1,2 +0,0 @@ -pseudo_dojo -=========== diff --git a/README.rst b/README.rst new file mode 100644 index 00000000..bfb345bd --- /dev/null +++ b/README.rst @@ -0,0 +1,125 @@ + +Pseudo_dojo is a open-source Python framework for generaring and validating pseudo potentials. + +pseudo_dojo is free to use. However, we also welcome your help to improve this +library by making your own contributions. These contributions can be in the +form of additional tools or modules you develop, or even simple things such +as bug reports. Please report any bugs and issues at pseudo_dojo's `Github page +`_. + +Getting pseudo_dojo +================ + +Stable version +-------------- + +The version at the Python Package Index (PyPI) is always the latest stable +release that will be hopefully, be relatively bug-free. The easiest way to +install pseudo_dojo on any system is to use easy_install or pip, as follows:: + + easy_install pseudo_dojo + +or:: + + pip install pseudo_dojo + +Developmental version +--------------------- + +The developmental version is at the pseudo_dojo's `Github repo `_. +After cloning the source, you can type:: + + python setup.py install + +or to install the package in developmental mode:: + + python setup.py develop + +Requirements +============ + +All required dependencies should be automatically taken care of if you +install pseudo_dojo using easy_install or pip. Otherwise, these packages should +be available on `PyPI `_. + +1. Python 2.7+ required. +2. numpy +3. scipy 0.10+ +4. matplotlib 1.1+ +5. periodictable + +Using pseudo_dojo +================= + +Basic usage +----------- + +Useful aliases for commonly used objects are now provided, similar in style to +numpy. Supported objects include Element, Composition, Structure, Molecule, +Spin and Orbital. Here are some quick examples of the core capabilities and +objects: + +.. code-block:: pycon + + >>> import pseudo_dojo as mg + >>> + >>> si = mg.Element("Si") + >>> si.atomic_mass + 28.0855 + >>> si.melting_point + u'1687 K' + >>> + >>> comp = mg.Composition("Fe2O3") + >>> comp.weight + 159.6882 + >>> #Note that Composition conveniently allows strings to be treated just + >>> #like an Element object. + >>> comp["Fe"] + 2.0 + >>> comp.get_atomic_fraction("Fe") + 0.4 + >>> lattice = mg.Lattice.cubic(4.2) + >>> structure = mg.Structure(lattice, ["Cs", "Cl"], + ... [[0, 0, 0], [0.5, 0.5, 0.5]]) + >>> structure.volume + 74.088000000000008 + >>> structure[0] + PeriodicSite: Cs (0.0000, 0.0000, 0.0000) [0.0000, 0.0000, 0.0000] + >>> + >>> #Integrated symmetry tools from spglib. + >>> from pseudo_dojo.symmetry.finder import SymmetryFinder + >>> finder = SymmetryFinder(structure) + >>> finder.get_spacegroup_symbol() + 'Pm-3m' + >>> + >>> #Writing out a POSCAR file for VASP calculations. + >>> poscar = Poscar(structure) + >>> mg.write_structure(structure, "POSCAR") + >>> + >>> #Reading a structure from a file. + >>> structure = mg.read_structure("POSCAR") + +The above illustrates only the most basic capabilities of pseudo_dojo. + +.. note:: Examples + + A good way to explore the functionality of pseudo_dojo is to look at examples. + We have created a `Github wiki page + `_ to allow users to + share their Github gists (essentially mini git repos of scripts) + performing various kinds of functions with pseudo_dojo. Please feel free to + check them out and we welcome your contributions as well! + +=============================== +Project PseudoDojo: global view +=============================== + +.. include:: PROJECT.rst + +License +======= + +pseudo_dojo is released under the GPL License. The terms of the license are as follows: + +.. literalinclude:: LICENSE + diff --git a/pseudo_dojo/core/atom.py b/pseudo_dojo/core/atom.py index 2275bee5..42245dd5 100644 --- a/pseudo_dojo/core/atom.py +++ b/pseudo_dojo/core/atom.py @@ -1,3 +1,4 @@ +"""This module provides objects and helper functions for atomic calculations.""" from __future__ import division, print_function import os @@ -7,6 +8,7 @@ from pseudo_dojo.refdata.nist import nist_database from scipy.interpolate import UnivariateSpline +from scipy.integrate import cumtrapz __version__ = "0.1" __author__ = "Matteo Giantomassi" @@ -45,7 +47,7 @@ def _asl(obj): def states_from_string(confstr): """ - Parse a string with an atomic configuration and buil a list of `QState` instance. + Parse a string with an atomic configuration and build a list of `QState` instance. """ states = [] tokens = confstr.split() @@ -261,7 +263,7 @@ def to_dict(self): @classmethod def from_dict(cls, d): - """Reconstitute the object from a dict representation created using to_dict.""" + """Initialize the object from a dict generated with to_dict.""" return cls(d["Z"], d["states"]) ########################################################################################## @@ -294,17 +296,20 @@ def __iter__(self): def __getitem__(self, rslice): return self.rmesh[rslice], self.values[rslice] - def __str__(self): + def __repr__(self): return "<%s, name = %s>" % (self.__class__.__name__, self.name) - #def __add__(self): - #def __sub__(self): - #def __mul__(self): + def __str__(self): + return str(self) + + #def __add__(self, other): + #def __sub__(self, other): + #def __mul__(self, other): @classmethod def from_filename(cls, filename, rfunc_name=None, cols=(0,1)): """ - Initialize the object by reading data from filename (txt format) + Initialize the object reading data from filename (txt format). Args: filename: @@ -368,7 +373,6 @@ def derivatives(self, r): def integral(self, a=None, b=None): """ Return definite integral of the spline of (r**2 values**2) between two given points a and b - Args: a: First point. rmesh[0] if a is None @@ -413,7 +417,7 @@ def ir_small(self, abs_tol=0.01): class RadialWaveFunction(RadialFunction): """ Extends RadialFunction adding info on the set of quantum numbers. - and special methods for elecronic wavefunctions. + and methods specialized for electronic wavefunctions. """ TOL_BOUND = 1.e-10 @@ -425,12 +429,20 @@ def __init__(self, state, name, rmesh, values): def isbound(self): """True if self is a bound state.""" back = min(10, len(self)) - return all(abs(self.values[-back:]) < self.TOL_BOUND) + return np.all(np.abs(self.values[-back:]) < self.TOL_BOUND) + def r2f2_integral(self): + """ + Cumulatively integrate r**2 f**2(r) using the composite trapezoidal rule. + """ + integ = cumtrapz(self.rmesh**2 * self.values**2, x=self.rmesh) + pad_intg = np.zeros(len(self)) + pad_intg[1:] = integ + return pad_intg ########################################################################################## -def plot_aepp(ae_funcs, pp_funcs=None, rmax=None, **kwargs): +def plot_aepp(ae_funcs, pp_funcs=None, **kwargs): """ Uses Matplotlib to plot the radial wavefunction (AE only or AE vs PP) @@ -439,15 +451,13 @@ def plot_aepp(ae_funcs, pp_funcs=None, rmax=None, **kwargs): All-electron radial functions. pp_funcs: Pseudo radial functions. - rmax: - float or dictionary {state: rmax} where rmax is the maximum r (Bohr) that will be be plotted. ============== ============================================================== kwargs Meaning ============== ============================================================== title Title of the plot (Default: None). show True to show the figure (Default: True). - savefig: 'abc.png' or 'abc.eps'* to save the figure to a file. + savefig 'abc.png' or 'abc.eps'* to save the figure to a file. ============== ============================================================== Returns: @@ -471,20 +481,19 @@ def plot_aepp(ae_funcs, pp_funcs=None, rmax=None, **kwargs): spl_idx += 1 ax = fig.add_subplot(num_funcs, 1, spl_idx) - if spl_idx==1 and title: + if spl_idx == 1 and title: ax.set_title(title) lines, legends = [], [] - if rmax is None: - ir = len(ae_phi) + 1 - else: - try: - rm = rmax[state] - except TypeError: - rm = float(rmax) - - ir = ae_phi.ifromr(rm) + ir = len(ae_phi) + 1 + # TODO + #if ae_phi.isbound: + # norm = ae_phi.r2f2_integral() + # for ir, v in enumerate(norm): + # if v > 0.95: break + # else: + # raise ValueError("AE wavefunction is not normalized") line, = ax.plot(ae_phi.rmesh[:ir], ae_phi.values[:ir], "-b", linewidth=2.0, markersize=1) @@ -497,7 +506,7 @@ def plot_aepp(ae_funcs, pp_funcs=None, rmax=None, **kwargs): lines.append(line) legends.append("PP: %s" % state) - ax.legend(lines, legends, 'lower right', shadow=True) + ax.legend(lines, legends, loc="best", shadow=True) # Set yticks and labels. ylabel = kwargs.get("ylabel", None) @@ -534,18 +543,18 @@ def plot_logders(ae_logders, pp_logders, **kwargs): ============== ============================================================== title Title of the plot (Default: None). show True to show the figure (Default: True). - savefig: 'abc.png' or 'abc.eps'* to save the figure to a file. + savefig 'abc.png' or 'abc.eps'* to save the figure to a file. ============== ============================================================== Returns: `matplotlib` figure. """ + assert len(ae_logders) == len(pp_logders) title = kwargs.pop("title", None) show = kwargs.pop("show", True) savefig = kwargs.pop("savefig", None) import matplotlib.pyplot as plt - assert len(ae_logders) == len(pp_logders) fig = plt.figure() num_logds, spl_idx = len(ae_logders), 0 @@ -553,6 +562,9 @@ def plot_logders(ae_logders, pp_logders, **kwargs): for (state, pp_logd) in pp_logders.items(): spl_idx += 1 ax = fig.add_subplot(num_logds, 1, spl_idx) + + if spl_idx == 1 and title: + ax.set_title(title) lines, legends = [], [] @@ -566,7 +578,7 @@ def plot_logders(ae_logders, pp_logders, **kwargs): lines.append(line) legends.append("PP logder %s" % state) - ax.legend(lines, legends, 'lower left', shadow=True) + ax.legend(lines, legends, loc="best", shadow=True) if spl_idx == num_logds: ax.set_xlabel("Energy [Ha]") @@ -588,12 +600,12 @@ class Dipole(object): """This object stores the dipole matrix elements.""" TOL_LRULE = 0.002 - def __init__(self, istate, ostate, aeres, ppres): + def __init__(self, istate, ostate, ae_res, pp_res): self.istate = istate self.ostate = ostate - self.aeres = float(aeres) - self.ppres = float(ppres) - self.aempp = self.aeres - self.ppres + self.ae_res = float(ae_res) + self.pp_res = float(pp_res) + self.aempp = self.ae_res - self.pp_res @property def fulfills_lrule(self): diff --git a/pseudo_dojo/core/ncparams.py b/pseudo_dojo/core/ncparams.py new file mode 100644 index 00000000..0fd22e87 --- /dev/null +++ b/pseudo_dojo/core/ncparams.py @@ -0,0 +1,84 @@ +from __future__ import print_function, division + +from collections import namedtuple + +class NcProjector(namedtuple("NCPROJ", "n l rcut scheme")): + """ + Descriptor for norm-conserving projectors. + + .. attributes: + + n: + Principal quantum number (the one associated to the + AE wavefunction that has been pseudized. + l: + Angular momentum. + rcut: + Cutoff radius in Bohr. + scheme: + String defining the pseudization scheme. + """ + @classmethod + def asprojector(cls, obj): + if isinstance(obj, cls): + return obj + else: + return cls(**obj) + +#class NLCC(dict): + +class NcParams(object): + """ + This object gathers the parameters used to generate a norm-conserving pseudo. + + .. attributes: + + reference_conf: + String defining the reference configuration. + Z: + Nuclear Charge. + Z_val: + Number of valence electrons. + l_max: + Maximum angular momentum used in the pseudization procedures. + l_local: + Angular momemtum used for the local part. + projectors: + List of `NcProjector` instances. + nlcc: + Defines the treatment of the non-linear core-correction. + wave_equation: + String defining the type of Hamiltonian. Possible values: + (schrodinger, scalar-relativistic, relativistic). + xc_functional + String defining the exchange-correlation functional (APE notation). + """ + def __init__(self, reference_conf, Z, Z_val, l_max, l_local, projectors, + nlcc, wave_equation=None, xc_functional=None, **extra_kwargs): + + self.reference_conf = reference_conf + self.Z = Z + self.Z_val = Z_val + self.l_max = l_max + self.l_local = l_local + self.projectors = projectors + self.nlcc = nlcc + self.wave_equation = wave_equation + self.xc_functional = xc_functional + self._extra_kwargs = extra_kwargs + + @property + def has_nlcc(self): + return bool(self.nlcc) + + def set_wave_equation(self, wave_equation): + self.wave_equation = wave_equation + + def set_xc_functional(self, xc_functional): + self.xc_functional = xc_functional + + #@classmethod + #def from_dict(cls, d): + + #@property + #def to_dict(self): diff --git a/pseudo_dojo/ppcodes/ape/apecore.py b/pseudo_dojo/ppcodes/ape/apecore.py index a8e31cbb..810d209b 100644 --- a/pseudo_dojo/ppcodes/ape/apecore.py +++ b/pseudo_dojo/ppcodes/ape/apecore.py @@ -11,7 +11,8 @@ from pseudo_dojo.core import plot_aepp, plot_logders from pseudo_dojo.ppcodes.ape.apeio import (ape_read_waves, ape_read_potentials, ape_read_densities, - ape_read_logders, ape_read_dipoles, ape_check_ppeigen, ape_check_ghosts) + ape_read_logders, ape_read_dipoles, ape_check_ppeigen, ape_check_ghosts + ) __version__ = "0.1" @@ -47,7 +48,7 @@ def __init__(self, workdir, inpgen, verbose=0): #@classmethod #def from_aeconf(cls, workdir, aeconf, verbose=0): - #inpgen = ApeInputGenerator.aeconf.to_input() + #inpgen = ApeInputGenerator.aeconf.to_apeinput() #return cls(workdir, inpgen, verbose=verbose) @classmethod @@ -112,13 +113,13 @@ def ae_densities(self): self._ae_densities = ape_read_densities(self.aedir) return self._ae_densities - def show_input(self, stream=sys.stdout): + def show_apeinput(self, stream=sys.stdout): lines = ["AE INPUT".center(80,"*")] - lines += self.to_input() + lines += self.to_apeinput() lines += ["END AE INPUT".center(80,"*")] stream.writelines("\n".join(lines)+"\n") - def to_input(self): + def to_apeinput(self): return self.inpgen.get_strlist() def solve(self, **kwargs): @@ -129,6 +130,9 @@ def solve(self, **kwargs): The exist status of the subprocess. """ print("Solving the all-electron problem...") + if self.verbose: + self.show_apeinput() + start = time.time() # Write the input file. @@ -141,7 +145,7 @@ def solve(self, **kwargs): os.makedirs(self.workdir) with open(self.input_path, "w") as fh: - fh.writelines("\n".join(self.to_input())) + fh.writelines("\n".join(self.to_apeinput())) # Launch APE in a subprocess. command = "%s < %s > %s 2> %s" % (self.exe, self.input_path, self.output_path, self.stderr_path) @@ -324,13 +328,13 @@ def path_in_workdir(self, filename): """Generates the absolute path in the working directory.""" return os.path.join(self.workdir, filename) - def show_input(self, stream=sys.stdout): + def show_apeinput(self, stream=sys.stdout): lines = ["PP INPUT".center(80,"*")] - lines += self.to_input() + lines += self.to_apeinput() lines += ["END PP INPUT".center(80,"*")] stream.writelines("\n".join(lines)+"\n") - def to_input(self): + def to_apeinput(self): """Returns the APE input (list of strings)""" return self.inpgen.get_strlist() @@ -341,7 +345,10 @@ def pseudize(self, **kwargs): Returns: The exist status of the subprocess. """ - if self.verbose: print("Starting pseudization...") + print("Starting pseudization...") + if self.verbose: + self.show_apeinput() + start = time.time() # Write the input file. @@ -358,7 +365,7 @@ def pseudize(self, **kwargs): shutil.copytree(self.ae_solver.aedir, os.path.join(self.workdir, "ae")) with open(self.input_path, "w") as fh: - fh.writelines("\n".join(self.to_input())) + fh.writelines("\n".join(self.to_apeinput())) # Launch APE in a subprocess. command = "%s < %s > %s 2> %s" % (self.exe, self.input_path, self.output_path, self.stderr_path) @@ -383,18 +390,22 @@ def pseudize(self, **kwargs): # Check ghost-states self._check_ghosts() + if self.ghosts: + print("Detected ghosts for states %s" % self.ghosts.keys()) # Check PP eigenvalues self._check_ppeigen() # Check logarithmic derivative. - return exit_code + # Generate pdf files. + #ape_plot_waves(dirpath, rmax=None, savefig=None) + #ape_plot_logders(dirpath, savefig=None) + #ape_plot_potentials(dirpath, rmax=rmax, savefig=None) - #def json_dump(self, protocol=0): - # path = os.path.join(self.workdir, self.json_fname) - # with open(path, "w") as fh: - # json.dump(self, fh, protocol=protocol) + # Generate final HTML report. + + return exit_code def _check_ppeigen(self): return ape_check_ppeigen(self.output) @@ -596,3 +607,17 @@ def get_strlist(self): lines += ["%s" % l] lines += ["%"] return lines + + +#def ape_pseudo_build(workdir, verbose=0): +# ae_solver = ApeAeSolver(workdir, inpgen, verbose=verbose) +# +# retcode = ae.solver(solve) +# if retcode != 0: +# raise ApeError("ae_solver returned %s" % retcode) +# +# pp_gen = ApePseudoGenerator(workdir, inpgen, ae_solver, verbose=verbose) +# +# retcode = pp_gen.pseudize() +# if retcode != 0: +# raise ApeError("pp_generator returned %s" % retcode) diff --git a/pseudo_dojo/ppcodes/ape/apeio.py b/pseudo_dojo/ppcodes/ape/apeio.py index 70c8ec24..bb0fae57 100644 --- a/pseudo_dojo/ppcodes/ape/apeio.py +++ b/pseudo_dojo/ppcodes/ape/apeio.py @@ -140,7 +140,7 @@ def __str__(self): lines = [self.filename+":",] for event in self: lines.append(str(event)) - return "\n".join(l for l in lines) + return "\n".join(lines) def append(self, event): """Add an event to the list.""" @@ -257,6 +257,33 @@ def ape_parse_events(out_lines): ########################################################################################## +#class ApePlotter(object): +# +# def build_figs(self, dirpath, **kwargs): +# """ +# Args: +# +# ============== ============================================================== +# kwargs Meaning +# ============== ============================================================== +# show True to show the figure (Default). +# ============== ============================================================== +# """ +# if "show" not in kwargs: +# kwargs["show"] = False +# +# self.figs = figs = {} +# +# figs["wave"] = ape_plot_waves(dirpath, **kwargs) +# +# figs["logd"] = ape_plot_logders(dirpath, **kwargs) +# +# figs["pot"] = ape_plot_potentials(dirpath, **kwargs) +# +# figs["den"] = ape_plot_densities(dirpath, **kwargs) + +########################################################################################## + def ape_read_waves(dirpath): """Read the APE radial wavefunctions located in directory dirpath.""" waves = {} @@ -275,14 +302,12 @@ def ape_read_waves(dirpath): return waves -def ape_plot_waves(dirpath, rmax=None, **kwargs): +def ape_plot_waves(dirpath, **kwargs): """ Uses Matplotlib to plot the radial wavefunction (AE vs PP) Args: dirpath: - rmax: - float or dictionary {state: rmax} where rmax is the maximum r (Bohr) that will be be plotted. ============== ============================================================== kwargs Meaning @@ -295,10 +320,6 @@ def ape_plot_waves(dirpath, rmax=None, **kwargs): Returns: `matplotlib` figure. """ - title = kwargs.pop("title", "Wavefunctions") - show = kwargs.pop("show", True) - savefig = kwargs.pop("savefig", None) - dirs = os.listdir(os.path.abspath(dirpath)) ae_waves, pp_waves = None, None @@ -310,10 +331,10 @@ def ape_plot_waves(dirpath, rmax=None, **kwargs): pp_waves = ape_read_waves(os.path.join(dirpath, "pp")) if ae_waves is not None: - fig = plot_aepp(ae_waves, pp_funcs=pp_waves, rmax=rmax, title=title, show=show, savefig=savefig) + fig = plot_aepp(ae_waves, pp_funcs=pp_waves, **kwargs) else: - fig = None print("Cannot find AE waves in dirpath %s" % dirpath) + fig = None return fig @@ -336,16 +357,13 @@ def ape_read_potentials(dirpath): return pots -def ape_plot_potentials(dirpath, rmax=None, **kwargs): +def ape_plot_potentials(dirpath, **kwargs): """ Uses Matplotlib to plot the potentials (AE vs PP) Args: dirpath: - rmax: - float or dictionary {state: rmax} where rmax is the maximum r (Bohr) that will be be plotted. - ============== ============================================================== kwargs Meaning ============== ============================================================== @@ -358,10 +376,6 @@ def ape_plot_potentials(dirpath, rmax=None, **kwargs): `matplotlib` figure. """ #raise NotImplementedError("") - title = kwargs.pop("title", "Potentials") - show = kwargs.pop("show", True) - savefig = kwargs.pop("savefig", None) - dirs = os.listdir(os.path.abspath(dirpath)) ae_pots, pp_pots = None, None @@ -373,10 +387,10 @@ def ape_plot_potentials(dirpath, rmax=None, **kwargs): pp_pots = ape_read_potentials(os.path.join(dirpath, "pp")) if ae_pots is not None: - fig = plot_aepp(ae_pots, pp_funcs=pp_pots, rmax=rmax, title=title, show=show, savefig=savefig) + fig = plot_aepp(ae_pots, pp_funcs=pp_pots, **kwargs) else: - fig = None print("Cannot find potentials in dirpath %s" % dirpath) + fig = None return fig @@ -399,15 +413,13 @@ def ape_read_densities(dirpath): return dens -def ape_plot_densities(dirpath, rmax=None, **kwargs): +def ape_plot_densities(dirpath, **kwargs): """ Uses Matplotlib to plot the densities (AE vs PP) Args: dirpath: Directory path. - rmax: - float or dictionary {state: rmax} where rmax is the maximum r (Bohr) that will be be plotted. ============== ============================================================== kwargs Meaning @@ -420,10 +432,6 @@ def ape_plot_densities(dirpath, rmax=None, **kwargs): Returns: `matplotlib` figure. """ - title = kwargs.pop("title", "Densities") - show = kwargs.pop("show", True) - savefig = kwargs.pop("savefig", None) - dirs = os.listdir(os.path.abspath(dirpath)) ae_dens, pp_dens = None, None @@ -435,10 +443,10 @@ def ape_plot_densities(dirpath, rmax=None, **kwargs): pp_dens = ape_read_densities(os.path.join(dirpath, "pp")) if ae_dens is not None: - fig = plot_aepp(ae_dens, pp_funcs=pp_dens, rmax=rmax, title=title, show=show, savefig=savefig) + fig = plot_aepp(ae_dens, pp_funcs=pp_dens, **kwargs) else: - fig = None print("Cannot find densities in dirpath %s" % dirpath) + fig = None return fig @@ -462,6 +470,7 @@ def ape_read_logders(dirpath): data = np.loadtxt(path) ae_logders[l_name] = RadialFunction(l_name, data[:,0], data[:,1]) pp_logders[l_name] = RadialFunction(l_name, data[:,0], data[:,2]) + return ae_logders, pp_logders @@ -484,18 +493,13 @@ def ape_plot_logders(dirpath, **kwargs): Returns: `matplotlib` figure. """ - title = kwargs.pop("title", "Logders") - show = kwargs.pop("show", True) - savefig = kwargs.pop("savefig", None) - dirpath = os.path.abspath(dirpath) if "tests" in os.listdir(dirpath): dirpath = os.path.join(dirpath, "tests") ae_logders, pp_logders = ape_read_logders(dirpath) - fig = plot_logders(ae_logders, pp_logders, show=show, title=title, savefig=savefig) - + fig = plot_logders(ae_logders, pp_logders, **kwargs) return fig ########################################################################################## diff --git a/pseudo_dojo/ppcodes/ape/apeobjects.py b/pseudo_dojo/ppcodes/ape/apeobjects.py index 9e06dafd..c747d264 100644 --- a/pseudo_dojo/ppcodes/ape/apeobjects.py +++ b/pseudo_dojo/ppcodes/ape/apeobjects.py @@ -19,9 +19,62 @@ ########################################################################################## +class ApeVariable(object): + + def __init__(self, name, type, doc, allowed_values=None): + self.name = name + self.type = type + self.doc = doc + self.allowed_values = allowed_values if allowed_values is not None else [] + + def set_value(self, value): + self._value = value + + @property + def value(self): + try: + return self._value + except AttributeError: + return None + + def isvalid(self): + if self.allowed_values: + return self.value in self.allowed_values + else: + return True + + def to_apeinput(self): + return "%s = %s" % (self.name, self.value) + + +class ApeBlock(object): + + def __init__(self, name, types, doc, allowed_values=None): + self.name = name + self.types = types + self.doc = doc + self.allowed_values = allowed_values if allowed_values is not None else [] + + self._table = [] + + def add_entry(self, entry): + row = [] + for i, item in enumerate(entry): + row.append(self.types[i](item)) + self._table.append(row) + + def to_apeinput(self): + lines = ["\%%s" % self.name] + for (state, rcut, scheme) in zip(self.states, self.core_radii, self.schemes): + lines += ["%s | %s | %s | %s" % (state.n, state.l, rcut, scheme)] + lines += ["%"] + return "\n".join(lines) + +########################################################################################## + class ApeAtomicConfiguration(AtomicConfiguration): - def to_input(self): + def to_apeinput(self): lines = ["NuclearCharge = %s " % self.Z] lines += ["SpinMode = %s" % self.spin_mode] lines += ["%Orbitals"] @@ -31,11 +84,11 @@ def to_input(self): return lines @classmethod - def from_input(cls, lines): + def from_apeinput(cls, lines): + if isinstance(lines, str): with open(lines, "r") as fh: lines = fh.readlines() - # TODO # Example (spin unpolarized, no spinor) #%Orbitals @@ -62,7 +115,7 @@ def from_input(cls, lines): print("s:",s) if "|" not in s: - # Handle noble configuration + # Handle noble gas. assert i == 0 s = s.translate(None, "'\"") noble_gas = ApeAtomicConfiguration.neutral_from_symbol(s) @@ -82,7 +135,7 @@ class ApeRadialMesh(dict): The radial mesh used by APE to represent radial functions. """ # Supported variables - _KEYS = [ + _VARS = [ "MeshType", "MeshStartingPoint", "MeshOutmostPoint", @@ -95,11 +148,11 @@ def __init__(self, *args, **kwargs): super(ApeRadialMesh, self).__init__(*args, **kwargs) for k in self: - if k not in self._KEYS: + if k not in self._VARS: raise ValueError("%s is not a registered key" % k) - def to_input(self): - return["%s = %s" % kv for kv in self.items()] + def to_apeinput(self): + return ["%s = %s" % kv for kv in self.items()] @property def to_dict(self): @@ -117,7 +170,7 @@ class ApeControl(dict): At each iteration the new guess potential is built mixing the input and output potentials. """ # Supported variables - _KEYS = [ + _VARS = [ # SCF "MaximumIter", "ConvAbsDens", @@ -143,11 +196,11 @@ def __init__(self, **kwargs): super(ApeControl, self).__init__(**kwargs) for k in self: - if k not in self._KEYS: + if k not in self._VARS: raise ValueError("%s is not a registered key" % k) - def to_input(self): - return["%s = %s" % kv for kv in self.items()] + def to_apeinput(self): + return ["%s = %s" % kv for kv in self.items()] @property def to_dict(self): @@ -197,10 +250,10 @@ def __init__(self, states, core_radii, schemes): self.core_radii = core_radii self.schemes = schemes - def to_input(self): + def to_apeinput(self): lines = ["%PPComponents"] for (state, rcut, scheme) in zip(self.states, self.core_radii, self.schemes): - lines += [" %s | %s | %s | %s " % (state.n, state.l, rcut, scheme)] + lines += ["%s | %s | %s | %s" % (state.n, state.l, rcut, scheme)] lines += ["%"] return lines @@ -213,8 +266,8 @@ def __init__(self, pp_components, core_correction=0, llocal=-1): self.core_correction = core_correction self.llocal = llocal - def to_input(self): - lines = ["# PseudoPotentials"] + def to_apeinput(self): + lines = ["# PseudoPotential Setup"] lines += ["CoreCorrection = %s" % self.core_correction] lines += ["Llocal = %s" % self.llocal] lines += self.pp_components.to_input() diff --git a/pseudo_dojo/pseudos/khein_lda_params.py b/pseudo_dojo/pseudos/khein_lda_params.py index 6752b7c1..f076e915 100644 --- a/pseudo_dojo/pseudos/khein_lda_params.py +++ b/pseudo_dojo/pseudos/khein_lda_params.py @@ -1240,98 +1240,8 @@ def get_params(symbol): return dict_params #return NcParams.from_dict(dict_params) -from collections import namedtuple - -class NcProjector(namedtuple("NCPROJ", "n l rcut scheme")): - """ - Descriptor for norm-conserving projectors. - - .. attributes: - - n: - Principal quantum number (the one associated to the - AE wavefunction that has been pseudized. - l: - Angular momentum. - rcut: - Cutoff radius in Bohr. - scheme: - String defining the pseudization scheme. - """ - @classmethod - def asprojector(cls, obj): - if isinstance(obj, cls): - return obj - else: - return cls(**obj) - -#class NLCC(dict): - -class NcParams(object): - """ - This object gathers the parameters used to generate a norm-conserving pseudo. - - .. attributes: - - reference_conf: - String defining the reference configuration. - Z: - Nuclear Charge. - Z_val: - Number of valence electrons. - l_max: - Maximum angular momentum used in the pseudization procedures. - l_local: - Angular momemtum used for the local part. - projectors: - List of `NcProjector` instances. - nlcc: - Defines the treatment of the non-linear core-correction. - wave_equation: - String defining the type of Hamiltonian. Possible values: - (schrodinger, scalar-relativistic, relativistic). - xc_functional - String defining the exchange-correlation functional (APE notation). - """ - def __init__(self, reference_conf, Z, Z_val, l_max, l_local, projectors, - nlcc, wave_equation=None, xc_functional=None, **extra_kwargs): - - self.reference_conf = reference_conf - self.Z = Z - self.Z_val = Z_val - self.l_max = l_max - self.l_local = l_local - self.projectors = projectors - self.nlcc = nlcc - self.wave_equation = wave_equation - self.xc_functional = xc_functional - self._extra_kwargs = extra_kwargs - - @property - def has_nlcc(self): - return bool(self.nlcc) - - def set_wave_equation(self, wave_equation): - self.wave_equation = wave_equation - - def set_xc_functional(self, xc_functional): - self.xc_functional = xc_functional - - #@classmethod - #def from_dict(cls, d): - - #@property - #def to_dict(self): - - -def ape_pseudize(symbol, wave_equation, xc_functional): - from pseudo_dojo.ppcodes.ape import ApeAtomicConfiguration, ApeAeSolver - params = get_params(symbol) - ae_conf = ApeAtomicConfiguration.from_string(params["Z"], params["reference_conf"]) - print(ae_conf.to_input()) - if __name__ == "__main__": - ape_pseudize("Si", "None", "None") + pass #from pseudo_dojo.ppcodes.nist import nist_database #from pprint import pprint #for k, d in _PARAMS.items(): diff --git a/pseudo_dojo/scripts/ape.py b/pseudo_dojo/scripts/ape.py index 8e76b0e3..efab39b9 100755 --- a/pseudo_dojo/scripts/ape.py +++ b/pseudo_dojo/scripts/ape.py @@ -3,18 +3,15 @@ import os import sys +import glob import argparse -from pseudo_dojo.ppcodes.nist import nist_database +from pseudo_dojo.refdata.nist import nist_database from pseudo_dojo.ppcodes.ape import * __version__ = "0.1" -__status__ = "Development" -__date__ = "$April 26, 2013M$" ########################################################################################## -# Helper functions. - def str_examples(): examples = """Example usage:\n @@ -24,7 +21,7 @@ def str_examples(): def show_examples_and_exit(err_msg=None, error_code=0): - "Display the usage of the script." + """Display the usage of the script.""" sys.stderr.write(str_examples()) if err_msg: sys.stderr.write("Fatal Error\n" + err_msg + "\n") @@ -32,13 +29,12 @@ def show_examples_and_exit(err_msg=None, error_code=0): def issymbol(string): - "True if string is a known element symbol" + """True if string is a known element symbol.""" return string in nist_database.allsymbols def show_nistdata(options): - "Handle nist command" - + """Handle nist command.""" for symbol in options.symbols: if symbol.endswith("+"): iontype = "Cation" @@ -56,19 +52,61 @@ def plot(options): "Handle plot command" dirpath = options.dirpath - rmax = options.rmax - if "w" in options.plots: - ape_plot_waves(dirpath, rmax=rmax, savefig=None) + ape_plot_waves(dirpath, savefig=None) if "l" in options.plots: ape_plot_logders(dirpath, savefig=None) if "p" in options.plots: - ape_plot_potentials(dirpath, rmax=rmax, savefig=None) + ape_plot_potentials(dirpath, savefig=None) if "d" in options.plots: - ape_plot_densities(dirpath, rmax=rmax, savefig=None) + ape_plot_densities(dirpath, savefig=None) + + +def analyze_llocal(options): + """Handle llocal command.""" + # Create an InputGenerator to facilitate the modification of the input file. + inpgen = ApeInputGenerator.from_template(options.template) + + # Init the AE solver from the template + ae_solver = ApeAeSolver(workdir, inpgen, verbose=verbose) + + if verbose: + ae_solver.show_input() + + if not dry_run: + ae_solver.solve(remove_wd=remove_wd) + + # Define the parameters for the pseudization. + # For each possible local L: + # 1) Pseudize. + # 2) Check ghosts + # 3) Plot wfs and logders + + pp_generators = [] + #for llocal in range(-1, 5, 1): + for llocal in range(0, 1, 1): + inpgen.reset() + inpgen.set_llocal(llocal) + #inpgen.set_ppcomponents() + #inpgen.set_correction() + + #pp_components = ApePPComponents.from_strings("3s|1.2|tm", "3p|1.27|tm", "3d|1.5|tm", "4f|1.9|tm") + #pp_setup = ApePPSetup(pp_components, core_correction=0, llocal=llocal) + + pp_workdir = os.path.join(workdir, "ppgen_loc%d" % llocal) + + pp_gen = ApePseudoGenerator(pp_workdir, inpgen, ae_solver, verbose=verbose) + + pp_generators.append(pp_gen) + + for pp_gen in pp_generators: + pp_gen.pseudize(remove_wd=remove_wd) + #if pp_gen.ghosts: + # print("Detected ghosts for states %s" % pp_gen.ghosts.keys()) + ########################################################################################## @@ -113,11 +151,14 @@ def main(): p_plot.add_argument('-p', '--plots', metavar='STRING', default= "w", help="Quantities to plot: " + "w for wavefunctions, l for logarithmic derivatives, p for potentials, d for densities (DEFAULT: w)") - p_plot.add_argument('-r', '--rmax', type=float, default=None, help="Max radius [Bohr] for plots.") - p_plot.add_argument('dirpath', metavar='DIRPATH', default = ".", help="Path to the directory containing APE output file") - # Access to the NIST database + # Analysis of LLocal option. + #p_llocal = subparsers.add_parser('llocal', help='Analyze the choice of the angular momentum for the local part') + + #p_llocal.add_argument('template', metavar='STRING', default= "", help="Path to the template file") + + # Access the NIST database p_nist = subparsers.add_parser('nist', help='Retrieve AE results from the NIST database') p_nist.add_argument('symbols', nargs="+", help="List of element symbols") @@ -131,7 +172,7 @@ def main(): dry_run = options.dry_run remove_wd = options.remove_wd - workdir = "helloape" + #if options.command == "autogen": # aconf_str = options.aconf_str @@ -171,12 +212,15 @@ def main(): # Create an InputGenerator to facilitate the modification of the input file. inpgen = ApeInputGenerator.from_template(options.template) + # Construct the name of the working directory. Ex: APERUN_1, APERUN_2 ... + prefix = "APERUN_" + nums = [int(dirname[len(prefix):]) for dirname in glob.glob(prefix + "*")] + max_num = max(nums) if nums else 0 + workdir = prefix + str(max_num+1) + # Init the AE solver from the template ae_solver = ApeAeSolver(workdir, inpgen, verbose=verbose) - - if verbose: - ae_solver.show_input() - + if not dry_run: ae_solver.solve(remove_wd=remove_wd) @@ -187,7 +231,8 @@ def main(): # 3) Plot wfs and logders pp_generators = [] - for llocal in range(-1, 5, 1): + #for llocal in range(-1, 5, 1): + for llocal in range(0, 1, 1): inpgen.reset() inpgen.set_llocal(llocal) @@ -204,19 +249,11 @@ def main(): pp_generators.append(pp_gen) for pp_gen in pp_generators: - if verbose: pp_gen.show_input() - pp_gen.pseudize(remove_wd=remove_wd) - if pp_gen.ghosts: - print("Detected ghosts for states %s" % pp_gen.ghosts.keys()) - - #print("LOGDER merit factors:") - #pprint(pp_gen.check_logders()) - #pprint(pp_gen.dipoles) - #pp_gen.plot_waves(savefig=) - #pp_gen.plot_logders(savefig=) + if options.command == "llocal": + analyze_llocal(options) if options.command == "plot": plot(options) diff --git a/pseudo_dojo/scripts/ppdojo_run.py b/pseudo_dojo/scripts/ppdojo_run.py index 5d84cda6..0db4cf0b 100755 --- a/pseudo_dojo/scripts/ppdojo_run.py +++ b/pseudo_dojo/scripts/ppdojo_run.py @@ -17,7 +17,7 @@ def str_examples(): examples = """ - Example usage: + Usage Example: \n ppdojo_run.py Si.fhi -m10 => Cutoff converge test for Si.fhi using up to 10 CPUs. ppdojo_run.py Si.fhi -l1 -n10 -m20 => Deltafactor test using up to 20 CPUs, each task uses 10 MPI nodes. diff --git a/setup.py b/setup.py index 5cd7b170..e8e23950 100755 --- a/setup.py +++ b/setup.py @@ -90,7 +90,8 @@ def cleanup(): "scipy>=0.10", "matplotlib>=1.1", "pymatgen>=2.7.3", -#"periodictable", +"periodictable", +#"django", ] #---------------------------------------------------------------------------