source: trunk/GSASIImath.py @ 664

Last change on this file since 664 was 664, checked in by vondreele, 11 years ago

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