Add multielement support for SPRKKR potential

src/msspec/__init__.py

@@ -35,6 +35,7 @@ def init_msspec():
     ase.atom.names['forward_angle'] = ('forward_angles', 20.)
     ase.atom.names['backward_angle'] = ('backward_angles', 20.)
     ase.atom.names['RA_cut_off'] = ('RA_cuts_off', 1)
+    ase.atom.names['atom_type'] = ('atom_types', None)
     ase.atoms.Atoms.absorber = None
 
 init_msspec()

src/msspec/calcio.py

@@ -15,6 +15,7 @@ from ase import Atom, Atoms
 from ase.data import chemical_symbols
 
 from msspec.misc import UREG, LOGGER
+from msspec.utils import get_atom_index
 
 
 class PhagenIO(object):
@@ -107,15 +108,14 @@ class PhagenIO(object):
         ionicity = self.parameters.ionicity
         ionicity_section = ''
         ionicity_format = '{:8d}{:9.2f}\n'
-        symbols = set(atoms.get_chemical_symbols())
+        numbers = np.unique(atoms.numbers)
-        for symbol in symbols:
+        for Z in numbers:
             try:
+                symbol = chemical_symbols[Z]
                 charge = ionicity[symbol]
             except KeyError:
-                charge = 0.
+                charge = 0
-            ionicity_section += ionicity_format.format(
+            ionicity_section += ionicity_format.format(Z, charge)
-                ase.data.atomic_numbers[symbol],
-                charge)
         ionicity_section += '{:8d}\n'.format(-1)
 
         content = header + cluster_section + ionicity_section
@@ -279,8 +279,9 @@ class PhagenIO(object):
         with open(filename, 'r') as fd:
             content = fd.readlines()
 
-        content = content[-len(self.parameters.atoms):]
+        atoms = self.parameters.atoms
-        cluster = Atoms()
+        content = content[-len(atoms):]
+        proto = Atoms()
         for line in content:
             datatxt = np.array(re.split(r'\s+', line.strip(' \n')))
             neq = int(datatxt[-1])
@@ -289,9 +290,19 @@ class PhagenIO(object):
                 xyz = (datatxt[[2, 3, 4]].astype(float) * UREG.a0).to(
                        'angstrom').magnitude
                 sym = chemical_symbols[Z]
-                cluster += Atom(sym, position=xyz)
+                proto += Atom(sym, position=xyz)
-        self.prototypical_atoms = cluster
+        # now translate the cluster so that the protocluster is at the
-        return cluster
+        # same origin than the original cluster
+        origin0 = atoms[atoms.absorber].position
+        origin1 = proto[0].position
+        proto.translate(origin0 - origin1)
+        # for each atom in this prototypical cluster, find the corresponding
+        # atom index in the original cluster
+        indices = []
+        for atom in proto:
+            indices.append(get_atom_index(atoms, *atom.position))
+        self.prototypical_atoms = atoms[indices]
+        return self.prototypical_atoms
 
     def load_potential_file(self, filename='plot_vc.dat'):
         a_index = 0

src/msspec/calculator.py

@@ -314,20 +314,22 @@ class _MSCALCULATOR(Calculator):
             shutil.copy(inpot_fname, dest)
 
         # Run if input files changed
+        os.chdir(os.path.join(self.tmp_folder, 'output'))
+        modified = True
         if modified:
             # Run phagen
-            os.chdir(os.path.join(self.tmp_folder, 'output'))
             do_phagen()
             # Copy some output files to be more explicit
             shutil.copy('clus/clus.out', 'cluster.clu')
             shutil.copy('fort.35', 'tmatrix.tl')
             shutil.copy('fort.55', 'tmatrix.rad')
-            # Load data
+
-            self.phagenio.load_tl_file('tmatrix.tl')
+        # Load data
-            self.phagenio.load_cluster_file('cluster.clu')
+        self.phagenio.load_tl_file('tmatrix.tl')
-            self.phagenio.load_prototypical_atoms('div/molinpot3.out')
+        self.phagenio.load_cluster_file('cluster.clu')
-            # Change back to the init_folder
+        self.phagenio.load_prototypical_atoms('div/molinpot3.out')
-            os.chdir(self.init_folder)
+        # Change back to the init_folder
+        os.chdir(self.init_folder)
 
 
     def run_spec(self, malloc={}):
