From 8e79e90fb5245c643f65eacaa4267253c0110dfd Mon Sep 17 00:00:00 2001 From: Sylvain Tricot Date: Tue, 7 Jul 2020 17:27:26 +0200 Subject: [PATCH] Add multielement support for SPRKKR potential --- src/msspec/__init__.py | 1 + src/msspec/calcio.py | 33 +++-- src/msspec/calculator.py | 21 +-- src/msspec/parameters.py | 7 +- src/msspec/phagen/fortran/phagen_scf_2.1_dp.f | 32 ++++- src/msspec/utils.py | 121 +++++++++++++----- tests/sprkkr/copper/Cu.py | 1 + 7 files changed, 161 insertions(+), 55 deletions(-) diff --git a/src/msspec/__init__.py b/src/msspec/__init__.py index 9a852c4..acffe9f 100644 --- a/src/msspec/__init__.py +++ b/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() diff --git a/src/msspec/calcio.py b/src/msspec/calcio.py index 324712a..965dfeb 100644 --- a/src/msspec/calcio.py +++ b/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()) - for symbol in symbols: + numbers = np.unique(atoms.numbers) + for Z in numbers: try: + symbol = chemical_symbols[Z] charge = ionicity[symbol] except KeyError: - charge = 0. - ionicity_section += ionicity_format.format( - ase.data.atomic_numbers[symbol], - charge) + charge = 0 + ionicity_section += ionicity_format.format(Z, 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):] - cluster = Atoms() + atoms = self.parameters.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) - self.prototypical_atoms = cluster - return cluster + proto += Atom(sym, position=xyz) + # now translate the cluster so that the protocluster is at the + # 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 diff --git a/src/msspec/calculator.py b/src/msspec/calculator.py index da4a2a0..2e1ffa2 100644 --- a/src/msspec/calculator.py +++ b/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') - self.phagenio.load_cluster_file('cluster.clu') - self.phagenio.load_prototypical_atoms('div/molinpot3.out') - # Change back to the init_folder - os.chdir(self.init_folder) + + # Load data + self.phagenio.load_tl_file('tmatrix.tl') + self.phagenio.load_cluster_file('cluster.clu') + self.phagenio.load_prototypical_atoms('div/molinpot3.out') + # 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) diff --git a/src/msspec/parameters.py b/src/msspec/parameters.py index 27db224..48eea44 100644 --- a/src/msspec/parameters.py +++ b/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, diff --git a/src/msspec/phagen/fortran/phagen_scf_2.1_dp.f b/src/msspec/phagen/fortran/phagen_scf_2.1_dp.f index 355768b..cfc0e01 100644 --- a/src/msspec/phagen/fortran/phagen_scf_2.1_dp.f +++ b/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 diff --git a/src/msspec/utils.py b/src/msspec/utils.py index f2ee65d..1dee866 100644 --- a/src/msspec/utils.py +++ b/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' : , # 'types': [ # { - # 'Z' : , - # 'RWS' : , - # 'data': + # 'atom_type': , + # 'Z' : , + # 'RWS' : , + # 'data' : # }, # ... # ... @@ -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\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.*?)\s+' + r'VMTZ_Coulomb\s*=\s*(?P.*?)\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*(?PIQ.*)\n' r'(?P(.*\n)+?))\*+'), [int] * 5 + [float]) - """ - for f in exported_files: - # get the IT from the filename - # m=re.match('SPRKKR-IT_(?P\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 = fd.readline().strip('\n') - second_line = fd.readline().strip('\n') + LOGGER.debug("SITES:") + for _ in self.sites_data: + LOGGER.debug(_) - # load Coulomb and Total interstitial potential - pattern = (r'#\s*VMTZ_TOTAL\s*=\s*(?P.*?)\s+' - r'VMTZ_Coulomb\s*=\s*(?P.*?)\s+.*$') - m = re.match(pattern, first_line) - self.phagen_data.update(VintCoulomb=float(m.group('COULOMB')), - VintTotal=float(m.group('TOTAL')), - RHOint=0.) + LOGGER.debug("TYPES:") + for _ in self.types_data: + LOGGER.debug(_) - # load Z, Wigner-Seitz radius from 2nd line of header - type_data = {} - _ = re.split(r'\s+', second_line.strip("# ")) - type_data.update(Z=int(_[0]), RWS=float(_[3])) + LOGGER.debug("OCCUPATION:") + for _ in self.occ_data: + LOGGER.debug(_) - # load the data - data = np.loadtxt(f, comments='#') - type_data.update(data=data) + for site in self.sites_data: + IQ = site['IQ'] + # 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 diff --git a/tests/sprkkr/copper/Cu.py b/tests/sprkkr/copper/Cu.py index b1f1996..5bd7ea2 100644 --- a/tests/sprkkr/copper/Cu.py +++ b/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()