source: trunk/GSASIImath.py @ 915

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

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