source: trunk/GSASIImath.py @ 853

Last change on this file since 853 was 853, checked in by vondreele, 9 years ago

more rigid body stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 45.2 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2013-02-19 16:27:54 +0000 (Tue, 19 Feb 2013) $
5# $Author: vondreele $
6# $Revision: 853 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 853 2013-02-19 16:27:54Z 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: 853 $")
23import GSASIIElem as G2el
24import GSASIIlattice as G2lat
25import GSASIIspc as G2spc
26import numpy.fft as fft
27
28sind = lambda x: np.sin(x*np.pi/180.)
29cosd = lambda x: np.cos(x*np.pi/180.)
30tand = lambda x: np.tan(x*np.pi/180.)
31asind = lambda x: 180.*np.arcsin(x)/np.pi
32acosd = lambda x: 180.*np.arccos(x)/np.pi
33atand = lambda x: 180.*np.arctan(x)/np.pi
34atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
35
36def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
37   
38    """
39    Minimize the sum of squares of a set of equations.
40
41    ::
42   
43                    Nobs
44        x = arg min(sum(func(y)**2,axis=0))
45                    y=0
46
47    Parameters
48    ----------
49    func : callable
50        should take at least one (possibly length N vector) argument and
51        returns M floating point numbers.
52    x0 : ndarray
53        The starting estimate for the minimization of length N
54    Hess : callable
55        A required function or method to compute the weighted vector and Hessian for func.
56        It must be a symmetric NxN array
57    args : tuple
58        Any extra arguments to func are placed in this tuple.
59    ftol : float
60        Relative error desired in the sum of squares.
61    xtol : float
62        Relative error desired in the approximate solution.
63    maxcyc : int
64        The maximum number of cycles of refinement to execute, if -1 refine
65        until other limits are met (ftol, xtol)
66
67    Returns
68    -------
69    x : ndarray
70        The solution (or the result of the last iteration for an unsuccessful
71        call).
72    cov_x : ndarray
73        Uses the fjac and ipvt optional outputs to construct an
74        estimate of the jacobian around the solution.  ``None`` if a
75        singular matrix encountered (indicates very flat curvature in
76        some direction).  This matrix must be multiplied by the
77        residual standard deviation to get the covariance of the
78        parameter estimates -- see curve_fit.
79    infodict : dict
80        a dictionary of optional outputs with the key s::
81
82            - 'fvec' : the function evaluated at the output
83
84
85    Notes
86    -----
87
88    """
89               
90    x0 = np.array(x0, ndmin=1)      #might be redundant?
91    n = len(x0)
92    if type(args) != type(()):
93        args = (args,)
94       
95    icycle = 0
96    One = np.ones((n,n))
97    lam = 0.001
98    lamMax = lam
99    nfev = 0
100    while icycle < maxcyc:
101        lamMax = max(lamMax,lam)
102        M = func(x0,*args)
103        nfev += 1
104        chisq0 = np.sum(M**2)
105        Yvec,Amat = Hess(x0,*args)
106        Adiag = np.sqrt(np.diag(Amat))
107        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
108        if np.any(psing):                #hard singularity in matrix
109            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
110        Anorm = np.outer(Adiag,Adiag)
111        Yvec /= Adiag
112        Amat /= Anorm       
113        while True:
114            Lam = np.eye(Amat.shape[0])*lam
115            Amatlam = Amat*(One+Lam)
116            try:
117                Xvec = nl.solve(Amatlam,Yvec)
118            except nl.LinAlgError:
119                print 'ouch #1'
120                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
121                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
122            Xvec /= Adiag
123            M2 = func(x0+Xvec,*args)
124            nfev += 1
125            chisq1 = np.sum(M2**2)
126            if chisq1 > chisq0:
127                lam *= 10.
128            else:
129                x0 += Xvec
130                lam /= 10.
131                break
132        if (chisq0-chisq1)/chisq0 < ftol:
133            break
134        icycle += 1
135    M = func(x0,*args)
136    nfev += 1
137    Yvec,Amat = Hess(x0,*args)
138    try:
139        Bmat = nl.inv(Amat)
140        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
141    except nl.LinAlgError:
142        print 'ouch #2 linear algebra error in LS'
143        psing = []
144        if maxcyc:
145            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
146        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
147   
148def getVCov(varyNames,varyList,covMatrix):
149    vcov = np.zeros((len(varyNames),len(varyNames)))
150    for i1,name1 in enumerate(varyNames):
151        for i2,name2 in enumerate(varyNames):
152            try:
153                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
154            except ValueError:
155                vcov[i1][i2] = 0.0
156    return vcov
157
158def FindAtomIndexByIDs(atomData,IDs,Draw=True):
159    indx = []
160    for i,atom in enumerate(atomData):
161        if Draw and atom[-3] in IDs:
162            indx.append(i)
163        elif atom[-1] in IDs:
164            indx.append(i)
165    return indx
166
167def FillAtomLookUp(atomData):
168    atomLookUp = {}
169    for iatm,atom in enumerate(atomData):
170        atomLookUp[atom[-1]] = iatm
171    return atomLookUp
172
173def GetAtomsById(atomData,atomLookUp,IdList):
174    atoms = []
175    for id in IdList:
176        atoms.append(atomData[atomLookUp[id]])
177    return atoms
178   
179def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
180    Items = []
181    if not isinstance(IdList,list):
182        IdList = [IdList,]
183    for id in IdList:
184        if numItems == 1:
185            Items.append(atomData[atomLookUp[id]][itemLoc])
186        else:
187            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
188    return Items
189   
190def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
191    pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
192    dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
193    XYZ = []
194    for ind in indx:
195        names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
196        dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
197        XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
198    return XYZ
199
200def UpdateResRBAtoms(Bmat,RBObj,RBData):
201    RBIds = GetResRBIds(RBData)
202    RBRes = RBData[RBIds[RBObj['RBname']]]
203    XYZ = np.array(RBRes['rbXYZ'])
204    for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
205        QuatA = AVdeg2Q(tor[0],XYZ[seq[0]]-XYZ[seq[1]])
206        for ride in seq[3]:
207            XYZ[ride] = prodQVQ(QuatA,XYZ[ride]-XYZ[seq[1]])+XYZ[seq[1]]
208    for i,xyz in enumerate(XYZ):
209        xyz = prodQVQ(RBObj['Orient'][0],xyz)
210        xyz = np.inner(Bmat,xyz)
211        xyz += RBObj['Orig'][0]
212        XYZ[i] = xyz
213    return XYZ
214   
215def UpdateRBAtoms(Bmat,RBObj,RBData,RBType):
216    RBIds = GetResRBIds(RBData[RBType])
217    RBRes = RBData[RBType][RBIds[RBObj['RBname']]]
218    if RBType == 'Vector':
219        vecs = RBRes['rbVect']
220        mags = RBRes['VectMag']
221        XYZ = np.zeros_like(vecs[0])
222        for vec,mag in zip(vecs,mags):
223            XYZ += vec*mag
224    elif RBType == 'Residue':
225        XYZ = np.array(RBRes['rbXYZ'])
226        for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
227            QuatA = AVdeg2Q(tor[0],XYZ[seq[0]]-XYZ[seq[1]])
228            for ride in seq[3]:
229                XYZ[ride] = prodQVQ(QuatA,XYZ[ride]-XYZ[seq[1]])+XYZ[seq[1]]
230    for i,xyz in enumerate(XYZ):
231        xyz = prodQVQ(RBObj['Orient'][0],xyz)
232        xyz = np.inner(Bmat,xyz)
233        xyz += RBObj['Orig'][0]
234        XYZ[i] = xyz
235    return XYZ   
236   
237def GetResRBIds(RBData):   
238    rbKeys = RBData.keys()
239    rbKeys.remove('AtInfo')
240    if not len(rbKeys):
241        return {}
242    RBNames = [RBData[k]['RBname'] for k in rbKeys]
243    return dict(zip(RBNames,rbKeys))
244   
245def GetSHCoeff(pId,parmDict,SHkeys):
246    SHCoeff = {}
247    for shkey in SHkeys:
248        shname = str(pId)+'::'+shkey
249        SHCoeff[shkey] = parmDict[shname]
250    return SHCoeff
251       
252def getMass(generalData):
253    mass = 0.
254    for i,elem in enumerate(generalData['AtomTypes']):
255        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
256    return mass   
257
258def getDensity(generalData):
259   
260    mass = getMass(generalData)
261    Volume = generalData['Cell'][7]
262    density = mass/(0.6022137*Volume)
263    return density,Volume/mass
264   
265def getSyXYZ(XYZ,ops,SGData):
266    XYZout = np.zeros_like(XYZ)
267    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
268        if op == '1':
269            XYZout[i] = xyz
270        else:
271            oprs = op.split('+')
272            unit = [0,0,0]
273            if oprs[1]:
274                unit = np.array(list(eval(oprs[1])))
275            syop =int(oprs[0])
276            inv = syop/abs(syop)
277            syop *= inv
278            cent = syop/100
279            syop %= 100
280            syop -= 1
281            M,T = SGData['SGOps'][syop]
282            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
283    return XYZout
284   
285def getRestDist(XYZ,Amat):
286    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
287   
288def getRestDeriv(Func,XYZ,Amat,ops,SGData):
289    deriv = np.zeros((len(XYZ),3))
290    dx = 0.00001
291    for j,xyz in enumerate(XYZ):
292        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
293            XYZ[j] += x
294            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
295            XYZ[j] -= 2*x
296            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
297            XYZ[j] += x
298            deriv[j][i] = (d1-d2)/(2*dx)
299    return deriv.flatten()
300
301def getRestAngle(XYZ,Amat):
302   
303    def calcVec(Ox,Tx,Amat):
304        return np.inner(Amat,(Tx-Ox))
305
306    VecA = calcVec(XYZ[1],XYZ[0],Amat)
307    VecA /= np.sqrt(np.sum(VecA**2))
308    VecB = calcVec(XYZ[1],XYZ[2],Amat)
309    VecB /= np.sqrt(np.sum(VecB**2))
310    edge = VecB-VecA
311    edge = np.sum(edge**2)
312    angle = (2.-edge)/2.
313    angle = max(angle,-1.)
314    return acosd(angle)
315   
316def getRestPlane(XYZ,Amat):
317    sumXYZ = np.zeros(3)
318    for xyz in XYZ:
319        sumXYZ += xyz
320    sumXYZ /= len(XYZ)
321    XYZ = np.array(XYZ)-sumXYZ
322    XYZ = np.inner(Amat,XYZ).T
323    Zmat = np.zeros((3,3))
324    for i,xyz in enumerate(XYZ):
325        Zmat += np.outer(xyz.T,xyz)
326    Evec,Emat = nl.eig(Zmat)
327    Evec = np.sqrt(Evec)/(len(XYZ)-3)
328    Order = np.argsort(Evec)
329    return Evec[Order[0]]
330   
331def getRestChiral(XYZ,Amat):   
332    VecA = np.empty((3,3))   
333    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
334    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
335    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
336    return nl.det(VecA)
337   
338def getRestTorsion(XYZ,Amat):
339    VecA = np.empty((3,3))
340    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
341    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
342    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
343    D = nl.det(VecA)
344    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
345    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
346    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
347    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
348    Ang = 1.0
349    if abs(P12) < 1.0 and abs(P23) < 1.0:
350        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
351    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
352    return TOR
353   
354def calcTorsionEnergy(TOR,Coeff=[]):
355    sum = 0.
356    if len(Coeff):
357        cof = np.reshape(Coeff,(3,3)).T
358        delt = TOR-cof[1]
359        delt = np.where(delt<-180.,delt+360.,delt)
360        delt = np.where(delt>180.,delt-360.,delt)
361        term = -cof[2]*delt**2
362        val = cof[0]*np.exp(term/1000.0)
363        pMax = cof[0][np.argmin(val)]
364        Eval = np.sum(val)
365        sum = Eval-pMax
366    return sum,Eval
367
368def getTorsionDeriv(XYZ,Amat,Coeff):
369    deriv = np.zeros((len(XYZ),3))
370    dx = 0.00001
371    for j,xyz in enumerate(XYZ):
372        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
373            XYZ[j] += x
374            tor = getRestTorsion(XYZ,Amat)
375            p,d1 = calcTorsionEnergy(tor,Coeff)
376            XYZ[j] -= 2*x
377            tor = getRestTorsion(XYZ,Amat)
378            p,d2 = calcTorsionEnergy(tor,Coeff)           
379            XYZ[j] += x
380            deriv[j][i] = (d1-d2)/(2*dx)
381    return deriv.flatten()
382
383def getRestRama(XYZ,Amat):
384    phi = getRestTorsion(XYZ[:5],Amat)
385    psi = getRestTorsion(XYZ[1:],Amat)
386    return phi,psi
387   
388def calcRamaEnergy(phi,psi,Coeff=[]):
389    sum = 0.
390    if len(Coeff):
391        cof = Coeff.T
392        dPhi = phi-cof[1]
393        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
394        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
395        dPsi = psi-cof[2]
396        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
397        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
398        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
399        val = cof[0]*np.exp(val/1000.)
400        pMax = cof[0][np.argmin(val)]
401        Eval = np.sum(val)
402        sum = Eval-pMax
403    return sum,Eval
404
405def getRamaDeriv(XYZ,Amat,Coeff):
406    deriv = np.zeros((len(XYZ),3))
407    dx = 0.00001
408    for j,xyz in enumerate(XYZ):
409        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
410            XYZ[j] += x
411            phi,psi = getRestRama(XYZ,Amat)
412            p,d1 = calcRamaEnergy(phi,psi,Coeff)
413            XYZ[j] -= 2*x
414            phi,psi = getRestRama(XYZ,Amat)
415            p,d2 = calcRamaEnergy(phi,psi,Coeff)
416            XYZ[j] += x
417            deriv[j][i] = (d1-d2)/(2*dx)
418    return deriv.flatten()
419
420def getRestPolefig(ODFln,SamSym,Grid):
421    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
422    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
423    R = np.where(R <= 1.,2.*atand(R),0.0)
424    Z = np.zeros_like(R)
425    Z = G2lat.polfcal(ODFln,SamSym,R,P)
426    Z = np.reshape(Z,(Grid,Grid))
427    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
428
429def getRestPolefigDerv(HKL,Grid,SHCoeff):
430    pass
431       
432def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
433   
434    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
435        TxT = inv*(np.inner(M,Tx)+T)+C+U
436        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
437       
438    inv = Top/abs(Top)
439    cent = abs(Top)/100
440    op = abs(Top)%100-1
441    M,T = SGData['SGOps'][op]
442    C = SGData['SGCen'][cent]
443    dx = .00001
444    deriv = np.zeros(6)
445    for i in [0,1,2]:
446        Oxyz[i] += dx
447        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
448        Oxyz[i] -= 2*dx
449        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
450        Oxyz[i] += dx
451        Txyz[i] += dx
452        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
453        Txyz[i] -= 2*dx
454        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
455        Txyz[i] += dx
456    return deriv
457   
458def getAngSig(VA,VB,Amat,SGData,covData={}):
459   
460    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
461        TxT = inv*(np.inner(M,Tx)+T)+C
462        TxT = G2spc.MoveToUnitCell(TxT)+U
463        return np.inner(Amat,(TxT-Ox))
464       
465    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
466        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
467        VecA /= np.sqrt(np.sum(VecA**2))
468        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
469        VecB /= np.sqrt(np.sum(VecB**2))
470        edge = VecB-VecA
471        edge = np.sum(edge**2)
472        angle = (2.-edge)/2.
473        angle = max(angle,-1.)
474        return acosd(angle)
475       
476    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
477    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
478    invA = invB = 1
479    invA = TopA/abs(TopA)
480    invB = TopB/abs(TopB)
481    centA = abs(TopA)/100
482    centB = abs(TopB)/100
483    opA = abs(TopA)%100-1
484    opB = abs(TopB)%100-1
485    MA,TA = SGData['SGOps'][opA]
486    MB,TB = SGData['SGOps'][opB]
487    CA = SGData['SGCen'][centA]
488    CB = SGData['SGCen'][centB]
489    if 'covMatrix' in covData:
490        covMatrix = covData['covMatrix']
491        varyList = covData['varyList']
492        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
493        dx = .00001
494        dadx = np.zeros(9)
495        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
496        for i in [0,1,2]:
497            OxA[i] += dx
498            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
499            OxA[i] -= 2*dx
500            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
501            OxA[i] += dx
502           
503            TxA[i] += dx
504            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
505            TxA[i] -= 2*dx
506            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
507            TxA[i] += dx
508           
509            TxB[i] += dx
510            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
511            TxB[i] -= 2*dx
512            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
513            TxB[i] += dx
514           
515        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
516        if sigAng < 0.01:
517            sigAng = 0.0
518        return Ang,sigAng
519    else:
520        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
521
522def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
523
524    def calcDist(Atoms,SyOps,Amat):
525        XYZ = []
526        for i,atom in enumerate(Atoms):
527            Inv,M,T,C,U = SyOps[i]
528            XYZ.append(np.array(atom[1:4]))
529            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
530            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
531        V1 = XYZ[1]-XYZ[0]
532        return np.sqrt(np.sum(V1**2))
533       
534    Inv = []
535    SyOps = []
536    names = []
537    for i,atom in enumerate(Oatoms):
538        names += atom[-1]
539        Op,unit = Atoms[i][-1]
540        inv = Op/abs(Op)
541        m,t = SGData['SGOps'][abs(Op)%100-1]
542        c = SGData['SGCen'][abs(Op)/100]
543        SyOps.append([inv,m,t,c,unit])
544    Dist = calcDist(Oatoms,SyOps,Amat)
545   
546    sig = -0.001
547    if 'covMatrix' in covData:
548        parmNames = []
549        dx = .00001
550        dadx = np.zeros(6)
551        for i in range(6):
552            ia = i/3
553            ix = i%3
554            Oatoms[ia][ix+1] += dx
555            a0 = calcDist(Oatoms,SyOps,Amat)
556            Oatoms[ia][ix+1] -= 2*dx
557            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
558        covMatrix = covData['covMatrix']
559        varyList = covData['varyList']
560        DistVcov = getVCov(names,varyList,covMatrix)
561        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
562        if sig < 0.001:
563            sig = -0.001
564   
565    return Dist,sig
566
567def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
568       
569    def calcAngle(Atoms,SyOps,Amat):
570        XYZ = []
571        for i,atom in enumerate(Atoms):
572            Inv,M,T,C,U = SyOps[i]
573            XYZ.append(np.array(atom[1:4]))
574            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
575            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
576        V1 = XYZ[1]-XYZ[0]
577        V1 /= np.sqrt(np.sum(V1**2))
578        V2 = XYZ[1]-XYZ[2]
579        V2 /= np.sqrt(np.sum(V2**2))
580        V3 = V2-V1
581        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
582        return acosd(cang)
583
584    Inv = []
585    SyOps = []
586    names = []
587    for i,atom in enumerate(Oatoms):
588        names += atom[-1]
589        Op,unit = Atoms[i][-1]
590        inv = Op/abs(Op)
591        m,t = SGData['SGOps'][abs(Op)%100-1]
592        c = SGData['SGCen'][abs(Op)/100]
593        SyOps.append([inv,m,t,c,unit])
594    Angle = calcAngle(Oatoms,SyOps,Amat)
595   
596    sig = -0.01
597    if 'covMatrix' in covData:
598        parmNames = []
599        dx = .00001
600        dadx = np.zeros(9)
601        for i in range(9):
602            ia = i/3
603            ix = i%3
604            Oatoms[ia][ix+1] += dx
605            a0 = calcAngle(Oatoms,SyOps,Amat)
606            Oatoms[ia][ix+1] -= 2*dx
607            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
608        covMatrix = covData['covMatrix']
609        varyList = covData['varyList']
610        AngVcov = getVCov(names,varyList,covMatrix)
611        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
612        if sig < 0.01:
613            sig = -0.01
614   
615    return Angle,sig
616
617def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
618   
619    def calcTorsion(Atoms,SyOps,Amat):
620       
621        XYZ = []
622        for i,atom in enumerate(Atoms):
623            Inv,M,T,C,U = SyOps[i]
624            XYZ.append(np.array(atom[1:4]))
625            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
626            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
627        V1 = XYZ[1]-XYZ[0]
628        V2 = XYZ[2]-XYZ[1]
629        V3 = XYZ[3]-XYZ[2]
630        V1 /= np.sqrt(np.sum(V1**2))
631        V2 /= np.sqrt(np.sum(V2**2))
632        V3 /= np.sqrt(np.sum(V3**2))
633        M = np.array([V1,V2,V3])
634        D = nl.det(M)
635        Ang = 1.0
636        P12 = np.dot(V1,V2)
637        P13 = np.dot(V1,V3)
638        P23 = np.dot(V2,V3)
639        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
640        return Tors
641           
642    Inv = []
643    SyOps = []
644    names = []
645    for i,atom in enumerate(Oatoms):
646        names += atom[-1]
647        Op,unit = Atoms[i][-1]
648        inv = Op/abs(Op)
649        m,t = SGData['SGOps'][abs(Op)%100-1]
650        c = SGData['SGCen'][abs(Op)/100]
651        SyOps.append([inv,m,t,c,unit])
652    Tors = calcTorsion(Oatoms,SyOps,Amat)
653   
654    sig = -0.01
655    if 'covMatrix' in covData:
656        parmNames = []
657        dx = .00001
658        dadx = np.zeros(12)
659        for i in range(12):
660            ia = i/3
661            ix = i%3
662            Oatoms[ia][ix+1] += dx
663            a0 = calcTorsion(Oatoms,SyOps,Amat)
664            Oatoms[ia][ix+1] -= 2*dx
665            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
666        covMatrix = covData['covMatrix']
667        varyList = covData['varyList']
668        TorVcov = getVCov(names,varyList,covMatrix)
669        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
670        if sig < 0.01:
671            sig = -0.01
672   
673    return Tors,sig
674       
675def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
676   
677    def calcDist(Atoms,SyOps,Amat):
678        XYZ = []
679        for i,atom in enumerate(Atoms):
680            Inv,M,T,C,U = SyOps[i]
681            XYZ.append(np.array(atom[1:4]))
682            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
683            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
684        V1 = XYZ[1]-XYZ[0]
685        return np.sqrt(np.sum(V1**2))
686       
687    def calcAngle(Atoms,SyOps,Amat):
688        XYZ = []
689        for i,atom in enumerate(Atoms):
690            Inv,M,T,C,U = SyOps[i]
691            XYZ.append(np.array(atom[1:4]))
692            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
693            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
694        V1 = XYZ[1]-XYZ[0]
695        V1 /= np.sqrt(np.sum(V1**2))
696        V2 = XYZ[1]-XYZ[2]
697        V2 /= np.sqrt(np.sum(V2**2))
698        V3 = V2-V1
699        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
700        return acosd(cang)
701
702    def calcTorsion(Atoms,SyOps,Amat):
703       
704        XYZ = []
705        for i,atom in enumerate(Atoms):
706            Inv,M,T,C,U = SyOps[i]
707            XYZ.append(np.array(atom[1:4]))
708            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
709            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
710        V1 = XYZ[1]-XYZ[0]
711        V2 = XYZ[2]-XYZ[1]
712        V3 = XYZ[3]-XYZ[2]
713        V1 /= np.sqrt(np.sum(V1**2))
714        V2 /= np.sqrt(np.sum(V2**2))
715        V3 /= np.sqrt(np.sum(V3**2))
716        M = np.array([V1,V2,V3])
717        D = nl.det(M)
718        Ang = 1.0
719        P12 = np.dot(V1,V2)
720        P13 = np.dot(V1,V3)
721        P23 = np.dot(V2,V3)
722        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
723        return Tors
724           
725    Inv = []
726    SyOps = []
727    names = []
728    for i,atom in enumerate(Oatoms):
729        names += atom[-1]
730        Op,unit = Atoms[i][-1]
731        inv = Op/abs(Op)
732        m,t = SGData['SGOps'][abs(Op)%100-1]
733        c = SGData['SGCen'][abs(Op)/100]
734        SyOps.append([inv,m,t,c,unit])
735    M = len(Oatoms)
736    if M == 2:
737        Val = calcDist(Oatoms,SyOps,Amat)
738    elif M == 3:
739        Val = calcAngle(Oatoms,SyOps,Amat)
740    else:
741        Val = calcTorsion(Oatoms,SyOps,Amat)
742   
743    sigVals = [-0.001,-0.01,-0.01]
744    sig = sigVals[M-3]
745    if 'covMatrix' in covData:
746        parmNames = []
747        dx = .00001
748        N = M*3
749        dadx = np.zeros(N)
750        for i in range(N):
751            ia = i/3
752            ix = i%3
753            Oatoms[ia][ix+1] += dx
754            if M == 2:
755                a0 = calcDist(Oatoms,SyOps,Amat)
756            elif M == 3:
757                a0 = calcAngle(Oatoms,SyOps,Amat)
758            else:
759                a0 = calcTorsion(Oatoms,SyOps,Amat)
760            Oatoms[ia][ix+1] -= 2*dx
761            if M == 2:
762                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
763            elif M == 3:
764                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
765            else:
766                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
767        covMatrix = covData['covMatrix']
768        varyList = covData['varyList']
769        Vcov = getVCov(names,varyList,covMatrix)
770        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
771        if sig < sigVals[M-3]:
772            sig = sigVals[M-3]
773   
774    return Val,sig
775       
776   
777def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
778    # returns value(esd) string; nTZ=True for no trailing zeros
779    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
780    #get the 2 significant digits in the esd
781    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
782    #get the number of digits to represent them
783    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
784   
785    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
786    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
787    if esd > 0:
788        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
789        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
790    elif esd < 0:
791         return str(round(value,mdec(esd)-1))
792    else:
793        text = str("%f"%(value))
794        if nTZ:
795            return text.rstrip('0')
796        else:
797            return text
798           
799def adjHKLmax(SGData,Hmax):
800    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
801        Hmax[0] = ((Hmax[0]+3)/6)*6
802        Hmax[1] = ((Hmax[1]+3)/6)*6
803        Hmax[2] = ((Hmax[2]+1)/4)*4
804    else:
805        Hmax[0] = ((Hmax[0]+2)/4)*4
806        Hmax[1] = ((Hmax[1]+2)/4)*4
807        Hmax[2] = ((Hmax[2]+1)/4)*4
808
809def FourierMap(data,reflData):
810   
811    generalData = data['General']
812    if not generalData['Map']['MapType']:
813        print '**** ERROR - Fourier map not defined'
814        return
815    mapData = generalData['Map']
816    dmin = mapData['Resolution']
817    SGData = generalData['SGData']
818    cell = generalData['Cell'][1:8]       
819    A = G2lat.cell2A(cell[:6])
820    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
821    adjHKLmax(SGData,Hmax)
822    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
823#    Fhkl[0,0,0] = generalData['F000X']
824    time0 = time.time()
825    for ref in reflData:
826        if ref[4] >= dmin:
827            Fosq,Fcsq,ph = ref[8:11]
828            for i,hkl in enumerate(ref[11]):
829                hkl = np.asarray(hkl,dtype='i')
830                dp = 360.*ref[12][i]
831                a = cosd(ph+dp)
832                b = sind(ph+dp)
833                phasep = complex(a,b)
834                phasem = complex(a,-b)
835                if 'Fobs' in mapData['MapType']:
836                    F = np.sqrt(Fosq)
837                    h,k,l = hkl+Hmax
838                    Fhkl[h,k,l] = F*phasep
839                    h,k,l = -hkl+Hmax
840                    Fhkl[h,k,l] = F*phasem
841                elif 'Fcalc' in mapData['MapType']:
842                    F = np.sqrt(Fcsq)
843                    h,k,l = hkl+Hmax
844                    Fhkl[h,k,l] = F*phasep
845                    h,k,l = -hkl+Hmax
846                    Fhkl[h,k,l] = F*phasem
847                elif 'delt-F' in mapData['MapType']:
848                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
849                    h,k,l = hkl+Hmax
850                    Fhkl[h,k,l] = dF*phasep
851                    h,k,l = -hkl+Hmax
852                    Fhkl[h,k,l] = dF*phasem
853                elif '2*Fo-Fc' in mapData['MapType']:
854                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
855                    h,k,l = hkl+Hmax
856                    Fhkl[h,k,l] = F*phasep
857                    h,k,l = -hkl+Hmax
858                    Fhkl[h,k,l] = F*phasem
859                elif 'Patterson' in mapData['MapType']:
860                    h,k,l = hkl+Hmax
861                    Fhkl[h,k,l] = complex(Fosq,0.)
862                    h,k,l = -hkl+Hmax
863                    Fhkl[h,k,l] = complex(Fosq,0.)
864    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
865    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
866    mapData['rho'] = np.real(rho)
867    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
868    return mapData
869   
870# map printing for testing purposes
871def printRho(SGLaue,rho,rhoMax):                         
872    dim = len(rho.shape)
873    if dim == 2:
874        ix,jy = rho.shape
875        for j in range(jy):
876            line = ''
877            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
878                line += (jy-j)*'  '
879            for i in range(ix):
880                r = int(100*rho[i,j]/rhoMax)
881                line += '%4d'%(r)
882            print line+'\n'
883    else:
884        ix,jy,kz = rho.shape
885        for k in range(kz):
886            print 'k = ',k
887            for j in range(jy):
888                line = ''
889                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
890                    line += (jy-j)*'  '
891                for i in range(ix):
892                    r = int(100*rho[i,j,k]/rhoMax)
893                    line += '%4d'%(r)
894                print line+'\n'
895## keep this
896               
897def findOffset(SGData,A,Fhkl):   
898    if SGData['SpGrp'] == 'P 1':
899        return [0,0,0]   
900    hklShape = Fhkl.shape
901    hklHalf = np.array(hklShape)/2
902    sortHKL = np.argsort(Fhkl.flatten())
903    Fdict = {}
904    for hkl in sortHKL:
905        HKL = np.unravel_index(hkl,hklShape)
906        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
907        if F == 0.:
908            break
909        Fdict['%.6f'%(np.absolute(F))] = hkl
910    Flist = np.flipud(np.sort(Fdict.keys()))
911    F = str(1.e6)
912    i = 0
913    DH = []
914    Dphi = []
915    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
916    for F in Flist:
917        hkl = np.unravel_index(Fdict[F],hklShape)
918        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
919        Uniq = np.array(Uniq,dtype='i')
920        Phi = np.array(Phi)
921        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
922        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
923        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
924        ang0 = np.angle(Fh0,deg=True)/360.
925        for H,phi in zip(Uniq,Phi)[1:]:
926            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
927            dH = H-hkl
928            dang = ang-ang0
929            if np.any(np.abs(dH)-Hmax > 0):    #keep low order DHs
930                continue
931            DH.append(dH)
932            Dphi.append((dang+.5) % 1.0)
933        if i > 20 or len(DH) > 30:
934            break
935        i += 1
936    DH = np.array(DH)
937    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
938    Dphi = np.array(Dphi)
939    steps = np.array(hklShape)
940    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
941    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
942    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
943    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
944    hist,bins = np.histogram(Mmap,bins=1000)
945#    for i,item in enumerate(hist[:10]):
946#        print item,bins[i]
947    chisq = np.min(Mmap)
948    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
949    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
950#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
951    return DX
952   
953def ChargeFlip(data,reflData,pgbar):
954    generalData = data['General']
955    mapData = generalData['Map']
956    flipData = generalData['Flip']
957    FFtable = {}
958    if 'None' not in flipData['Norm element']:
959        normElem = flipData['Norm element'].upper()
960        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
961        for ff in FFs:
962            if ff['Symbol'] == normElem:
963                FFtable.update(ff)
964    dmin = flipData['Resolution']
965    SGData = generalData['SGData']
966    cell = generalData['Cell'][1:8]       
967    A = G2lat.cell2A(cell[:6])
968    Vol = cell[6]
969    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
970    adjHKLmax(SGData,Hmax)
971    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
972    time0 = time.time()
973    for ref in reflData:
974        dsp = ref[4]
975        if dsp >= dmin:
976            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
977            if FFtable:
978                SQ = 0.25/dsp**2
979                ff *= G2el.ScatFac(FFtable,SQ)[0]
980            if ref[8] > 0.:         #use only +ve Fobs**2
981                E = np.sqrt(ref[8])/ff
982            else:
983                E = 0.
984            ph = ref[10]
985            ph = rn.uniform(0.,360.)
986            for i,hkl in enumerate(ref[11]):
987                hkl = np.asarray(hkl,dtype='i')
988                dp = 360.*ref[12][i]
989                a = cosd(ph+dp)
990                b = sind(ph+dp)
991                phasep = complex(a,b)
992                phasem = complex(a,-b)
993                h,k,l = hkl+Hmax
994                Ehkl[h,k,l] = E*phasep
995                h,k,l = -hkl+Hmax       #Friedel pair refl.
996                Ehkl[h,k,l] = E*phasem
997#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
998    CEhkl = copy.copy(Ehkl)
999    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1000    Emask = ma.getmask(MEhkl)
1001    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1002    Ncyc = 0
1003    old = np.seterr(all='raise')
1004    while True:       
1005        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1006        CEsig = np.std(CErho)
1007        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1008        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem! make 20. adjustible
1009        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1010        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
1011        phase = CFhkl/np.absolute(CFhkl)
1012        CEhkl = np.absolute(Ehkl)*phase
1013        Ncyc += 1
1014        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1015        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1016        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1017        if Rcf < 5.:
1018            break
1019        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
1020        if not GoOn or Ncyc > 10000:
1021            break
1022    np.seterr(**old)
1023    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
1024    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
1025    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
1026    roll = findOffset(SGData,A,CEhkl)
1027       
1028    mapData['Rcf'] = Rcf
1029    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1030    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1031    mapData['rollMap'] = [0,0,0]
1032    return mapData
1033   
1034def SearchMap(data):
1035    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1036   
1037    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
1038   
1039    def noDuplicate(xyz,peaks,Amat):
1040        XYZ = np.inner(Amat,xyz)
1041        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1042            print ' Peak',xyz,' <0.5A from another peak'
1043            return False
1044        return True
1045                           
1046    def fixSpecialPos(xyz,SGData,Amat):
1047        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
1048        X = []
1049        xyzs = [equiv[0] for equiv in equivs]
1050        for x in xyzs:
1051            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
1052                X.append(x)
1053        if len(X) > 1:
1054            return np.average(X,axis=0)
1055        else:
1056            return xyz
1057       
1058    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
1059        Mag,x0,y0,z0,sig = parms
1060        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
1061#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
1062        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
1063       
1064    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
1065        Mag,x0,y0,z0,sig = parms
1066        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1067        return M
1068       
1069    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
1070        Mag,x0,y0,z0,sig = parms
1071        dMdv = np.zeros(([5,]+list(rX.shape)))
1072        delt = .01
1073        for i in range(5):
1074            parms[i] -= delt
1075            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1076            parms[i] += 2.*delt
1077            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1078            parms[i] -= delt
1079            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
1080        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1081        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
1082        dMdv = np.reshape(dMdv,(5,rX.size))
1083        Hess = np.inner(dMdv,dMdv)
1084       
1085        return Vec,Hess
1086       
1087    generalData = data['General']
1088    phaseName = generalData['Name']
1089    SGData = generalData['SGData']
1090    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1091    drawingData = data['Drawing']
1092    peaks = []
1093    mags = []
1094    dzeros = []
1095    try:
1096        mapData = generalData['Map']
1097        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
1098        rho = copy.copy(mapData['rho'])     #don't mess up original
1099        mapHalf = np.array(rho.shape)/2
1100        res = mapData['Resolution']
1101        incre = np.array(rho.shape,dtype=np.float)
1102        step = max(1.0,1./res)+1
1103        steps = np.array(3*[step,])
1104    except KeyError:
1105        print '**** ERROR - Fourier map not defined'
1106        return peaks,mags
1107    rhoMask = ma.array(rho,mask=(rho<contLevel))
1108    indices = (-1,0,1)
1109    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1110    for roll in rolls:
1111        if np.any(roll):
1112            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
1113    indx = np.transpose(rhoMask.nonzero())
1114    peaks = indx/incre
1115    mags = rhoMask[rhoMask.nonzero()]
1116    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
1117        rho = rollMap(rho,ind)
1118        rMM = mapHalf-steps
1119        rMP = mapHalf+steps+1
1120        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1121        peakInt = np.sum(rhoPeak)*res**3
1122        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1123        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
1124        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
1125            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
1126        x1 = result[0]
1127        if not np.any(x1 < 0):
1128            mag = x1[0]
1129            peak = (np.array(x1[1:4])-ind)/incre
1130        peak = fixSpecialPos(peak,SGData,Amat)
1131        rho = rollMap(rho,-ind)       
1132    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
1133    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
1134   
1135def sortArray(data,pos,reverse=False):
1136    #data is a list of items
1137    #sort by pos in list; reverse if True
1138    T = []
1139    for i,M in enumerate(data):
1140        T.append((M[pos],i))
1141    D = dict(zip(T,data))
1142    T.sort()
1143    if reverse:
1144        T.reverse()
1145    X = []
1146    for key in T:
1147        X.append(D[key])
1148    return X
1149
1150def PeaksEquiv(data,Ind):
1151
1152    def Duplicate(xyz,peaks,Amat):
1153        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1154            return True
1155        return False
1156                           
1157    generalData = data['General']
1158    cell = generalData['Cell'][1:7]
1159    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1160    A = G2lat.cell2A(cell)
1161    SGData = generalData['SGData']
1162    mapPeaks = data['Map Peaks']
1163    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
1164    Indx = {}
1165    for ind in Ind:
1166        xyz = np.array(mapPeaks[ind][1:4])
1167        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
1168#        for x in xyzs: print x
1169        for jnd,xyz in enumerate(XYZ):       
1170            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
1171    Ind = []
1172    for ind in Indx:
1173        if Indx[ind]:
1174            Ind.append(ind)
1175    return Ind
1176               
1177def PeaksUnique(data,Ind):
1178#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
1179
1180    def noDuplicate(xyz,peaks,Amat):
1181        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1182            return False
1183        return True
1184                           
1185    generalData = data['General']
1186    cell = generalData['Cell'][1:7]
1187    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1188    A = G2lat.cell2A(cell)
1189    SGData = generalData['SGData']
1190    mapPeaks = data['Map Peaks']
1191    Indx = {}
1192    XYZ = {}
1193    for ind in Ind:
1194        XYZ[ind] = np.array(mapPeaks[ind][1:4])
1195        Indx[ind] = True
1196    for ind in Ind:
1197        if Indx[ind]:
1198            xyz = XYZ[ind]
1199            for jnd in Ind:
1200                if ind != jnd and Indx[jnd]:                       
1201                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
1202                    xyzs = np.array([equiv[0] for equiv in Equiv])
1203                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
1204    Ind = []
1205    for ind in Indx:
1206        if Indx[ind]:
1207            Ind.append(ind)
1208    return Ind
1209   
1210def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
1211    ind = 0
1212    if useFit:
1213        ind = 1
1214    ins = {}
1215    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
1216        for x in ['U','V','W','X','Y']:
1217            ins[x] = Parms[x][ind]
1218        if ifQ:                              #qplot - convert back to 2-theta
1219            pos = 2.0*asind(pos*wave/(4*math.pi))
1220        sig = ins['U']*tand(pos/2.0)**2+ins['V']*tand(pos/2.0)+ins['W']
1221        gam = ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)           
1222        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
1223    else:
1224        if ifQ:
1225            dsp = 2.*np.pi/pos
1226            pos = Parms['difC']*dsp
1227        else:
1228            dsp = pos/Parms['difC'][1]
1229        if 'Pdabc' in Parms2:
1230            for x in ['sig-0','sig-1','X','Y']:
1231                ins[x] = Parms[x][ind]
1232            Pdabc = Parms2['Pdabc'].T
1233            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1234            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1235        else:
1236            for x in ['alpha','beta-0','beta-1','sig-0','sig-1','X','Y']:
1237                ins[x] = Parms[x][ind]
1238            alp = ins['alpha']/dsp
1239            bet = ins['beta-0']+ins['beta-1']/dsp**4
1240        sig = ins['sig-0']+ins['sig-1']*dsp**2
1241        gam = ins['X']*dsp+ins['Y']*dsp**2
1242        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
1243    return XY
1244       
1245def getWave(Parms):
1246    try:
1247        return Parms['Lam'][1]
1248    except KeyError:
1249        return Parms['Lam1'][1]
1250   
1251def prodQQ(QA,QB):
1252    ''' Grassman quaternion product
1253        QA,QB quaternions; q=r+ai+bj+ck
1254    '''
1255    D = np.zeros(4)
1256    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
1257    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
1258    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
1259    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
1260    return D
1261   
1262def normQ(QA):
1263    ''' get length of quaternion & normalize it
1264        q=r+ai+bj+ck
1265    '''
1266    n = np.sqrt(np.sum(np.array(QA)**2))
1267    return QA/n
1268   
1269def invQ(Q):
1270    '''
1271        get inverse of quaternion
1272        q=r+ai+bj+ck; q* = r-ai-bj-ck
1273    '''
1274    return Q*np.array([1,-1,-1,-1])
1275   
1276def prodQVQ(Q,V):
1277    ''' compute the quaternion vector rotation qvq-1 = v'
1278        q=r+ai+bj+ck
1279    '''
1280    VP = np.zeros(3)
1281    T2 = Q[0]*Q[1]
1282    T3 = Q[0]*Q[2]
1283    T4 = Q[0]*Q[3]
1284    T5 = -Q[1]*Q[1]
1285    T6 = Q[1]*Q[2]
1286    T7 = Q[1]*Q[3]
1287    T8 = -Q[2]*Q[2]
1288    T9 = Q[2]*Q[3]
1289    T10 = -Q[3]*Q[3]
1290    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
1291    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
1292    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2] 
1293    return VP   
1294   
1295def Q2Mat(Q):
1296    ''' make rotation matrix from quaternion
1297        q=r+ai+bj+ck
1298    '''
1299    aa = Q[0]**2
1300    ab = Q[0]*Q[1]
1301    ac = Q[0]*Q[2]
1302    ad = Q[0]*Q[3]
1303    bb = Q[1]**2
1304    bc = Q[1]*Q[2]
1305    bd = Q[1]*Q[3]
1306    cc = Q[2]**2
1307    cd = Q[2]*Q[3]
1308    dd = Q[3]**2
1309    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
1310        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
1311        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
1312    return np.array(M)
1313   
1314def AV2Q(A,V):
1315    ''' convert angle (radians -pi to pi) & vector to quaternion
1316        q=r+ai+bj+ck
1317    '''
1318    Q = np.zeros(4)
1319    d = np.sqrt(np.sum(np.array(V)**2))
1320    if d:
1321        V /= d
1322    else:
1323        return [1.,0.,0.,0.]    #identity
1324    p = A/2.
1325    Q[0] = np.cos(p)
1326    Q[1:4] = V*np.sin(p)
1327    return Q
1328   
1329def AVdeg2Q(A,V):
1330    ''' convert angle (degrees -180 to 180) & vector to quaternion
1331        q=r+ai+bj+ck
1332    '''
1333    Q = np.zeros(4)
1334    d = np.sqrt(np.sum(np.array(V)**2))
1335    if d:
1336        V /= d
1337    else:
1338        return [1.,0.,0.,0.]    #identity
1339    p = A/2.
1340    Q[0] = cosd(p)
1341    Q[1:4] = V*sind(p)
1342    return Q
1343
1344def makeQuat(A,B,C):
1345    ''' Make quaternion from rotation of A vector to B vector about C axis
1346        A,B,C are np.array Cartesian 3-vectors
1347    Returns quaternion & rotation angle in radians
1348        q=r+ai+bj+ck
1349    '''
1350
1351    V1 = np.cross(A,C)
1352    V2 = np.cross(B,C)
1353    if nl.norm(V1)*nl.norm(V2):
1354        V1 /= nl.norm(V1)
1355        V2 /= nl.norm(V2)
1356        V3 = np.cross(V1,V2)
1357    else:
1358        V3 = np.zeros(3)
1359    Q = np.array([1.0,0.0,0.0,0.0])
1360    D = 0.
1361    if nl.norm(V3):
1362        V3 /= nl.norm(V3)
1363        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
1364        D = np.arccos(D1)/2.0
1365        V1 = C-V3
1366        V2 = C+V3
1367        DM = nl.norm(V1)
1368        DP = nl.norm(V2)
1369        S = np.sin(D)
1370        Q[0] = np.cos(D)
1371        Q[1:] = V3*S
1372        D *= 2.
1373        if DM > DP:
1374            D *= -1.
1375    return Q,D
1376   
Note: See TracBrowser for help on using the repository browser.