source: trunk/GSASIImath.py @ 860

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

more fixes to rigid body stuff

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