source: trunk/GSASIImath.py @ 568

Last change on this file since 568 was 568, checked in by vondreele, 10 years ago

remove gsas data window placement
moved peaks are C atoms
start on a relocate map routine???

File size: 27.2 KB
Line 
1#GSASIImath - major mathematics routines
2########### SVN repository information ###################
3# $Date: 2012-01-13 11:48:53 -0600 (Fri, 13 Jan 2012) $
4# $Author: vondreele & toby $
5# $Revision: 451 $
6# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIImath.py $
7# $Id: GSASIImath.py 451 2012-01-13 17:48:53Z vondreele $
8########### SVN repository information ###################
9import sys
10import os
11import os.path as ospath
12import random as rn
13import numpy as np
14import numpy.linalg as nl
15import numpy.ma as ma
16import cPickle
17import time
18import math
19import copy
20import GSASIIpath
21import GSASIIElem as G2el
22import GSASIIlattice as G2lat
23import GSASIIspc as G2spc
24import scipy.optimize as so
25import scipy.linalg as sl
26
27sind = lambda x: np.sin(x*np.pi/180.)
28cosd = lambda x: np.cos(x*np.pi/180.)
29tand = lambda x: np.tan(x*np.pi/180.)
30asind = lambda x: 180.*np.arcsin(x)/np.pi
31acosd = lambda x: 180.*np.arccos(x)/np.pi
32atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
33
34def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
35   
36    """
37    Minimize the sum of squares of a set of equations.
38
39    ::
40   
41                    Nobs
42        x = arg min(sum(func(y)**2,axis=0))
43                    y=0
44
45    Parameters
46    ----------
47    func : callable
48        should take at least one (possibly length N vector) argument and
49        returns M floating point numbers.
50    x0 : ndarray
51        The starting estimate for the minimization of length N
52    Hess : callable
53        A required function or method to compute the weighted vector and Hessian for func.
54        It must be a symmetric NxN array
55    args : tuple
56        Any extra arguments to func are placed in this tuple.
57    ftol : float
58        Relative error desired in the sum of squares.
59    xtol : float
60        Relative error desired in the approximate solution.
61    maxcyc : int
62        The maximum number of cycles of refinement to execute, if -1 refine
63        until other limits are met (ftol, xtol)
64
65    Returns
66    -------
67    x : ndarray
68        The solution (or the result of the last iteration for an unsuccessful
69        call).
70    cov_x : ndarray
71        Uses the fjac and ipvt optional outputs to construct an
72        estimate of the jacobian around the solution.  ``None`` if a
73        singular matrix encountered (indicates very flat curvature in
74        some direction).  This matrix must be multiplied by the
75        residual standard deviation to get the covariance of the
76        parameter estimates -- see curve_fit.
77    infodict : dict
78        a dictionary of optional outputs with the key s::
79
80            - 'fvec' : the function evaluated at the output
81
82
83    Notes
84    -----
85
86    """
87               
88    x0 = np.array(x0, ndmin=1)      #might be redundant?
89    n = len(x0)
90    if type(args) != type(()):
91        args = (args,)
92       
93    icycle = 0
94    One = np.ones((n,n))
95    lam = 0.001
96    lamMax = lam
97    nfev = 0
98    while icycle < maxcyc:
99        lamMax = max(lamMax,lam)
100        M = func(x0,*args)
101        nfev += 1
102        chisq0 = np.sum(M**2)
103        Yvec,Amat = Hess(x0,*args)
104        Adiag = np.sqrt(np.diag(Amat))
105        if 0.0 in Adiag:                #hard singularity in matrix
106            psing = list(np.where(Adiag == 0.)[0])
107            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
108        Anorm = np.outer(Adiag,Adiag)
109        Yvec /= Adiag
110        Amat /= Anorm       
111        while True:
112            Lam = np.eye(Amat.shape[0])*lam
113            Amatlam = Amat*(One+Lam)
114            try:
115                Xvec = nl.solve(Amatlam,Yvec)
116            except LinAlgError:
117                psing = list(np.where(np.diag(nl.gr(Amatlam)[1]) < 1.e-14)[0])
118                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
119            Xvec /= Adiag
120            M2 = func(x0+Xvec,*args)
121            nfev += 1
122            chisq1 = np.sum(M2**2)
123            if chisq1 > chisq0:
124                lam *= 10.
125            else:
126                x0 += Xvec
127                lam /= 10.
128                break
129        if (chisq0-chisq1)/chisq0 < ftol:
130            break
131        icycle += 1
132    M = func(x0,*args)
133    nfev += 1
134    Yvec,Amat = Hess(x0,*args)
135    try:
136        Bmat = nl.inv(Amat)
137        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
138    except nl.LinAlgError:
139        psing = []
140        if maxcyc:
141            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
142        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
143   
144def getVCov(varyNames,varyList,covMatrix):
145    vcov = np.zeros((len(varyNames),len(varyNames)))
146    for i1,name1 in enumerate(varyNames):
147        for i2,name2 in enumerate(varyNames):
148            try:
149                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
150            except ValueError:
151                vcov[i1][i2] = 0.0
152    return vcov
153   
154def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
155   
156    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
157        TxT = inv*(np.inner(M,Tx)+T)+C+U
158        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
159       
160    inv = Top/abs(Top)
161    cent = abs(Top)/100
162    op = abs(Top)%100-1
163    M,T = SGData['SGOps'][op]
164    C = SGData['SGCen'][cent]
165    dx = .00001
166    deriv = np.zeros(6)
167    for i in [0,1,2]:
168        Oxyz[i] += dx
169        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
170        Oxyz[i] -= 2*dx
171        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
172        Oxyz[i] += dx
173        Txyz[i] += dx
174        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
175        Txyz[i] -= 2*dx
176        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
177        Txyz[i] += dx
178    return deriv
179   
180def getAngSig(VA,VB,Amat,SGData,covData={}):
181   
182    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
183        TxT = inv*(np.inner(M,Tx)+T)+C
184        TxT = G2spc.MoveToUnitCell(TxT)+U
185        return np.inner(Amat,(TxT-Ox))
186       
187    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
188        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
189        VecA /= np.sqrt(np.sum(VecA**2))
190        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
191        VecB /= np.sqrt(np.sum(VecB**2))
192        edge = VecB-VecA
193        edge = np.sum(edge**2)
194        angle = (2.-edge)/2.
195        angle = max(angle,-1.)
196        return acosd(angle)
197       
198    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
199    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
200    invA = invB = 1
201    invA = TopA/abs(TopA)
202    invB = TopB/abs(TopB)
203    centA = abs(TopA)/100
204    centB = abs(TopB)/100
205    opA = abs(TopA)%100-1
206    opB = abs(TopB)%100-1
207    MA,TA = SGData['SGOps'][opA]
208    MB,TB = SGData['SGOps'][opB]
209    CA = SGData['SGCen'][centA]
210    CB = SGData['SGCen'][centB]
211    if 'covMatrix' in covData:
212        covMatrix = covData['covMatrix']
213        varyList = covData['varyList']
214        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
215        dx = .00001
216        dadx = np.zeros(9)
217        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
218        for i in [0,1,2]:
219            OxA[i] += dx
220            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
221            OxA[i] -= 2*dx
222            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
223            OxA[i] += dx
224           
225            TxA[i] += dx
226            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
227            TxA[i] -= 2*dx
228            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
229            TxA[i] += dx
230           
231            TxB[i] += dx
232            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
233            TxB[i] -= 2*dx
234            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
235            TxB[i] += dx
236           
237        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
238        if sigAng < 0.01:
239            sigAng = 0.0
240        return Ang,sigAng
241    else:
242        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
243
244def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
245
246    def calcDist(Atoms,SyOps,Amat):
247        XYZ = []
248        for i,atom in enumerate(Atoms):
249            Inv,M,T,C,U = SyOps[i]
250            XYZ.append(np.array(atom[1:4]))
251            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
252            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
253        V1 = XYZ[1]-XYZ[0]
254        return np.sqrt(np.sum(V1**2))
255       
256    Inv = []
257    SyOps = []
258    names = []
259    for i,atom in enumerate(Oatoms):
260        names += atom[-1]
261        Op,unit = Atoms[i][-1]
262        inv = Op/abs(Op)
263        m,t = SGData['SGOps'][abs(Op)%100-1]
264        c = SGData['SGCen'][abs(Op)/100]
265        SyOps.append([inv,m,t,c,unit])
266    Dist = calcDist(Oatoms,SyOps,Amat)
267   
268    sig = -0.001
269    if 'covMatrix' in covData:
270        parmNames = []
271        dx = .00001
272        dadx = np.zeros(6)
273        for i in range(6):
274            ia = i/3
275            ix = i%3
276            Oatoms[ia][ix+1] += dx
277            a0 = calcDist(Oatoms,SyOps,Amat)
278            Oatoms[ia][ix+1] -= 2*dx
279            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
280        covMatrix = covData['covMatrix']
281        varyList = covData['varyList']
282        DistVcov = getVCov(names,varyList,covMatrix)
283        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
284        if sig < 0.001:
285            sig = -0.001
286   
287    return Dist,sig
288
289def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
290       
291    def calcAngle(Atoms,SyOps,Amat):
292        XYZ = []
293        for i,atom in enumerate(Atoms):
294            Inv,M,T,C,U = SyOps[i]
295            XYZ.append(np.array(atom[1:4]))
296            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
297            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
298        V1 = XYZ[1]-XYZ[0]
299        V1 /= np.sqrt(np.sum(V1**2))
300        V2 = XYZ[1]-XYZ[2]
301        V2 /= np.sqrt(np.sum(V2**2))
302        V3 = V2-V1
303        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
304        return acosd(cang)
305
306    Inv = []
307    SyOps = []
308    names = []
309    for i,atom in enumerate(Oatoms):
310        names += atom[-1]
311        Op,unit = Atoms[i][-1]
312        inv = Op/abs(Op)
313        m,t = SGData['SGOps'][abs(Op)%100-1]
314        c = SGData['SGCen'][abs(Op)/100]
315        SyOps.append([inv,m,t,c,unit])
316    Angle = calcAngle(Oatoms,SyOps,Amat)
317   
318    sig = -0.01
319    if 'covMatrix' in covData:
320        parmNames = []
321        dx = .00001
322        dadx = np.zeros(9)
323        for i in range(9):
324            ia = i/3
325            ix = i%3
326            Oatoms[ia][ix+1] += dx
327            a0 = calcAngle(Oatoms,SyOps,Amat)
328            Oatoms[ia][ix+1] -= 2*dx
329            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
330        covMatrix = covData['covMatrix']
331        varyList = covData['varyList']
332        AngVcov = getVCov(names,varyList,covMatrix)
333        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
334        if sig < 0.01:
335            sig = -0.01
336   
337    return Angle,sig
338
339def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
340   
341    def calcTorsion(Atoms,SyOps,Amat):
342       
343        XYZ = []
344        for i,atom in enumerate(Atoms):
345            Inv,M,T,C,U = SyOps[i]
346            XYZ.append(np.array(atom[1:4]))
347            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
348            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
349        V1 = XYZ[1]-XYZ[0]
350        V2 = XYZ[2]-XYZ[1]
351        V3 = XYZ[3]-XYZ[2]
352        V1 /= np.sqrt(np.sum(V1**2))
353        V2 /= np.sqrt(np.sum(V2**2))
354        V3 /= np.sqrt(np.sum(V3**2))
355        M = np.array([V1,V2,V3])
356        D = nl.det(M)
357        Ang = 1.0
358        P12 = np.dot(V1,V2)
359        P13 = np.dot(V1,V3)
360        P23 = np.dot(V2,V3)
361        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
362        return Tors
363           
364    Inv = []
365    SyOps = []
366    names = []
367    for i,atom in enumerate(Oatoms):
368        names += atom[-1]
369        Op,unit = Atoms[i][-1]
370        inv = Op/abs(Op)
371        m,t = SGData['SGOps'][abs(Op)%100-1]
372        c = SGData['SGCen'][abs(Op)/100]
373        SyOps.append([inv,m,t,c,unit])
374    Tors = calcTorsion(Oatoms,SyOps,Amat)
375   
376    sig = -0.01
377    if 'covMatrix' in covData:
378        parmNames = []
379        dx = .00001
380        dadx = np.zeros(12)
381        for i in range(12):
382            ia = i/3
383            ix = i%3
384            Oatoms[ia][ix+1] += dx
385            a0 = calcTorsion(Oatoms,SyOps,Amat)
386            Oatoms[ia][ix+1] -= 2*dx
387            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
388        covMatrix = covData['covMatrix']
389        varyList = covData['varyList']
390        TorVcov = getVCov(names,varyList,covMatrix)
391        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
392        if sig < 0.01:
393            sig = -0.01
394   
395    return Tors,sig
396       
397def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
398   
399    def calcDist(Atoms,SyOps,Amat):
400        XYZ = []
401        for i,atom in enumerate(Atoms):
402            Inv,M,T,C,U = SyOps[i]
403            XYZ.append(np.array(atom[1:4]))
404            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
405            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
406        V1 = XYZ[1]-XYZ[0]
407        return np.sqrt(np.sum(V1**2))
408       
409    def calcAngle(Atoms,SyOps,Amat):
410        XYZ = []
411        for i,atom in enumerate(Atoms):
412            Inv,M,T,C,U = SyOps[i]
413            XYZ.append(np.array(atom[1:4]))
414            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
415            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
416        V1 = XYZ[1]-XYZ[0]
417        V1 /= np.sqrt(np.sum(V1**2))
418        V2 = XYZ[1]-XYZ[2]
419        V2 /= np.sqrt(np.sum(V2**2))
420        V3 = V2-V1
421        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
422        return acosd(cang)
423
424    def calcTorsion(Atoms,SyOps,Amat):
425       
426        XYZ = []
427        for i,atom in enumerate(Atoms):
428            Inv,M,T,C,U = SyOps[i]
429            XYZ.append(np.array(atom[1:4]))
430            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
431            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
432        V1 = XYZ[1]-XYZ[0]
433        V2 = XYZ[2]-XYZ[1]
434        V3 = XYZ[3]-XYZ[2]
435        V1 /= np.sqrt(np.sum(V1**2))
436        V2 /= np.sqrt(np.sum(V2**2))
437        V3 /= np.sqrt(np.sum(V3**2))
438        M = np.array([V1,V2,V3])
439        D = nl.det(M)
440        Ang = 1.0
441        P12 = np.dot(V1,V2)
442        P13 = np.dot(V1,V3)
443        P23 = np.dot(V2,V3)
444        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
445        return Tors
446           
447    Inv = []
448    SyOps = []
449    names = []
450    for i,atom in enumerate(Oatoms):
451        names += atom[-1]
452        Op,unit = Atoms[i][-1]
453        inv = Op/abs(Op)
454        m,t = SGData['SGOps'][abs(Op)%100-1]
455        c = SGData['SGCen'][abs(Op)/100]
456        SyOps.append([inv,m,t,c,unit])
457    M = len(Oatoms)
458    if M == 2:
459        Val = calcDist(Oatoms,SyOps,Amat)
460    elif M == 3:
461        Val = calcAngle(Oatoms,SyOps,Amat)
462    else:
463        Val = calcTorsion(Oatoms,SyOps,Amat)
464   
465    sigVals = [-0.001,-0.01,-0.01]
466    sig = sigVals[M-3]
467    if 'covMatrix' in covData:
468        parmNames = []
469        dx = .00001
470        N = M*3
471        dadx = np.zeros(N)
472        for i in range(N):
473            ia = i/3
474            ix = i%3
475            Oatoms[ia][ix+1] += dx
476            if M == 2:
477                a0 = calcDist(Oatoms,SyOps,Amat)
478            elif M == 3:
479                a0 = calcAngle(Oatoms,SyOps,Amat)
480            else:
481                a0 = calcTorsion(Oatoms,SyOps,Amat)
482            Oatoms[ia][ix+1] -= 2*dx
483            if M == 2:
484                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
485            elif M == 3:
486                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
487            else:
488                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
489        covMatrix = covData['covMatrix']
490        varyList = covData['varyList']
491        Vcov = getVCov(names,varyList,covMatrix)
492        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
493        if sig < sigVals[M-3]:
494            sig = sigVals[M-3]
495   
496    return Val,sig
497       
498   
499def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
500    # returns value(esd) string; nTZ=True for no trailing zeros
501    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
502    #get the 2 significant digits in the esd
503    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
504    #get the number of digits to represent them
505    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
506   
507    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
508    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
509    if esd > 0:
510        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
511        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
512    elif esd < 0:
513         return str(round(value,mdec(esd)-1))
514    else:
515        text = str("%f"%(value))
516        if nTZ:
517            return text.rstrip('0')
518        else:
519            return text
520
521def FourierMap(data,reflData):
522   
523    import scipy.fftpack as fft
524    generalData = data['General']
525    if not generalData['Map']['MapType']:
526        print '**** ERROR - Fourier map not defined'
527        return
528    mapData = generalData['Map']
529    dmin = mapData['Resolution']
530    SGData = generalData['SGData']
531    cell = generalData['Cell'][1:8]       
532    A = G2lat.cell2A(cell[:6])
533    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
534    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
535#    Fhkl[0,0,0] = generalData['F000X']
536    time0 = time.time()
537    for ref in reflData:
538        if ref[4] >= dmin:
539            Fosq,Fcsq,ph = ref[8:11]
540            for i,hkl in enumerate(ref[11]):
541                hkl = np.asarray(hkl,dtype='i')
542                dp = 360.*ref[12][i]
543                a = cosd(ph+dp)
544                b = sind(ph+dp)
545                phasep = complex(a,b)
546                phasem = complex(a,-b)
547                if 'Fobs' in mapData['MapType']:
548                    F = np.sqrt(Fosq)
549                    h,k,l = hkl+Hmax
550                    Fhkl[h,k,l] = F*phasep
551                    h,k,l = -hkl+Hmax
552                    Fhkl[h,k,l] = F*phasem
553                elif 'Fcalc' in mapData['MapType']:
554                    F = np.sqrt(Fcsq)
555                    h,k,l = hkl+Hmax
556                    Fhkl[h,k,l] = F*phasep
557                    h,k,l = -hkl+Hmax
558                    Fhkl[h,k,l] = F*phasem
559                elif 'delt-F' in mapData['MapType']:
560                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
561                    h,k,l = hkl+Hmax
562                    Fhkl[h,k,l] = dF*phasep
563                    h,k,l = -hkl+Hmax
564                    Fhkl[h,k,l] = dF*phasem
565                elif '2*Fo-Fc' in mapData['MapType']:
566                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
567                    h,k,l = hkl+Hmax
568                    Fhkl[h,k,l] = F*phasep
569                    h,k,l = -hkl+Hmax
570                    Fhkl[h,k,l] = F*phasem
571                elif 'Patterson' in mapData['MapType']:
572                    h,k,l = hkl+Hmax
573                    Fhkl[h,k,l] = complex(Fosq,0.)
574                    h,k,l = -hkl+Hmax
575                    Fhkl[h,k,l] = complex(Fosq,0.)
576    Fhkl = fft.fftshift(Fhkl)
577    rho = fft.fftn(Fhkl)/cell[6]
578    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
579    mapData['rho'] = np.real(rho)
580    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
581    return mapData
582   
583def findRoll(SGData,rho,Fhkl):
584    mapShape = rho.shape
585    hklShape = Fhkl.shape
586    mapHalf = np.array(mapShape)/2
587    hklHalf = np.array(hklShape)/2
588    mapMax = np.unravel_index(np.argmax(rho),mapShape)
589    print mapMax,rho[mapMax],mapHalf
590    hklMax = np.unravel_index(np.argmax(Fhkl),hklShape)
591    print hklMax,Fhkl[hklMax],hklMax-hklHalf,hklHalf
592   
593    return [0,0,0]
594
595def ChargeFlip(data,reflData,pgbar):
596#    import scipy.fftpack as fft
597    import numpy.fft as fft
598    generalData = data['General']
599    mapData = generalData['Map']
600    flipData = generalData['Flip']
601    FFtable = {}
602    if 'None' not in flipData['Norm element']:
603        normElem = flipData['Norm element'].upper()
604        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
605        for ff in FFs:
606            if ff['Symbol'] == normElem:
607                FFtable.update(ff)
608    dmin = flipData['Resolution']
609    SGData = generalData['SGData']
610    cell = generalData['Cell'][1:8]       
611    A = G2lat.cell2A(cell[:6])
612    Vol = cell[6]
613    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
614    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
615    time0 = time.time()
616    for ref in reflData:
617        dsp = ref[4]
618        if dsp >= dmin:
619            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
620            if FFtable:
621                SQ = 0.25/dsp**2
622                ff *= G2el.ScatFac(FFtable,SQ)[0]
623            E = np.sqrt(ref[8])/ff
624            ph = ref[10]
625            if SGData['SGInv']:
626                ph = rn.randint(0,1)*180.
627            else:
628                ph = rn.uniform(0.,360.)
629            for i,hkl in enumerate(ref[11]):
630                hkl = np.asarray(hkl,dtype='i')
631                dp = 360.*ref[12][i]
632                a = cosd(ph+dp)
633                b = sind(ph+dp)
634                phasep = complex(a,b)
635                phasem = complex(a,-b)
636                h,k,l = hkl+Hmax
637                Ehkl[h,k,l] = E*phasep
638                h,k,l = -hkl+Hmax       #Friedel pair refl.
639                Ehkl[h,k,l] = E*phasem
640#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
641    CEhkl = copy.copy(Ehkl)
642    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
643    Emask = ma.getmask(MEhkl)
644    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
645    Ncyc = 0
646    while True:       
647        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
648        CEsig = np.std(CErho)
649        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
650        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
651        phase = CFhkl/np.absolute(CFhkl)
652        CEhkl = np.absolute(Ehkl)*phase
653        Ncyc += 1
654        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
655        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
656        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
657        if Rcf < 5.:
658            break
659        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
660        if not GoOn:
661            break
662    print 'Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
663    print 'No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')
664    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
665    roll = findRoll(SGData,CErho,CEhkl)
666    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
667    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
668    return mapData
669   
670def SearchMap(data,keepDup=False):
671   
672    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
673   
674    def noDuplicate(xyz,peaks,SGData):                  #be careful where this is used - it's slow
675        equivs = G2spc.GenAtom(xyz,SGData)
676        xyzs = [equiv[0] for equiv in equivs]
677        for x in xyzs:
678            if True in [np.allclose(x,peak,atol=0.002) for peak in peaks]:
679                return False
680        return True
681           
682    def findRoll(rhoMask,mapHalf):
683        indx = np.array(ma.nonzero(rhoMask)).T
684        rhoList = np.array([rho[i,j,k] for i,j,k in indx])
685        rhoMax = np.max(rhoList)
686        return mapHalf-indx[np.argmax(rhoList)]
687       
688    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
689        Mag,x0,y0,z0,sig = parms
690        if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
691            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(x0-rX)*(y0-rY)+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
692        else:
693            return norm*Mag*np.exp(-((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2))/(sig*res**3)
694       
695    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
696        Mag,x0,y0,z0,sig = parms
697        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
698        return M
699       
700    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
701        Mag,x0,y0,z0,sig = parms
702        dMdv = np.zeros(([5,]+list(rX.shape)))
703        delt = .01
704        for i in range(5):
705            parms[i] -= delt
706            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
707            parms[i] += 2.*delt
708            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
709            parms[i] -= delt
710            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
711        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
712        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
713        dMdv = np.reshape(dMdv,(5,rX.size))
714        Hess = np.inner(dMdv,dMdv)
715       
716        return Vec,Hess
717       
718    generalData = data['General']
719    phaseName = generalData['Name']
720    SGData = generalData['SGData']
721    cell = generalData['Cell'][1:8]       
722    A = G2lat.cell2A(cell[:6])
723    drawingData = data['Drawing']
724    peaks = []
725    mags = []
726    try:
727        mapData = generalData['Map']
728        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
729        rho = copy.copy(mapData['rho'])     #don't mess up original
730        mapHalf = np.array(rho.shape)/2
731        res = mapData['Resolution']
732        incre = 1./np.array(rho.shape)
733        step = max(1.0,1./res)+1
734        steps = np.array(3*[step,])
735    except KeyError:
736        print '**** ERROR - Fourier map not defined'
737        return peaks,mags
738    while True:
739        rhoMask = ma.array(rho,mask=(rho<contLevel))
740        if not ma.count(rhoMask):
741            break
742        rMI = findRoll(rhoMask,mapHalf)
743        rho = np.roll(np.roll(np.roll(rho,rMI[0],axis=0),rMI[1],axis=1),rMI[2],axis=2)
744        rMM = mapHalf-steps
745        rMP = mapHalf+steps+1
746        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
747        peakInt = np.sum(rhoPeak)*res**3
748        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
749        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
750        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
751            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.0001,maxcyc=100)
752        x1 = result[0]
753        if np.any(x1 < 0):
754            break
755        peak = (np.array(x1[1:4])-rMI)*incre
756        if not len(peaks):
757            peaks.append(peak)
758            mags.append(x1[0])
759        else:
760            if keepDup:
761                peaks.append(peak)
762                mags.append(x1[0])
763            elif noDuplicate(peak,peaks,SGData) and x1[0] > 0.:
764                peaks.append(peak)
765                mags.append(x1[0])
766        rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]] = peakFunc(result[0],rX,rY,rZ,rhoPeak,res,SGData['SGLaue'])
767        rho = np.roll(np.roll(np.roll(rho,-rMI[2],axis=2),-rMI[1],axis=1),-rMI[0],axis=0)
768    return np.array(peaks),np.array([mags,]).T
Note: See TracBrowser for help on using the repository browser.