1################################################################## 2## (c) Copyright 2015- by Jaron T. Krogel ## 3################################################################## 4 5 6#====================================================================# 7# physical_system.py # 8# Representations of matter, particles, and particles collected # 9# together in complete systems. # 10# # 11# Content summary: # 12# PhysicalSystem # 13# Class representing electrons+ions for a simulation. # 14# # 15# generate_physical_system # 16# User function to create arbitrary physical systems. # 17# # 18# Matter # 19# Base class for all forms of matter. # 20# Class contains a list of all matter known to Nexus. # 21# # 22# Particle # 23# Class representing a particular particle species. # 24# # 25# Ion, PseudoIon # 26# Specialized Particle classes for ion species and ions # 27# represented by pseudopotentials. # 28# # 29# Particles # 30# A collection of particles. # 31# PhysicalSystem objects contain a Particles instance. # 32# # 33#====================================================================# 34 35 36import numpy as np 37from numpy import dot,array 38from numpy.linalg import inv 39from generic import obj 40from developer import DevBase 41from unit_converter import convert 42from periodic_table import is_element 43from structure import Structure 44from debug import * 45 46 47class Matter(DevBase): 48 particle_collection = None 49 50 @classmethod 51 def set_elements(cls,elements): 52 cls.elements = set(elements) 53 #end def set_elements 54 55 @classmethod 56 def set_particle_collection(cls,pc): 57 cls.particle_collection = pc 58 #end def set_particle_collection 59 60 @classmethod 61 def new_particles(cls,*particles,**named_particles): 62 cls.particle_collection.add_particles(*particles,**named_particles) 63 #end def new_particles 64 65 def is_element(self,name,symbol=False): 66 return is_element(name,symbol=symbol) 67 #end def is_element 68#end class Matter 69 70 71class Particle(Matter): 72 def __init__(self,name=None,mass=None,charge=None,spin=None): 73 self.name = name 74 self.mass = mass 75 self.charge = charge 76 self.spin = spin 77 #end def __init__ 78 79 def set_count(self,count): 80 self.count = count 81 #end def set_count 82#end class Particle 83 84 85class Ion(Particle): 86 def __init__(self,name=None,mass=None,charge=None,spin=None, 87 protons=None,neutrons=None): 88 Particle.__init__(self,name,mass,charge,spin) 89 self.protons = protons 90 self.neutrons = neutrons 91 #end def __init__ 92 93 def pseudize(self,valence): 94 ps = PseudoIon() 95 ps.transfer_from(self) 96 ps.charge = valence 97 ps.core_electrons = ps.protons - valence 98 return ps 99 #end def pseudize 100#end class Ion 101 102 103class PseudoIon(Ion): 104 def __init__(self,name=None,mass=None,charge=None,spin=None, 105 protons=None,neutrons=None,core_electrons=None): 106 Ion.__init__(self,name,mass,charge,spin,protons,neutrons) 107 self.core_electrons = core_electrons 108 #end def __init__ 109#end class PseudoIon 110 111 112class Particles(Matter): 113 def __init__(self,*particles,**named_particles): 114 self.add_particles(*particles,**named_particles) 115 #end def __init__ 116 117 def add_particles(self,*particles,**named_particles): 118 if len(particles)==1 and isinstance(particles[0],list): 119 particles = particles[0] 120 #end if 121 for particle in particles: 122 self[particle.name] = particle 123 #end for 124 for name,particle in named_particles.items(): 125 self[name] = particle 126 #end for 127 #end def add_particles 128 129 def get_particle(self,name): 130 p = None 131 if name in self: 132 p = self[name] 133 else: 134 iselem,symbol = self.is_element(name,symbol=True) 135 if iselem and symbol in self: 136 p = self[symbol].copy() 137 p.name = name 138 self[name] = p 139 #end if 140 #end if 141 return p 142 #end def get_particle 143 144 # test needed 145 def get(self,quantity): 146 q = obj() 147 for name,particle in self.items(): 148 q[name] = particle[quantity] 149 #end for 150 return q 151 #end def get 152 153 def rename(self,**name_pairs): 154 for old,new in name_pairs.items(): 155 if old in self: 156 o = self[old] 157 o.name = new 158 del self[old] 159 if new in self: 160 self[new].count += o.count 161 else: 162 self[new] = o 163 #end if 164 #end if 165 #end for 166 #end def rename 167 168 def get_ions(self): 169 ions = obj() 170 for name,particle in self.items(): 171 if self.is_element(name): 172 ions[name] = particle 173 #end if 174 #end for 175 return ions 176 #end def get_ions 177 178 def count_ions(self,species=False): 179 nions = 0 180 nspecies = 0 181 for name,particle in self.items(): 182 if self.is_element(name): 183 nspecies += 1 184 nions += particle.count 185 #end if 186 #end for 187 if species: 188 return nions,nspecies 189 else: 190 return nions 191 #end if 192 #end def count_ions 193 194 def get_electrons(self): 195 electrons = obj() 196 for electron in ('up_electron','down_electron'): 197 if electron in self: 198 electrons[electron] = self[electron] 199 #end if 200 #end for 201 return electrons 202 #end def get_electrons 203 204 def count_electrons(self): 205 nelectrons = 0 206 for electron in ('up_electron','down_electron'): 207 if electron in self: 208 nelectrons += self[electron].count 209 #end if 210 #end for 211 return nelectrons 212 #end def count_electrons 213 214 def electron_counts(self): 215 counts = [] 216 for electron in ('up_electron','down_electron'): 217 if electron in self: 218 counts.append(self[electron].count) 219 else: 220 counts.append(0) 221 #end if 222 #end for 223 return counts 224 #end def electron_counts 225#end class Particles 226 227me_amu = convert(1.,'me','amu') 228 229plist = [ 230 Particle('up_electron' ,1.0,-1, 1), 231 Particle('down_electron',1.0,-1,-1), 232 ] 233from periodic_table import ptable 234for name,a in ptable.elements.items(): 235 spin = 0 # don't have this data 236 protons = a.atomic_number 237 neutrons = int(round(a.atomic_weight['amu']-a.atomic_number)) 238 p = Ion(a.symbol,a.atomic_weight['me'],a.atomic_number,spin,protons,neutrons) 239 plist.append(p) 240#end for 241for name,iso in ptable.isotopes.items(): 242 for mass_number,a in iso.items(): 243 spin = 0 # don't have this data 244 protons = a.atomic_number 245 neutrons = int(round(a.atomic_weight['amu']-a.atomic_number)) 246 p = Ion(a.symbol+'_'+str(mass_number),a.atomic_weight['me'],a.atomic_number,spin,protons,neutrons) 247 plist.append(p) 248 #end for 249#end for 250 251Matter.set_elements(ptable.elements.keys()) 252Matter.set_particle_collection(Particles(plist)) 253 254del plist 255 256 257 258class PhysicalSystem(Matter): 259 260 def __init__(self,structure=None,net_charge=0,net_spin=0,particles=None,**valency): 261 self.pseudized = False 262 263 if structure is None: 264 self.structure = Structure() 265 else: 266 self.structure = structure 267 #end if 268 if particles is None: 269 self.particles = Particles() 270 else: 271 self.particles = particles.copy() 272 #end if 273 274 self.folded_system = None 275 if self.structure.has_folded(): 276 if self.structure.is_tiled(): 277 vratio = structure.volume()/structure.folded_structure.volume() 278 ncells = int(round(vratio)) 279 if abs(vratio-ncells)>1e-4: 280 self.error('volume of system does not divide evenly into folded system') 281 #end if 282 if net_charge%ncells!=0: 283 self.error('net charge of system does not divide evenly into folded system') 284 #end if 285 if isinstance(net_spin,str): 286 net_spin_fold = net_spin 287 elif net_spin%ncells!=0: 288 self.error('net_spin of system does not divide evenly into folded system') 289 else: 290 net_spin_fold = net_spin//ncells 291 #end if 292 net_charge_fold = net_charge//ncells 293 elif not self.structure.has_axes(): # folded molecule 294 # net charge/spin are not physically meaningful 295 # for a point group folded molecule 296 # set them to safe values; they will not be used later 297 net_charge_fold = 0 298 net_spin_fold = 'low' 299 else: 300 self.error('folded structure is not correctly integrated with full structure\nfolded physical system cannot be constructed') 301 #end if 302 303 self.folded_system = PhysicalSystem( 304 structure = structure.folded_structure, 305 net_charge = net_charge_fold, 306 net_spin = net_spin_fold, 307 particles = particles, 308 **valency 309 ) 310 #end if 311 312 self.valency_in = obj(**valency) 313 self.net_charge_in = net_charge 314 self.net_spin_in = net_spin 315 316 self.update_particles(clear=False) 317 318 self.check_folded_system() 319 #end def __init__ 320 321 322 def update_particles(self,clear=True): 323 #add ions 324 pc = dict() 325 elem = list(self.structure.elem) 326 for ion in set(elem): 327 pc[ion] = elem.count(ion) 328 #end for 329 missing = set(pc.keys())-set(self.particles.keys()) 330 if len(missing)>0 or len(elem)==0: 331 if clear: 332 self.particles.clear() 333 #end if 334 self.add_particles(**pc) 335 336 #pseudize 337 if len(self.valency_in)>0: 338 self.pseudize(**self.valency_in) 339 #end if 340 341 #add electrons 342 self.generate_electrons(self.net_charge_in,self.net_spin_in) 343 #end if 344 #end def update_particles 345 346 347 def update(self): 348 self.net_charge = self.structure.background_charge 349 self.net_spin = 0 350 for p in self.particles: 351 self.net_charge += p.count*p.charge 352 self.net_spin += p.count*p.spin 353 #end for 354 self.net_charge = int(round(float(self.net_charge))) 355 self.net_spin = int(round(float(self.net_spin))) 356 #end def update 357 358 359 def add_particles(self,**particle_counts): 360 pc = self.particle_collection # all known particles 361 plist = [] 362 for name,count in particle_counts.items(): 363 particle = pc.get_particle(name) 364 if particle is None: 365 self.error('particle {0} is unknown'.format(name)) 366 else: 367 particle = particle.copy() 368 #end if 369 particle.set_count(count) 370 plist.append(particle) 371 #end for 372 self.particles.add_particles(plist) 373 self.update() 374 #end def add_particles 375 376 377 def generate_electrons(self,net_charge=0,net_spin=0): 378 nelectrons = -net_charge + self.net_charge 379 if net_spin=='low': 380 net_spin = nelectrons%2 381 #end if 382 nup = float(nelectrons + net_spin - self.net_spin)/2 383 ndown = float(nelectrons - net_spin + self.net_spin)/2 384 if abs(nup-int(nup))>1e-3: 385 self.error('requested spin state {0} incompatible with {1} electrons'.format(net_spin,nelectrons)) 386 #end if 387 nup = int(nup) 388 ndown = int(ndown) 389 self.add_particles(up_electron=nup,down_electron=ndown) 390 #end def generate_electrons 391 392 393 def pseudize(self,**valency): 394 errors = False 395 for ion,valence_charge in valency.items(): 396 if ion in self.particles: 397 ionp = self.particles[ion] 398 if isinstance(ionp,Ion): 399 self.particles[ion] = ionp.pseudize(valence_charge) 400 self.pseudized = True 401 else: 402 self.error(ion+' cannot be pseudized',exit=False) 403 #end if 404 else: 405 self.error(ion+' is not in the physical system',exit=False) 406 errors = True 407 #end if 408 #end for 409 if errors: 410 self.error('system cannot be generated') 411 #end if 412 self.valency = obj(**valency) 413 self.update() 414 #end def pseudize 415 416 417 def check_folded_system(self,exit=True,message=False): 418 msg = '' 419 sys_folded = self.folded_system!=None 420 struct_folded = self.structure.folded_structure!=None 421 if sys_folded!=struct_folded: 422 msg+='folding of physical system and structure is not consistent\nsystem folded: {0}\nstructure folded: {1}\n'.format(sys_folded,struct_folded) 423 #end if 424 if sys_folded and id(self.structure.folded_structure)!=id(self.folded_system.structure): 425 msg+='structure of folded system and folded structure are distinct\nthis is not allowed and may be a developer error' 426 #end if 427 success = len(msg)==0 428 if not success and exit: 429 self.error(msg) 430 #end if 431 if not message: 432 return success 433 else: 434 return success,msg 435 #end if 436 #end def check_folded_system 437 438 439 def check_consistent(self,tol=1e-8,exit=True,message=False): 440 fs,fm = self.check_folded_system(exit=False,message=True) 441 cs,cm = self.structure.check_consistent(tol,exit=False,message=True) 442 msg = '' 443 if not fs: 444 msg += fm+'\n' 445 #end if 446 if not cs: 447 msg += cm+'\n' 448 #end if 449 consistent = len(msg)==0 450 if not consistent and exit: 451 self.error(msg) 452 #end if 453 if not message: 454 return consistent 455 else: 456 return consistent,msg 457 #end if 458 #end def check_consistent 459 460 461 def is_valid(self): 462 return self.check_consistent(exit=False) 463 #end def is_valid 464 465 466 def change_units(self,units): 467 self.structure.change_units(units,folded=False) 468 if self.folded_system!=None: 469 self.folded_system.change_units(units) 470 #end if 471 #end def change_units 472 473 474 def group_atoms(self): 475 self.structure.group_atoms(folded=False) 476 if self.folded_system!=None: 477 self.folded_system.group_atoms() 478 #end if 479 #end def group_atoms 480 481 482 def rename(self,folded=True,**name_pairs): 483 self.particles.rename(**name_pairs) 484 self.structure.rename(folded=False,**name_pairs) 485 if self.pseudized: 486 for old,new in name_pairs.items(): 487 if old in self.valency: 488 if new not in self.valency: 489 self.valency[new] = self.valency[old] 490 #end if 491 del self.valency[old] 492 #end if 493 #end for 494 self.valency_in = self.valency 495 #end if 496 if self.folded_system!=None and folded: 497 self.folded_system.rename(folded=folded,**name_pairs) 498 #end if 499 #end def rename 500 501 502 def copy(self): 503 cp = DevBase.copy(self) 504 if self.folded_system!=None and self.structure.folded_structure!=None: 505 del cp.folded_system.structure 506 cp.folded_system.structure = cp.structure.folded_structure 507 #end if 508 return cp 509 #end def copy 510 511 512 def load(self,filepath): 513 DevBase.load(self,filepath) 514 if self.folded_system!=None and self.structure.folded_structure!=None: 515 del self.folded_system.structure 516 self.folded_system.structure = self.structure.folded_structure 517 #end if 518 #end def load 519 520 521 def tile(self,*td,**kwargs): 522 extensive = True 523 net_spin = None 524 if 'extensive' in kwargs: 525 extensive = kwargs['extensive'] 526 #end if 527 if 'net_spin' in kwargs: 528 net_spin = kwargs['net_spin'] 529 #end if 530 supercell = self.structure.tile(*td) 531 supercell.remove_folded() 532 if extensive: 533 ncells = int(round(supercell.volume()/self.structure.volume())) 534 net_charge = ncells*self.net_charge 535 if net_spin is None: 536 net_spin = ncells*self.net_spin 537 #end if 538 else: 539 net_charge = self.net_charge 540 if net_spin is None: 541 net_spin = self.net_spin 542 #end if 543 #end if 544 system = self.copy() 545 supersystem = PhysicalSystem( 546 structure = supercell, 547 net_charge = net_charge, 548 net_spin = net_spin, 549 **self.valency 550 ) 551 supersystem.folded_system = system 552 supersystem.structure.set_folded(system.structure) 553 return supersystem 554 #end def tile 555 556 557 def has_folded(self): 558 return self.folded_system is not None 559 #end def has_folded 560 561 562 def remove_folded_system(self): 563 self.folded_system = None 564 self.structure.remove_folded_structure() 565 #end def remove_folded_system 566 567 568 def remove_folded(self): 569 self.remove_folded_system() 570 #end def remove_folded 571 572 573 def get_smallest(self): 574 if self.has_folded(): 575 return self.folded_system 576 else: 577 return self 578 #end if 579 #end def get_smallest 580 581 582 def is_magnetic(self): 583 return self.net_spin!=0 or self.structure.is_magnetic() 584 #end def is_magnetic 585 586 587 def spin_polarized_orbitals(self): 588 return self.is_magnetic() 589 #end def spin_polarized_orbitals 590 591 592 # test needed 593 def large_Zeff_elem(self,Zmin): 594 elem = [] 595 for atom,Zeff in self.valency.items(): 596 if Zeff>Zmin: 597 elem.append(atom) 598 #end if 599 #end for 600 return elem 601 #end def large_Zeff_elem 602 603 604 # test needed 605 def ae_pp_species(self): 606 species = set(self.structure.elem) 607 if self.pseudized: 608 pp_species = set(self.valency.keys()) 609 ae_species = species-pp_species 610 else: 611 pp_species = set() 612 ae_species = species 613 #end if 614 return ae_species,pp_species 615 #end def ae_pp_species 616 617 618 def kf_rpa(self): 619 nelecs = self.particles.electron_counts() 620 volume = self.structure.volume() 621 kvol1 = (2*np.pi)**3/volume # k-space volume per particle 622 kfs = [(3*nelec*kvol1/(4*np.pi))**(1./3) for nelec in nelecs] 623 return np.array(kfs) 624 #end def kf_rpa 625#end class PhysicalSystem 626 627 628 629import os 630from structure import generate_structure,read_structure 631from copy import deepcopy 632ps_defaults = dict( 633 type='crystal', 634 kshift = (0,0,0), 635 net_charge=0, 636 net_spin=0, 637 pretile=None, 638 tiling=None, 639 tiled_spin=None, 640 extensive=True 641 ) 642def generate_physical_system(**kwargs): 643 for var,val in ps_defaults.items(): 644 if not var in kwargs: 645 kwargs[var] = val 646 #end if 647 #end for 648 type = kwargs['type'] 649 if type=='atom' or type=='dimer' or type=='trimer': 650 del kwargs['kshift'] 651 del kwargs['tiling'] 652 #if not 'units' in kwargs: 653 # kwargs['units'] = 'B' 654 ##end if 655 tiling = None 656 else: 657 tiling = kwargs['tiling'] 658 #end if 659 660 if 'structure' in kwargs: 661 s = kwargs['structure'] 662 is_str = isinstance(s,str) 663 if is_str: 664 if os.path.exists(s): 665 if 'elem' in kwargs: 666 s = read_structure(s,elem=kwargs['elem']) 667 else: 668 s = read_structure(s) 669 #end if 670 if 'axes' in kwargs: 671 s.reset_axes(kwargs['axes']) 672 #end if 673 kwargs['structure'] = s 674 else: 675 slow = s.lower() 676 format = None 677 if '.' in slow: 678 format = slow.rsplit('.')[1] 679 elif 'poscar' in slow: 680 format = 'poscar' 681 #end if 682 is_path = '/' in s 683 is_file = format in set('xyz xsf poscar cif fhi-aims'.split()) 684 if is_path or is_file: 685 PhysicalSystem.class_error('user provided structure file does not exist\nstructure file path: '+s,'generate_physical_system') 686 #end if 687 #end if 688 #end if 689 #end if 690 691 generation_info = obj() 692 generation_info.transfer_from(deepcopy(kwargs)) 693 694 net_charge = kwargs['net_charge'] 695 net_spin = kwargs['net_spin'] 696 tiled_spin = kwargs['tiled_spin'] 697 extensive = kwargs['extensive'] 698 del kwargs['net_spin'] 699 del kwargs['net_charge'] 700 del kwargs['tiled_spin'] 701 del kwargs['extensive'] 702 if 'particles' in kwargs: 703 particles = kwargs['particles'] 704 del kwargs['particles'] 705 else: 706 generation_info.particles = None 707 #end if 708 pretile = kwargs['pretile'] 709 del kwargs['pretile'] 710 valency = dict() 711 remove = [] 712 for var in kwargs: 713 #if var in Matter.elements: 714 if is_element(var): 715 valency[var] = kwargs[var] 716 remove.append(var) 717 #end if 718 #end if 719 generation_info.valency = deepcopy(valency) 720 for var in remove: 721 del kwargs[var] 722 #end for 723 724 if pretile is None: 725 structure = generate_structure(**kwargs) 726 else: 727 for d in range(len(pretile)): 728 if tiling[d]%pretile[d]!=0: 729 PhysicalSystem.class_error('pretile does not divide evenly into tiling\n tiling provided: {0}\n pretile provided: {1}'.format(tiling,pretile),'generate_physical_system') 730 #end if 731 #end for 732 tiling = tuple(array(tiling)//array(pretile)) 733 kwargs['tiling'] = pretile 734 pre = generate_structure(**kwargs) 735 pre.remove_folded_structure() 736 structure = pre.tile(tiling) 737 #end if 738 739 if tiling!=None and tiling!=(1,1,1) and structure.has_folded(): 740 fps = PhysicalSystem( 741 structure = structure.folded_structure, 742 net_charge = net_charge, 743 net_spin = net_spin, 744 **valency 745 ) 746 structure.remove_folded() 747 folded_structure = fps.structure 748 if extensive: 749 ncells = int(round(structure.volume()/folded_structure.volume())) 750 net_charge = ncells*net_charge 751 if not isinstance(net_spin,str): 752 net_spin = ncells*net_spin 753 #end if 754 #end if 755 if tiled_spin!=None: 756 net_spin = tiled_spin 757 #end if 758 ps = PhysicalSystem( 759 structure = structure, 760 net_charge = net_charge, 761 net_spin = net_spin, 762 **valency 763 ) 764 structure.set_folded(folded_structure) 765 ps.folded_system = fps 766 else: 767 ps = PhysicalSystem( 768 structure = structure, 769 net_charge = net_charge, 770 net_spin = net_spin, 771 **valency 772 ) 773 #end if 774 775 ps.generation_info = generation_info 776 777 return ps 778#end def generate_physical_system 779 780 781 782# test needed 783def ghost_atoms(*particles): 784 for particle in particles: 785 Matter.particle_collection.add_particles(Ion(name=particle,mass=0,charge=0,spin=0,protons=0,neutrons=0)) 786 #end for 787#end def ghost_atoms 788 789