diff --git a/nexus/bin/nxs-test b/nexus/bin/nxs-test index bc95972acb..806a6f18d2 100755 --- a/nexus/bin/nxs-test +++ b/nexus/bin/nxs-test @@ -2410,12 +2410,283 @@ def machines(): +def structure(): + # divert logging function + nlog_divert() + + nlabel('imports') + from generic import obj + from structure import Structure,Crystal + from structure import generate_structure + from structure import read_structure + + + nlabel('definitions') + def structure_same(s1,s2): + keys = ('units','elem','pos','axes','kpoints','kweights','kaxes') + o1 = s1.obj(keys) + o2 = s2.obj(keys) + return object_eq(o1,o2) + #end def structure_same + + def structure_diff(s1,s2): + keys = ('units','elem','pos','axes','kpoints','kweights','kaxes') + o1 = s1.obj(keys) + o2 = s2.obj(keys) + return object_diff(o1,o2,full=True) + #end def structure_diff + + + nlabel('empty_init') + s1 = Structure() + s2 = generate_structure('empty') + nassert(object_eq(s1,s2)) + + + nlabel('ref_inputs') + ref_in = obj() + ref_in.diamond_prim = obj( + units = 'A', + axes = [[1.785, 1.785, 0. ], + [0. , 1.785, 1.785], + [1.785, 0. , 1.785]], + elem = ['C','C'], + pos = [[0. , 0. , 0. ], + [0.8925, 0.8925, 0.8925]], + ) + ref_in.diamond_conv = obj( + units = 'A', + axes = [[3.57, 0 , 0. ], + [0. , 3.57, 0. ], + [0 , 0. , 3.57]], + elem = ['C','C','C','C','C','C','C','C'], + pos = [[0.0000, 0.0000, 0.0000], + [0.8925, 0.8925, 0.8925], + [0.0000, 1.7850, 1.7850], + [0.8925, 2.6775, 2.6775], + [1.7850, 0.0000, 1.7850], + [2.6775, 0.8925, 2.6775], + [1.7850, 1.7850, 0.0000], + [2.6775, 2.6775, 0.8925]], + ) + ref_in.wurtzite_prim = obj( + units = 'A', + axes = [[ 3.350, 0.00000000, 0.00], + [-1.675, 2.90118510, 0.00], + [ 0.000, 0.00000000, 5.22]], + elem = ['Zn','O','Zn','O'], + pos = [[0.00000000e+00, 0.00000000e+00, 3.26250000e+00], + [0.00000000e+00, 0.00000000e+00, 0.00000000e+00], + [1.67500000e+00, 9.67061701e-01, 6.52500000e-01], + [1.67500000e+00, 9.67061701e-01, 2.61000000e+00]], + ) + ref_in.oxygen_prim = obj( + units = 'A', + axes = [[ 2.70150000e+00, -1.71450000e+00, 0.00000000e+00], + [ 2.70150000e+00, 1.71450000e+00, 0.00000000e+00], + [-3.43801471e+00, 0.00000000e+00, 3.74799291e+00]], + elem = ['O','O'], + pos = [[ 0., 0., 0.575], + [ 0., 0., -0.575]], + ) + ref_in.CuO_prim = obj( + units = 'A', + axes = [[ 2.34150000e+00, -1.71100000e+00, 0.00000000e+00], + [ 2.34150000e+00, 1.71100000e+00, 0.00000000e+00], + [-8.49894838e-01, 0.00000000e+00, 5.05708046e+00]], + elem = ['Cu','O','Cu','O'], + pos = [[ 1.17075 , 0.8555 , 0. ], + [-0.21247371, 1.430396, 1.26427011], + [ 0.74580258, 2.5665 , 2.52854023], + [ 1.70407887, 0.280604, 3.79281034]], + ) + ref_in.Ca2CuO3_prim = obj( + units = 'A', + axes = [[ 1.885, 1.625, -6.115], + [-1.885, 1.625, 6.115], + [ 1.885, -1.625, 6.115]], + elem = ['Cu','O','O','O','Ca','Ca'], + pos = [[0.000, 0.000, 0.000], + [1.885, 0.000, 0.000], + [0.000, 0.000, 1.960], + [0.000, 0.000, 10.270], + [0.000, 0.000, 4.290], + [0.000, 0.000, 7.940]], + ) + ref_in.La2CuO4_prim = obj( + units = 'A', + axes = [[ 1.9045, 1.9045, -6.5845], + [-1.9045, 1.9045, 6.5845], + [ 1.9045, -1.9045, 6.5845]], + elem = ['Cu','O','O','O','O','La','La'], + pos = [[ 0. , 0. , 0. ], + [ 0.95225 , 0.95225 , -3.29225 ], + [-0.95225 , 0.95225 , 3.29225 ], + [ 0.346619, -0.346619, 1.198379], + [-0.346619, 0.346619, -1.198379], + [ 0.689429, -0.689429, 2.383589], + [-0.689429, 0.689429, -2.383589]], + ) + ref_in.graphene_prim = obj( + units = 'A', + axes = [[ 2.462, 0.00000000, 0.00], + [-1.231, 2.13215454, 0.00], + [ 0.000, 0.00000000, 15.00]], + elem = ['C','C'], + pos = [[0. , 0. , 0. ], + [1.231, 0.71071818, 0. ]], + ) + + + nlabel('direct_init') + ref = obj() + for name,inputs in ref_in.iteritems(): + ref[name] = Structure(**inputs) + #end for + + + nlabel('generate_init') + gen = obj() + for name,inputs in ref_in.iteritems(): + gen[name] = generate_structure(**inputs) + #end for + + + nlabel('crystal_init') + crys = obj() + for (latt,cell),inputs in Crystal.known_crystals.iteritems(): + s = generate_structure(structure=latt,cell=cell) + crys[latt+'_'+cell] = s + #end for + + + nlabel('check_generate_init') + for name,s in ref.iteritems(): + nassert(structure_same(s,gen[name])) + #end for + + + nlabel('check_crystal_init') + for name,s in ref.iteritems(): + nassert(structure_same(s,crys[name])) + #end for + + + nlabel('diagonal_tiling') + diag_tilings = [ + (1, 1, 1), + (1, 2, 3), + (1, 2, 4), + (1, 3, 1), + (1, 3, 2), + (1, 3, 3), + (1, 5, 3), + (1, 5, 5), + (1, 5, 6), + (2, 1, 2), + (2, 1, 3), + (2, 2, 2), + (2, 3, 1), + (2, 3, 2), + (2, 3, 3), + (2, 4, 4), + (2, 5, 1), + (2, 5, 2), + (2, 5, 3), + (3, 1, 1), + (3, 1, 2), + (3, 1, 3), + (3, 2, 1), + (3, 2, 3), + (3, 3, 6), + (4, 6, 4), + (5, 5, 1), + (5, 5, 5), + (6, 2, 6), + (6, 3, 1), + (6, 3, 6), + (6, 4, 6), + (6, 6, 4), + ] + for name,s in ref.iteritems(): + for tvec in diag_tilings: + st = s.tile(tvec) + st.check_tiling() + #end for + #end for + + + nlabel('matrix_tiling') + matrix_tilings = [ + (-3,-2, 0,-2,-2,-3, 3, 2,-1), + (-3,-1, 3, 0, 1, 1,-3, 3,-3), + (-3, 0,-2,-3, 0, 3, 2,-3, 0), + (-3, 1, 0,-3, 3,-3, 0,-2, 1), + (-2,-1, 2, 3,-3,-3,-2, 2,-2), + (-2, 0, 3,-2, 1,-1,-3, 0, 1), + (-2, 1,-3,-3, 0, 3, 2, 0, 3), + (-2, 2,-3,-1, 0, 1,-1,-2, 3), + (-1, 3, 2,-3, 2,-2, 3,-3,-1), + (-1, 3, 3,-3, 3, 0,-3,-1, 1), + ( 0,-1,-2, 2, 3, 1, 3,-2, 0), + ( 0, 1, 3,-2, 3, 1, 2, 1, 2), + ( 0, 3,-1,-1,-1, 2, 2, 3, 2), + ( 1,-3,-3,-3, 3, 0,-2, 2,-2), + ( 1,-2, 2, 1, 1, 1, 3, 2,-1), + ( 1, 0, 1,-3, 2, 1, 1, 2,-2), + ( 1, 2, 1,-2, 1,-1,-2, 3, 2), + ( 1, 2, 1, 2, 1,-1, 3, 2, 2), + ( 1, 2, 2,-2, 1, 0, 1, 0,-1), + ( 1, 2, 2, 1, 0, 1,-3, 0,-2), + ( 1, 3,-1, 0, 1,-2, 1, 0, 2), + ( 1, 3, 0,-1, 0,-3, 2, 0, 2), + ( 1, 3, 0, 0,-2,-2,-1,-1, 1), + ( 1, 3, 1,-3, 2, 0, 0, 0, 2), + ( 2,-3, 0, 2, 2,-3,-1,-2, 0), + ( 2,-2, 1, 3, 3,-2,-2, 1,-3), + ( 2,-1, 2,-1, 0,-3, 2, 0, 1), + ( 2,-1, 2, 0,-2, 0,-2, 0, 0), + ( 2,-1, 3, 2, 1, 2, 1,-1,-3), + ( 2, 1, 2, 0,-3,-2,-2,-3,-1), + ( 2, 2, 1,-3,-2,-3, 2,-1,-2), + ( 2, 3,-3, 1, 3, 1, 2,-2,-1), + ( 3,-3, 0, 1,-2, 3, 2, 1, 1), + ( 3,-2, 3, 3,-3,-3, 0, 0,-1), + ( 3,-1, 2,-3, 2,-2, 3, 2,-1), + ( 3, 0,-3, 0, 2,-1, 3,-2,-1), + ( 3, 0, 2, 3,-2,-3,-3,-2, 1), + ( 3, 1, 1,-2,-1,-3, 3, 1, 0), + ( 3, 2,-3,-2, 2, 1, 1, 0,-1), + ( 3, 2,-2,-1,-3,-1, 3, 3, 2), + ( 3, 2, 3,-1, 2, 2, 0,-1, 0), + ( 3, 3,-1, 0,-3, 2,-3,-1, 0), + ( 3, 3, 1,-3,-2,-3, 2,-3, 3), + ( 3, 3, 1, 1,-2, 2, 2,-2, 0), + ] + #for name in sorted(ref.keys()): + for name in ['diamond_prim']: + s = ref[name] + npass = 0 + for tmat in matrix_tilings: + tmat = array(tmat,dtype=int) + tmat.shape = 3,3 + st = s.tile(tmat) + st.check_tiling() + #end for + #end for + + # restore logging function + nlog_restore() +#end def structure + + + def pwscf_input(): # divert logging function nlog_divert() nlabel('directories') - loc_test_dir = os.path.join(test_dir,'pwscf_input') + loc_test_dir = os.path.join(test_dir,'pwscf_input') sample_inputs = os.path.join(loc_test_dir,'sample_inputs') sample_structures = os.path.join(loc_test_dir,'sample_structures') @@ -2974,6 +3245,7 @@ NexusTest( generic_extensions ) NexusTest( nexus_imports ) NexusTest( settings_operation , 'settings' ) NexusTest( machines ) +NexusTest( structure ) NexusTest( pwscf_input ) for label in example_information.keys(): diff --git a/nexus/lib/structure.py b/nexus/lib/structure.py index a374ca7289..b26249dfc2 100755 --- a/nexus/lib/structure.py +++ b/nexus/lib/structure.py @@ -266,6 +266,9 @@ def reduce_tilematrix(tiling): dim = len(t) matrix_tiling = t.shape == (dim,dim) if matrix_tiling: + if abs(det(t))==0: + Structure.class_error('requested tiling matrix is singular\ntiling requested: {0}'.format(t)) + #end if #find a tiling tuple from the tiling matrix # do this by shearing the tiling matrix (or equivalently the tiled cell) # until it is orthogonal (in the untiled cell axes) @@ -276,7 +279,6 @@ def reduce_tilematrix(tiling): T = t #tiling matrix tilematrix = T.copy() del t - Tnew = array(T,dtype=float) #sheared/orthogonal tiling matrix tbar = identity(dim) #basis for shearing dr = range(dim) #dr = [1,0,2] @@ -287,8 +289,10 @@ def reduce_tilematrix(tiling): #move each axis to be parallel to barred directions # these are volume preserving shears of the supercell # each shear keeps two cell face planes fixed while moving the others + tvecs = [] for dp in [(0,1,2),(2,0,1),(1,2,0),(2,1,0),(0,2,1),(1,0,2)]: success = True + Tnew = array(T,dtype=float) #sheared/orthogonal tiling matrix for d in dr: tb = tbar[dp[d]] t = T[d] @@ -304,29 +308,32 @@ def reduce_tilematrix(tiling): Tnew[d] = tn #end for if success: - break + # apply inverse permutation, if needed + Tn = Tnew.copy() + for d in dr: + d2 = dp[d] + Tnew[d2] = Tn[d] + #end for + #the resulting tiling matrix should be diagonal and integer + tr = diag(Tnew) + nondiagonal = abs(Tnew-diag(tr)).sum()>1e-6 + if nondiagonal: + Structure.class_error('could not find a diagonal tiling matrix for generating tiled coordinates') + #end if + tvecs.append(abs(tr)) #end if #end for - # apply inverse permutation, if needed - Tn = Tnew.copy() - for d in dr: - d2 = dp[d] - Tnew[d2] = Tn[d] + tvecs_old = tvecs + tvecs = [] + tvset = set() + for tv in tvecs_old: + tvk = tuple(array(around(1e7*tv),dtype=uint64)) + if tvk not in tvset: + tvset.add(tvk) + tvecs.append(tv) + #end if #end for - #the resulting tiling matrix should be diagonal and integer - tr = diag(Tnew) - t = array(around(tr),dtype=int) - nondiagonal = abs(Tnew-diag(tr)).sum()>1e-6 - noninteger = abs(tr-t).sum()>1e-6 - if nondiagonal: - Structure.class_error('could not find a diagonal tiling matrix for generating tiled coordinates') - #end if - #non-integer tile vectors are handled directly by tile_points now - #if noninteger: - # Structure.class_error('calculated diagonal tiling matrix is non-integer\n tiled coordinates cannot be determined') - ##end if - #tilevector = abs(t) - tilevector = abs(tr) + tilevector = array(tvecs) else: tilevector = t tilematrix = diag(t) @@ -336,6 +343,7 @@ def reduce_tilematrix(tiling): #end def reduce_tilematrix + def tile_magnetization(mag,tilevec,mag_order,mag_prim): # jtk mark current # implement true magnetic tiling based on the magnetic order @@ -348,7 +356,7 @@ def tile_magnetization(mag,tilevec,mag_order,mag_prim): # if magnetic primitive is requested, a small tiling vector should first be used # (ie 221,212,122,222 periods all involve a 211 prim tiling w/ a simple reshaping/reassignment of the cell axes # Structure should have a function providing labels to each magnetic species - mag = array(int(round( tilevec.prod() ))*list(mag),dtype=object) + #mag = array(int(round( tilevec.prod() ))*list(mag),dtype=object) return mag #end def tile_magnetization @@ -2950,6 +2958,7 @@ def tile(self,*td,**kwargs): in_place = kwargs.pop('in_place',False) magnetic_order = kwargs.pop('magnetic_order',None) magnetic_primitive = kwargs.pop('magnetic_primitive',True) + check = kwargs.pop('check',False) dim = self.dim if len(td)==1: @@ -3016,6 +3025,10 @@ def tile(self,*td,**kwargs): ts = self #end if + if check: + ts.check_tiling() + #end if + return ts #end def tile @@ -3024,6 +3037,48 @@ def tile_points(self,points,axes,tilemat,tilevec=None): if tilevec is None: tilemat,tilevec = reduce_tilematrix(tilemat) #end if + if not isinstance(tilemat,ndarray): + tilemat = array(tilemat) + #end if + matrix_tiling = abs(tilemat-diag(diag(tilemat))).sum()>0.1 + if not matrix_tiling: + return self.tile_points_simple(points,axes,diag(tilemat)) + else: + if not isinstance(axes,ndarray): + axes = array(axes) + #end if + if not isinstance(tilevec,ndarray): + tilevec = array(tilevec) + #end if + dim = len(axes) + npoints = len(points) + ntpoints = npoints*int(round(abs(det(tilemat)))) + if tilevec.size==dim: + tilevec.shape = 1,dim + #end if + taxes = dot(tilemat,axes) + success = False + for tvec in tilevec: + tpoints = self.tile_points_simple(points,axes,tvec) + tpoints,weights,pmap = self.unique_points_fast(tpoints,taxes) + if len(tpoints)==ntpoints: + success = True + break + #end if + #end for + if not success: + tpoints = self.tile_points_brute(points,axes,tilemat) + tpoints,weights,pmap = self.unique_points_fast(tpoints,taxes) + if len(tpoints)!=ntpoints: + self.error('brute force tiling failed') + #end if + #end if + #end if + return tpoints + #end def tile_points + + + def tile_points_simple(self,points,axes,tilevec): if not isinstance(points,ndarray): points = array(points) #end if @@ -3033,57 +3088,64 @@ def tile_points(self,points,axes,tilemat,tilevec=None): if not isinstance(axes,ndarray): axes = array(axes) #end if - t = tilevec - ti = array(around(t),dtype=int) - noninteger = abs(t-ti).sum()>1e-6 if len(points.shape)==1: npoints,dim = len(points),1 else: npoints,dim = points.shape #end if - ntpoints = npoints*int(round( t.prod() )) - if not noninteger: + t = tilevec + ti = array(around(t),dtype=int) + noninteger = abs(t-ti).sum()>1e-6 + if noninteger: + tp = t.prod() + if abs(tp-int(tp))>1e-6: + self.error('tiling vector does not correspond to an integer volume change\ntiling vector: {0}\nvolume change: {1} {2} {3}'.format(tilevec,tilevec.prod(),ntpoints,int(ntpoints))) + #end if + t = array(ceil(t),dtype=int)+1 + else: t = ti - if ntpoints==0: - tpoints = array([]) - else: - tpoints = empty((ntpoints,dim)) - ns=0 - ne=npoints - for k in range(t[2]): - for j in range(t[1]): - for i in range(t[0]): - v = dot(array([[i,j,k]]),axes) - for d in range(dim): - tpoints[ns:ne,d] = points[:,d]+v[0,d] - #end for - ns+=npoints - ne+=npoints + #end if + ntpoints = npoints*int(round( t.prod() )) + if ntpoints==0: + tpoints = array([]) + else: + tpoints = empty((ntpoints,dim)) + ns=0 + ne=npoints + for k in range(t[2]): + for j in range(t[1]): + for i in range(t[0]): + v = dot(array([[i,j,k]]),axes) + for d in range(dim): + tpoints[ns:ne,d] = points[:,d]+v[0,d] #end for + ns+=npoints + ne+=npoints #end for #end for - #end if - else: - if abs(ntpoints-int(around(ntpoints)))>1e-6: - self.error('tiling vector does not correspond to an integer volume change\ntiling vector: {0}\nvolume change: {1} {2} {3}'.format(tilevec,tilevec.prod(),ntpoints,int(ntpoints))) - #end if - ntpoints = int(around(ntpoints)) - # round up to larger tiling - # +1 added for greater cell coverage - # add more if error below is tripped w/ fewer points than expected - t = array(ceil(t),dtype=int)+1 - # get the tiled points - tpoints = self.tile_points(points,axes,tilemat,t) - # remove any that are not unique - taxes = dot(tilemat,axes) - #tpoints,weights,pmap = self.unique_points(tpoints,taxes) - tpoints,weights,pmap = self.unique_points_fast(tpoints,taxes) - if len(tpoints)!=ntpoints: - self.error('tiling by non-integer tiling vector failed\npoints expected after tiling: {0}\npoints resulted from tiling: {1}'.format(ntpoints,len(tpoints))) - #end if + #end for #end if return tpoints - #end def tile_points + #end def tile_points_simple + + + def tile_points_brute(self,points,axes,tilemat): + tcorners = [[0,0,0], + [1,0,0], + [0,1,0], + [0,0,1], + [0,1,1], + [1,0,1], + [1,1,0], + [1,1,1]] + tcorners = dot(tcorners,tilemat) + tmin = tcorners.min(axis=0) + tmax = tcorners.max(axis=0) + tilevec = tmax-tmin + tpoints = self.tile_points_simple(points,axes,tilevec) + return tpoints + #end def tile_points_brute + def opt_tilematrix(self,*args,**kwargs): @@ -3097,6 +3159,43 @@ def tile_opt(self,*args,**kwargs): #end def tile_opt + def check_tiling(self,tol=1e-6,exit=True): + msg = '' + if not self.is_tiled(): + return msg + #end if + msgs = [] + st = self + s = self.folded_structure + nt = len(st.pos) + n = len(s.pos) + if nt%n!=0: + msgs.append('tiled atom count does is not divisible by untiled atom count') + #end if + vratio = st.volume()/s.volume() + if abs(vratio-float(nt)/n)>tol: + msgs.append('tiled/untiled volume ratio does not match tiled/untiled atom count ratio') + #end if + if abs(vratio-abs(det(st.tmatrix)))>tol: + msgs.append('tiled/untiled volume ratio does not match tiling matrix determinant') + #end if + p,w,pmap = self.unique_points_fast(st.pos,st.axes) + if len(p)!=nt: + msgs.append('tiled positions are not unique') + #end if + if len(msgs)>0: + msg = 'tiling check failed' + for m in msgs: + msg += '\n'+m + #end for + if exit: + self.error(msg) + #end if + #end if + return msg + #end def check_tiling + + def kfold(self,tiling,kpoints,kweights): if isinstance(tiling,int): tiling = self.dim*[tiling]