From 27b5745d0de0c7383095d21e3613056255633cfc Mon Sep 17 00:00:00 2001 From: brockdyer03 Date: Wed, 8 Jul 2026 10:16:50 -0400 Subject: [PATCH 1/4] Nexus: First step to ridding us of `Particles` --- nexus/nexus/gamess.py | 2 +- nexus/nexus/physical_system.py | 22 ++++++++ nexus/nexus/pwscf_input.py | 24 ++++----- nexus/nexus/qmcpack.py | 16 +++--- nexus/nexus/qmcpack_input.py | 76 +++++++++++++--------------- nexus/nexus/quantum_package_input.py | 2 +- nexus/nexus/rmg_input.py | 2 +- nexus/nexus/vasp_input.py | 2 +- 8 files changed, 79 insertions(+), 67 deletions(-) diff --git a/nexus/nexus/gamess.py b/nexus/nexus/gamess.py index 1fcc81f8f6..e7940a468a 100644 --- a/nexus/nexus/gamess.py +++ b/nexus/nexus/gamess.py @@ -154,7 +154,7 @@ def incorporate_result(self,result_name,result,sim): self.block() #end if guess_inputs = obj() - ecounts = self.system.particles.electron_counts() + ecounts = self.system.electron_counts() orbs = result.orbitals order_map = obj(up='iorder',down='jorder') nelec_map = obj(up=ecounts[0],down=ecounts[1]) diff --git a/nexus/nexus/physical_system.py b/nexus/nexus/physical_system.py index 4d4af3ea59..856c70aa58 100644 --- a/nexus/nexus/physical_system.py +++ b/nexus/nexus/physical_system.py @@ -634,6 +634,28 @@ def kf_rpa(self): kfs = [(3*nelec*kvol1/(4*np.pi))**(1./3) for nelec in nelecs] return np.array(kfs) #end def kf_rpa + + + def electron_counts(self): + return self.particles.electron_counts() + + def count_electrons(self): + return self.particles.count_electrons() + + def count_ions(self, species=False): + return self.particles.count_ions(species) + + def get_ions(self): + return self.particles.get_ions() + + @property + def up_electron(self): + return self.particles.up_electron + + @property + def down_electron(self): + return self.particles.down_electron + #end class PhysicalSystem diff --git a/nexus/nexus/pwscf_input.py b/nexus/nexus/pwscf_input.py index 544d29f7dc..c3bde26283 100644 --- a/nexus/nexus/pwscf_input.py +++ b/nexus/nexus/pwscf_input.py @@ -1675,17 +1675,16 @@ def incorporate_system(self,system,elem_order=None): system.check_folded_system() system.update_particles() system.change_units('B') - p = system.particles s = system.structure nc = system.net_charge ns = system.net_spin - nup = p.up_electron.count - ndn = p.down_electron.count + nup = system.up_electron.count + ndn = system.down_electron.count self.system.ibrav = 0 # self.system['celldm(1)'] = 1.0e0 - nions,nspecies = p.count_ions(species=True) + nions,nspecies = system.count_ions(species=True) self.system.nat = nions self.system.ntyp = nspecies self.system.tot_charge = nc @@ -1717,7 +1716,7 @@ def incorporate_system(self,system,elem_order=None): #end if #end if - atoms = p.get_ions() + atoms = system.get_ions() if 'masses' not in self.atomic_species: self.atomic_species.masses = obj() #end if @@ -1739,7 +1738,7 @@ def incorporate_system(self,system,elem_order=None): pp = self.atomic_species.pseudopotentials for atom in self.atomic_species.atoms: if atom not in pp: - iselem,symbol = p.is_element(atom,symbol=True) + iselem,symbol = system.is_element(atom,symbol=True) if iselem and symbol in pp: pp[atom] = str(pp[symbol]) #end if @@ -1770,17 +1769,16 @@ def incorporate_system_old(self,system,spin_polarized=None): system.check_folded_system() system.update_particles() system.change_units('B') - p = system.particles s = system.structure nc = system.net_charge ns = system.net_spin - nup = p.up_electron.count - ndn = p.down_electron.count + nup = system.up_electron.count + ndn = system.down_electron.count self.system.ibrav = 0 # self.system['celldm(1)'] = 1.0e0 - nions,nspecies = p.count_ions(species=True) + nions,nspecies = system.count_ions(species=True) self.system.nat = nions self.system.ntyp = nspecies #self.system.nelec = nup+ndn @@ -1819,7 +1817,7 @@ def incorporate_system_old(self,system,spin_polarized=None): #end if #end if - atoms = p.get_ions() + atoms = system.get_ions() masses = obj() for name,a in atoms.items(): masses[name] = convert(a.mass,'me','amu') @@ -1830,7 +1828,7 @@ def incorporate_system_old(self,system,spin_polarized=None): pp = self.atomic_species.pseudopotentials for atom in self.atomic_species.atoms: if atom not in pp: - iselem,symbol = p.is_element(atom,symbol=True) + iselem,symbol = system.is_element(atom,symbol=True) if iselem and symbol in pp: pp[atom] = str(pp[symbol]) #end if @@ -2265,7 +2263,7 @@ def generate_any_pwscf_input(**kwargs): # set nbnd using bandfac, if provided if nbnd is None and bandfac is not None: - nocc = max(system.particles.electron_counts()) + nocc = max(system.electron_counts()) pw.system.nbnd = int(np.ceil(nocc*bandfac)) #end if diff --git a/nexus/nexus/qmcpack.py b/nexus/nexus/qmcpack.py index 5e7348987c..0c7f49bef0 100644 --- a/nexus/nexus/qmcpack.py +++ b/nexus/nexus/qmcpack.py @@ -188,11 +188,11 @@ def unfolded_nelecs(self): Returns the number of electrons in the primitive cell ''' if self.system.folded_system is None: - n_up = self.system.particles.up_electron.count - n_dn = self.system.particles.down_electron.count + n_up = self.system.up_electron.count + n_dn = self.system.down_electron.count else: - n_up = self.system.folded_system.particles.up_electron.count - n_dn = self.system.folded_system.particles.down_electron.count + n_up = self.system.folded_system.up_electron.count + n_dn = self.system.folded_system.down_electron.count #end if nelecs = n_up + n_dn return nelecs @@ -463,8 +463,8 @@ def sum_charge_twists(self): ''' Returns the net charge of a system with multiple twists (not averaged) ''' - n_up = self.system.particles.up_electron.count - n_dn = self.system.particles.down_electron.count + n_up = self.system.up_electron.count + n_dn = self.system.down_electron.count n_total = n_up + n_dn nelecs_at_twist = self.nelecs_at_twist kweights = np.array(self.system.structure.kweights) @@ -536,8 +536,8 @@ def write_gcta_report(self, locdir, fermi_level, scf_magnet = None): spin_sum_twists = self.sum_spin_twists() qmc_magnet = spin_sum_twists / nosym_kpoints #end if - n_up = self.system.particles.up_electron.count - n_dn = self.system.particles.down_electron.count + n_up = self.system.up_electron.count + n_dn = self.system.down_electron.count n_total = n_up + n_dn nelecs_at_twist = self.nelecs_at_twist fermi_level = np.array(fermi_level) diff --git a/nexus/nexus/qmcpack_input.py b/nexus/nexus/qmcpack_input.py index 258f92b2af..9f91a223c9 100644 --- a/nexus/nexus/qmcpack_input.py +++ b/nexus/nexus/qmcpack_input.py @@ -3825,7 +3825,6 @@ def incorporate_system(self,system): system.change_units('B') #system.structure.group_atoms() system.structure.order_by_species() - particles = system.particles structure = system.structure net_charge = system.net_charge net_spin = system.net_spin @@ -3898,10 +3897,9 @@ def incorporate_system(self,system): sc.lattice = axes #end if - elns = particles.get_electrons() - ions = particles.get_ions() - eup = elns.up_electron - edn = elns.down_electron + ions = system.get_ions() + eup = system.up_electron + edn = system.down_electron particlesets = [] eps = particleset( @@ -3960,10 +3958,10 @@ def incorporate_system(self,system): udet,ddet = self.get('updet','downdet') if udet is not None: - udet.size = elns.up_electron.count + udet.size = eup.count #end if if ddet is not None: - ddet.size = elns.down_electron.count + ddet.size = edn.count #end if if np.abs(net_spin) > 1e-1: @@ -5387,7 +5385,7 @@ def generate_simulationcell(bconds='ppp',lr_dim_cutoff=15,lr_tol=None,lr_handler structure = system.structure if isinstance(structure,Jellium): sc.rs = structure.rs() - sc.nparticles = system.particles.count_electrons() + sc.nparticles = system.count_electrons() else: #setting the 'lattice' (cell axes) requires some delicate care # qmcpack will fail if this is even 1e-10 off of what is in @@ -5440,15 +5438,13 @@ def generate_particlesets(electrons = 'e', system.check_folded_system() system.change_units('B') - particles = system.particles structure = system.structure net_charge = system.net_charge net_spin = system.net_spin - elns = particles.get_electrons() - ions = particles.get_ions() - eup = elns.up_electron - edn = elns.down_electron + ions = system.get_ions() + eup = system.up_electron + edn = system.down_electron use_spinor = spinor is not None and spinor @@ -5556,9 +5552,8 @@ def generate_sposets(type = None, if system is None and not have_counts: QmcpackInput.class_error('cannot generate sposets in occupation mode {0}\n arguments nup & ndown or system must be given to generate_sposets'.format(occupation)) elif not have_counts: - elns = system.particles.get_electrons() - nup = elns.up_electron.count - ndn = elns.down_electron.count + nup = system.up_electron.count + ndn = system.down_electron.count else: ndn = ndown #end if @@ -5805,9 +5800,8 @@ def generate_determinantset(up = 'u', if system is None: QmcpackInput.class_error('generate_determinantset argument system must not be None') #end if - elns = system.particles.get_electrons() - nup = elns.up_electron.count - ndn = elns.down_electron.count + nup = system.up_electron.count + ndn = system.down_electron.count use_spinor = spinor is not None and spinor if not spin_polarized and nup==ndn and not use_spinor: spo_u = 'spo_ud' @@ -5986,7 +5980,6 @@ def generate_determinantset_old(type = 'bspline', if system is None: QmcpackInput.class_error('generate_determinantset argument system must not be None') #end if - elns = system.particles.get_electrons() down_spin = 0 if spin_polarized: down_spin=1 @@ -5996,8 +5989,8 @@ def generate_determinantset_old(type = 'bspline', tilematrix = system.structure.tilematrix() #end if use_spinor = spinor is not None and spinor - nup = elns.up_electron.count - ndn = elns.down_electron.count + nup = system.up_electron.count + ndn = system.down_electron.count determinants_list = [] if not use_spinor: if nup > 0: @@ -6075,8 +6068,8 @@ def generate_determinantset_old(type = 'bspline', # Are there an equal number of up and down electrons? # If no, then exit. Currently, singlet and triplet # excitations are assumed to have ms = 0. - if elns.down_electron.count != elns.up_electron.count: - QmcpackInput.class_error('The \'singlet\' and \'triplet\' excitation types currently assume number of up and down electrons is the same for the reference ground state. Otherwise, one should use \'up\' or \'down\' types.\nFor your system: Nup={} and Ndown={}.\nWe plan to expand to additional cases in the future.'.format(elns.up_electron.count,elns.down_electron.count)) + if system.down_electron.count != system.up_electron.count: + QmcpackInput.class_error('The \'singlet\' and \'triplet\' excitation types currently assume number of up and down electrons is the same for the reference ground state. Otherwise, one should use \'up\' or \'down\' types.\nFor your system: Nup={} and Ndown={}.\nWe plan to expand to additional cases in the future.'.format(system.up_electron.count,system.down_electron.count)) #end if coeff_sign = '' @@ -6087,14 +6080,14 @@ def generate_determinantset_old(type = 'bspline', if down_spin: sposet_list = [sposet(name = 'spo_u', spindataset = 0, - size = elns.up_electron.count+1, + size = system.up_electron.count+1, occupation = section(mode='ground'), coefficient = section(size=90,spindataset=0), spos = '' ), sposet(name = 'spo_d', spindataset = 1, - size = elns.up_electron.count+1, + size = system.up_electron.count+1, occupation = section(mode='ground'), coefficient = section(spindataset=1), spos = '' @@ -6102,7 +6095,7 @@ def generate_determinantset_old(type = 'bspline', else: sposet_list = [sposet(name = 'spo_ud', spindataset = 0, - size = elns.up_electron.count+1, + size = system.up_electron.count+1, occupation = section(mode='ground'), coefficient = section(spindataset=0), spos = '' @@ -6127,8 +6120,8 @@ def generate_determinantset_old(type = 'bspline', type = 'CSF', nca = '0', ncb = '0', - nea = elns.up_electron.count, - neb = elns.down_electron.count, + nea = system.up_electron.count, + neb = system.down_electron.count, cutoff = '0.001', csf = csf( id = 'CSF_0', @@ -6153,7 +6146,7 @@ def generate_determinantset_old(type = 'bspline', if exc_type in (exc_types.energy,exc_types.lowest): - nup = elns.up_electron.count + nup = system.up_electron.count if exc_type==exc_types.lowest: exc_orbs = [nup,nup+1] else: @@ -6277,9 +6270,9 @@ def generate_determinantset_old(type = 'bspline', elif exc_type == exc_types.lowest: # Type 4 occ.format = 'energy' if exc_spin == exc_spins.up: - nel = elns.up_electron.count + nel = system.up_electron.count else: - nel = elns.down_electron.count + nel = system.down_electron.count #end if excitation = '-{} +{}'.format(nel,nel+1) occ.contents = '\n'+excitation+'\n' @@ -6319,17 +6312,16 @@ def generate_hamiltonian(name = 'h0', del pseudos del wavefunction - particles = system.particles - if particles.count_electrons()==0: + if system.count_electrons()==0: QmcpackInput.class_error('cannot generate hamiltonian, no electrons present') #end if pairpots = [] if interactions is not None: pairpots.append(coulomb(name='ElecElec',type='coulomb',source=ename,target=ename)) - if particles.count_ions()>0: + if system.count_ions()>0: pairpots.append(coulomb(name='IonIon',type='coulomb',source=iname,target=iname)) - ions = particles.get_ions() + ions = system.get_ions() if not system.pseudized: pairpots.append(coulomb(name='ElecIon',type='coulomb',source=iname,target=ename)) else: @@ -6577,7 +6569,7 @@ def generate_jastrows(jastrows,system=None,return_list=False,check_ions=False): jin = [] have_ions = True if check_ions and system is not None: - have_ions = system.particles.count_ions()>0 + have_ions = system.count_ions()>0 #end if if isinstance(jastrows,str): jorders = set(jastrows.replace('generate','')) @@ -6591,7 +6583,7 @@ def generate_jastrows(jastrows,system=None,return_list=False,check_ions=False): jterm = generate_jastrow('J3','polynomial',3,3,4.0,system=system) #end if if 'k' in jorders: - kcut = max(system.rpa_kf()) + kcut = max(system.kf_rpa()) nksh = system.structure.count_kshells(kcut) jterm = generate_kspace_jastrow(kc1=0, kc2=kcut, nk1=0, nk2=nksh) #end if @@ -6682,8 +6674,8 @@ def generate_jastrows_alt( openbc = system.structure.is_open() - natoms = system.particles.count_ions() - nelec = system.particles.count_electrons() + natoms = system.count_ions() + nelec = system.count_electrons() jastrows = [] J2 |= J3 @@ -6950,7 +6942,7 @@ def generate_bspline_jastrow2(size=8,rcut=None,coeff=None,spins=('u','d'),densit if rcut is None and isperiodic: rcut = rwigner #end if - nelectrons = system.particles.count_electrons() + nelectrons = system.count_electrons() density = nelectrons/volume #end if elif init=='rpa': @@ -7053,7 +7045,7 @@ def generate_jastrow2(function='bspline',*args,**kwargs): QmcpackInput.class_error('function is invalid\n you provided: {0}\n valid options are: bspline or pade'.format(function),'generate_jastrow2') #end if if 'system' in kwargs and kwargs['system'] is not None: - nup,ndn = kwargs['system'].particles.electron_counts() + nup,ndn = kwargs['system'].electron_counts() if nup<2: del j2.correlations.uu #end if diff --git a/nexus/nexus/quantum_package_input.py b/nexus/nexus/quantum_package_input.py index b24ece13d1..5ac4df8c98 100644 --- a/nexus/nexus/quantum_package_input.py +++ b/nexus/nexus/quantum_package_input.py @@ -750,7 +750,7 @@ def generate_quantum_package_input(**kwargs): if 'electrons' not in qpi: qpi.electrons = Section() #end if - nup,ndn = system.particles.electron_counts() + nup,ndn = system.electron_counts() qpi.electrons.elec_alpha_num = nup qpi.electrons.elec_beta_num = ndn diff --git a/nexus/nexus/rmg_input.py b/nexus/nexus/rmg_input.py index f0a00903a5..c5ed4d8031 100644 --- a/nexus/nexus/rmg_input.py +++ b/nexus/nexus/rmg_input.py @@ -3525,7 +3525,7 @@ def generate_any_rmg_input(**kwargs): del ri[k] #end if #end for - nup,ndn = system.particles.electron_counts() + nup,ndn = system.electron_counts() nvirt = int(np.ceil(virtual_frac*max(nup,ndn))) nptot = max(nup,ndn) + nvirt nup_virt = nptot-nup diff --git a/nexus/nexus/vasp_input.py b/nexus/nexus/vasp_input.py index 9a0b4f0e1e..2e9d7b19ac 100644 --- a/nexus/nexus/vasp_input.py +++ b/nexus/nexus/vasp_input.py @@ -1525,7 +1525,7 @@ def incorporate_system(self,system,incorp_kpoints=True,coord='cartesian',set_nel # handle charged systems if set_nelect or system.net_charge!=0: # warning: spin polarization is handled by the user! - self.incar.nelect = system.particles.count_electrons() + self.incar.nelect = system.count_electrons() #end if return species From 2356ef7d03c0449a38bd0d5763e00655e70bb239 Mon Sep 17 00:00:00 2001 From: brockdyer03 Date: Wed, 8 Jul 2026 10:59:49 -0400 Subject: [PATCH 2/4] Nexus: Add test for particle equiv --- nexus/nexus/tests/test_physical_system.py | 27 +++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/nexus/nexus/tests/test_physical_system.py b/nexus/nexus/tests/test_physical_system.py index fd438f1545..387700083e 100644 --- a/nexus/nexus/tests/test_physical_system.py +++ b/nexus/nexus/tests/test_physical_system.py @@ -8,6 +8,7 @@ import numpy as np from .. import testing from ..testing import value_eq,object_eq +from nexus.physical_system import generate_physical_system, Particles from .test_structure import structure_same @@ -511,3 +512,29 @@ 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.electron_counts()) + assert(ref.particles.count_electrons() == ref.count_electrons()) + assert(ref.particles.count_ions() == ref.count_ions()) + assert(ref.particles.count_ions(species=True) == ref.count_ions(species=True)) + assert(object_eq(ref.particles.get_ions(), ref.get_ions())) + assert(object_eq(ref.particles.up_electron, ref.up_electron)) + assert(object_eq(ref.particles.down_electron, ref.down_electron)) +#end def test_particle_equiv From e1c7996b2b58b831a6e97a323c978a25f3b7c1f5 Mon Sep 17 00:00:00 2001 From: brockdyer03 Date: Wed, 8 Jul 2026 17:24:33 -0400 Subject: [PATCH 3/4] Nexus: Update to address review comments --- nexus/nexus/gamess.py | 2 +- nexus/nexus/physical_system.py | 54 +++++++--- nexus/nexus/pwscf_input.py | 44 ++++---- nexus/nexus/qmcpack.py | 16 +-- nexus/nexus/qmcpack_input.py | 124 +++++++++++----------- nexus/nexus/quantum_package_input.py | 2 +- nexus/nexus/rmg_input.py | 2 +- nexus/nexus/tests/test_physical_system.py | 24 +++-- nexus/nexus/vasp_input.py | 2 +- 9 files changed, 153 insertions(+), 117 deletions(-) diff --git a/nexus/nexus/gamess.py b/nexus/nexus/gamess.py index e7940a468a..aad2c3b07a 100644 --- a/nexus/nexus/gamess.py +++ b/nexus/nexus/gamess.py @@ -154,7 +154,7 @@ def incorporate_result(self,result_name,result,sim): self.block() #end if guess_inputs = obj() - ecounts = self.system.electron_counts() + ecounts = self.system.n_up, self.system.n_down orbs = result.orbitals order_map = obj(up='iorder',down='jorder') nelec_map = obj(up=ecounts[0],down=ecounts[1]) diff --git a/nexus/nexus/physical_system.py b/nexus/nexus/physical_system.py index 856c70aa58..4bd4271801 100644 --- a/nexus/nexus/physical_system.py +++ b/nexus/nexus/physical_system.py @@ -636,26 +636,54 @@ def kf_rpa(self): #end def kf_rpa - def electron_counts(self): - return self.particles.electron_counts() - - def count_electrons(self): - return self.particles.count_electrons() + @property + def n_elec(self): + ions = self.structure.elem.tolist() + tot_charge = 0 + for ion in ions: + if hasattr(self, "valency"): + if ion in self.valency: + tot_charge += self.valency[ion] + else: + _, element = Elements.is_element(ion, return_element=True) + tot_charge += element.atomic_number + else: + _, element = Elements.is_element(ion, return_element=True) + tot_charge += element.atomic_number - def count_ions(self, species=False): - return self.particles.count_ions(species) + return tot_charge - self.net_charge - def get_ions(self): - return self.particles.get_ions() + @property + def n_up(self): + return (self.n_elec + self.net_spin) // 2 @property - def up_electron(self): - return self.particles.up_electron + def n_down(self): + return (self.n_elec - self.net_spin) // 2 @property - def down_electron(self): - return self.particles.down_electron + def n_species(self): + return len(set(self.structure.elem)) + + @property + def n_ions(self): + return len(self.structure.elem) + + @property + def ion_labels(self): + return set(self.structure.elem) + + @property + def Zeff(self): + if hasattr(self, "valency"): + return self.valency + + Zeff = {} + for ion in self.ion_labels: + _, element = Elements.is_element(ion, return_element=True) + Zeff[ion] = element.atomic_number + return Zeff #end class PhysicalSystem diff --git a/nexus/nexus/pwscf_input.py b/nexus/nexus/pwscf_input.py index c3bde26283..f1bc2f3410 100644 --- a/nexus/nexus/pwscf_input.py +++ b/nexus/nexus/pwscf_input.py @@ -1679,12 +1679,13 @@ def incorporate_system(self,system,elem_order=None): nc = system.net_charge ns = system.net_spin - nup = system.up_electron.count - ndn = system.down_electron.count + nup = system.n_up + ndn = system.n_down self.system.ibrav = 0 # self.system['celldm(1)'] = 1.0e0 - nions,nspecies = system.count_ions(species=True) + nions = system.n_ions + nspecies = system.n_species self.system.nat = nions self.system.ntyp = nspecies self.system.tot_charge = nc @@ -1716,19 +1717,19 @@ def incorporate_system(self,system,elem_order=None): #end if #end if - atoms = system.get_ions() if 'masses' not in self.atomic_species: self.atomic_species.masses = obj() #end if - for name,a in atoms.items(): - self.atomic_species.masses[name] = convert(a.mass,'me','amu') + for name in system.ion_labels: + is_elem, element = Elements.is_element(name, return_element=True) + self.atomic_species.masses[name] = element.atomic_weight #end for if elem_order is None: - self.atomic_species.atoms = list(sorted(atoms.keys())) + self.atomic_species.atoms = list(sorted(system.ion_labels)) else: - if set(elem_order)!=set(atoms.keys()): - self.error('elem_order is missing some atomic species\natomic species present: {0}\nelem_order: {1}'.format(sorted(atoms.keys()),elem_order)) - elif len(elem_order)!=len(atoms): + if set(elem_order)!=set(system.ion_labels): + self.error('elem_order is missing some atomic species\natomic species present: {0}\nelem_order: {1}'.format(sorted(masses.keys()),elem_order)) + elif len(elem_order)!=system.n_ions: self.error('elem_order has repeated elements\nelem_order: {0}'.format(elem_order)) #end if self.atomic_species.atoms = list(elem_order) @@ -1738,7 +1739,8 @@ def incorporate_system(self,system,elem_order=None): pp = self.atomic_species.pseudopotentials for atom in self.atomic_species.atoms: if atom not in pp: - iselem,symbol = system.is_element(atom,symbol=True) + iselem,element = Elements.is_element(atom,return_element=True) + symbol = element.symbol if iselem and symbol in pp: pp[atom] = str(pp[symbol]) #end if @@ -1773,12 +1775,13 @@ def incorporate_system_old(self,system,spin_polarized=None): nc = system.net_charge ns = system.net_spin - nup = system.up_electron.count - ndn = system.down_electron.count + nup = system.n_up + ndn = system.n_down self.system.ibrav = 0 # self.system['celldm(1)'] = 1.0e0 - nions,nspecies = system.count_ions(species=True) + nions = system.n_ions + nspecies = system.n_species self.system.nat = nions self.system.ntyp = nspecies #self.system.nelec = nup+ndn @@ -1817,18 +1820,19 @@ def incorporate_system_old(self,system,spin_polarized=None): #end if #end if - atoms = system.get_ions() masses = obj() - for name,a in atoms.items(): - masses[name] = convert(a.mass,'me','amu') + for name in system.ion_labels: + is_elem, element = Elements.is_element(name, return_element=True) + masses[name] = element.atomic_weight #end for - self.atomic_species.atoms = list(sorted(atoms.keys())) + self.atomic_species.atoms = list(sorted(system.ion_labels)) self.atomic_species.masses = masses # set pseudopotentials for renamed atoms (e.g. Cu3 is same as Cu) pp = self.atomic_species.pseudopotentials for atom in self.atomic_species.atoms: if atom not in pp: - iselem,symbol = system.is_element(atom,symbol=True) + iselem,element = Elements.is_element(atom,return_element=True) + symbol = element.symbol if iselem and symbol in pp: pp[atom] = str(pp[symbol]) #end if @@ -2263,7 +2267,7 @@ def generate_any_pwscf_input(**kwargs): # set nbnd using bandfac, if provided if nbnd is None and bandfac is not None: - nocc = max(system.electron_counts()) + nocc = max(system.n_up, system.n_down) pw.system.nbnd = int(np.ceil(nocc*bandfac)) #end if diff --git a/nexus/nexus/qmcpack.py b/nexus/nexus/qmcpack.py index 0c7f49bef0..b0c4f57823 100644 --- a/nexus/nexus/qmcpack.py +++ b/nexus/nexus/qmcpack.py @@ -188,11 +188,11 @@ def unfolded_nelecs(self): Returns the number of electrons in the primitive cell ''' if self.system.folded_system is None: - n_up = self.system.up_electron.count - n_dn = self.system.down_electron.count + n_up = self.system.n_up + n_dn = self.system.n_down else: - n_up = self.system.folded_system.up_electron.count - n_dn = self.system.folded_system.down_electron.count + n_up = self.system.folded_system.n_up + n_dn = self.system.folded_system.n_down #end if nelecs = n_up + n_dn return nelecs @@ -463,8 +463,8 @@ def sum_charge_twists(self): ''' Returns the net charge of a system with multiple twists (not averaged) ''' - n_up = self.system.up_electron.count - n_dn = self.system.down_electron.count + n_up = self.system.n_up + n_dn = self.system.n_down n_total = n_up + n_dn nelecs_at_twist = self.nelecs_at_twist kweights = np.array(self.system.structure.kweights) @@ -536,8 +536,8 @@ def write_gcta_report(self, locdir, fermi_level, scf_magnet = None): spin_sum_twists = self.sum_spin_twists() qmc_magnet = spin_sum_twists / nosym_kpoints #end if - n_up = self.system.up_electron.count - n_dn = self.system.down_electron.count + n_up = self.system.n_up + n_dn = self.system.n_down n_total = n_up + n_dn nelecs_at_twist = self.nelecs_at_twist fermi_level = np.array(fermi_level) diff --git a/nexus/nexus/qmcpack_input.py b/nexus/nexus/qmcpack_input.py index 9f91a223c9..cade1fa5a2 100644 --- a/nexus/nexus/qmcpack_input.py +++ b/nexus/nexus/qmcpack_input.py @@ -147,6 +147,7 @@ from .simulation import SimulationInput, SimulationInputTemplate from .pwscf_input import array_to_string as pwscf_array_string from .utilities import path_string +from .unit_converter import convert from . import numpy_extensions as npe yesno_dict = {True:'yes' ,False:'no'} @@ -3897,20 +3898,17 @@ def incorporate_system(self,system): sc.lattice = axes #end if - ions = system.get_ions() - eup = system.up_electron - edn = system.down_electron particlesets = [] eps = particleset( name='e',random=True, groups = [ - group(name='u',charge=-1,mass=eup.mass,size=eup.count), - group(name='d',charge=-1,mass=edn.mass,size=edn.count) + group(name='u',charge=-1,mass=1.0,size=system.n_up), + group(name='d',charge=-1,mass=1.0,size=system.n_down) ] ) particlesets.append(eps) - if len(ions)>0: + if system.n_ions>0: if sc is not None and 'bconds' in sc and tuple(sc.bconds)!=('p','p','p'): eps.randomsrc = 'ion0' #end if @@ -3927,20 +3925,21 @@ def incorporate_system(self,system): pp.pseudos = pseudos #end if #end if - for ion in ions: - gpos = pos[elem==ion.name] + for ion in system.ion_labels: + gpos = pos[elem==ion] + is_elem, element = Elements.is_element(ion, return_element=True) g = group( - name = ion.name, - charge = ion.charge, - valence = ion.charge, - atomicnumber = ion.protons, - mass = ion.mass, + name = ion, + charge = system.Zeff[ion], + valence = system.Zeff[ion], + atomicnumber = element.atomic_number, + mass = convert(element.atomic_weight, "amu", "me"), position = gpos, size = len(gpos) ) groups.append(g) - if pseudos is not None and ion.name not in pseudos: - pseudos[ion.name] = pseudo(elementtype=ion.name,href='MISSING.xml') + if pseudos is not None and ion not in pseudos: + pseudos[ion] = pseudo(elementtype=ion,href='MISSING.xml') #end if #end for ips.groups = make_collection(groups) @@ -3951,17 +3950,17 @@ def incorporate_system(self,system): if old_eps_name is not None: self.replace(old_eps_name,'e') #end if - if old_ips_name is not None and len(ions)>0: + if old_ips_name is not None and system.n_ions>0: self.replace(old_ips_name,'ion0') #end if udet,ddet = self.get('updet','downdet') if udet is not None: - udet.size = eup.count + udet.size = system.n_up #end if if ddet is not None: - ddet.size = edn.count + ddet.size = system.n_down #end if if np.abs(net_spin) > 1e-1: @@ -5385,7 +5384,7 @@ def generate_simulationcell(bconds='ppp',lr_dim_cutoff=15,lr_tol=None,lr_handler structure = system.structure if isinstance(structure,Jellium): sc.rs = structure.rs() - sc.nparticles = system.count_electrons() + sc.nparticles = system.n_elec else: #setting the 'lattice' (cell axes) requires some delicate care # qmcpack will fail if this is even 1e-10 off of what is in @@ -5439,27 +5438,24 @@ def generate_particlesets(electrons = 'e', system.change_units('B') structure = system.structure - net_charge = system.net_charge - net_spin = system.net_spin - ions = system.get_ions() - eup = system.up_electron - edn = system.down_electron + eup = system.n_up + edn = system.n_down use_spinor = spinor is not None and spinor particleset_groups = [] if not use_spinor: - if eup.count > 0: - particleset_groups.append(group(name=uname,charge=-1,mass=eup.mass,size=eup.count)) + if eup > 0: + particleset_groups.append(group(name=uname,charge=-1,mass=1.0,size=eup)) #end if - if edn.count > 0: - particleset_groups.append(group(name=dname,charge=-1,mass=edn.mass,size=edn.count)) + if edn > 0: + particleset_groups.append(group(name=dname,charge=-1,mass=1.0,size=edn)) #end if else: - ecount = eup.count+edn.count + ecount = eup+edn if ecount>0: - particleset_groups.append(group(name=uname,charge=-1,mass=eup.mass,size=ecount)) + particleset_groups.append(group(name=uname,charge=-1,mass=1.0,size=ecount)) #end if #end if @@ -5473,7 +5469,7 @@ def generate_particlesets(electrons = 'e', eps.spinor = True #end if particlesets.append(eps) - if len(ions)>0: + if system.n_ions>0: # maintain consistent order ion_species,ion_counts = structure.order_by_species() elem = structure.elem @@ -5501,14 +5497,14 @@ def generate_particlesets(electrons = 'e', # make groups groups = [] for ion_spec in ion_species: - ion = ions[ion_spec] - gpos = pos[elem==ion.name] + is_elem, element = Elements.is_element(ion_spec, return_element=True) + gpos = pos[elem==ion_spec] g = group( - name = ion.name, - charge = ion.charge, - valence = ion.charge, - atomicnumber = ion.protons, - mass = ion.mass, + name = ion_spec, + charge = system.Zeff[ion_spec], + valence = system.Zeff[ion_spec], + atomicnumber = element.atomic_number, + mass = convert(element.atomic_weight, "amu", "me"), position = gpos, size = len(gpos) ) @@ -5552,8 +5548,8 @@ def generate_sposets(type = None, if system is None and not have_counts: QmcpackInput.class_error('cannot generate sposets in occupation mode {0}\n arguments nup & ndown or system must be given to generate_sposets'.format(occupation)) elif not have_counts: - nup = system.up_electron.count - ndn = system.down_electron.count + nup = system.n_up + ndn = system.n_down else: ndn = ndown #end if @@ -5800,8 +5796,8 @@ def generate_determinantset(up = 'u', if system is None: QmcpackInput.class_error('generate_determinantset argument system must not be None') #end if - nup = system.up_electron.count - ndn = system.down_electron.count + nup = system.n_up + ndn = system.n_down use_spinor = spinor is not None and spinor if not spin_polarized and nup==ndn and not use_spinor: spo_u = 'spo_ud' @@ -5989,8 +5985,8 @@ def generate_determinantset_old(type = 'bspline', tilematrix = system.structure.tilematrix() #end if use_spinor = spinor is not None and spinor - nup = system.up_electron.count - ndn = system.down_electron.count + nup = system.n_up + ndn = system.n_down determinants_list = [] if not use_spinor: if nup > 0: @@ -6068,7 +6064,7 @@ def generate_determinantset_old(type = 'bspline', # Are there an equal number of up and down electrons? # If no, then exit. Currently, singlet and triplet # excitations are assumed to have ms = 0. - if system.down_electron.count != system.up_electron.count: + if system.n_down != system.n_up: QmcpackInput.class_error('The \'singlet\' and \'triplet\' excitation types currently assume number of up and down electrons is the same for the reference ground state. Otherwise, one should use \'up\' or \'down\' types.\nFor your system: Nup={} and Ndown={}.\nWe plan to expand to additional cases in the future.'.format(system.up_electron.count,system.down_electron.count)) #end if @@ -6080,14 +6076,14 @@ def generate_determinantset_old(type = 'bspline', if down_spin: sposet_list = [sposet(name = 'spo_u', spindataset = 0, - size = system.up_electron.count+1, + size = system.n_up+1, occupation = section(mode='ground'), coefficient = section(size=90,spindataset=0), spos = '' ), sposet(name = 'spo_d', spindataset = 1, - size = system.up_electron.count+1, + size = system.n_up+1, occupation = section(mode='ground'), coefficient = section(spindataset=1), spos = '' @@ -6095,7 +6091,7 @@ def generate_determinantset_old(type = 'bspline', else: sposet_list = [sposet(name = 'spo_ud', spindataset = 0, - size = system.up_electron.count+1, + size = system.n_up+1, occupation = section(mode='ground'), coefficient = section(spindataset=0), spos = '' @@ -6120,8 +6116,8 @@ def generate_determinantset_old(type = 'bspline', type = 'CSF', nca = '0', ncb = '0', - nea = system.up_electron.count, - neb = system.down_electron.count, + nea = system.n_up, + neb = system.n_down, cutoff = '0.001', csf = csf( id = 'CSF_0', @@ -6146,7 +6142,7 @@ def generate_determinantset_old(type = 'bspline', if exc_type in (exc_types.energy,exc_types.lowest): - nup = system.up_electron.count + nup = system.n_up if exc_type==exc_types.lowest: exc_orbs = [nup,nup+1] else: @@ -6270,9 +6266,9 @@ def generate_determinantset_old(type = 'bspline', elif exc_type == exc_types.lowest: # Type 4 occ.format = 'energy' if exc_spin == exc_spins.up: - nel = system.up_electron.count + nel = system.n_up else: - nel = system.down_electron.count + nel = system.n_down #end if excitation = '-{} +{}'.format(nel,nel+1) occ.contents = '\n'+excitation+'\n' @@ -6312,16 +6308,16 @@ def generate_hamiltonian(name = 'h0', del pseudos del wavefunction - if system.count_electrons()==0: + if system.n_elec==0: QmcpackInput.class_error('cannot generate hamiltonian, no electrons present') #end if pairpots = [] if interactions is not None: pairpots.append(coulomb(name='ElecElec',type='coulomb',source=ename,target=ename)) - if system.count_ions()>0: + if system.n_ions>0: pairpots.append(coulomb(name='IonIon',type='coulomb',source=iname,target=iname)) - ions = system.get_ions() + ions = system.ion_labels if not system.pseudized: pairpots.append(coulomb(name='ElecIon',type='coulomb',source=iname,target=ename)) else: @@ -6346,16 +6342,15 @@ def generate_hamiltonian(name = 'h0', #end if pseudos = collection() for ion in ions: - label = ion.name - iselem, element = Elements.is_element(ion.name, return_element=True) - if label in ppfiles: - ppfile = ppfiles[label] + iselem, element = Elements.is_element(ion, return_element=True) + if ion in ppfiles: + ppfile = ppfiles[ion] elif element.symbol in ppfiles: ppfile = ppfiles[element.symbol] else: QmcpackInput.class_error('pseudos provided to generate_hamiltonian are incomplete\n a pseudopotential for ion of type {0} is missing\n pseudos provided:\n{1}'.format(ion.name,str(ppfiles))) #end if - pseudos.add(pseudo(elementtype=label,href=ppfile)) + pseudos.add(pseudo(elementtype=ion,href=ppfile)) #end for pp = pseudopotential(name='PseudoPot',type='pseudo',source=iname,wavefunction=wfname,format=format,pseudos=pseudos) if algorithm is not None: @@ -6569,7 +6564,7 @@ def generate_jastrows(jastrows,system=None,return_list=False,check_ions=False): jin = [] have_ions = True if check_ions and system is not None: - have_ions = system.count_ions()>0 + have_ions = system.n_ions>0 #end if if isinstance(jastrows,str): jorders = set(jastrows.replace('generate','')) @@ -6942,7 +6937,7 @@ def generate_bspline_jastrow2(size=8,rcut=None,coeff=None,spins=('u','d'),densit if rcut is None and isperiodic: rcut = rwigner #end if - nelectrons = system.count_electrons() + nelectrons = system.n_elec density = nelectrons/volume #end if elif init=='rpa': @@ -7045,7 +7040,8 @@ def generate_jastrow2(function='bspline',*args,**kwargs): QmcpackInput.class_error('function is invalid\n you provided: {0}\n valid options are: bspline or pade'.format(function),'generate_jastrow2') #end if if 'system' in kwargs and kwargs['system'] is not None: - nup,ndn = kwargs['system'].electron_counts() + system = kwargs['system'] + nup,ndn = system.n_up, system.n_down if nup<2: del j2.correlations.uu #end if diff --git a/nexus/nexus/quantum_package_input.py b/nexus/nexus/quantum_package_input.py index 5ac4df8c98..87351cb9c6 100644 --- a/nexus/nexus/quantum_package_input.py +++ b/nexus/nexus/quantum_package_input.py @@ -750,7 +750,7 @@ def generate_quantum_package_input(**kwargs): if 'electrons' not in qpi: qpi.electrons = Section() #end if - nup,ndn = system.electron_counts() + nup,ndn = system.n_up, system.n_down qpi.electrons.elec_alpha_num = nup qpi.electrons.elec_beta_num = ndn diff --git a/nexus/nexus/rmg_input.py b/nexus/nexus/rmg_input.py index c5ed4d8031..f6ec3bdb26 100644 --- a/nexus/nexus/rmg_input.py +++ b/nexus/nexus/rmg_input.py @@ -3525,7 +3525,7 @@ def generate_any_rmg_input(**kwargs): del ri[k] #end if #end for - nup,ndn = system.electron_counts() + nup,ndn = system.n_up, system.n_down nvirt = int(np.ceil(virtual_frac*max(nup,ndn))) nptot = max(nup,ndn) + nvirt nup_virt = nptot-nup diff --git a/nexus/nexus/tests/test_physical_system.py b/nexus/nexus/tests/test_physical_system.py index 387700083e..2a015ffe71 100644 --- a/nexus/nexus/tests/test_physical_system.py +++ b/nexus/nexus/tests/test_physical_system.py @@ -8,7 +8,9 @@ import numpy as np from .. import testing from ..testing import value_eq,object_eq -from nexus.physical_system import generate_physical_system, Particles +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 @@ -530,11 +532,17 @@ def test_particle_equiv(): C = 4, ) - assert(ref.particles.electron_counts() == ref.electron_counts()) - assert(ref.particles.count_electrons() == ref.count_electrons()) - assert(ref.particles.count_ions() == ref.count_ions()) - assert(ref.particles.count_ions(species=True) == ref.count_ions(species=True)) - assert(object_eq(ref.particles.get_ions(), ref.get_ions())) - assert(object_eq(ref.particles.up_electron, ref.up_electron)) - assert(object_eq(ref.particles.down_electron, ref.down_electron)) + 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 diff --git a/nexus/nexus/vasp_input.py b/nexus/nexus/vasp_input.py index 2e9d7b19ac..b67554496f 100644 --- a/nexus/nexus/vasp_input.py +++ b/nexus/nexus/vasp_input.py @@ -1525,7 +1525,7 @@ def incorporate_system(self,system,incorp_kpoints=True,coord='cartesian',set_nel # handle charged systems if set_nelect or system.net_charge!=0: # warning: spin polarization is handled by the user! - self.incar.nelect = system.count_electrons() + self.incar.nelect = system.n_elec #end if return species From 8d4179d21068d01a5ec20a34e2c38cb54888e563 Mon Sep 17 00:00:00 2001 From: brockdyer03 Date: Thu, 9 Jul 2026 11:28:13 -0400 Subject: [PATCH 4/4] Nexus: Address review comments --- nexus/nexus/pwscf_input.py | 2 +- nexus/nexus/qmcpack_input.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/nexus/nexus/pwscf_input.py b/nexus/nexus/pwscf_input.py index f1bc2f3410..86d43076d0 100644 --- a/nexus/nexus/pwscf_input.py +++ b/nexus/nexus/pwscf_input.py @@ -1728,7 +1728,7 @@ def incorporate_system(self,system,elem_order=None): self.atomic_species.atoms = list(sorted(system.ion_labels)) else: if set(elem_order)!=set(system.ion_labels): - self.error('elem_order is missing some atomic species\natomic species present: {0}\nelem_order: {1}'.format(sorted(masses.keys()),elem_order)) + self.error('elem_order is missing some atomic species\natomic species present: {0}\nelem_order: {1}'.format(sorted(system.ion_labels),elem_order)) elif len(elem_order)!=system.n_ions: self.error('elem_order has repeated elements\nelem_order: {0}'.format(elem_order)) #end if diff --git a/nexus/nexus/qmcpack_input.py b/nexus/nexus/qmcpack_input.py index cade1fa5a2..f456b1a2dd 100644 --- a/nexus/nexus/qmcpack_input.py +++ b/nexus/nexus/qmcpack_input.py @@ -6065,7 +6065,7 @@ def generate_determinantset_old(type = 'bspline', # If no, then exit. Currently, singlet and triplet # excitations are assumed to have ms = 0. if system.n_down != system.n_up: - QmcpackInput.class_error('The \'singlet\' and \'triplet\' excitation types currently assume number of up and down electrons is the same for the reference ground state. Otherwise, one should use \'up\' or \'down\' types.\nFor your system: Nup={} and Ndown={}.\nWe plan to expand to additional cases in the future.'.format(system.up_electron.count,system.down_electron.count)) + QmcpackInput.class_error('The \'singlet\' and \'triplet\' excitation types currently assume number of up and down electrons is the same for the reference ground state. Otherwise, one should use \'up\' or \'down\' types.\nFor your system: Nup={} and Ndown={}.\nWe plan to expand to additional cases in the future.'.format(system.n_up,system.n_down)) #end if coeff_sign = '' @@ -6669,8 +6669,8 @@ def generate_jastrows_alt( openbc = system.structure.is_open() - natoms = system.count_ions() - nelec = system.count_electrons() + natoms = system.n_ions + nelec = system.n_elec jastrows = [] J2 |= J3