source: trunk/GSASIImath.py @ 859

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

still more rigid body fixes
add thermal motion models
allow copy of RB parms
import of rb models from gpx files
seq refinement results fixes

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