From d1963130a127c2280f3332bd06ab6b14f06fa761 Mon Sep 17 00:00:00 2001 From: Brock Dyer Date: Mon, 13 Jul 2026 18:02:49 -0400 Subject: [PATCH] Nexus: Remove `Matter`, `Particle`, `Particles`, `Ion`, and `PseudoIon`, along with the related functions --- nexus/nexus/physical_system.py | 381 ++-------------------- nexus/nexus/pwscf_input.py | 8 - nexus/nexus/qmcpack_input.py | 4 - nexus/nexus/rmg_input.py | 1 - nexus/nexus/tests/test_physical_system.py | 169 +--------- 5 files changed, 36 insertions(+), 527 deletions(-) diff --git a/nexus/nexus/physical_system.py b/nexus/nexus/physical_system.py index 4bd4271801..916e48859b 100644 --- a/nexus/nexus/physical_system.py +++ b/nexus/nexus/physical_system.py @@ -5,8 +5,8 @@ #====================================================================# # physical_system.py # -# Representations of matter, particles, and particles collected # -# together in complete systems. # +# Representations of particles collected together in complete # +# systems. # # # # Content summary: # # PhysicalSystem # @@ -15,273 +15,28 @@ # generate_physical_system # # User function to create arbitrary physical systems. # # # -# Matter # -# Base class for all forms of matter. # -# Class contains a list of all matter known to Nexus. # -# # -# Particle # -# Class representing a particular particle species. # -# # -# Ion, PseudoIon # -# Specialized Particle classes for ion species and ions # -# represented by pseudopotentials. # -# # -# Particles # -# A collection of particles. # -# PhysicalSystem objects contain a Particles instance. # -# # #====================================================================# import os from pathlib import Path from copy import deepcopy import numpy as np -from .developer import DevBase, obj -from .unit_converter import convert +from .developer import DevBase, obj, warn from .periodic_table import Elements from .structure import Structure, generate_structure, read_structure -class Matter(DevBase): - particle_collection = None - - @classmethod - def set_elements(cls,elements): - cls.elements = set(elements) - #end def set_elements - - @classmethod - def set_particle_collection(cls,pc): - cls.particle_collection = pc - #end def set_particle_collection - - @classmethod - def new_particles(cls,*particles,**named_particles): - cls.particle_collection.add_particles(*particles,**named_particles) - #end def new_particles - - def is_element(self,name,symbol=False): - if symbol: - is_elem, element = Elements.is_element(name, return_element=symbol) - return is_elem, element.symbol - else: - is_elem = Elements.is_element(name) - return is_elem - #end def is_element -#end class Matter - - -class Particle(Matter): - def __init__(self,name=None,mass=None,charge=None,spin=None): - self.name = name - self.mass = mass - self.charge = charge - self.spin = spin - #end def __init__ - - def set_count(self,count): - self.count = count - #end def set_count -#end class Particle - - -class Ion(Particle): - def __init__(self,name=None,mass=None,charge=None,spin=None, - protons=None,neutrons=None): - Particle.__init__(self,name,mass,charge,spin) - self.protons = protons - self.neutrons = neutrons - #end def __init__ - - def pseudize(self,valence): - ps = PseudoIon() - ps.transfer_from(self) - ps.charge = valence - ps.core_electrons = ps.protons - valence - return ps - #end def pseudize -#end class Ion - - -class PseudoIon(Ion): - def __init__(self,name=None,mass=None,charge=None,spin=None, - protons=None,neutrons=None,core_electrons=None): - Ion.__init__(self,name,mass,charge,spin,protons,neutrons) - self.core_electrons = core_electrons - #end def __init__ -#end class PseudoIon - - -class Particles(Matter): - def __init__(self,*particles,**named_particles): - self.add_particles(*particles,**named_particles) - #end def __init__ - - def add_particles(self,*particles,**named_particles): - if len(particles)==1 and isinstance(particles[0],list): - particles = particles[0] - #end if - for particle in particles: - self[particle.name] = particle - #end for - for name,particle in named_particles.items(): - self[name] = particle - #end for - #end def add_particles - - def get_particle(self,name): - p = None - if name in self: - p = self[name] - else: - iselem,symbol = self.is_element(name,symbol=True) - if iselem and symbol in self: - p = self[symbol].copy() - p.name = name - self[name] = p - #end if - #end if - return p - #end def get_particle - - # test needed - def get(self,quantity): - q = obj() - for name,particle in self.items(): - q[name] = particle[quantity] - #end for - return q - #end def get - - def rename(self,**name_pairs): - for old,new in name_pairs.items(): - if old in self: - o = self[old] - o.name = new - del self[old] - if new in self: - self[new].count += o.count - else: - self[new] = o - #end if - #end if - #end for - #end def rename - - def get_ions(self): - ions = obj() - for name,particle in self.items(): - if self.is_element(name): - ions[name] = particle - #end if - #end for - return ions - #end def get_ions - - def count_ions(self,species=False): - nions = 0 - nspecies = 0 - for name,particle in self.items(): - if self.is_element(name): - nspecies += 1 - nions += particle.count - #end if - #end for - if species: - return nions,nspecies - else: - return nions - #end if - #end def count_ions - - def get_electrons(self): - electrons = obj() - for electron in ('up_electron','down_electron'): - if electron in self: - electrons[electron] = self[electron] - #end if - #end for - return electrons - #end def get_electrons - - def count_electrons(self): - nelectrons = 0 - for electron in ('up_electron','down_electron'): - if electron in self: - nelectrons += self[electron].count - #end if - #end for - return nelectrons - #end def count_electrons - - def electron_counts(self): - counts = [] - for electron in ('up_electron','down_electron'): - if electron in self: - counts.append(self[electron].count) - else: - counts.append(0) - #end if - #end for - return counts - #end def electron_counts -#end class Particles - -amu_me = convert(1.,'amu','me') - -plist = [ - Particle('up_electron' ,1.0,-1, 1), - Particle('down_electron',1.0,-1,-1), - ] - -for elem in Elements: - plist.append( - Ion( - name = elem.symbol, - mass = elem.atomic_weight * amu_me, - charge = elem.atomic_number, - spin = 0, # Don't have this data - protons = elem.atomic_number, - neutrons = round(elem.atomic_weight-elem.atomic_number), - ) - ) -#end for -for elem in Elements: - for mass_number, rel_atomic_mass in elem.isotopes.items(): - plist.append( - Ion( - name = f"{elem.symbol}_{mass_number}", - mass = rel_atomic_mass * amu_me, - charge = elem.atomic_number, - spin = 0, # Don't have this data - protons = elem.atomic_number, - neutrons = round(rel_atomic_mass - elem.atomic_number), - ) - ) - #end for -#end for - -Matter.set_elements([e.symbol for e in Elements]) -Matter.set_particle_collection(Particles(plist)) - -del plist - +class PhysicalSystem(DevBase): + ghost_aliases = ["Xx"] -class PhysicalSystem(Matter): - - def __init__(self,structure=None,net_charge=0,net_spin=0,particles=None,**valency): + def __init__(self,structure=None,net_charge=0,net_spin=0,**valency): self.pseudized = False if structure is None: self.structure = Structure() else: self.structure = structure #end if - if particles is None: - self.particles = Particles() - else: - self.particles = particles.copy() - #end if self.folded_system = None if self.structure.has_folded(): @@ -316,113 +71,27 @@ def __init__(self,structure=None,net_charge=0,net_spin=0,particles=None,**valenc structure = structure.folded_structure, net_charge = net_charge_fold, net_spin = net_spin_fold, - particles = particles, **valency ) #end if - - self.valency_in = obj(**valency) - self.net_charge_in = net_charge - self.net_spin_in = net_spin - - self.update_particles(clear=False) + if valency is not None and len(valency) > 0: + self.pseudize(**valency) + else: + self.valency = None + self.net_charge = net_charge + self.net_spin = net_spin self.check_folded_system() #end def __init__ - def update_particles(self,clear=True): - #add ions - pc = dict() - elem = list(self.structure.elem) - for ion in set(elem): - pc[ion] = elem.count(ion) - #end for - missing = set(pc.keys())-set(self.particles.keys()) - if len(missing)>0 or len(elem)==0: - if clear: - self.particles.clear() - #end if - self.add_particles(**pc) - - #pseudize - if len(self.valency_in)>0: - self.pseudize(**self.valency_in) - #end if - - #add electrons - self.generate_electrons(self.net_charge_in,self.net_spin_in) - #end if - #end def update_particles - - - def update(self): - self.net_charge = self.structure.background_charge - self.net_spin = 0 - for p in self.particles: - self.net_charge += p.count*p.charge - self.net_spin += p.count*p.spin - #end for - self.net_charge = int(round(float(self.net_charge))) - self.net_spin = int(round(float(self.net_spin))) - #end def update - - - def add_particles(self,**particle_counts): - pc = self.particle_collection # all known particles - plist = [] - for name,count in particle_counts.items(): - particle = pc.get_particle(name) - if particle is None: - self.error('particle {0} is unknown'.format(name)) - else: - particle = particle.copy() - #end if - particle.set_count(count) - plist.append(particle) - #end for - self.particles.add_particles(plist) - self.update() - #end def add_particles - - - def generate_electrons(self,net_charge=0,net_spin=0): - nelectrons = -net_charge + self.net_charge - if net_spin=='low': - net_spin = nelectrons%2 - #end if - nup = float(nelectrons + net_spin - self.net_spin)/2 - ndown = float(nelectrons - net_spin + self.net_spin)/2 - if abs(nup-int(nup))>1e-3: - self.error('requested spin state {0} incompatible with {1} electrons'.format(net_spin,nelectrons)) - #end if - nup = int(nup) - ndown = int(ndown) - self.add_particles(up_electron=nup,down_electron=ndown) - #end def generate_electrons - - def pseudize(self,**valency): - errors = False - for ion,valence_charge in valency.items(): - if ion in self.particles: - ionp = self.particles[ion] - if isinstance(ionp,Ion): - self.particles[ion] = ionp.pseudize(valence_charge) - self.pseudized = True - else: - self.error(ion+' cannot be pseudized',exit=False) - #end if - else: + for ion in valency.keys(): + if ion not in self.ion_labels: self.error(ion+' is not in the physical system',exit=False) - errors = True - #end if - #end for - if errors: - self.error('system cannot be generated') - #end if + self.valency = obj(**valency) - self.update() + self.pseudized = True #end def pseudize @@ -492,7 +161,6 @@ def group_atoms(self): def rename(self,folded=True,**name_pairs): - self.particles.rename(**name_pairs) self.structure.rename(folded=False,**name_pairs) if self.pseudized: for old,new in name_pairs.items(): @@ -503,7 +171,6 @@ def rename(self,folded=True,**name_pairs): del self.valency[old] #end if #end for - self.valency_in = self.valency #end if if self.folded_system is not None and folded: self.folded_system.rename(folded=folded,**name_pairs) @@ -628,7 +295,7 @@ def ae_pp_species(self): def kf_rpa(self): - nelecs = self.particles.electron_counts() + nelecs = (self.n_up, self.n_down) volume = self.structure.volume() kvol1 = (2*np.pi)**3/volume # k-space volume per particle kfs = [(3*nelec*kvol1/(4*np.pi))**(1./3) for nelec in nelecs] @@ -641,7 +308,7 @@ def n_elec(self): ions = self.structure.elem.tolist() tot_charge = 0 for ion in ions: - if hasattr(self, "valency"): + if self.valency is not None: if ion in self.valency: tot_charge += self.valency[ion] else: @@ -675,7 +342,7 @@ def ion_labels(self): @property def Zeff(self): - if hasattr(self, "valency"): + if self.valency is not None: return self.valency Zeff = {} @@ -758,11 +425,9 @@ def generate_physical_system(**kwargs): del kwargs['tiled_spin'] del kwargs['extensive'] if 'particles' in kwargs: - particles = kwargs['particles'] + warn("The keyword `particles` is no longer valid. Please remove from your scripts!") del kwargs['particles'] - else: - generation_info.particles = None - #end if + pretile = kwargs['pretile'] del kwargs['pretile'] valency = dict() @@ -848,7 +513,5 @@ def generate_physical_system(**kwargs): # test needed def ghost_atoms(*particles): for particle in particles: - Matter.particle_collection.add_particles(Ion(name=particle,mass=0,charge=0,spin=0,protons=0,neutrons=0)) - #end for + PhysicalSystem.ghost_aliases.append(particle) #end def ghost_atoms - diff --git a/nexus/nexus/pwscf_input.py b/nexus/nexus/pwscf_input.py index 86d43076d0..a911dcea12 100644 --- a/nexus/nexus/pwscf_input.py +++ b/nexus/nexus/pwscf_input.py @@ -1673,7 +1673,6 @@ def incorporate_hubbard(self, hubbard_result): def incorporate_system(self,system,elem_order=None): system.check_folded_system() - system.update_particles() system.change_units('B') s = system.structure nc = system.net_charge @@ -1769,7 +1768,6 @@ def incorporate_system(self,system,elem_order=None): def incorporate_system_old(self,system,spin_polarized=None): system.check_folded_system() - system.update_particles() system.change_units('B') s = system.structure nc = system.net_charge @@ -1967,12 +1965,6 @@ def standardize_types(self): def generate_pwscf_input(selector,**kwargs): - if 'system' in kwargs: - system = kwargs['system'] - if isinstance(system,PhysicalSystem): - system.update_particles() - #end if - #end if if selector=='generic': return generate_any_pwscf_input(**kwargs) if selector=='scf': diff --git a/nexus/nexus/qmcpack_input.py b/nexus/nexus/qmcpack_input.py index f456b1a2dd..3f6fb68533 100644 --- a/nexus/nexus/qmcpack_input.py +++ b/nexus/nexus/qmcpack_input.py @@ -8422,10 +8422,6 @@ def generate_batched_dmc_calculations( def generate_qmcpack_input(**kwargs): QIcollections.clear() - system = kwargs.get('system',None) - if isinstance(system,PhysicalSystem): - system.update_particles() - #end if selector = kwargs.pop('input_type','basic') if selector=='basic': inp = generate_basic_input(**kwargs) diff --git a/nexus/nexus/rmg_input.py b/nexus/nexus/rmg_input.py index f6ec3bdb26..0ff37ba488 100644 --- a/nexus/nexus/rmg_input.py +++ b/nexus/nexus/rmg_input.py @@ -3410,7 +3410,6 @@ def generate_any_rmg_input(**kwargs): system = system.get_smallest() #end if system.check_folded_system() - system.update_particles() # set atomic species, positions, magnetic moments and mobility if 'atomic_coordinate_type' not in ri: diff --git a/nexus/nexus/tests/test_physical_system.py b/nexus/nexus/tests/test_physical_system.py index 2a015ffe71..482650780a 100644 --- a/nexus/nexus/tests/test_physical_system.py +++ b/nexus/nexus/tests/test_physical_system.py @@ -6,18 +6,15 @@ generic_settings.raise_error = True import numpy as np -from .. import testing from ..testing import value_eq,object_eq from nexus.physical_system import generate_physical_system -from nexus.periodic_table import Elements -from nexus.unit_converter import convert from .test_structure import structure_same def system_same(s1,s2,pseudized=True,tiled=False): same = True - keys = ('net_charge','net_spin','pseudized','particles') + keys = ('net_charge','net_spin','pseudized') o1 = s1.obj(keys) o2 = s2.obj(keys) qsame = object_eq(o1,o2) @@ -35,75 +32,7 @@ def system_same(s1,s2,pseudized=True,tiled=False): #end def system_same -def test_particle_initialization(): - from ..physical_system import Matter,Particle,Ion,PseudoIon,Particles - - # empty initialization - Matter() - p = Particle() - i = Ion() - pi = PseudoIon() - Particles() - - def check_none(v,vfields): - for f in vfields: - assert(f in v) - assert(v[f] is None) - #end for - #end def check_none - - pfields = 'name mass charge spin'.split() - ifields = pfields + 'protons neutrons'.split() - pifields = ifields+['core_electrons'] - - check_none(p,pfields) - check_none(i,ifields) - check_none(pi,pifields) - - # matter - elements = Matter.elements - assert(len(elements)==119) - assert('Si' in elements) - - pc = Matter.particle_collection - for e in elements: - assert(e in pc) - #end for - assert('up_electron' in pc) - assert('down_electron' in pc) - - u = pc.up_electron - assert(u.name=='up_electron') - assert(value_eq(u.mass,1.0)) - assert(u.charge==-1) - assert(u.spin==1) - - d = pc.down_electron - assert(d.name=='down_electron') - assert(value_eq(d.mass,1.0)) - assert(d.charge==-1) - assert(d.spin==-1) - - si = pc.Si - assert(si.name=='Si') - print(si.mass) - assert(value_eq(si.mass,51195.82309476658)) - assert(si.charge==14) - assert(si.protons==14) - assert(si.neutrons==14) - - # test get_particle - assert(object_eq(pc.get_particle('Si'),si)) - si1 = si.copy() - si1.name = 'Si1' - assert(object_eq(pc.get_particle('Si1'),si1)) - -#end def test_particle_initialization - - - def test_physical_system_initialization(tmp_path): - import os from ..developer import obj from ..structure import generate_structure from ..physical_system import generate_physical_system @@ -242,33 +171,6 @@ def test_physical_system_initialization(tmp_path): C = 4, ) - pref = obj( - C = obj( - charge = 4, - core_electrons = 2, - count = 8, - mass = 21894.7135906, - name = 'C', - neutrons = 6, - protons = 6, - spin = 0, - ), - down_electron = obj( - charge = -1, - count = 16, - mass = 1.0, - name = 'down_electron', - spin = -1, - ), - up_electron = obj( - charge = -1, - count = 16, - mass = 1.0, - name = 'up_electron', - spin = 1, - ), - ) - # check direct system w/o tiling ref = direct_notile sref = ref.structure @@ -276,15 +178,12 @@ def test_physical_system_initialization(tmp_path): assert(ref.net_spin==0) assert(ref.pseudized) assert(object_eq(ref.valency,obj(C=4))) - assert(object_eq(ref.particles.to_obj(),pref)) assert(structure_same(sref,d8)) assert(value_eq(sref.axes,3.57*np.eye(3))) assert(tuple(sref.bconds)==tuple('ppp')) assert(list(sref.elem)==8*['C']) assert(value_eq(tuple(sref.pos[-1]),(2.6775,2.6775,0.8925))) assert(sref.units=='A') - assert(object_eq(ref.particles.get_ions().to_obj(),obj(C=pref.C))) - assert(object_eq(ref.particles.get_electrons().to_obj(),obj(down_electron=pref.down_electron,up_electron=pref.up_electron))) # check direct system w/ tiling ref = direct_tile @@ -293,7 +192,6 @@ def test_physical_system_initialization(tmp_path): assert(ref.net_spin==0) assert(ref.pseudized) assert(object_eq(ref.valency,obj(C=4))) - assert(object_eq(ref.particles.to_obj(),pref)) assert(structure_same(sref,d8_tile)) assert(value_eq(sref.axes,3.57*np.eye(3))) assert(tuple(sref.bconds)==tuple('ppp')) @@ -302,14 +200,10 @@ def test_physical_system_initialization(tmp_path): assert(sref.units=='A') ref = direct_tile.folded_system sref = ref.structure - pref.C.count = 2 - pref.down_electron.count = 4 - pref.up_electron.count = 4 assert(ref.net_charge==0) assert(ref.net_spin==0) assert(ref.pseudized) assert(object_eq(ref.valency,obj(C=4))) - assert(object_eq(ref.particles.to_obj(),pref)) assert(structure_same(sref,d2)) assert(value_eq(sref.axes,1.785*np.array([[1.,1,0],[0,1,1],[1,0,1]]))) assert(tuple(sref.bconds)==tuple('ppp')) @@ -377,16 +271,14 @@ def test_physical_system_initialization(tmp_path): #end for # test particle counts - p = direct_notile.particles - assert(p.count_ions()==8) - assert(p.count_ions(species=True)==(8,1)) - assert(p.count_electrons()==32) - assert(p.electron_counts()==[16,16]) - + assert(direct_notile.n_ions == 8) + assert(direct_notile.n_species == 1) + assert(direct_notile.n_elec == 32) + assert(direct_notile.n_up == 16) + assert(direct_notile.n_down == 16) #end def test_physical_system_initialization - def test_change_units(): from ..physical_system import generate_physical_system @@ -438,26 +330,26 @@ def test_rename(): ref = sys assert(object_eq(ref.valency,obj(C1=4,C2=4))) assert(list(ref.structure.elem)==4*['C1','C2']) - assert(ref.particles.count_ions()==8) - assert(ref.particles.count_ions(species=True)==(8,2)) + assert(ref.n_ions==8) + assert(ref.n_species==2) ref = sys.folded_system assert(object_eq(ref.valency,obj(C1=4,C2=4))) assert(list(ref.structure.elem)==['C1','C2']) - assert(ref.particles.count_ions()==2) - assert(ref.particles.count_ions(species=True)==(2,2)) + assert(ref.n_ions==2) + assert(ref.n_species==2) sys.rename(C1='C',C2='C') ref = sys assert(object_eq(ref.valency,obj(C=4))) assert(list(ref.structure.elem)==8*['C']) - assert(ref.particles.count_ions()==8) - assert(ref.particles.count_ions(species=True)==(8,1)) + assert(ref.n_ions==8) + assert(ref.n_species==1) ref = sys.folded_system assert(object_eq(ref.valency,obj(C=4))) assert(list(ref.structure.elem)==2*['C']) - assert(ref.particles.count_ions()==2) - assert(ref.particles.count_ions(species=True)==(2,1)) + assert(ref.n_ions==2) + assert(ref.n_species==1) #end def test_rename @@ -502,7 +394,6 @@ def test_tile(): def test_kf_rpa(): from .test_structure import example_structure_h4 - from ..physical_system import generate_physical_system s1 = example_structure_h4() ps = generate_physical_system( structure = s1, @@ -514,35 +405,3 @@ def test_kf_rpa(): assert np.isclose(kfs[0], 1.465, atol=1e-3) assert np.isclose(kfs[1], 1.465/2**(1./3), atol=1e-3) #end def test_kf_rpa - - -def test_particle_equiv(): - - ref = generate_physical_system( - units = 'A', - axes = [[1.785, 1.785, 0. ], - [0. , 1.785, 1.785], - [1.785, 0. , 1.785]], - elem = 2*['C'], - posu = [[0.00, 0.00, 0.00], - [0.25, 0.25, 0.25]], - tiling = [[ 1, -1, 1], - [ 1, 1, -1], - [-1, 1, 1]], - C = 4, - ) - - assert(ref.particles.electron_counts() == [ref.n_up, ref.n_down]) - assert(ref.particles.count_electrons() == ref.n_elec) - assert(ref.particles.count_ions() == ref.n_ions) - assert(ref.particles.count_ions(species=True) == (ref.n_ions, ref.n_species)) - assert(set(ref.particles.get_ions().keys()) == ref.ion_labels) - - ions = ref.particles.get_ions() - for ion in ref.ion_labels: - is_elem, element = Elements.is_element(ion, return_element=True) - assert(ions[ion].name == ion) - assert(ions[ion].charge == ref.Zeff[ion]) - assert(value_eq(ions[ion].mass, convert(element.atomic_weight, "amu", "me"))) - -#end def test_particle_equiv