source: trunk/GSASIImath.py @ 855

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

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