Changeset 3434


Ignore:
Timestamp:
Jun 13, 2018 5:48:33 PM (3 years ago)
Author:
vondreele
Message:

working magnetic extinction check routine w. reference & comments - works

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r3423 r3434  
    191191            nfev += 1
    192192            chisq1 = np.sum(M2**2)
    193             if chisq1 > chisq0*(1.+chitol):     #TODO put Alan Coehlo's criteria for lambda here
     193            if chisq1 > chisq0*(1.+chitol):     #TODO put Alan Coehlo's criteria for lambda here?
    194194                lam *= 10.
    195195                if Print:
     
    223223    try:
    224224        Bmat,Nzero = pinv(Amatlam,xtol)    #Moore-Penrose inversion (via SVD) & count of zeros
    225 #        if Print:
    226225        if Print: print ('Found %d SVD zeros'%(Nzero))
    227 #        Bmat = nl.inv(Amatlam); Nzeros = 0
    228226        Bmat = Bmat/Anorm
    229227        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[],'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}]
  • trunk/GSASIIspc.py

    r3433 r3434  
    19101910    phi = f[:Nuniq]
    19111911    return iabsnt,mulp,Uniq,phi
    1912 
    1913 #def MagHKLchk(HKL,SGData):
    1914 #    SpnFlp = SGData['SpnFlp']
    1915 #    print(HKL)
    1916 #    Uniq = GenHKL(HKL,SGData)
    19171912   
    19181913def checkSSLaue(HKL,SGData,SSGData):
     
    19641959       
    19651960def checkHKLextc(HKL,SGData):
     1961    '''
     1962    Checks if reflection extinct
     1963
     1964    :param HKL:  [h,k,l]
     1965    :param SGData: space group data obtained from SpcGroup
     1966    :returns: False if  extinct; True if allowed
     1967
     1968    '''
    19661969    Ops = SGData['SGOps']
    19671970    OpM = np.array([op[0] for op in Ops])
     
    19791982
    19801983def checkMagextc(HKL,SGData):
     1984    '''
     1985    Checks if reflection magnetically extinct;
     1986    uses algorthm from Gallego, et al., J. Appl. Cryst. 45, 1236-1247 (2012)
     1987
     1988    :param HKL:  [h,k,l]
     1989    :param SGData: space group data obtained from SpcGroup; must have magnetic symmetry SpnFlp data
     1990    :returns: False if magnetically extinct; True if allowed
     1991
     1992    '''
    19811993    Ops = SGData['SGOps']
     1994    Ncen = len(SGData['SGCen'])
    19821995    OpM = np.array([op[0] for op in Ops])
    19831996    OpT = np.array([op[1] for op in Ops])
     
    19851998        OpM = np.vstack((OpM,-OpM))
    19861999        OpT = np.vstack((OpT,-OpT))%1.
     2000    OpM = np.reshape(np.array(list(OpM)*Ncen),(-1,3,3))
     2001    OpT = np.reshape(np.array([OpT+cen for cen in SGData['SGCen']]),(-1,3))
    19872002    Spn = SGData['SpnFlp'][:len(OpM)]
    19882003    Mag = np.array([nl.det(opm) for opm in OpM])*Spn
    19892004    DHKL = np.reshape(np.inner(HKL,OpM),(-1,3))
    1990     PHKL = np.reshape(np.cos(twopi*np.inner(HKL,OpT))*Mag,(-1,))[:,nxs,nxs]*OpM
    1991     Ftest = np.random.rand(3)
     2005    PHKL = np.reshape(np.cos(twopi*np.inner(HKL,OpT))*Mag,(-1,))[:,nxs,nxs]*OpM     #compute T(R,theta) eq(7)
     2006    Ftest = np.random.rand(3)       #random magnetic moment
    19922007    Psum = np.zeros(3)
     2008    nsum = 0.
     2009    nA = 0
    19932010    for dhkl,phkl in zip(DHKL,PHKL):
    1994         if not np.allclose(dhkl,HKL):
     2011        if not np.allclose(dhkl,HKL):           #test for eq(5)
    19952012            continue
    19962013        else:
    1997             Psum += np.inner(Ftest,phkl)
     2014            nA += 1
     2015            nsum += np.trace(phkl)          #eq(8)
     2016            pterm = np.inner(Ftest,phkl)    #eq(9)
     2017            Psum += pterm
     2018    if nsum/nA > 1.:        #only need to look at nA=1 frok eq(8)
     2019        return True
    19982020    if np.allclose(Psum,np.zeros(3)):
    19992021        return False
Note: See TracChangeset for help on using the changeset viewer.