Changeset 763 for trunk/GSASIImath.py


Ignore:
Timestamp:
Sep 26, 2012 1:34:54 AM (10 years ago)
Author:
vondreele
Message:

change peak search algorithm - now much faster & more complete

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r762 r763  
     1<<<<<<< .mine
    12# -*- coding: utf-8 -*-
    23#GSASIImath - major mathematics routines
     
    2425import GSASIIlattice as G2lat
    2526import GSASIIspc as G2spc
    26 import scipy.optimize as so
    27 import scipy.linalg as sl
    2827import numpy.fft as fft
    2928
     
    807806    return mapData
    808807   
     808def SearchMap(data):
     809    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
     810   
     811    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
     812   
     813    def noDuplicate(xyz,peaks,Amat):
     814        XYZ = np.inner(Amat,xyz)
     815        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
     816            print ' Peak',xyz,' <0.5A from another peak'
     817            return False
     818        return True
     819                           
     820    def fixSpecialPos(xyz,SGData,Amat):
     821        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
     822        X = []
     823        xyzs = [equiv[0] for equiv in equivs]
     824        for x in xyzs:
     825            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
     826                X.append(x)
     827        if len(X) > 1:
     828            return np.average(X,axis=0)
     829        else:
     830            return xyz
     831       
     832    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
     833        Mag,x0,y0,z0,sig = parms
     834        return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
     835       
     836    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
     837        Mag,x0,y0,z0,sig = parms
     838        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     839        return M
     840       
     841    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
     842        Mag,x0,y0,z0,sig = parms
     843        dMdv = np.zeros(([5,]+list(rX.shape)))
     844        delt = .01
     845        for i in range(5):
     846            parms[i] -= delt
     847            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     848            parms[i] += 2.*delt
     849            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     850            parms[i] -= delt
     851            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
     852        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
     853        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
     854        dMdv = np.reshape(dMdv,(5,rX.size))
     855        Hess = np.inner(dMdv,dMdv)
     856       
     857        return Vec,Hess
     858       
     859    generalData = data['General']
     860    phaseName = generalData['Name']
     861    SGData = generalData['SGData']
     862    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     863    drawingData = data['Drawing']
     864    peaks = []
     865    mags = []
     866    dzeros = []
     867    try:
     868        mapData = generalData['Map']
     869        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
     870        rho = copy.copy(mapData['rho'])     #don't mess up original
     871        mapHalf = np.array(rho.shape)/2
     872        res = mapData['Resolution']
     873        incre = np.array(rho.shape,dtype=np.float)
     874        step = max(1.0,1./res)+1
     875        steps = np.array(3*[step,])
     876    except KeyError:
     877        print '**** ERROR - Fourier map not defined'
     878        return peaks,mags
     879    rhoMask = ma.array(rho,mask=(rho<contLevel))
     880    indices = (-1,0,1)
     881    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
     882    for roll in rolls:
     883        if np.any(roll):
     884            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
     885    indx = np.transpose(rhoMask.nonzero())
     886    peaks = indx/incre
     887    mags = rhoMask[rhoMask.nonzero()]
     888    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
     889        rho = rollMap(rho,ind)
     890        rMM = mapHalf-steps
     891        rMP = mapHalf+steps+1
     892        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
     893        peakInt = np.sum(rhoPeak)*res**3
     894        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
     895        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
     896        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
     897            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
     898        x1 = result[0]
     899        if not np.any(x1 < 0):
     900            mag = x1[0]
     901            peak = (np.array(x1[1:4])-ind)/incre
     902        peak = fixSpecialPos(peak,SGData,Amat)
     903        rho = rollMap(rho,-ind)       
     904    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
     905    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
     906   
     907def sortArray(data,pos,reverse=False):
     908    #data is a list of items
     909    #sort by pos in list; reverse if True
     910    T = []
     911    for i,M in enumerate(data):
     912        T.append((M[pos],i))
     913    D = dict(zip(T,data))
     914    T.sort()
     915    if reverse:
     916        T.reverse()
     917    X = []
     918    for key in T:
     919        X.append(D[key])
     920    return X
     921               
     922def PeaksUnique(data,Ind):
     923
     924    def noDuplicate(xyz,peaks,Amat):
     925        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
     926            return False
     927        return True
     928                           
     929    generalData = data['General']
     930    cell = generalData['Cell'][1:7]
     931    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     932    A = G2lat.cell2A(cell)
     933    SGData = generalData['SGData']
     934    mapPeaks = data['Map Peaks']
     935    Indx = {}
     936    XYZ = {}
     937    for ind in Ind:
     938        XYZ[ind] = np.array(mapPeaks[ind][1:4])
     939        Indx[ind] = True
     940    for ind in Ind:
     941        if Indx[ind]:
     942            xyz = XYZ[ind]
     943            for jnd in Ind:
     944                if ind != jnd and Indx[jnd]:                       
     945                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
     946                    xyzs = np.array([equiv[0] for equiv in Equiv])
     947                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
     948    Ind = []
     949    for ind in Indx:
     950        if Indx[ind]:
     951            Ind.append(ind)
     952    return Ind
     953   
     954def prodQQ(QA,QB):
     955    ''' Grassman quaternion product
     956        QA,QB quaternions; q=r+ai+bj+ck
     957    '''
     958    D = np.zeros(4)
     959    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
     960    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
     961    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
     962    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
     963    return D
     964   
     965def normQ(QA):
     966    ''' get length of quaternion & normalize it
     967        q=r+ai+bj+ck
     968    '''
     969    n = np.sqrt(np.sum(np.array(QA)**2))
     970    return QA/n
     971   
     972def invQ(Q):
     973    '''
     974        get inverse of quaternion
     975        q=r+ai+bj+ck; q* = r-ai-bj-ck
     976    '''
     977    return Q*np.array([1,-1,-1,-1])
     978   
     979def prodQVQ(Q,V):
     980    ''' compute the quaternion vector rotation qvq-1 = v'
     981        q=r+ai+bj+ck
     982    '''
     983    VP = np.zeros(3)
     984    T2 = Q[0]*Q[1]
     985    T3 = Q[0]*Q[2]
     986    T4 = Q[0]*Q[3]
     987    T5 = -Q[1]*Q[1]
     988    T6 = Q[1]*Q[2]
     989    T7 = Q[1]*Q[3]
     990    T8 = -Q[2]*Q[2]
     991    T9 = Q[2]*Q[3]
     992    T10 = -Q[3]*Q[3]
     993    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
     994    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
     995    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2]
     996    return VP   
     997   
     998def Q2Mat(Q):
     999    ''' make rotation matrix from quaternion
     1000        q=r+ai+bj+ck
     1001    '''
     1002    aa = Q[0]**2
     1003    ab = Q[0]*Q[1]
     1004    ac = Q[0]*Q[2]
     1005    ad = Q[0]*Q[3]
     1006    bb = Q[1]**2
     1007    bc = Q[1]*Q[2]
     1008    bd = Q[1]*Q[3]
     1009    cc = Q[2]**2
     1010    cd = Q[2]*Q[3]
     1011    dd = Q[3]**2
     1012    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
     1013        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
     1014        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
     1015    return np.array(M)
     1016   
     1017def AV2Q(A,V):
     1018    ''' convert angle (radians -pi to pi) & vector to quaternion
     1019        q=r+ai+bj+ck
     1020    '''
     1021    Q = np.zeros(4)
     1022    d = np.sqrt(np.sum(np.array(V)**2))
     1023    if d:
     1024        V /= d
     1025    else:
     1026        return [1.,0.,0.,0.]    #identity
     1027    p = A/2.
     1028    Q[0] = np.cos(p)
     1029    s = np.sin(p)
     1030    Q[1:4] = V*s
     1031    return Q
     1032   
     1033def AVdeg2Q(A,V):
     1034    ''' convert angle (degrees -180 to 180) & vector to quaternion
     1035        q=r+ai+bj+ck
     1036    '''
     1037    Q = np.zeros(4)
     1038    d = np.sqrt(np.sum(np.array(V)**2))
     1039    if d:
     1040        V /= d
     1041    else:
     1042        return [1.,0.,0.,0.]    #identity
     1043    p = A/2.
     1044    Q[0] = cosd(p)
     1045    S = sind(p)
     1046    Q[1:4] = V*S
     1047    return Q
     1048   
     1049=======
     1050# -*- coding: utf-8 -*-
     1051#GSASIImath - major mathematics routines
     1052########### SVN repository information ###################
     1053# $Date$
     1054# $Author$
     1055# $Revision$
     1056# $URL$
     1057# $Id$
     1058########### SVN repository information ###################
     1059import sys
     1060import os
     1061import os.path as ospath
     1062import random as rn
     1063import numpy as np
     1064import numpy.linalg as nl
     1065import numpy.ma as ma
     1066import cPickle
     1067import time
     1068import math
     1069import copy
     1070import GSASIIpath
     1071GSASIIpath.SetVersionNumber("$Revision$")
     1072import GSASIIElem as G2el
     1073import GSASIIlattice as G2lat
     1074import GSASIIspc as G2spc
     1075import scipy.optimize as so
     1076import scipy.linalg as sl
     1077import numpy.fft as fft
     1078
     1079sind = lambda x: np.sin(x*np.pi/180.)
     1080cosd = lambda x: np.cos(x*np.pi/180.)
     1081tand = lambda x: np.tan(x*np.pi/180.)
     1082asind = lambda x: 180.*np.arcsin(x)/np.pi
     1083acosd = lambda x: 180.*np.arccos(x)/np.pi
     1084atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
     1085
     1086def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
     1087   
     1088    """
     1089    Minimize the sum of squares of a set of equations.
     1090
     1091    ::
     1092   
     1093                    Nobs
     1094        x = arg min(sum(func(y)**2,axis=0))
     1095                    y=0
     1096
     1097    Parameters
     1098    ----------
     1099    func : callable
     1100        should take at least one (possibly length N vector) argument and
     1101        returns M floating point numbers.
     1102    x0 : ndarray
     1103        The starting estimate for the minimization of length N
     1104    Hess : callable
     1105        A required function or method to compute the weighted vector and Hessian for func.
     1106        It must be a symmetric NxN array
     1107    args : tuple
     1108        Any extra arguments to func are placed in this tuple.
     1109    ftol : float
     1110        Relative error desired in the sum of squares.
     1111    xtol : float
     1112        Relative error desired in the approximate solution.
     1113    maxcyc : int
     1114        The maximum number of cycles of refinement to execute, if -1 refine
     1115        until other limits are met (ftol, xtol)
     1116
     1117    Returns
     1118    -------
     1119    x : ndarray
     1120        The solution (or the result of the last iteration for an unsuccessful
     1121        call).
     1122    cov_x : ndarray
     1123        Uses the fjac and ipvt optional outputs to construct an
     1124        estimate of the jacobian around the solution.  ``None`` if a
     1125        singular matrix encountered (indicates very flat curvature in
     1126        some direction).  This matrix must be multiplied by the
     1127        residual standard deviation to get the covariance of the
     1128        parameter estimates -- see curve_fit.
     1129    infodict : dict
     1130        a dictionary of optional outputs with the key s::
     1131
     1132            - 'fvec' : the function evaluated at the output
     1133
     1134
     1135    Notes
     1136    -----
     1137
     1138    """
     1139               
     1140    x0 = np.array(x0, ndmin=1)      #might be redundant?
     1141    n = len(x0)
     1142    if type(args) != type(()):
     1143        args = (args,)
     1144       
     1145    icycle = 0
     1146    One = np.ones((n,n))
     1147    lam = 0.001
     1148    lamMax = lam
     1149    nfev = 0
     1150    while icycle < maxcyc:
     1151        lamMax = max(lamMax,lam)
     1152        M = func(x0,*args)
     1153        nfev += 1
     1154        chisq0 = np.sum(M**2)
     1155        Yvec,Amat = Hess(x0,*args)
     1156        Adiag = np.sqrt(np.diag(Amat))
     1157        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
     1158        if np.any(psing):                #hard singularity in matrix
     1159            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
     1160        Anorm = np.outer(Adiag,Adiag)
     1161        Yvec /= Adiag
     1162        Amat /= Anorm       
     1163        while True:
     1164            Lam = np.eye(Amat.shape[0])*lam
     1165            Amatlam = Amat*(One+Lam)
     1166            try:
     1167                Xvec = nl.solve(Amatlam,Yvec)
     1168            except nl.LinAlgError:
     1169                print 'ouch #1'
     1170                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
     1171                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
     1172            Xvec /= Adiag
     1173            M2 = func(x0+Xvec,*args)
     1174            nfev += 1
     1175            chisq1 = np.sum(M2**2)
     1176            if chisq1 > chisq0:
     1177                lam *= 10.
     1178            else:
     1179                x0 += Xvec
     1180                lam /= 10.
     1181                break
     1182        if (chisq0-chisq1)/chisq0 < ftol:
     1183            break
     1184        icycle += 1
     1185    M = func(x0,*args)
     1186    nfev += 1
     1187    Yvec,Amat = Hess(x0,*args)
     1188    try:
     1189        Bmat = nl.inv(Amat)
     1190        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
     1191    except nl.LinAlgError:
     1192        print 'ouch #2 linear algebra error in LS'
     1193        psing = []
     1194        if maxcyc:
     1195            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
     1196        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
     1197   
     1198def getVCov(varyNames,varyList,covMatrix):
     1199    vcov = np.zeros((len(varyNames),len(varyNames)))
     1200    for i1,name1 in enumerate(varyNames):
     1201        for i2,name2 in enumerate(varyNames):
     1202            try:
     1203                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
     1204            except ValueError:
     1205                vcov[i1][i2] = 0.0
     1206    return vcov
     1207
     1208def getMass(generalData):
     1209    mass = 0.
     1210    for i,elem in enumerate(generalData['AtomTypes']):
     1211        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
     1212    return mass   
     1213
     1214def getDensity(generalData):
     1215   
     1216    mass = getMass(generalData)
     1217    Volume = generalData['Cell'][7]
     1218    density = mass/(0.6022137*Volume)
     1219    return density,Volume/mass
     1220   
     1221def getRestDist(XYZ,Amat):
     1222    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
     1223
     1224def getRestAngle(XYZ,Amat):
     1225   
     1226    def calcVec(Ox,Tx,Amat):
     1227        return np.inner(Amat,(Tx-Ox))
     1228
     1229    VecA = calcVec(XYZ[1],XYZ[0],Amat)
     1230    VecA /= np.sqrt(np.sum(VecA**2))
     1231    VecB = calcVec(XYZ[1],XYZ[2],Amat)
     1232    VecB /= np.sqrt(np.sum(VecB**2))
     1233    edge = VecB-VecA
     1234    edge = np.sum(edge**2)
     1235    angle = (2.-edge)/2.
     1236    angle = max(angle,-1.)
     1237    return acosd(angle)
     1238   
     1239def getRestPlane(XYZ,Amat):
     1240    sumXYZ = np.zeros(3)
     1241    for xyz in XYZ:
     1242        sumXYZ += xyz
     1243    sumXYZ /= len(XYZ)
     1244    XYZ = np.array(XYZ)-sumXYZ
     1245    XYZ = np.inner(Amat,XYZ).T
     1246    Zmat = np.zeros((3,3))
     1247    for i,xyz in enumerate(XYZ):
     1248        Zmat += np.outer(xyz.T,xyz)
     1249    Evec,Emat = nl.eig(Zmat)
     1250    Evec = np.sqrt(Evec)/(len(XYZ)-3)
     1251    Order = np.argsort(Evec)
     1252    return Evec[Order[0]]
     1253   
     1254def getRestChiral(XYZ,Amat):
     1255   
     1256    VecA = np.empty((3,3))   
     1257    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
     1258    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
     1259    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
     1260    return nl.det(VecA)
     1261       
     1262def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
     1263   
     1264    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
     1265        TxT = inv*(np.inner(M,Tx)+T)+C+U
     1266        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
     1267       
     1268    inv = Top/abs(Top)
     1269    cent = abs(Top)/100
     1270    op = abs(Top)%100-1
     1271    M,T = SGData['SGOps'][op]
     1272    C = SGData['SGCen'][cent]
     1273    dx = .00001
     1274    deriv = np.zeros(6)
     1275    for i in [0,1,2]:
     1276        Oxyz[i] += dx
     1277        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
     1278        Oxyz[i] -= 2*dx
     1279        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
     1280        Oxyz[i] += dx
     1281        Txyz[i] += dx
     1282        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
     1283        Txyz[i] -= 2*dx
     1284        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
     1285        Txyz[i] += dx
     1286    return deriv
     1287   
     1288def getAngSig(VA,VB,Amat,SGData,covData={}):
     1289   
     1290    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
     1291        TxT = inv*(np.inner(M,Tx)+T)+C
     1292        TxT = G2spc.MoveToUnitCell(TxT)+U
     1293        return np.inner(Amat,(TxT-Ox))
     1294       
     1295    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
     1296        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
     1297        VecA /= np.sqrt(np.sum(VecA**2))
     1298        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
     1299        VecB /= np.sqrt(np.sum(VecB**2))
     1300        edge = VecB-VecA
     1301        edge = np.sum(edge**2)
     1302        angle = (2.-edge)/2.
     1303        angle = max(angle,-1.)
     1304        return acosd(angle)
     1305       
     1306    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
     1307    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
     1308    invA = invB = 1
     1309    invA = TopA/abs(TopA)
     1310    invB = TopB/abs(TopB)
     1311    centA = abs(TopA)/100
     1312    centB = abs(TopB)/100
     1313    opA = abs(TopA)%100-1
     1314    opB = abs(TopB)%100-1
     1315    MA,TA = SGData['SGOps'][opA]
     1316    MB,TB = SGData['SGOps'][opB]
     1317    CA = SGData['SGCen'][centA]
     1318    CB = SGData['SGCen'][centB]
     1319    if 'covMatrix' in covData:
     1320        covMatrix = covData['covMatrix']
     1321        varyList = covData['varyList']
     1322        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
     1323        dx = .00001
     1324        dadx = np.zeros(9)
     1325        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
     1326        for i in [0,1,2]:
     1327            OxA[i] += dx
     1328            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
     1329            OxA[i] -= 2*dx
     1330            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
     1331            OxA[i] += dx
     1332           
     1333            TxA[i] += dx
     1334            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
     1335            TxA[i] -= 2*dx
     1336            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
     1337            TxA[i] += dx
     1338           
     1339            TxB[i] += dx
     1340            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
     1341            TxB[i] -= 2*dx
     1342            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
     1343            TxB[i] += dx
     1344           
     1345        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
     1346        if sigAng < 0.01:
     1347            sigAng = 0.0
     1348        return Ang,sigAng
     1349    else:
     1350        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
     1351
     1352def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
     1353
     1354    def calcDist(Atoms,SyOps,Amat):
     1355        XYZ = []
     1356        for i,atom in enumerate(Atoms):
     1357            Inv,M,T,C,U = SyOps[i]
     1358            XYZ.append(np.array(atom[1:4]))
     1359            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
     1360            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
     1361        V1 = XYZ[1]-XYZ[0]
     1362        return np.sqrt(np.sum(V1**2))
     1363       
     1364    Inv = []
     1365    SyOps = []
     1366    names = []
     1367    for i,atom in enumerate(Oatoms):
     1368        names += atom[-1]
     1369        Op,unit = Atoms[i][-1]
     1370        inv = Op/abs(Op)
     1371        m,t = SGData['SGOps'][abs(Op)%100-1]
     1372        c = SGData['SGCen'][abs(Op)/100]
     1373        SyOps.append([inv,m,t,c,unit])
     1374    Dist = calcDist(Oatoms,SyOps,Amat)
     1375   
     1376    sig = -0.001
     1377    if 'covMatrix' in covData:
     1378        parmNames = []
     1379        dx = .00001
     1380        dadx = np.zeros(6)
     1381        for i in range(6):
     1382            ia = i/3
     1383            ix = i%3
     1384            Oatoms[ia][ix+1] += dx
     1385            a0 = calcDist(Oatoms,SyOps,Amat)
     1386            Oatoms[ia][ix+1] -= 2*dx
     1387            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
     1388        covMatrix = covData['covMatrix']
     1389        varyList = covData['varyList']
     1390        DistVcov = getVCov(names,varyList,covMatrix)
     1391        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
     1392        if sig < 0.001:
     1393            sig = -0.001
     1394   
     1395    return Dist,sig
     1396
     1397def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
     1398       
     1399    def calcAngle(Atoms,SyOps,Amat):
     1400        XYZ = []
     1401        for i,atom in enumerate(Atoms):
     1402            Inv,M,T,C,U = SyOps[i]
     1403            XYZ.append(np.array(atom[1:4]))
     1404            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
     1405            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
     1406        V1 = XYZ[1]-XYZ[0]
     1407        V1 /= np.sqrt(np.sum(V1**2))
     1408        V2 = XYZ[1]-XYZ[2]
     1409        V2 /= np.sqrt(np.sum(V2**2))
     1410        V3 = V2-V1
     1411        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
     1412        return acosd(cang)
     1413
     1414    Inv = []
     1415    SyOps = []
     1416    names = []
     1417    for i,atom in enumerate(Oatoms):
     1418        names += atom[-1]
     1419        Op,unit = Atoms[i][-1]
     1420        inv = Op/abs(Op)
     1421        m,t = SGData['SGOps'][abs(Op)%100-1]
     1422        c = SGData['SGCen'][abs(Op)/100]
     1423        SyOps.append([inv,m,t,c,unit])
     1424    Angle = calcAngle(Oatoms,SyOps,Amat)
     1425   
     1426    sig = -0.01
     1427    if 'covMatrix' in covData:
     1428        parmNames = []
     1429        dx = .00001
     1430        dadx = np.zeros(9)
     1431        for i in range(9):
     1432            ia = i/3
     1433            ix = i%3
     1434            Oatoms[ia][ix+1] += dx
     1435            a0 = calcAngle(Oatoms,SyOps,Amat)
     1436            Oatoms[ia][ix+1] -= 2*dx
     1437            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
     1438        covMatrix = covData['covMatrix']
     1439        varyList = covData['varyList']
     1440        AngVcov = getVCov(names,varyList,covMatrix)
     1441        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
     1442        if sig < 0.01:
     1443            sig = -0.01
     1444   
     1445    return Angle,sig
     1446
     1447def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
     1448   
     1449    def calcTorsion(Atoms,SyOps,Amat):
     1450       
     1451        XYZ = []
     1452        for i,atom in enumerate(Atoms):
     1453            Inv,M,T,C,U = SyOps[i]
     1454            XYZ.append(np.array(atom[1:4]))
     1455            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
     1456            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
     1457        V1 = XYZ[1]-XYZ[0]
     1458        V2 = XYZ[2]-XYZ[1]
     1459        V3 = XYZ[3]-XYZ[2]
     1460        V1 /= np.sqrt(np.sum(V1**2))
     1461        V2 /= np.sqrt(np.sum(V2**2))
     1462        V3 /= np.sqrt(np.sum(V3**2))
     1463        M = np.array([V1,V2,V3])
     1464        D = nl.det(M)
     1465        Ang = 1.0
     1466        P12 = np.dot(V1,V2)
     1467        P13 = np.dot(V1,V3)
     1468        P23 = np.dot(V2,V3)
     1469        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
     1470        return Tors
     1471           
     1472    Inv = []
     1473    SyOps = []
     1474    names = []
     1475    for i,atom in enumerate(Oatoms):
     1476        names += atom[-1]
     1477        Op,unit = Atoms[i][-1]
     1478        inv = Op/abs(Op)
     1479        m,t = SGData['SGOps'][abs(Op)%100-1]
     1480        c = SGData['SGCen'][abs(Op)/100]
     1481        SyOps.append([inv,m,t,c,unit])
     1482    Tors = calcTorsion(Oatoms,SyOps,Amat)
     1483   
     1484    sig = -0.01
     1485    if 'covMatrix' in covData:
     1486        parmNames = []
     1487        dx = .00001
     1488        dadx = np.zeros(12)
     1489        for i in range(12):
     1490            ia = i/3
     1491            ix = i%3
     1492            Oatoms[ia][ix+1] += dx
     1493            a0 = calcTorsion(Oatoms,SyOps,Amat)
     1494            Oatoms[ia][ix+1] -= 2*dx
     1495            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
     1496        covMatrix = covData['covMatrix']
     1497        varyList = covData['varyList']
     1498        TorVcov = getVCov(names,varyList,covMatrix)
     1499        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
     1500        if sig < 0.01:
     1501            sig = -0.01
     1502   
     1503    return Tors,sig
     1504       
     1505def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
     1506   
     1507    def calcDist(Atoms,SyOps,Amat):
     1508        XYZ = []
     1509        for i,atom in enumerate(Atoms):
     1510            Inv,M,T,C,U = SyOps[i]
     1511            XYZ.append(np.array(atom[1:4]))
     1512            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
     1513            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
     1514        V1 = XYZ[1]-XYZ[0]
     1515        return np.sqrt(np.sum(V1**2))
     1516       
     1517    def calcAngle(Atoms,SyOps,Amat):
     1518        XYZ = []
     1519        for i,atom in enumerate(Atoms):
     1520            Inv,M,T,C,U = SyOps[i]
     1521            XYZ.append(np.array(atom[1:4]))
     1522            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
     1523            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
     1524        V1 = XYZ[1]-XYZ[0]
     1525        V1 /= np.sqrt(np.sum(V1**2))
     1526        V2 = XYZ[1]-XYZ[2]
     1527        V2 /= np.sqrt(np.sum(V2**2))
     1528        V3 = V2-V1
     1529        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
     1530        return acosd(cang)
     1531
     1532    def calcTorsion(Atoms,SyOps,Amat):
     1533       
     1534        XYZ = []
     1535        for i,atom in enumerate(Atoms):
     1536            Inv,M,T,C,U = SyOps[i]
     1537            XYZ.append(np.array(atom[1:4]))
     1538            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
     1539            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
     1540        V1 = XYZ[1]-XYZ[0]
     1541        V2 = XYZ[2]-XYZ[1]
     1542        V3 = XYZ[3]-XYZ[2]
     1543        V1 /= np.sqrt(np.sum(V1**2))
     1544        V2 /= np.sqrt(np.sum(V2**2))
     1545        V3 /= np.sqrt(np.sum(V3**2))
     1546        M = np.array([V1,V2,V3])
     1547        D = nl.det(M)
     1548        Ang = 1.0
     1549        P12 = np.dot(V1,V2)
     1550        P13 = np.dot(V1,V3)
     1551        P23 = np.dot(V2,V3)
     1552        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
     1553        return Tors
     1554           
     1555    Inv = []
     1556    SyOps = []
     1557    names = []
     1558    for i,atom in enumerate(Oatoms):
     1559        names += atom[-1]
     1560        Op,unit = Atoms[i][-1]
     1561        inv = Op/abs(Op)
     1562        m,t = SGData['SGOps'][abs(Op)%100-1]
     1563        c = SGData['SGCen'][abs(Op)/100]
     1564        SyOps.append([inv,m,t,c,unit])
     1565    M = len(Oatoms)
     1566    if M == 2:
     1567        Val = calcDist(Oatoms,SyOps,Amat)
     1568    elif M == 3:
     1569        Val = calcAngle(Oatoms,SyOps,Amat)
     1570    else:
     1571        Val = calcTorsion(Oatoms,SyOps,Amat)
     1572   
     1573    sigVals = [-0.001,-0.01,-0.01]
     1574    sig = sigVals[M-3]
     1575    if 'covMatrix' in covData:
     1576        parmNames = []
     1577        dx = .00001
     1578        N = M*3
     1579        dadx = np.zeros(N)
     1580        for i in range(N):
     1581            ia = i/3
     1582            ix = i%3
     1583            Oatoms[ia][ix+1] += dx
     1584            if M == 2:
     1585                a0 = calcDist(Oatoms,SyOps,Amat)
     1586            elif M == 3:
     1587                a0 = calcAngle(Oatoms,SyOps,Amat)
     1588            else:
     1589                a0 = calcTorsion(Oatoms,SyOps,Amat)
     1590            Oatoms[ia][ix+1] -= 2*dx
     1591            if M == 2:
     1592                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
     1593            elif M == 3:
     1594                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
     1595            else:
     1596                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
     1597        covMatrix = covData['covMatrix']
     1598        varyList = covData['varyList']
     1599        Vcov = getVCov(names,varyList,covMatrix)
     1600        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
     1601        if sig < sigVals[M-3]:
     1602            sig = sigVals[M-3]
     1603   
     1604    return Val,sig
     1605       
     1606   
     1607def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
     1608    # returns value(esd) string; nTZ=True for no trailing zeros
     1609    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
     1610    #get the 2 significant digits in the esd
     1611    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
     1612    #get the number of digits to represent them
     1613    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
     1614   
     1615    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
     1616    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
     1617    if esd > 0:
     1618        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
     1619        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
     1620    elif esd < 0:
     1621         return str(round(value,mdec(esd)-1))
     1622    else:
     1623        text = str("%f"%(value))
     1624        if nTZ:
     1625            return text.rstrip('0')
     1626        else:
     1627            return text
     1628           
     1629def adjHKLmax(SGData,Hmax):
     1630    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
     1631        Hmax[0] = ((Hmax[0]+3)/6)*6
     1632        Hmax[1] = ((Hmax[1]+3)/6)*6
     1633        Hmax[2] = ((Hmax[2]+1)/4)*4
     1634    else:
     1635        Hmax[0] = ((Hmax[0]+2)/4)*4
     1636        Hmax[1] = ((Hmax[1]+2)/4)*4
     1637        Hmax[2] = ((Hmax[2]+1)/4)*4
     1638
     1639def FourierMap(data,reflData):
     1640   
     1641    generalData = data['General']
     1642    if not generalData['Map']['MapType']:
     1643        print '**** ERROR - Fourier map not defined'
     1644        return
     1645    mapData = generalData['Map']
     1646    dmin = mapData['Resolution']
     1647    SGData = generalData['SGData']
     1648    cell = generalData['Cell'][1:8]       
     1649    A = G2lat.cell2A(cell[:6])
     1650    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
     1651    adjHKLmax(SGData,Hmax)
     1652    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
     1653#    Fhkl[0,0,0] = generalData['F000X']
     1654    time0 = time.time()
     1655    for ref in reflData:
     1656        if ref[4] >= dmin:
     1657            Fosq,Fcsq,ph = ref[8:11]
     1658            for i,hkl in enumerate(ref[11]):
     1659                hkl = np.asarray(hkl,dtype='i')
     1660                dp = 360.*ref[12][i]
     1661                a = cosd(ph+dp)
     1662                b = sind(ph+dp)
     1663                phasep = complex(a,b)
     1664                phasem = complex(a,-b)
     1665                if 'Fobs' in mapData['MapType']:
     1666                    F = np.sqrt(Fosq)
     1667                    h,k,l = hkl+Hmax
     1668                    Fhkl[h,k,l] = F*phasep
     1669                    h,k,l = -hkl+Hmax
     1670                    Fhkl[h,k,l] = F*phasem
     1671                elif 'Fcalc' in mapData['MapType']:
     1672                    F = np.sqrt(Fcsq)
     1673                    h,k,l = hkl+Hmax
     1674                    Fhkl[h,k,l] = F*phasep
     1675                    h,k,l = -hkl+Hmax
     1676                    Fhkl[h,k,l] = F*phasem
     1677                elif 'delt-F' in mapData['MapType']:
     1678                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
     1679                    h,k,l = hkl+Hmax
     1680                    Fhkl[h,k,l] = dF*phasep
     1681                    h,k,l = -hkl+Hmax
     1682                    Fhkl[h,k,l] = dF*phasem
     1683                elif '2*Fo-Fc' in mapData['MapType']:
     1684                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
     1685                    h,k,l = hkl+Hmax
     1686                    Fhkl[h,k,l] = F*phasep
     1687                    h,k,l = -hkl+Hmax
     1688                    Fhkl[h,k,l] = F*phasem
     1689                elif 'Patterson' in mapData['MapType']:
     1690                    h,k,l = hkl+Hmax
     1691                    Fhkl[h,k,l] = complex(Fosq,0.)
     1692                    h,k,l = -hkl+Hmax
     1693                    Fhkl[h,k,l] = complex(Fosq,0.)
     1694    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
     1695    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
     1696    mapData['rho'] = np.real(rho)
     1697    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
     1698    return mapData
     1699   
     1700# map printing for testing purposes
     1701def printRho(SGLaue,rho,rhoMax):                         
     1702    dim = len(rho.shape)
     1703    if dim == 2:
     1704        ix,jy = rho.shape
     1705        for j in range(jy):
     1706            line = ''
     1707            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
     1708                line += (jy-j)*'  '
     1709            for i in range(ix):
     1710                r = int(100*rho[i,j]/rhoMax)
     1711                line += '%4d'%(r)
     1712            print line+'\n'
     1713    else:
     1714        ix,jy,kz = rho.shape
     1715        for k in range(kz):
     1716            print 'k = ',k
     1717            for j in range(jy):
     1718                line = ''
     1719                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
     1720                    line += (jy-j)*'  '
     1721                for i in range(ix):
     1722                    r = int(100*rho[i,j,k]/rhoMax)
     1723                    line += '%4d'%(r)
     1724                print line+'\n'
     1725## keep this
     1726               
     1727def findOffset(SGData,A,Fhkl):   
     1728    if SGData['SpGrp'] == 'P 1':
     1729        return [0,0,0]   
     1730    hklShape = Fhkl.shape
     1731    steps = np.array(hklShape)
     1732    Hmax = 2*np.asarray(G2lat.getHKLmax(4.5,SGData,A),dtype='i')
     1733    Fmax = np.max(np.absolute(Fhkl))
     1734    hklHalf = np.array(hklShape)/2
     1735    sortHKL = np.argsort(Fhkl.flatten())
     1736    Fdict = {}
     1737    for hkl in sortHKL:
     1738        HKL = np.unravel_index(hkl,hklShape)
     1739        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
     1740        if F == 0.:
     1741            break
     1742        Fdict['%.6f'%(np.absolute(F))] = hkl
     1743    Flist = np.flipud(np.sort(Fdict.keys()))
     1744    F = str(1.e6)
     1745    i = 0
     1746    DH = []
     1747    Dphi = []
     1748    while i < 20 and len(DH) < 30:
     1749        F = Flist[i]
     1750        hkl = np.unravel_index(Fdict[F],hklShape)
     1751        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
     1752        Uniq = np.array(Uniq,dtype='i')
     1753        Phi = np.array(Phi)
     1754        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
     1755        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
     1756        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
     1757        ang0 = np.angle(Fh0,deg=True)/360.
     1758        for j,H in enumerate(Uniq[1:]):
     1759            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-Phi[j+1])
     1760            dH = H-hkl
     1761            dang = ang-ang0
     1762            if np.any(np.abs(dH)-Hmax > 0):    #keep low order DHs
     1763                continue
     1764            DH.append(dH)
     1765            Dphi.append((dang+0.5) % 1.0)
     1766        i += 1
     1767    DH = np.array(DH)
     1768    print ' map offset no.of terms: %d'%(len(DH))
     1769    Dphi = np.array(Dphi)
     1770    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
     1771    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
     1772    Mmap = np.reshape(np.sum(((np.dot(XYZ,DH.T)+.5)%1.-Dphi)**2,axis=1),newshape=steps)
     1773    chisq = np.min(Mmap)
     1774    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
     1775    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
     1776    return DX
     1777   
     1778def ChargeFlip(data,reflData,pgbar):
     1779    generalData = data['General']
     1780    mapData = generalData['Map']
     1781    flipData = generalData['Flip']
     1782    FFtable = {}
     1783    if 'None' not in flipData['Norm element']:
     1784        normElem = flipData['Norm element'].upper()
     1785        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
     1786        for ff in FFs:
     1787            if ff['Symbol'] == normElem:
     1788                FFtable.update(ff)
     1789    dmin = flipData['Resolution']
     1790    SGData = generalData['SGData']
     1791    cell = generalData['Cell'][1:8]       
     1792    A = G2lat.cell2A(cell[:6])
     1793    Vol = cell[6]
     1794    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
     1795    adjHKLmax(SGData,Hmax)
     1796    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
     1797    time0 = time.time()
     1798    for ref in reflData:
     1799        dsp = ref[4]
     1800        if dsp >= dmin:
     1801            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
     1802            if FFtable:
     1803                SQ = 0.25/dsp**2
     1804                ff *= G2el.ScatFac(FFtable,SQ)[0]
     1805            if ref[8] > 0.:
     1806                E = np.sqrt(ref[8])/ff
     1807            else:
     1808                E = 0.
     1809            ph = ref[10]
     1810            ph = rn.uniform(0.,360.)
     1811            for i,hkl in enumerate(ref[11]):
     1812                hkl = np.asarray(hkl,dtype='i')
     1813                dp = 360.*ref[12][i]
     1814                a = cosd(ph+dp)
     1815                b = sind(ph+dp)
     1816                phasep = complex(a,b)
     1817                phasem = complex(a,-b)
     1818                h,k,l = hkl+Hmax
     1819                Ehkl[h,k,l] = E*phasep
     1820                h,k,l = -hkl+Hmax       #Friedel pair refl.
     1821                Ehkl[h,k,l] = E*phasem
     1822#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
     1823    CEhkl = copy.copy(Ehkl)
     1824    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
     1825    Emask = ma.getmask(MEhkl)
     1826    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
     1827    Ncyc = 0
     1828    old = np.seterr(all='raise')
     1829    while True:       
     1830        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
     1831        CEsig = np.std(CErho)
     1832        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
     1833        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
     1834        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
     1835        phase = CFhkl/np.absolute(CFhkl)
     1836        CEhkl = np.absolute(Ehkl)*phase
     1837        Ncyc += 1
     1838        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
     1839        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
     1840        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
     1841        if Rcf < 5.:
     1842            break
     1843        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
     1844        if not GoOn or Ncyc > 10000:
     1845            break
     1846    np.seterr(**old)
     1847    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
     1848    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
     1849    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
     1850    roll = findOffset(SGData,A,CEhkl)
     1851       
     1852    mapData['Rcf'] = Rcf
     1853    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
     1854    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
     1855    mapData['rollMap'] = [0,0,0]
     1856    return mapData
     1857   
    8091858def SearchMap(data,keepDup=False,Pgbar=None):
    8101859   
     
    10912140    return Q
    10922141   
     2142>>>>>>> .r762
Note: See TracChangeset for help on using the changeset viewer.