@@ -353,6 +355,7 @@ class _MSCALCULATOR(Calculator):
         #self._make('compute')
         #t1 = time.time()
         #self.resources['spec_time'] = t1 - t0
+        self.modified = True
         if self.modified:
             #self.get_tmatrix()
             t0 = time.time()
@@ -479,7 +482,9 @@ class _MSCALCULATOR(Calculator):
 
         max_tl = self.tmatrix_parameters.get_parameter('max_tl').value
         cluster = self.phagen_parameters.get_parameter('atoms').value
-        proto_indices = cluster.get_array('proto_indices')
+        #proto_indices = cluster.get_array('proto_indices')
+        proto_indices = np.array(
+            [atom.get('proto_index') for atom in self.atoms])
 
         if max_tl != None:
             LOGGER.debug("    applying max_tl: %s", max_tl)

src/msspec/parameters.py

@@ -201,6 +201,9 @@ class Parameter(object):
         value = args[-1]
         self.value = value
 
+    def reset(self):
+        self.value = self.default
+
     def __str__(self):
         fmt = '{:' + self.fmt + '}'
         try:
@@ -358,6 +361,8 @@ class PhagenParameters(BaseParameters):
             Parameter('l2h', types=(int,), fmt='d', default=4),
 
             Parameter('ionicity', types=dict, default={}),
+            Parameter('noproto', allowed_values=('.true.', '.false.'),
+                      types=(str,), fmt='>7s', default='.true.'),
             #Parameter('absorber', types=(int,), limits=(1, None), fmt='d', default=1),
             #Parameter('nosym', types=(str,), allowed_values=('.true.', '.false.'), fmt='s', default='.true.'),
             #Parameter('outersph', types=(str,), allowed_values=('.true.', '.false.'), fmt='s', default='.false.'),
@@ -625,7 +630,7 @@ class SpecParameters(BaseParameters):
                       fmt='d'),
             Parameter('calc_npath', types=int, limits=[0, None], default=100,
                       fmt='d'),
