Changeset 646


Ignore:
Timestamp:
Jun 4, 2012 2:43:15 PM (9 years ago)
Author:
vondreele
Message:

eliminate Friedel argument in genHKLf - didn't do anything
work on molecule rotation in structure drawing stuff
a working version of findOffset, but slow

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r645 r646  
    584584   
    585585def findOffset(SGData,rho,Fhkl):
     586   
     587    def calcPhase(DX,DH,Dphi):
     588        H,K,L = DH.T
     589        Phi = (H*DX[0]+K*DX[1]+L*DX[2]+0.5) % 1.
     590        return Dphi-Phi
     591   
     592    if SGData['SpGrp'] == 'P 1':
     593        return [0,0,0]   
     594# will need to consider 'SGPolax': one of '','x','y','x y','z','x z','y z','xyz','111'
    586595    mapShape = rho.shape
     596    steps = np.array(mapShape)
    587597    hklShape = Fhkl.shape
    588598    mapHalf = np.array(mapShape)/2
     
    600610    F = str(1.e6)
    601611    i = 0
     612    DH = []
     613    Dphi = []
    602614    while float(F) > 0.5*Fmax:
    603615        F = Flist[i]
    604616        hkl = np.unravel_index(Fdict[F],hklShape)
    605         iabsnt,mulp,Uniq,phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData,Friedel=True)
    606         Uniq = np.array(Uniq,dtype='i')+hklHalf
    607         print hkl-hklHalf
     617        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
     618        Uniq = np.array(Uniq,dtype='i')
     619        Phi = np.array(Phi)
     620        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
     621        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
    608622        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
    609         ph0 = np.angle(Fh0,deg=True)/360.
    610         for j,H in enumerate(Uniq):
    611             Fh = Fhkl[H[0],H[1],H[2]]
    612             h,k,l = H-hklHalf
    613             ang = (np.angle(Fh,deg=True)/360.-phi[j]) % 1.
    614             print '(%3d,%3d,%3d) %9.5f'%(h,k,l,ang)       
     623        ang0 = np.angle(Fh0,deg=True)/360.
     624        for j,H in enumerate(Uniq[1:]):
     625            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-Phi[j+1])
     626            dH = H-hkl
     627            dang = ang-ang0
     628            DH.append(dH)
     629            Dphi.append((dang+0.5) % 1.0)
    615630        i += 1
    616        
    617        
    618    
    619     return [0,0,0]
     631#    for item in zip(DH,Dphi):
     632#        print item
     633    DH = np.array(DH)
     634    Dphi = np.array(Dphi)
     635    DX = np.zeros(3)
     636    X,Y,Z = np.mgrid[0:1:steps[0]*(0+1j),0:1:steps[1]*(0+1j),0:1:steps[2]*(0+1j)]
     637    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
     638    Mmin = 1.e10
     639   
     640    for xyz in XYZ:
     641        M = np.sum(calcPhase(xyz,DH,Dphi)**2)
     642        if M < Mmin:
     643            DX = xyz
     644            Mmin = M
     645   
     646    print DX,Mmin
     647    result = so.leastsq(calcPhase,DX,full_output=True,args=(DH,Dphi))
     648    chisq = np.sum(result[2]['fvec']**2)
     649    DX = np.array(np.rint(-result[0]*steps),dtype='i')
     650    print chisq,DX
     651    print result[0],result[0]*steps
     652    return DX
    620653
    621654def ChargeFlip(data,reflData,pgbar):
  • trunk/GSASIIphsGUI.py

    r636 r646  
    34973497            H = ref[:3]
    34983498            ref[4] = np.sqrt(1./G2lat.calc_rDsq2(H,G))
    3499             iabsnt,mulp,Uniq,phi = G2spc.GenHKLf(H,SGData,Friedel=False)
     3499            iabsnt,mulp,Uniq,phi = G2spc.GenHKLf(H,SGData)
    35003500            ref[3] = mulp/2             #convert from powder mulp.
    35013501            ref[11] = Uniq
  • trunk/GSASIIplot.py

    r644 r646  
    28632863        glMatrixMode(GL_MODELVIEW)
    28642864        glLoadIdentity()
    2865         glRotate(anglez,0,0,1)
    2866         glRotate(angley,sind(anglez),cosd(anglez),0)
    2867         glRotate(anglex,cosd(anglez),-sind(anglez),0)
     2865        matX = G2lat.rotdMat(anglex,axis=0)
     2866        matY = G2lat.rotdMat(angley,axis=1)
     2867        matZ = G2lat.rotdMat(anglez,axis=2)
     2868        matRot = np.inner(matX,np.inner(matY,matZ))
     2869        matRot = np.concatenate((np.concatenate((matRot,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
     2870        glMultMatrixf(matRot.T)
    28682871        glMultMatrixf(A4mat.T)
    28692872        glTranslate(-Tx,-Ty,-Tz)
  • trunk/GSASIIspc.py

    r584 r646  
    3535             'SGOps': symmetry operations as [M,T] so that M*x+T = x'
    3636             'SGSys': one of 'triclinic','monoclinic','orthorhombic','tetragonal','rhombohedral','trigonal','hexagonal','cubic'
     37             'SGPolax': one of '','x','y','x y','z','x z','y z','xyz','111' for arbitrary axes
    3738       '''
    3839    LaueSym = ('-1','2/m','mmm','4/m','4/mmm','3R','3mR','3','3m1','31m','6/m','6/mmm','m3','m3m')
     
    312313        return zip(XYZEquiv,Idup,Cell)
    313314
    314 def GenHKLf(HKL,SGData,Friedel=False):
     315def GenHKLf(HKL,SGData):
    315316    '''
    316317    Uses old GSAS Fortran routine genhkl.for
     
    318319        HKL - [h,k,l]
    319320        SGData - space group data obtained from SpcGroup
    320         Friedel = True to retain Friedel pairs for centrosymmetric case
    321321    returns:
    322322        iabsnt = True is reflection is forbidden by symmetry
    323         mulp = reflection multiplicity including Fridel pairs
     323        mulp = reflection multiplicity including Friedel pairs
    324324        Uniq = numpy array of equivalent hkl in descending order of h,k,l
    325325    '''
     
    335335    Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq]))
    336336    phi = f[:Nuniq]
    337     Uniq = np.array(Uniq)
    338    
    339     return iabsnt,2*mulp,Uniq,phi       #include Friedel pairs in powder mulp
     337   
     338    return iabsnt,mulp,Uniq,phi
    340339                                 
    341340def GetOprPtrName(key):           
  • trunk/GSASIIstruct.py

    r637 r646  
    12111211                for h,k,l,d in HKLd:
    12121212                    ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
     1213                    mul *= 2      # for powder overlap of Friedel pairs
    12131214                    if ext:
    12141215                        continue
     
    19831984        biso = -SQfactor*Uisodata
    19841985        Tiso = np.where(biso<1.,np.exp(biso),1.0)
    1985 #        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
    19861986        HbH = -np.inner(H,np.inner(bij,H))
    19871987        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
Note: See TracChangeset for help on using the changeset viewer.