scnc
/
qmcpack
mirror of https://github.com/QMCPACK/qmcpack.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
qmcpack/nexus/tests/unit/test_physical_system.py

521 lines
15 KiB
Python
Raw Permalink Blame History

#! /usr/bin/env python3
import numpy as np
import testing
from testing import value_eq,object_eq,object_diff,print_diff
from test_structure import structure_same
def system_same(s1,s2,pseudized=True,tiled=False):
same = True
keys = ('net_charge','net_spin','pseudized','particles')
o1 = s1.obj(keys)
o2 = s2.obj(keys)
qsame = object_eq(o1,o2)
vsame = True
if pseudized:
vsame = s1.valency==s2.valency
#end if
ssame = structure_same(s1.structure,s2.structure)
fsame = True
if tiled:
fsame = system_same(s1.folded_system,s2.folded_system)
#end if
same = qsame and vsame and ssame and fsame
return same
#end def system_same
def test_import():
import physical_system
from physical_system import Matter,Particle,Ion,PseudoIon,Particles
from physical_system import PhysicalSystem,generate_physical_system
#end def test_import
def test_particle_initialization():
from generic import obj
from physical_system import Matter,Particle,Ion,PseudoIon,Particles
from physical_system import PhysicalSystem,generate_physical_system
# 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)==103)
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')
assert(value_eq(si.mass,51197.6459833))
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():
import os
from generic import obj
from structure import generate_structure
from physical_system import generate_physical_system
from physical_system import PhysicalSystem
tpath = testing.setup_unit_test_output_directory('physical_system','test_physical_system_initialization')
d2 = generate_structure(
structure = 'diamond',
cell = 'prim',
)
d2_path = os.path.join(tpath,'diamond2.xsf')
d2.write(d2_path)
d8 = generate_structure(
structure = 'diamond',
cell = 'conv',
)
d8_path = os.path.join(tpath,'diamond8.xsf')
d8.write(d8_path)
d8_tile = d2.tile([[ 1, -1, 1],
[ 1, 1, -1],
[-1, 1, 1]])
d8_tile_pos_ref = np.array([
[0. , 0. , 0. ],
[0.25, 0.25, 0.25],
[0.5 , 0.5 , 0. ],
[0.75, 0.75, 0.25],
[0. , 0.5 , 0.5 ],
[0.25, 0.75, 0.75],
[0.5 , 0. , 0.5 ],
[0.75, 0.25, 0.75]])
assert(value_eq(d8_tile.pos_unit(),d8_tile_pos_ref,atol=1e-8))
direct_notile = generate_physical_system(
units = 'A',
axes = [[3.57, 0.00, 0.00],
[0.00, 3.57, 0.00],
[0.00, 0.00, 3.57]],
elem = 8*['C'],
posu = [[0.00, 0.00, 0.00],
[0.25, 0.25, 0.25],
[0.00, 0.50, 0.50],
[0.25, 0.75, 0.75],
[0.50, 0.00, 0.50],
[0.75, 0.25, 0.75],
[0.50, 0.50, 0.00],
[0.75, 0.75, 0.25]],
C = 4,
)
direct_tile = 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,
)
struct_notile = generate_physical_system(
structure = d8,
C = 4,
)
struct_tile = generate_physical_system(
structure = d2,
tiling = [[ 1, -1, 1],
[ 1, 1, -1],
[-1, 1, 1]],
C = 4,
)
read_notile = generate_physical_system(
structure = d8_path,
C = 4,
)
read_tile = generate_physical_system(
structure = d2_path,
tiling = [[ 1, -1, 1],
[ 1, 1, -1],
[-1, 1, 1]],
C = 4,
)
gen_notile = generate_physical_system(
lattice = 'cubic', # cubic tetragonal orthorhombic rhombohedral
# hexagonal triclinic monoclinic
cell = 'conventional', # primitive, conventional
centering = 'F', # P A B C I R F
constants = 3.57, # a,b,c,alpha,beta,gamma
units = 'A', # A or B
atoms = 'C', # species in primitive cell
basis = [[0,0,0], # basis vectors (optional)
[.25,.25,.25]],
C = 4,
)
gen_tile = generate_physical_system(
lattice = 'cubic', # cubic tetragonal orthorhombic rhombohedral
# hexagonal triclinic monoclinic
cell = 'primitive', # primitive, conventional
centering = 'F', # P A B C I R F
constants = 3.57, # a,b,c,alpha,beta,gamma
units = 'A', # A or B
atoms = 'C', # species in primitive cell
basis = [[0,0,0], # basis vectors (optional)
[.25,.25,.25]],
tiling = [[ 1, -1, 1],
[ 1, 1, -1],
[-1, 1, 1]],
C = 4,
)
lookup_notile = generate_physical_system(
structure = 'diamond',
cell = 'conv',
C = 4,
)
lookup_tile = generate_physical_system(
structure = 'diamond',
cell = 'prim',
tiling = [[ 1, -1, 1],
[ 1, 1, -1],
[-1, 1, 1]],
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
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,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
sref = ref.structure
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,d8_tile))
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,0.8925,2.6775)))
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'))
assert(list(sref.elem)==2*['C'])
assert(value_eq(tuple(sref.pos[-1]),(0.8925,0.8925,0.8925)))
assert(sref.units=='A')
ref_notile = direct_notile
ref_tile = direct_tile
assert(system_same(struct_notile,ref_notile))
assert(system_same(read_notile ,ref_notile))
assert(system_same(gen_notile ,ref_notile))
assert(system_same(lookup_notile,ref_notile))
assert(system_same(struct_tile,ref_tile,tiled=True))
assert(system_same(read_tile ,ref_tile,tiled=True))
assert(system_same(gen_tile ,ref_tile,tiled=True))
assert(system_same(lookup_tile,ref_tile,tiled=True))
systems_notile = [
direct_notile,
struct_notile,
read_notile ,
gen_notile ,
lookup_notile,
]
systems_tile = [
direct_tile,
struct_tile,
read_tile ,
gen_tile ,
lookup_tile,
]
systems = systems_notile+systems_tile
for sys in systems:
assert(sys.is_valid())
#end for
# test has_folded
for sys in systems_notile:
assert(not sys.has_folded())
#end for
for sys in systems_tile:
assert(sys.has_folded())
#end for
# test copy
for sys in systems:
c = sys.copy()
assert(id(c)!=id(sys))
assert(c.is_valid())
assert(system_same(c,sys,tiled=sys.has_folded()))
#end for
# test load
for i,sys in enumerate(systems):
path = os.path.join(tpath,'system_{}'.format(i))
sys.save(path)
sys2 = PhysicalSystem()
sys2.load(path)
assert(sys2.is_valid())
assert(system_same(sys2,sys,tiled=sys.has_folded()))
#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])
#end def test_physical_system_initialization
def test_change_units():
from physical_system import generate_physical_system
sys = generate_physical_system(
units = 'A',
axes = [[3.57, 0.00, 0.00],
[0.00, 3.57, 0.00],
[0.00, 0.00, 3.57]],
elem = 8*['C'],
posu = [[0.00, 0.00, 0.00],
[0.25, 0.25, 0.25],
[0.00, 0.50, 0.50],
[0.25, 0.75, 0.75],
[0.50, 0.00, 0.50],
[0.75, 0.25, 0.75],
[0.50, 0.50, 0.00],
[0.75, 0.75, 0.25]],
C = 4,
)
s = sys.structure
assert(value_eq(s.pos[-1],np.array([2.6775,2.6775,0.8925])))
sys.change_units('B')
assert(value_eq(s.pos[-1],np.array([5.05974172,5.05974172,1.68658057])))
#end def test_change_units
def test_rename():
from generic import obj
from physical_system import generate_physical_system
sys = generate_physical_system(
units = 'A',
axes = [[1.785, 1.785, 0. ],
[0. , 1.785, 1.785],
[1.785, 0. , 1.785]],
elem = ['C1','C2'],
posu = [[0.00, 0.00, 0.00],
[0.25, 0.25, 0.25]],
tiling = [[ 1, -1, 1],
[ 1, 1, -1],
[-1, 1, 1]],
C1 = 4,
C2 = 4,
)
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))
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))
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))
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))
#end def test_rename
def test_tile():
from physical_system import generate_physical_system
d2_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]],
C = 4,
)
d8_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,
)
d8 = d2_ref.tile([[ 1, -1, 1],
[ 1, 1, -1],
[-1, 1, 1]])
assert(system_same(d8,d8_ref,tiled=True))
#end 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,
net_charge = 1,
net_spin = 1,
H = 1
)
kfs = ps.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
Reference in New Issue View Git Blame Copy Permalink