-            Parameter('calc_vint', types=float, default=0., fmt='.2f'),
+            Parameter('calc_vint', types=(float,int), default=0., fmt='.2f'),
             Parameter('calc_ifwd', types=int, limits=[0, 1], default=0,
                       fmt='d'),
             Parameter('calc_nthout', types=int, limits=[0, None], default=1,

src/msspec/phagen/fortran/phagen_scf_2.1_dp.f

@@ -335,6 +335,9 @@ CST   Phagen to python shared object modification <==
 C
       complex*16 eelsme,p1,p2,p3,ramfsr1,ramfsr2,ramfsr3
       complex*16 p3irreg,p2irreg
+CST
+      logical noproto
+CST
 C
 C.....................................................................
 C
@@ -347,6 +350,9 @@ C
      5                   lmxels(3,ua_),p3irreg(rdx_,7),p2irreg(rdx_,7)
       common /typot/     ipot
       common /v_type/    version
+CST
+      common /misc/ noproto
+CST
 C
 C.....................................................................
 C
@@ -354,7 +360,11 @@ C
      1               emin,emax,delta,gamma,eftri,cip,vc0,rs0,vinput,
      2               eikappr,rsh,db,lmaxt,ovlpfac,ionzst,charelx,
      3               calctype,potgen,lmax_mode,relc,einc,esct,scangl,
-     4               optrsh,enunit,lambda,expmode,nosym,version
+     4               optrsh,enunit,lambda,expmode,nosym,version,
+CST   ==> Phagen to python shared object modification
+CST   I added an option to control the number of prototypical atoms
+     5               noproto
+CST   Phagen to python shared object modification <==
 C
 C.....................................................................
 C
@@ -385,6 +395,11 @@ C
       enunit='Ryd'          ! energy unit
 C
       version='2.0'         ! MsSpec version number
+
+CST   ==> Phagen to python shared object modification
+      noproto = .false.
+CST   Phagen to python shared object modification
+
 C
       vc0 = -0.7d0
       rs0 = 3.d0
@@ -2391,11 +2406,24 @@ c
       common/charge_center/cc_dif(3,1),z_shift,i_z_shift,shift_cc
       common/transform/trans
       logical shift_cc
+CST
+      common /misc/ noproto
+      logical noproto
+CST
 c
-      data zero,thr/0.0d0,0.001d0/ !if thr is negative, all cluster atoms are considered prototypical 
+CST   ==> Phagen to python shared object modification
+CST   data zero,thr/0.0d0,0.001d0/ !if thr is negative, all cluster atoms are considered prototypical 
+      data zero/0.0d0/
+CST   Phagen to python shared object modification <==
 c
       data jtape/21/
       data lunout/7/
+
+CST   ==> Phagen to python shared object modification
+CST
+      thr=0.001d0
+      if (noproto) thr=-1*thr
+CST   Phagen to python shared object modification
 c
 c-----------------------------------------------------------------------
 c  find the center of charge of the nuclear framework and

src/msspec/utils.py

@@ -30,9 +30,11 @@ Module utils
 """
 
 
+import os
 import re
 
 import numpy as np
+import ase.atom
 from ase import Atom
 from ase import Atoms
 
@@ -51,9 +53,10 @@ class ForeignPotential(object):
         #     'RHOint'     : <float>,
         #     'types': [
         #        {
-        #            'Z'   : <int>,
+        #            'atom_type': <int>,
-        #            'RWS' : <float>,
+        #            'Z'        : <int>,
-        #            'data': <np.array(..., 4, dtype=float)>
+        #            'RWS'      : <float>,
+        #            'data'     : <np.array(..., 4, dtype=float)>
         #        },
         #        ...
         #        ...
@@ -65,7 +68,7 @@ class ForeignPotential(object):
     def write(self, filename, prototypical_atoms):
         LOGGER.debug(f"Writing Phagen input potential file: {filename}")
 
-        def append_atom_potential(atom):
+        def DEPRECATEDappend_atom_potential(atom):
             Z = atom.number
             # Find the right type (Z) in the phagen_data
             itypes = []
@@ -85,6 +88,18 @@ class ForeignPotential(object):
                 s += line_fmt.format(*row)
             return s
 
+        def append_atom_potential(atom):
+            line_fmt = "{:+1.8e} " * 4 + "\n"
+            atom_type = atom.get('atom_type')
+            assert atom_type != None, f"Unable get the atom type!"
+            for t in self.phagen_data['types']:
+                if t['atom_type'] == atom_type:
+                    s = "{:<7d}{:<10d}{:1.4f}\n".format(
+                        t['Z'], len(t['data']), t['RWS'])
+                    for row in t['data']:
+                        s += line_fmt.format(*row)
+                    return s
+
         content = ""
         # Start by writing the header line
         content += "{:.2f}  {:.4f}   {:.2f}\n".format(
@@ -115,7 +130,43 @@ class ForeignPotential(object):
 class SPRKKRPotential(ForeignPotential):
     def __init__(self, atoms, potfile, *exported_files):
         super().__init__()
-        """
+        self.atoms = atoms
+        self.potfile = potfile
+        self.load_sprkkr_atom_types()
+        for f in exported_files:
+            LOGGER.info(f"Loading file {f}...")
+            # get the IT from the filename
+            m=re.match('SPRKKR-IT_(?P<IT>\d+)-PHAGEN.*', os.path.basename(f))
+            it = int(m.group('IT'))
+
+            # load the content of the header (2 lines)
+            with open(f, 'r') as fd:
+                first_line, second_line = [fd.readline().strip()
+                                           for _ in range(2)]
+
+            # load Coulomb and Total interstitial potential
+            pattern = (r'#\s*VMTZ_TOTAL\s*=\s*(?P<TOTAL>.*?)\s+'
+                       r'VMTZ_Coulomb\s*=\s*(?P<COULOMB>.*?)\s+.*$')
+            m = re.match(pattern, first_line)
+            self.phagen_data.update(VintCoulomb=float(m.group('COULOMB')),
+                                    VintTotal=float(m.group('TOTAL')),
+                                    RHOint=0.)
+
+            # load Z, Wigner-Seitz radius from 2nd line of header
+            type_data = {}
+            _ = re.split(r'\s+', second_line.strip("# "))
+            type_data.update(atom_type=int(it), Z=int(_[0]), RWS=float(_[3]))
+
+            # load the data
+            data = np.loadtxt(f, comments='#')
+            type_data.update(data=data)
+
+            self.phagen_data['types'].append(type_data)
+
+        # store the sprkkr type number in the sprkkr_info property of each
+        # atom in the provided self.atoms
+
+    def load_sprkkr_atom_types(self):
         def read(content, pattern, types):
             # compile the pattern for regex matching
             pat = re.compile(pattern, re.MULTILINE)
@@ -137,7 +188,8 @@ class SPRKKRPotential(ForeignPotential):
             return data
 
         # load info in *.pot file
-        with open(potfile, 'r') as fd:
+        LOGGER.info(f"Loading SPRKKR *.pot file {self.potfile}...")
+        with open(self.potfile, 'r') as fd:
             content = fd.read()
 
         self.sites_data = read(content,
@@ -152,35 +204,33 @@ class SPRKKRPotential(ForeignPotential):
                              (r'^\s*OCCUPATION\s*\n(\s*(?P<KEYS>IQ.*)\n'
                               r'(?P<DATA>(.*\n)+?))\*+'),
                              [int] * 5 + [float])
-        """
-        for f in exported_files:
-            # get the IT from the filename
-            # m=re.match('SPRKKR-IT_(?P<IT>\d+)-PHAGEN.*', os.path.basename(f))
-            # it = int(m.group('IT'))
 
-            # load the content of the header (2 lines)
+        LOGGER.debug("SITES:")
-            with open(f, 'r') as fd:
+        for _ in self.sites_data:
-                first_line = fd.readline().strip('\n')
+            LOGGER.debug(_)
-                second_line = fd.readline().strip('\n')
 
-            # load Coulomb and Total interstitial potential
+        LOGGER.debug("TYPES:")
-            pattern = (r'#\s*VMTZ_TOTAL\s*=\s*(?P<TOTAL>.*?)\s+'
+        for _ in self.types_data:
-                       r'VMTZ_Coulomb\s*=\s*(?P<COULOMB>.*?)\s+.*$')
+            LOGGER.debug(_)
-            m = re.match(pattern, first_line)
-            self.phagen_data.update(VintCoulomb=float(m.group('COULOMB')),
-                                    VintTotal=float(m.group('TOTAL')),
-                                    RHOint=0.)
 
-            # load Z, Wigner-Seitz radius from 2nd line of header
+        LOGGER.debug("OCCUPATION:")
-            type_data = {}
+        for _ in self.occ_data:
-            _ = re.split(r'\s+', second_line.strip("# "))
+            LOGGER.debug(_)
-            type_data.update(Z=int(_[0]), RWS=float(_[3]))
 
-            # load the data
+        for site in self.sites_data:
-            data = np.loadtxt(f, comments='#')
+            IQ = site['IQ']
-            type_data.update(data=data)
+            # Get the Atom at x, y, z
+            x = site['QBAS(X)']
+            y = site['QBAS(Y)']
+            z = site['QBAS(Z)']
+            i = get_atom_index(self.atoms, x, y, z, scaled=True)
+            for occupation in self.occ_data:
+                if occupation['IQ'] == IQ:
+                    IT = occupation['ITOQ']
+                    atom = self.atoms[i]
+                    atom.set('atom_type', IT)
+            LOGGER.debug(f"Site #{IQ} is type #{IT}, atom {atom}")
 
-            self.phagen_data['types'].append(type_data)
 
 
 class EmptySphere(Atom):
@@ -189,7 +239,7 @@ class EmptySphere(Atom):
         self.symbol = 'X'
 
 
-def get_atom_index(atoms, x, y, z):
+def get_atom_index(atoms, x, y, z, scaled=False):
     """ Return the index of the atom that is the closest to the coordiantes
     given as parameters.
 
@@ -197,15 +247,20 @@ def get_atom_index(atoms, x, y, z):
     :param float x: the x position in angstroms
     :param float y: the y position in angstroms
     :param float z: the z position in angstroms
+    :param bool scaled: whether the x,y,z coordinates are scaled
 
     :return: the index of the atom as an integer
     :rtype: int
     """
     # get all distances
-    d = np.linalg.norm(atoms.get_positions() - np.array([x, y, z]), axis=1)
+    if scaled:
+        positions = atoms.get_scaled_positions()
+    else:
+        positions = atoms.get_positions()
+    d = np.linalg.norm(positions - np.array([x, y, z]), axis=1)
     # get the index of the min distance
     i = np.argmin(d)
-    # return the index and the corresponding distance
+    # return the index
     return i
 
 

tests/sprkkr/copper/Cu.py

@@ -67,6 +67,7 @@ if 'msspec' in sys.argv:
 
     pot = SPRKKRPotential(Cu, "Cu/Cu.pot", *glob.glob("Cu/*PHAGEN.pot"))
     calc.tmatrix_parameters.potential = pot
+    calc.phagen_parameters.noproto = '.true.'
 
     data = calc.get_theta_scan(level='2p3/2')
     data.view()