pyiron_atomistics icon indicating copy to clipboard operation
pyiron_atomistics copied to clipboard

No JTH potential for Si

Open samwaseda opened this issue 1 year ago • 27 comments

from pyiron_atomistics import Project
pr = Project("TEST")
job = pr.create.job.Sphinx("test")
job.structure = pr.create.structure.bulk("Si")  
job.input["VaspPot"] = False
job.run()

This raises KeyError 'Si'. The reason is simply because there is no JTH Si potential in the resources. Who created the potentials? @jan-janssen

samwaseda avatar Nov 09 '23 17:11 samwaseda

The potentials are packaged in https://github.com/eisenforschung/sphinx-data

jan-janssen avatar Nov 09 '23 17:11 jan-janssen

The potentials are packaged in https://github.com/eisenforschung/sphinx-data

Ok so I see that Si is indeed missing

samwaseda avatar Nov 09 '23 17:11 samwaseda

@freyso Do you remember if there was a specific reason not to include the Si potential?

jan-janssen avatar Nov 09 '23 17:11 jan-janssen

If I remember correctly the potentials were originally downloaded from https://www.abinit.org/downloads/PAW2

jan-janssen avatar Nov 10 '23 12:11 jan-janssen

I took a look but the data is given in XML (while our data isn't). Was there a converter? Otherwise I'm gonna do it manually

samwaseda avatar Nov 10 '23 14:11 samwaseda

It should work with http://users.wfu.edu/natalie/papers/pwpaw/man.html if I remember correctly, this can be installed from conda https://anaconda.org/conda-forge/atompaw

jan-janssen avatar Nov 10 '23 14:11 jan-janssen

If I remember correctly, Si was removed because it was pathological (PAW-corrected norm not garantueed to be positive). I do not remember which version of the PAW set this was, maybe there is an updated version by now. If this is still the same potential, it is highly likely that this is broken.

Regarding xml format: sphinx is able to read xml PAW format since version 3.0.5 (March 2022). However, it may need a fix in the pyiron-sphinx interface to enable that.

freyso avatar Nov 10 '23 15:11 freyso

And the way I got the PAW potentials in atomicdata format (which is the intrinsic format of atompaw): I downloaded the xml files, extracted the input file for the atompaw generator which is embedded as a comment in the end, and reran atompaw.

freyso avatar Nov 10 '23 15:11 freyso

@freyso Can you share the syntax to load the xml files? Currently, we use:

species {
    name = "Fe";
    potType = "AtomPAW";
    element = "Fe";
    potential = "Fe_GGA.atomicdata";
}

jan-janssen avatar Nov 10 '23 15:11 jan-janssen

potType = "xml"; see https://sxrepo.mpie.de/attachments/download/68/sphinx-manual-3.0.5.pdf Section 3.3, Table at bottom of page 13

freyso avatar Nov 10 '23 15:11 freyso

but does it mean we should not simply download the Si file and start using it?

samwaseda avatar Nov 10 '23 15:11 samwaseda

Is there actually an established workflow that can be used to test pseudo potentials?

samwaseda avatar Nov 10 '23 15:11 samwaseda

No. Mostafa Joulaian did some work in his master thesis analogous to the Delta project, but it turned out to be quite a mess in the end, in part because the delta project lacked a useful database for structures and parameters (cutoff/k-points/etc.), in part because too many cases needed some finetuning in the setup (e.g. convergence parameters).

It had been my hope to learn from Mostafa's master project how to generalize, but in the end, I felt it was no good basis for a potential testing according to my expectations.

freyso avatar Nov 10 '23 16:11 freyso

but does it mean we should not simply download the Si file and start using it?

bare sphinx can. pyiron_atomistics will need a fix, or a manual fixup of the pawPot group (+ adding the pot file to the input files, no clue how this works).

freyso avatar Nov 10 '23 16:11 freyso

bare sphinx can. pyiron_atomistics will need a fix, or a manual fixup of the pawPot group (+ adding the pot file to the input files, no clue how this works).

Sorry I meant we should not blindly trust the latest Abinit Si potential?

samwaseda avatar Nov 10 '23 16:11 samwaseda

Better not. I am sure there was a reason why I excluded it at the time, but this is min. 3 years ago, so I don't remember the details.

freyso avatar Nov 10 '23 16:11 freyso

ok from the pyiron developer perspective that presents a huge problem, because that means there's no way to do Si calculations in the open source version. Do you see a possibility to suggest a quick workflow? We'll take care of the actual implementation.

samwaseda avatar Nov 10 '23 16:11 samwaseda

I found something - I generated in 2018 a directory named JTH-1.1, and it contains a list of PAW .atomicdata files that looks suspiciously similar to the one in pyiron.

The README reads:

This is the JTH PBE 1.1 data set as downloaded from the Abinit homepage 2019-09-13

The potentials have been regenerated from the atompaw input files with atompaw 4.1.0.6 (compiled locally from tarball) linked with libxc 3.0.0-1build1 (from Ubuntu) as follows (using bash shell)

cd JTH-PBE-atomicdata-1.1 mkdir work JTH-PBE-1.1-atomicdata cd work for file in cd ../UTILS/InputFiles && ls -1 ; do sed -e'/XMLOUT/i' -e'5' ../UTILS/InputFiles/$file | ../../atompaw-4.1.0.6/src/atompaw > sed -e's/_input/.out/' <<<"$file" ; mv *.atomicdata ../JTH-PBE-1.1-atomicdata/sed -e's/_input//;s/.*/&.atomicdata/' <<<"$file" ; done #note the sed command on the input file adds output format "5" (.atomicdata for pwpaw program)

freyso avatar Nov 10 '23 16:11 freyso

So I suspect this is the current abinit database, and the Si potential shouldn't be used.

freyso avatar Nov 10 '23 16:11 freyso

Ok then we got one problem solved. The other question is now whether there's a possibility to offer some pseudo potential that can be reliably used. Do you perhaps see a solution?

samwaseda avatar Nov 10 '23 16:11 samwaseda

And here is the end of output of the PAW generator when it generates the Si potential:

Eigenvalues of overlap operator (in the basis of projectors):
    1       -4.50102841E-04
    2        7.26235054E-02
    3        2.51794387E+01
    4        3.09917054E+01

So as suspected, the overlap operator is pathological. (must be positive)

freyso avatar Nov 10 '23 16:11 freyso

If you need a Si potential, what about the GPAW ones (I think, xml format).

freyso avatar Nov 10 '23 16:11 freyso

If you need a Si potential, what about the GPAW ones (I think, xml format).

And this can be used with SPHInX right? Should we test it beforehand? If yes, could you maybe tell us what to take a look at?

samwaseda avatar Nov 10 '23 16:11 samwaseda

Probably start with the basic stuff: Murnaghan plot for bulk Si? And then maybe energy-volume curve for your material of interest and compare with literature or VASP pot?

Silicon shouldn't be very problematic, that's why I now remember that I was surprised that JTH produced such a lousy one (or maybe exactly because: if you run with low enough cutoff, you may be able to deform the projector enough to rise the overlap eigenvalue above zero (in the end, it is almost zero), and then you never notice that it will fail pw cutoff convergence tests).

freyso avatar Nov 10 '23 16:11 freyso

Btw: Jörg is from the "we always have to test the pseudopotentials we use" age (=early bronze age of DFT), so he may have additional hints what to test.

freyso avatar Nov 10 '23 16:11 freyso

I found this page. I copy and paste the first 20 lines of Si.PBE:

<?xml version="1.0"?>
<paw_setup version="0.6">
  <!-- Silicon setup for the Projector Augmented Wave method. -->
  <!-- Units: Hartree and Bohr radii.                         -->
  <atom symbol="Si" Z="14" core="10.0" valence="4"/>
  <xc_functional type="GGA" name="PBE"/>
  <generator type="scalar-relativistic" name="gpaw-0.4.2039">
    Frozen core: [Ne]
  </generator>
  <ae_energy kinetic="290.505129" xc="-20.648317"
             electrostatic="-559.674620" total="-289.817808"/>
  <core_energy kinetic="285.331842"/>
  <valence_states>
    <state n="3" l="0" f="2"  rc="2.000" e="-0.39733" id="Si-3s"/>
    <state n="3" l="1" f="2"  rc="2.000" e="-0.14997" id="Si-3p"/>
    <state       l="0"        rc="2.000" e=" 0.60267" id="Si-s1"/>
    <state       l="1"        rc="2.000" e=" 0.85003" id="Si-p1"/>
    <state       l="2"        rc="2.000" e=" 0.00000" id="Si-d1"/>
  </valence_states>
  <radial_grid eq="r=a*i/(n-i)" a="0.400000" n="450" istart="0" iend="449" id="g1"/>

This is the kind of xml that can be read?

samwaseda avatar Nov 10 '23 16:11 samwaseda

I hope so. No garantuee, since I had tested a few examples, and there is a lot of flexibility in the format (e.g. grid types), that needs to be accommodated.

freyso avatar Nov 10 '23 17:11 freyso