source: trunk/GSASIImath.py @ 1792

Last change on this file since 1792 was 1792, checked in by vondreele, 7 years ago

remove Exceptions from seqRefine; now returns & makes Message dialog
make histograms go in order of tree items not random dict order
make sph. harm. stuff work on arrays of hkl - speed up of texture calculations
implement selected copy for Data parameters (phase fraction, etc.)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 131.2 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2015-04-17 16:50:01 +0000 (Fri, 17 Apr 2015) $
5# $Author: vondreele $
6# $Revision: 1792 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 1792 2015-04-17 16:50:01Z vondreele $
9########### SVN repository information ###################
10'''
11*GSASIImath: computation module*
12================================
13
14Routines for least-squares minimization and other stuff
15
16'''
17import sys
18import os
19import os.path as ospath
20import random as rn
21import numpy as np
22import numpy.linalg as nl
23import numpy.ma as ma
24import cPickle
25import time
26import math
27import copy
28import GSASIIpath
29GSASIIpath.SetVersionNumber("$Revision: 1792 $")
30import GSASIIElem as G2el
31import GSASIIlattice as G2lat
32import GSASIIspc as G2spc
33import GSASIIpwd as G2pwd
34import numpy.fft as fft
35import scipy.optimize as so
36import pypowder as pyd
37
38sind = lambda x: np.sin(x*np.pi/180.)
39cosd = lambda x: np.cos(x*np.pi/180.)
40tand = lambda x: np.tan(x*np.pi/180.)
41asind = lambda x: 180.*np.arcsin(x)/np.pi
42acosd = lambda x: 180.*np.arccos(x)/np.pi
43atand = lambda x: 180.*np.arctan(x)/np.pi
44atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
45twopi = 2.0*np.pi
46twopisq = 2.0*np.pi**2
47   
48################################################################################
49##### Hessian least-squares Levenberg-Marquardt routine
50################################################################################
51
52def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0,Print=False):
53   
54    """
55    Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
56    values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`
57   
58    ::
59   
60                    Nobs
61        x = arg min(sum(func(y)**2,axis=0))
62                    y=0
63
64    :param function func: callable method or function
65        should take at least one (possibly length N vector) argument and
66        returns M floating point numbers.
67    :param np.ndarray x0: The starting estimate for the minimization of length N
68    :param function Hess: callable method or function
69        A required function or method to compute the weighted vector and Hessian for func.
70        It must be a symmetric NxN array
71    :param tuple args: Any extra arguments to func are placed in this tuple.
72    :param float ftol: Relative error desired in the sum of squares.
73    :param float xtol: Relative error desired in the approximate solution.
74    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
75        until other limits are met (ftol, xtol)
76    :param bool Print: True for printing results (residuals & times) by cycle
77
78    :returns: (x,cov_x,infodict) where
79
80      * x : ndarray
81        The solution (or the result of the last iteration for an unsuccessful
82        call).
83      * cov_x : ndarray
84        Uses the fjac and ipvt optional outputs to construct an
85        estimate of the jacobian around the solution.  ``None`` if a
86        singular matrix encountered (indicates very flat curvature in
87        some direction).  This matrix must be multiplied by the
88        residual standard deviation to get the covariance of the
89        parameter estimates -- see curve_fit.
90      * infodict : dict
91        a dictionary of optional outputs with the keys:
92
93         * 'fvec' : the function evaluated at the output
94         * 'num cyc':
95         * 'nfev':
96         * 'lamMax':
97         * 'psing':
98           
99    """
100               
101    ifConverged = False
102    deltaChi2 = -10.
103    x0 = np.array(x0, ndmin=1)      #might be redundant?
104    n = len(x0)
105    if type(args) != type(()):
106        args = (args,)
107       
108    icycle = 0
109    One = np.ones((n,n))
110    lam = 0.001
111    lamMax = lam
112    nfev = 0
113    if Print:
114        print ' Hessian refinement on %d variables:'%(n)
115    Lam = np.zeros((n,n))
116    while icycle < maxcyc:
117        time0 = time.time()
118        M = func(x0,*args)
119        nfev += 1
120        chisq0 = np.sum(M**2)
121        Yvec,Amat = Hess(x0,*args)
122        Adiag = np.sqrt(np.diag(Amat))
123        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
124        if np.any(psing):                #hard singularity in matrix
125            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
126        Anorm = np.outer(Adiag,Adiag)
127        Yvec /= Adiag
128        Amat /= Anorm
129        while True:
130            Lam = np.eye(Amat.shape[0])*lam
131            Amatlam = Amat*(One+Lam)
132            try:
133                Xvec = nl.solve(Amatlam,Yvec)
134            except nl.LinAlgError:
135                print 'ouch #1'
136                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
137                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
138            Xvec /= Adiag
139            M2 = func(x0+Xvec,*args)
140            nfev += 1
141            chisq1 = np.sum(M2**2)
142            if chisq1 > chisq0:
143                lam *= 10.
144            else:
145                x0 += Xvec
146                lam /= 10.
147                break
148            if lam > 10.e5:
149                print 'ouch #3 chisq1 ',chisq1,' stuck > chisq0 ',chisq0
150                break
151        lamMax = max(lamMax,lam)
152        deltaChi2 = (chisq0-chisq1)/chisq0
153        if Print:
154            print ' Cycle: %d, Time: %.2fs, Chi**2: %.5g, Lambda: %.3g,  Delta: %.3g'%(
155                icycle,time.time()-time0,chisq1,lam,deltaChi2)
156        if deltaChi2 < ftol:
157            ifConverged = True
158            if Print: print "converged"
159            break
160        icycle += 1
161    M = func(x0,*args)
162    nfev += 1
163    Yvec,Amat = Hess(x0,*args)
164    Adiag = np.sqrt(np.diag(Amat))
165    Anorm = np.outer(Adiag,Adiag)
166    Lam = np.eye(Amat.shape[0])*lam
167    Amatlam = Amat/Anorm  #*(One+Lam)              #don't scale Amat to Marquardt array       
168    try:
169        Bmat = nl.inv(Amatlam)/Anorm  #*(One+Lam)      #don't rescale Bmat to Marquardt array
170        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[], 'Converged': ifConverged, 'DelChi2':deltaChi2}]
171    except nl.LinAlgError:
172        print 'ouch #2 linear algebra error in making v-cov matrix'
173        psing = []
174        if maxcyc:
175            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
176        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
177   
178def getVCov(varyNames,varyList,covMatrix):
179    '''obtain variance-covariance terms for a set of variables. NB: the varyList
180    and covMatrix were saved by the last least squares refinement so they must match.
181   
182    :param list varyNames: variable names to find v-cov matric for
183    :param list varyList: full list of all variables in v-cov matrix
184    :param nparray covMatrix: full variance-covariance matrix from the last
185     least squares refinement
186   
187    :returns: nparray vcov: variance-covariance matrix for the variables given
188     in varyNames
189   
190    '''
191    vcov = np.zeros((len(varyNames),len(varyNames)))
192    for i1,name1 in enumerate(varyNames):
193        for i2,name2 in enumerate(varyNames):
194            try:
195                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
196            except ValueError:
197                vcov[i1][i2] = 0.0
198#                if i1 == i2:
199#                    vcov[i1][i2] = 1e-20
200#                else:
201#                    vcov[i1][i2] = 0.0
202    return vcov
203   
204################################################################################
205##### Atom manipulations
206################################################################################
207
208def FindMolecule(ind,generalData,atomData):                    #uses numpy & masks - very fast even for proteins!
209
210    def getNeighbors(atom,radius):
211        neighList = [] 
212        Dx = UAtoms-np.array(atom[cx:cx+3])
213        dist = ma.masked_less(np.sqrt(np.sum(np.inner(Amat,Dx)**2,axis=0)),0.5) #gets rid of disorder "bonds" < 0.5A
214        sumR = Radii+radius
215        return set(ma.nonzero(ma.masked_greater(dist-factor*sumR,0.))[0])                #get indices of bonded atoms
216
217    import numpy.ma as ma
218    indices = (-1,0,1)
219    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices],dtype='f')
220    cx,ct,cs,ci = generalData['AtomPtrs']
221    DisAglCtls = generalData['DisAglCtls']
222    SGData = generalData['SGData']
223    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
224    radii = DisAglCtls['BondRadii']
225    atomTypes = DisAglCtls['AtomTypes']
226    factor = DisAglCtls['Factors'][0]
227    unit = np.zeros(3)
228    try:
229        indH = atomTypes.index('H')
230        radii[indH] = 0.5
231    except:
232        pass           
233    nAtom = len(atomData)
234    Indx = range(nAtom)
235    UAtoms = []
236    Radii = []
237    for atom in atomData:
238        UAtoms.append(np.array(atom[cx:cx+3]))
239        Radii.append(radii[atomTypes.index(atom[ct])])
240    UAtoms = np.array(UAtoms)
241    Radii = np.array(Radii)
242    for nOp,Op in enumerate(SGData['SGOps'][1:]):
243        UAtoms = np.concatenate((UAtoms,(np.inner(Op[0],UAtoms[:nAtom]).T+Op[1])))
244        Radii = np.concatenate((Radii,Radii[:nAtom]))
245        Indx += Indx[:nAtom]
246    for icen,cen in enumerate(SGData['SGCen'][1:]):
247        UAtoms = np.concatenate((UAtoms,(UAtoms+cen)))
248        Radii = np.concatenate((Radii,Radii))
249        Indx += Indx[:nAtom]
250    if SGData['SGInv']:
251        UAtoms = np.concatenate((UAtoms,-UAtoms))
252        Radii = np.concatenate((Radii,Radii))
253        Indx += Indx
254    UAtoms %= 1.
255    mAtoms = len(UAtoms)
256    for unit in Units:
257        if np.any(unit):    #skip origin cell
258            UAtoms = np.concatenate((UAtoms,UAtoms[:mAtoms]+unit))
259            Radii = np.concatenate((Radii,Radii[:mAtoms]))
260            Indx += Indx[:mAtoms]
261    UAtoms = np.array(UAtoms)
262    Radii = np.array(Radii)
263    newAtoms = [atomData[ind],]
264    atomData[ind] = None
265    radius = Radii[ind]
266    IndB = getNeighbors(newAtoms[-1],radius)
267    while True:
268        if not len(IndB):
269            break
270        indb = IndB.pop()
271        if atomData[Indx[indb]] == None:
272            continue
273        while True:
274            try:
275                jndb = Indb.index(indb)
276                Indb.remove(jndb)
277            except:
278                break
279        newAtom = copy.copy(atomData[Indx[indb]])
280        newAtom[cx:cx+3] = UAtoms[indb]     #NB: thermal Uij, etc. not transformed!
281        newAtoms.append(newAtom)
282        atomData[Indx[indb]] = None
283        IndB = set(list(IndB)+list(getNeighbors(newAtoms[-1],radius)))
284        if len(IndB) > nAtom:
285            return 'Assemble molecule cannot be used on extended structures'
286    for atom in atomData:
287        if atom != None:
288            newAtoms.append(atom)
289    return newAtoms
290       
291def FindAtomIndexByIDs(atomData,IDs,Draw=True):
292    '''finds the set of atom array indices for a list of atom IDs. Will search
293    either the Atom table or the drawAtom table.
294   
295    :param list atomData: Atom or drawAtom table containting coordinates, etc.
296    :param list IDs: atom IDs to be found
297    :param bool Draw: True if drawAtom table to be searched; False if Atom table
298      is searched
299   
300    :returns: list indx: atom (or drawAtom) indices
301   
302    '''
303    indx = []
304    for i,atom in enumerate(atomData):
305        if Draw and atom[-3] in IDs:
306            indx.append(i)
307        elif atom[-1] in IDs:
308            indx.append(i)
309    return indx
310
311def FillAtomLookUp(atomData,indx):
312    '''create a dictionary of atom indexes with atom IDs as keys
313   
314    :param list atomData: Atom table to be used
315   
316    :returns: dict atomLookUp: dictionary of atom indexes with atom IDs as keys
317   
318    '''
319    atomLookUp = {}
320    for iatm,atom in enumerate(atomData):
321        atomLookUp[atom[indx]] = iatm
322    return atomLookUp
323
324def GetAtomsById(atomData,atomLookUp,IdList):
325    '''gets a list of atoms from Atom table that match a set of atom IDs
326   
327    :param list atomData: Atom table to be used
328    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
329    :param list IdList: atom IDs to be found
330   
331    :returns: list atoms: list of atoms found
332   
333    '''
334    atoms = []
335    for id in IdList:
336        atoms.append(atomData[atomLookUp[id]])
337    return atoms
338   
339def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
340    '''gets atom parameters for atoms using atom IDs
341   
342    :param list atomData: Atom table to be used
343    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
344    :param list IdList: atom IDs to be found
345    :param int itemLoc: pointer to desired 1st item in an atom table entry
346    :param int numItems: number of items to be retrieved
347   
348    :returns: type name: description
349   
350    '''
351    Items = []
352    if not isinstance(IdList,list):
353        IdList = [IdList,]
354    for id in IdList:
355        if numItems == 1:
356            Items.append(atomData[atomLookUp[id]][itemLoc])
357        else:
358            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
359    return Items
360   
361def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
362    '''default doc string
363   
364    :param type name: description
365   
366    :returns: type name: description
367   
368    '''
369    pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
370    dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
371    XYZ = []
372    for ind in indx:
373        names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
374        dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
375        XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
376    return XYZ
377
378def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
379    '''default doc string
380   
381    :param type name: description
382   
383    :returns: type name: description
384   
385    '''
386    for atom in atomData:
387        XYZ = np.inner(Amat,atom[cx:cx+3])
388        if atom[cia] == 'A':
389            UIJ = atom[cia+2:cia+8]
390               
391def TLS2Uij(xyz,g,Amat,rbObj):    #not used anywhere, but could be?
392    '''default doc string
393   
394    :param type name: description
395   
396    :returns: type name: description
397   
398    '''
399    TLStype,TLS = rbObj['ThermalMotion'][:2]
400    Tmat = np.zeros((3,3))
401    Lmat = np.zeros((3,3))
402    Smat = np.zeros((3,3))
403    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
404        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
405    if 'T' in TLStype:
406        Tmat = G2lat.U6toUij(TLS[:6])
407    if 'L' in TLStype:
408        Lmat = G2lat.U6toUij(TLS[6:12])
409    if 'S' in TLStype:
410        Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
411    XYZ = np.inner(Amat,xyz)
412    Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
413    Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
414    beta = np.inner(np.inner(g,Umat),g)
415    return G2lat.UijtoU6(beta)*gvec   
416       
417def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj):    #not used anywhere, but could be?
418    '''default doc string
419   
420    :param type name: description
421   
422    :returns: type name: description
423   
424    '''
425    cx,ct,cs,cia = atPtrs
426    TLStype,TLS = rbObj['ThermalMotion'][:2]
427    Tmat = np.zeros((3,3))
428    Lmat = np.zeros((3,3))
429    Smat = np.zeros((3,3))
430    G,g = G2lat.A2Gmat(Amat)
431    gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
432    if 'T' in TLStype:
433        Tmat = G2lat.U6toUij(TLS[:6])
434    if 'L' in TLStype:
435        Lmat = G2lat.U6toUij(TLS[6:12])
436    if 'S' in TLStype:
437        Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
438    for atom in atomData:
439        XYZ = np.inner(Amat,atom[cx:cx+3])
440        Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
441        if 'U' in TSLtype:
442            atom[cia+1] = TLS[0]
443            atom[cia] = 'I'
444        else:
445            atom[cia] = 'A'
446            Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
447            beta = np.inner(np.inner(g,Umat),g)
448            atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
449
450def GetXYZDist(xyz,XYZ,Amat):
451    '''gets distance from position xyz to all XYZ, xyz & XYZ are np.array
452        and are in crystal coordinates; Amat is crystal to Cart matrix
453   
454    :param type name: description
455   
456    :returns: type name: description
457   
458    '''
459    return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
460
461def getAtomXYZ(atoms,cx):
462    '''default doc string
463   
464    :param type name: description
465   
466    :returns: type name: description
467   
468    '''
469    XYZ = []
470    for atom in atoms:
471        XYZ.append(atom[cx:cx+3])
472    return np.array(XYZ)
473
474def RotateRBXYZ(Bmat,Cart,oriQ):
475    '''rotate & transform cartesian coordinates to crystallographic ones
476    no translation applied. To be used for numerical derivatives
477   
478    :param type name: description
479   
480    :returns: type name: description
481   
482    '''
483    ''' returns crystal coordinates for atoms described by RBObj
484    '''
485    XYZ = np.zeros_like(Cart)
486    for i,xyz in enumerate(Cart):
487        XYZ[i] = np.inner(Bmat,prodQVQ(oriQ,xyz))
488    return XYZ
489
490def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
491    '''default doc string
492   
493    :param type name: description
494   
495    :returns: type name: description
496   
497    '''
498    ''' returns crystal coordinates for atoms described by RBObj
499    '''
500    RBRes = RBData[RBType][RBObj['RBId']]
501    if RBType == 'Vector':
502        vecs = RBRes['rbVect']
503        mags = RBRes['VectMag']
504        Cart = np.zeros_like(vecs[0])
505        for vec,mag in zip(vecs,mags):
506            Cart += vec*mag
507    elif RBType == 'Residue':
508        Cart = np.array(RBRes['rbXYZ'])
509        for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
510            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
511            Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
512    XYZ = np.zeros_like(Cart)
513    for i,xyz in enumerate(Cart):
514        XYZ[i] = np.inner(Bmat,prodQVQ(RBObj['Orient'][0],xyz))+RBObj['Orig'][0]
515    return XYZ,Cart
516
517def UpdateMCSAxyz(Bmat,MCSA):
518    '''default doc string
519   
520    :param type name: description
521   
522    :returns: type name: description
523   
524    '''
525    xyz = []
526    atTypes = []
527    iatm = 0
528    for model in MCSA['Models'][1:]:        #skip the MD model
529        if model['Type'] == 'Atom':
530            xyz.append(model['Pos'][0])
531            atTypes.append(model['atType'])
532            iatm += 1
533        else:
534            RBRes = MCSA['rbData'][model['Type']][model['RBId']]
535            Pos = np.array(model['Pos'][0])
536            Ori = np.array(model['Ori'][0])
537            Qori = AVdeg2Q(Ori[0],Ori[1:])
538            if model['Type'] == 'Vector':
539                vecs = RBRes['rbVect']
540                mags = RBRes['VectMag']
541                Cart = np.zeros_like(vecs[0])
542                for vec,mag in zip(vecs,mags):
543                    Cart += vec*mag
544            elif model['Type'] == 'Residue':
545                Cart = np.array(RBRes['rbXYZ'])
546                for itor,seq in enumerate(RBRes['rbSeq']):
547                    QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
548                    Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
549            if model['MolCent'][1]:
550                Cart -= model['MolCent'][0]
551            for i,x in enumerate(Cart):
552                xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
553                atType = RBRes['rbTypes'][i]
554                atTypes.append(atType)
555                iatm += 1
556    return np.array(xyz),atTypes
557   
558def SetMolCent(model,RBData):
559    '''default doc string
560   
561    :param type name: description
562   
563    :returns: type name: description
564   
565    '''
566    rideList = []
567    RBRes = RBData[model['Type']][model['RBId']]
568    if model['Type'] == 'Vector':
569        vecs = RBRes['rbVect']
570        mags = RBRes['VectMag']
571        Cart = np.zeros_like(vecs[0])
572        for vec,mag in zip(vecs,mags):
573            Cart += vec*mag
574    elif model['Type'] == 'Residue':
575        Cart = np.array(RBRes['rbXYZ'])
576        for seq in RBRes['rbSeq']:
577            rideList += seq[3]
578    centList = set(range(len(Cart)))-set(rideList)
579    cent = np.zeros(3)
580    for i in centList:
581        cent += Cart[i]
582    model['MolCent'][0] = cent/len(centList)   
583   
584def UpdateRBUIJ(Bmat,Cart,RBObj):
585    '''default doc string
586   
587    :param type name: description
588   
589    :returns: type name: description
590   
591    '''
592    ''' returns atom I/A, Uiso or UIJ for atoms at XYZ as described by RBObj
593    '''
594    TLStype,TLS = RBObj['ThermalMotion'][:2]
595    T = np.zeros(6)
596    L = np.zeros(6)
597    S = np.zeros(8)
598    if 'T' in TLStype:
599        T = TLS[:6]
600    if 'L' in TLStype:
601        L = np.array(TLS[6:12])*(np.pi/180.)**2
602    if 'S' in TLStype:
603        S = np.array(TLS[12:])*(np.pi/180.)
604    g = nl.inv(np.inner(Bmat,Bmat))
605    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
606        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
607    Uout = []
608    Q = RBObj['Orient'][0]
609    for X in Cart:
610        X = prodQVQ(Q,X)
611        if 'U' in TLStype:
612            Uout.append(['I',TLS[0],0,0,0,0,0,0])
613        elif not 'N' in TLStype:
614            U = [0,0,0,0,0,0]
615            U[0] = T[0]+L[1]*X[2]**2+L[2]*X[1]**2-2.0*L[5]*X[1]*X[2]+2.0*(S[2]*X[2]-S[4]*X[1])
616            U[1] = T[1]+L[0]*X[2]**2+L[2]*X[0]**2-2.0*L[4]*X[0]*X[2]+2.0*(S[5]*X[0]-S[0]*X[2])
617            U[2] = T[2]+L[1]*X[0]**2+L[0]*X[1]**2-2.0*L[3]*X[1]*X[0]+2.0*(S[1]*X[1]-S[3]*X[0])
618            U[3] = T[3]+L[4]*X[1]*X[2]+L[5]*X[0]*X[2]-L[3]*X[2]**2-L[2]*X[0]*X[1]+  \
619                S[4]*X[0]-S[5]*X[1]-(S[6]+S[7])*X[2]
620            U[4] = T[4]+L[3]*X[1]*X[2]+L[5]*X[0]*X[1]-L[4]*X[1]**2-L[1]*X[0]*X[2]+  \
621                S[3]*X[2]-S[2]*X[0]+S[6]*X[1]
622            U[5] = T[5]+L[3]*X[0]*X[2]+L[4]*X[0]*X[1]-L[5]*X[0]**2-L[0]*X[2]*X[1]+  \
623                S[0]*X[1]-S[1]*X[2]+S[7]*X[0]
624            Umat = G2lat.U6toUij(U)
625            beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
626            Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
627        else:
628            Uout.append(['N',])
629    return Uout
630   
631def GetSHCoeff(pId,parmDict,SHkeys):
632    '''default doc string
633   
634    :param type name: description
635   
636    :returns: type name: description
637   
638    '''
639    SHCoeff = {}
640    for shkey in SHkeys:
641        shname = str(pId)+'::'+shkey
642        SHCoeff[shkey] = parmDict[shname]
643    return SHCoeff
644       
645def getMass(generalData):
646    '''Computes mass of unit cell contents
647   
648    :param dict generalData: The General dictionary in Phase
649   
650    :returns: float mass: Crystal unit cell mass in AMU.
651   
652    '''
653    mass = 0.
654    for i,elem in enumerate(generalData['AtomTypes']):
655        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
656    return mass   
657
658def getDensity(generalData):
659    '''calculate crystal structure density
660   
661    :param dict generalData: The General dictionary in Phase
662   
663    :returns: float density: crystal density in gm/cm^3
664   
665    '''
666    mass = getMass(generalData)
667    Volume = generalData['Cell'][7]
668    density = mass/(0.6022137*Volume)
669    return density,Volume/mass
670   
671def getWave(Parms):
672    '''returns wavelength from Instrument parameters dictionary
673   
674    :param dict Parms: Instrument parameters;
675        must contain:
676        Lam: single wavelength
677        or
678        Lam1: Ka1 radiation wavelength
679   
680    :returns: float wave: wavelength
681   
682    '''
683    try:
684        return Parms['Lam'][1]
685    except KeyError:
686        return Parms['Lam1'][1]
687       
688def El2Mass(Elements):
689    '''compute molecular weight from Elements
690   
691    :param dict Elements: elements in molecular formula;
692        each must contain
693        Num: number of atoms in formula
694        Mass: at. wt.
695   
696    :returns: float mass: molecular weight.
697   
698    '''
699    mass = 0
700    for El in Elements:
701        mass += Elements[El]['Num']*Elements[El]['Mass']
702    return mass
703       
704def Den2Vol(Elements,density):
705    '''converts density to molecular volume
706   
707    :param dict Elements: elements in molecular formula;
708        each must contain
709        Num: number of atoms in formula
710        Mass: at. wt.
711    :param float density: material density in gm/cm^3
712   
713    :returns: float volume: molecular volume in A^3
714   
715    '''
716    return El2Mass(Elements)/(density*0.6022137)
717   
718def Vol2Den(Elements,volume):
719    '''converts volume to density
720   
721    :param dict Elements: elements in molecular formula;
722        each must contain
723        Num: number of atoms in formula
724        Mass: at. wt.
725    :param float volume: molecular volume in A^3
726   
727    :returns: float density: material density in gm/cm^3
728   
729    '''
730    return El2Mass(Elements)/(volume*0.6022137)
731   
732def El2EstVol(Elements):
733    '''Estimate volume from molecular formula; assumes atom volume = 10A^3
734   
735    :param dict Elements: elements in molecular formula;
736        each must contain
737        Num: number of atoms in formula
738   
739    :returns: float volume: estimate of molecular volume in A^3
740   
741    '''
742    vol = 0
743    for El in Elements:
744        vol += 10.*Elements[El]['Num']
745    return vol
746   
747def XScattDen(Elements,vol,wave=0.):
748    '''Estimate X-ray scattering density from molecular formula & volume;
749    ignores valence, but includes anomalous effects
750   
751    :param dict Elements: elements in molecular formula;
752        each element must contain
753        Num: number of atoms in formula
754        Z: atomic number
755    :param float vol: molecular volume in A^3
756    :param float wave: optional wavelength in A
757   
758    :returns: float rho: scattering density in 10^10cm^-2;
759        if wave > 0 the includes f' contribution
760    :returns: float mu: if wave>0 absorption coeff in cm^-1 ; otherwise 0
761   
762    '''
763    rho = 0
764    mu = 0
765    if wave:
766        Xanom = XAnomAbs(Elements,wave)
767    for El in Elements:
768        f0 = Elements[El]['Z']
769        if wave:
770            f0 += Xanom[El][0]
771            mu += Xanom[El][2]*Elements[El]['Num']
772        rho += Elements[El]['Num']*f0
773    return 28.179*rho/vol,0.1*mu/vol
774   
775def wavekE(wavekE):
776    '''Convert wavelength to energy & vise versa
777   
778    :param float waveKe:wavelength in A or energy in kE
779   
780    :returns float waveKe:the other one
781   
782    '''
783    return 12.397639/wavekE
784       
785def XAnomAbs(Elements,wave):
786    kE = wavekE(wave)
787    Xanom = {}
788    for El in Elements:
789        Orbs = G2el.GetXsectionCoeff(El)
790        Xanom[El] = G2el.FPcalc(Orbs, kE)
791    return Xanom
792   
793def Modulation(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata):
794    import scipy.special as sp
795    import scipy.integrate as si
796    m = SSUniq.T[3]
797    nh = np.zeros(1)
798    if XSSdata.ndim > 2:
799        nh = np.arange(XSSdata.shape[1])       
800    M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])
801    A = np.array(XSSdata[:3])
802    B = np.array(XSSdata[3:])
803    HdotA = (np.inner(A.T,SSUniq.T[:3].T)+SSPhi)
804    HdotB = (np.inner(B.T,SSUniq.T[:3].T)+SSPhi)
805    GpA = sp.jn(M[:,np.newaxis],twopi*HdotA)
806    GpB = sp.jn(M[:,np.newaxis],twopi*HdotB)*(1.j)**M
807    Gp = np.sum(GpA+GpB,axis=0)
808    return np.real(Gp).T,np.imag(Gp).T
809   
810def ModulationDerv(waveTypes,SSUniq,SSPhi,FSSdata,XSSdata,USSdata):
811    import scipy.special as sp
812    m = SSUniq.T[3]
813    nh = np.zeros(1)
814    if XSSdata.ndim > 2:
815        nh = np.arange(XSSdata.shape[1])       
816    M = np.where(m>0,m+nh[:,np.newaxis],m-nh[:,np.newaxis])
817    A = np.array([[a,b] for a,b in zip(XSSdata[:3],XSSdata[3:])])
818    HdotA = twopi*(np.inner(SSUniq.T[:3].T,A.T)+SSPhi)
819    Gpm = sp.jn(M[:,np.newaxis,:]-1,HdotA)
820    Gpp = sp.jn(M[:,np.newaxis,:]+1,HdotA)
821    if Gpm.ndim > 3: #sum over multiple harmonics
822        Gpm = np.sum(Gpm,axis=0)
823        Gpp = np.sum(Gpp,axis=0)
824    dGpdk = 0.5*(Gpm+Gpp)
825    return np.real(dGpdk),np.imag(dGpdk)
826   
827def posFourier(tau,psin,pcos):
828    A = np.array([ps[:,np.newaxis]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)])
829    B = np.array([pc[:,np.newaxis]*np.cos(2*np.pi*(i+1)*tau) for i,pc in enumerate(pcos)])
830    return np.sum(A,axis=0)+np.sum(B,axis=0)
831   
832def posSawtooth(tau,Toff,slopes):
833    Tau = (tau-Toff)%1.
834    A = slopes[:,np.newaxis]*Tau
835    return A
836
837def posZigZag(tau,Toff,slopes):
838    Tau = (tau-Toff)%1.
839    A = np.where(Tau <= 0.5,slopes[:,np.newaxis]*Tau,slopes[:,np.newaxis]*(1.-Tau))
840    return A
841   
842def fracCrenel(tau,Toff,Twid):
843    Tau = (tau-Toff)%1.
844    A = np.where(Tau<Twid,1.,0.)
845    return A
846   
847def fracFourier(tau,fsin,fcos):
848    A = np.array([fs[:,np.newaxis]*np.sin(2.*np.pi*(i+1)*tau) for i,fs in enumerate(fsin)])
849    B = np.array([fc[:,np.newaxis]*np.cos(2.*np.pi*(i+1)*tau) for i,fc in enumerate(fcos)])
850    return np.sum(A,axis=0)+np.sum(B,axis=0)
851   
852################################################################################
853##### distance, angle, planes, torsion stuff
854################################################################################
855
856def getSyXYZ(XYZ,ops,SGData):
857    '''default doc string
858   
859    :param type name: description
860   
861    :returns: type name: description
862   
863    '''
864    XYZout = np.zeros_like(XYZ)
865    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
866        if op == '1':
867            XYZout[i] = xyz
868        else:
869            oprs = op.split('+')
870            unit = [0,0,0]
871            if len(oprs)>1:
872                unit = np.array(list(eval(oprs[1])))
873            syop =int(oprs[0])
874            inv = syop/abs(syop)
875            syop *= inv
876            cent = syop/100
877            syop %= 100
878            syop -= 1
879            M,T = SGData['SGOps'][syop]
880            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
881    return XYZout
882   
883def getRestDist(XYZ,Amat):
884    '''default doc string
885   
886    :param type name: description
887   
888    :returns: type name: description
889   
890    '''
891    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
892   
893def getRestDeriv(Func,XYZ,Amat,ops,SGData):
894    '''default doc string
895   
896    :param type name: description
897   
898    :returns: type name: description
899   
900    '''
901    deriv = np.zeros((len(XYZ),3))
902    dx = 0.00001
903    for j,xyz in enumerate(XYZ):
904        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
905            XYZ[j] -= x
906            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
907            XYZ[j] += 2*x
908            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
909            XYZ[j] -= x
910            deriv[j][i] = (d1-d2)/(2*dx)
911    return deriv.flatten()
912
913def getRestAngle(XYZ,Amat):
914    '''default doc string
915   
916    :param type name: description
917   
918    :returns: type name: description
919   
920    '''
921   
922    def calcVec(Ox,Tx,Amat):
923        return np.inner(Amat,(Tx-Ox))
924
925    VecA = calcVec(XYZ[1],XYZ[0],Amat)
926    VecA /= np.sqrt(np.sum(VecA**2))
927    VecB = calcVec(XYZ[1],XYZ[2],Amat)
928    VecB /= np.sqrt(np.sum(VecB**2))
929    edge = VecB-VecA
930    edge = np.sum(edge**2)
931    angle = (2.-edge)/2.
932    angle = max(angle,-1.)
933    return acosd(angle)
934   
935def getRestPlane(XYZ,Amat):
936    '''default doc string
937   
938    :param type name: description
939   
940    :returns: type name: description
941   
942    '''
943    sumXYZ = np.zeros(3)
944    for xyz in XYZ:
945        sumXYZ += xyz
946    sumXYZ /= len(XYZ)
947    XYZ = np.array(XYZ)-sumXYZ
948    XYZ = np.inner(Amat,XYZ).T
949    Zmat = np.zeros((3,3))
950    for i,xyz in enumerate(XYZ):
951        Zmat += np.outer(xyz.T,xyz)
952    Evec,Emat = nl.eig(Zmat)
953    Evec = np.sqrt(Evec)/(len(XYZ)-3)
954    Order = np.argsort(Evec)
955    return Evec[Order[0]]
956   
957def getRestChiral(XYZ,Amat):   
958    '''default doc string
959   
960    :param type name: description
961   
962    :returns: type name: description
963   
964    '''
965    VecA = np.empty((3,3))   
966    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
967    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
968    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
969    return nl.det(VecA)
970   
971def getRestTorsion(XYZ,Amat):
972    '''default doc string
973   
974    :param type name: description
975   
976    :returns: type name: description
977   
978    '''
979    VecA = np.empty((3,3))
980    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
981    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
982    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
983    D = nl.det(VecA)
984    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
985    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
986    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
987    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
988    Ang = 1.0
989    if abs(P12) < 1.0 and abs(P23) < 1.0:
990        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
991    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
992    return TOR
993   
994def calcTorsionEnergy(TOR,Coeff=[]):
995    '''default doc string
996   
997    :param type name: description
998   
999    :returns: type name: description
1000   
1001    '''
1002    sum = 0.
1003    if len(Coeff):
1004        cof = np.reshape(Coeff,(3,3)).T
1005        delt = TOR-cof[1]
1006        delt = np.where(delt<-180.,delt+360.,delt)
1007        delt = np.where(delt>180.,delt-360.,delt)
1008        term = -cof[2]*delt**2
1009        val = cof[0]*np.exp(term/1000.0)
1010        pMax = cof[0][np.argmin(val)]
1011        Eval = np.sum(val)
1012        sum = Eval-pMax
1013    return sum,Eval
1014
1015def getTorsionDeriv(XYZ,Amat,Coeff):
1016    '''default doc string
1017   
1018    :param type name: description
1019   
1020    :returns: type name: description
1021   
1022    '''
1023    deriv = np.zeros((len(XYZ),3))
1024    dx = 0.00001
1025    for j,xyz in enumerate(XYZ):
1026        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1027            XYZ[j] -= x
1028            tor = getRestTorsion(XYZ,Amat)
1029            p1,d1 = calcTorsionEnergy(tor,Coeff)
1030            XYZ[j] += 2*x
1031            tor = getRestTorsion(XYZ,Amat)
1032            p2,d2 = calcTorsionEnergy(tor,Coeff)           
1033            XYZ[j] -= x
1034            deriv[j][i] = (p2-p1)/(2*dx)
1035    return deriv.flatten()
1036
1037def getRestRama(XYZ,Amat):
1038    '''Computes a pair of torsion angles in a 5 atom string
1039   
1040    :param nparray XYZ: crystallographic coordinates of 5 atoms
1041    :param nparray Amat: crystal to cartesian transformation matrix
1042   
1043    :returns: list (phi,psi) two torsion angles in degrees
1044   
1045    '''
1046    phi = getRestTorsion(XYZ[:5],Amat)
1047    psi = getRestTorsion(XYZ[1:],Amat)
1048    return phi,psi
1049   
1050def calcRamaEnergy(phi,psi,Coeff=[]):
1051    '''Computes pseudo potential energy from a pair of torsion angles and a
1052    numerical description of the potential energy surface. Used to create
1053    penalty function in LS refinement:     
1054    :math:`Eval(\\phi,\\psi) = C[0]*exp(-V/1000)`
1055
1056    where :math:`V = -C[3] * (\\phi-C[1])^2 - C[4]*(\\psi-C[2])^2 - 2*(\\phi-C[1])*(\\psi-C[2])`
1057   
1058    :param float phi: first torsion angle (:math:`\\phi`)
1059    :param float psi: second torsion angle (:math:`\\psi`)
1060    :param list Coeff: pseudo potential coefficients
1061   
1062    :returns: list (sum,Eval): pseudo-potential difference from minimum & value;
1063      sum is used for penalty function.
1064   
1065    '''
1066    sum = 0.
1067    Eval = 0.
1068    if len(Coeff):
1069        cof = Coeff.T
1070        dPhi = phi-cof[1]
1071        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
1072        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
1073        dPsi = psi-cof[2]
1074        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
1075        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
1076        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
1077        val = cof[0]*np.exp(val/1000.)
1078        pMax = cof[0][np.argmin(val)]
1079        Eval = np.sum(val)
1080        sum = Eval-pMax
1081    return sum,Eval
1082
1083def getRamaDeriv(XYZ,Amat,Coeff):
1084    '''Computes numerical derivatives of torsion angle pair pseudo potential
1085    with respect of crystallographic atom coordinates of the 5 atom sequence
1086   
1087    :param nparray XYZ: crystallographic coordinates of 5 atoms
1088    :param nparray Amat: crystal to cartesian transformation matrix
1089    :param list Coeff: pseudo potential coefficients
1090   
1091    :returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
1092     crystallographic xyz coordinates.
1093   
1094    '''
1095    deriv = np.zeros((len(XYZ),3))
1096    dx = 0.00001
1097    for j,xyz in enumerate(XYZ):
1098        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1099            XYZ[j] -= x
1100            phi,psi = getRestRama(XYZ,Amat)
1101            p1,d1 = calcRamaEnergy(phi,psi,Coeff)
1102            XYZ[j] += 2*x
1103            phi,psi = getRestRama(XYZ,Amat)
1104            p2,d2 = calcRamaEnergy(phi,psi,Coeff)
1105            XYZ[j] -= x
1106            deriv[j][i] = (p2-p1)/(2*dx)
1107    return deriv.flatten()
1108
1109def getRestPolefig(ODFln,SamSym,Grid):
1110    '''default doc string
1111   
1112    :param type name: description
1113   
1114    :returns: type name: description
1115   
1116    '''
1117    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
1118    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
1119    R = np.where(R <= 1.,2.*atand(R),0.0)
1120    Z = np.zeros_like(R)
1121    Z = G2lat.polfcal(ODFln,SamSym,R,P)
1122    Z = np.reshape(Z,(Grid,Grid))
1123    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
1124
1125def getRestPolefigDerv(HKL,Grid,SHCoeff):
1126    '''default doc string
1127   
1128    :param type name: description
1129   
1130    :returns: type name: description
1131   
1132    '''
1133    pass
1134       
1135def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
1136    '''default doc string
1137   
1138    :param type name: description
1139   
1140    :returns: type name: description
1141   
1142    '''
1143    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
1144        TxT = inv*(np.inner(M,Tx)+T)+C+U
1145        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
1146       
1147    inv = Top/abs(Top)
1148    cent = abs(Top)/100
1149    op = abs(Top)%100-1
1150    M,T = SGData['SGOps'][op]
1151    C = SGData['SGCen'][cent]
1152    dx = .00001
1153    deriv = np.zeros(6)
1154    for i in [0,1,2]:
1155        Oxyz[i] -= dx
1156        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1157        Oxyz[i] += 2*dx
1158        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1159        Oxyz[i] -= dx
1160        Txyz[i] -= dx
1161        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1162        Txyz[i] += 2*dx
1163        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1164        Txyz[i] -= dx
1165    return deriv
1166   
1167def getAngSig(VA,VB,Amat,SGData,covData={}):
1168    '''default doc string
1169   
1170    :param type name: description
1171   
1172    :returns: type name: description
1173   
1174    '''
1175    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
1176        TxT = inv*(np.inner(M,Tx)+T)+C+U
1177        return np.inner(Amat,(TxT-Ox))
1178       
1179    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
1180        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
1181        VecA /= np.sqrt(np.sum(VecA**2))
1182        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
1183        VecB /= np.sqrt(np.sum(VecB**2))
1184        edge = VecB-VecA
1185        edge = np.sum(edge**2)
1186        angle = (2.-edge)/2.
1187        angle = max(angle,-1.)
1188        return acosd(angle)
1189       
1190    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
1191    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
1192    invA = invB = 1
1193    invA = TopA/abs(TopA)
1194    invB = TopB/abs(TopB)
1195    centA = abs(TopA)/100
1196    centB = abs(TopB)/100
1197    opA = abs(TopA)%100-1
1198    opB = abs(TopB)%100-1
1199    MA,TA = SGData['SGOps'][opA]
1200    MB,TB = SGData['SGOps'][opB]
1201    CA = SGData['SGCen'][centA]
1202    CB = SGData['SGCen'][centB]
1203    if 'covMatrix' in covData:
1204        covMatrix = covData['covMatrix']
1205        varyList = covData['varyList']
1206        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
1207        dx = .00001
1208        dadx = np.zeros(9)
1209        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1210        for i in [0,1,2]:
1211            OxA[i] -= dx
1212            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1213            OxA[i] += 2*dx
1214            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1215            OxA[i] -= dx
1216           
1217            TxA[i] -= dx
1218            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1219            TxA[i] += 2*dx
1220            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1221            TxA[i] -= dx
1222           
1223            TxB[i] -= dx
1224            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1225            TxB[i] += 2*dx
1226            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1227            TxB[i] -= dx
1228           
1229        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1230        if sigAng < 0.01:
1231            sigAng = 0.0
1232        return Ang,sigAng
1233    else:
1234        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
1235
1236def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
1237    '''default doc string
1238   
1239    :param type name: description
1240   
1241    :returns: type name: description
1242   
1243    '''
1244    def calcDist(Atoms,SyOps,Amat):
1245        XYZ = []
1246        for i,atom in enumerate(Atoms):
1247            Inv,M,T,C,U = SyOps[i]
1248            XYZ.append(np.array(atom[1:4]))
1249            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1250            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1251        V1 = XYZ[1]-XYZ[0]
1252        return np.sqrt(np.sum(V1**2))
1253       
1254    Inv = []
1255    SyOps = []
1256    names = []
1257    for i,atom in enumerate(Oatoms):
1258        names += atom[-1]
1259        Op,unit = Atoms[i][-1]
1260        inv = Op/abs(Op)
1261        m,t = SGData['SGOps'][abs(Op)%100-1]
1262        c = SGData['SGCen'][abs(Op)/100]
1263        SyOps.append([inv,m,t,c,unit])
1264    Dist = calcDist(Oatoms,SyOps,Amat)
1265   
1266    sig = -0.001
1267    if 'covMatrix' in covData:
1268        parmNames = []
1269        dx = .00001
1270        dadx = np.zeros(6)
1271        for i in range(6):
1272            ia = i/3
1273            ix = i%3
1274            Oatoms[ia][ix+1] += dx
1275            a0 = calcDist(Oatoms,SyOps,Amat)
1276            Oatoms[ia][ix+1] -= 2*dx
1277            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1278        covMatrix = covData['covMatrix']
1279        varyList = covData['varyList']
1280        DistVcov = getVCov(names,varyList,covMatrix)
1281        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
1282        if sig < 0.001:
1283            sig = -0.001
1284   
1285    return Dist,sig
1286
1287def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
1288    '''default doc string
1289   
1290    :param type name: description
1291   
1292    :returns: type name: description
1293   
1294    '''
1295
1296    def calcAngle(Atoms,SyOps,Amat):
1297        XYZ = []
1298        for i,atom in enumerate(Atoms):
1299            Inv,M,T,C,U = SyOps[i]
1300            XYZ.append(np.array(atom[1:4]))
1301            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1302            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1303        V1 = XYZ[1]-XYZ[0]
1304        V1 /= np.sqrt(np.sum(V1**2))
1305        V2 = XYZ[1]-XYZ[2]
1306        V2 /= np.sqrt(np.sum(V2**2))
1307        V3 = V2-V1
1308        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1309        return acosd(cang)
1310
1311    Inv = []
1312    SyOps = []
1313    names = []
1314    for i,atom in enumerate(Oatoms):
1315        names += atom[-1]
1316        Op,unit = Atoms[i][-1]
1317        inv = Op/abs(Op)
1318        m,t = SGData['SGOps'][abs(Op)%100-1]
1319        c = SGData['SGCen'][abs(Op)/100]
1320        SyOps.append([inv,m,t,c,unit])
1321    Angle = calcAngle(Oatoms,SyOps,Amat)
1322   
1323    sig = -0.01
1324    if 'covMatrix' in covData:
1325        parmNames = []
1326        dx = .00001
1327        dadx = np.zeros(9)
1328        for i in range(9):
1329            ia = i/3
1330            ix = i%3
1331            Oatoms[ia][ix+1] += dx
1332            a0 = calcAngle(Oatoms,SyOps,Amat)
1333            Oatoms[ia][ix+1] -= 2*dx
1334            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1335        covMatrix = covData['covMatrix']
1336        varyList = covData['varyList']
1337        AngVcov = getVCov(names,varyList,covMatrix)
1338        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1339        if sig < 0.01:
1340            sig = -0.01
1341   
1342    return Angle,sig
1343
1344def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
1345    '''default doc string
1346   
1347    :param type name: description
1348   
1349    :returns: type name: description
1350   
1351    '''
1352
1353    def calcTorsion(Atoms,SyOps,Amat):
1354       
1355        XYZ = []
1356        for i,atom in enumerate(Atoms):
1357            Inv,M,T,C,U = SyOps[i]
1358            XYZ.append(np.array(atom[1:4]))
1359            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1360            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1361        V1 = XYZ[1]-XYZ[0]
1362        V2 = XYZ[2]-XYZ[1]
1363        V3 = XYZ[3]-XYZ[2]
1364        V1 /= np.sqrt(np.sum(V1**2))
1365        V2 /= np.sqrt(np.sum(V2**2))
1366        V3 /= np.sqrt(np.sum(V3**2))
1367        M = np.array([V1,V2,V3])
1368        D = nl.det(M)
1369        Ang = 1.0
1370        P12 = np.dot(V1,V2)
1371        P13 = np.dot(V1,V3)
1372        P23 = np.dot(V2,V3)
1373        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1374        return Tors
1375           
1376    Inv = []
1377    SyOps = []
1378    names = []
1379    for i,atom in enumerate(Oatoms):
1380        names += atom[-1]
1381        Op,unit = Atoms[i][-1]
1382        inv = Op/abs(Op)
1383        m,t = SGData['SGOps'][abs(Op)%100-1]
1384        c = SGData['SGCen'][abs(Op)/100]
1385        SyOps.append([inv,m,t,c,unit])
1386    Tors = calcTorsion(Oatoms,SyOps,Amat)
1387   
1388    sig = -0.01
1389    if 'covMatrix' in covData:
1390        parmNames = []
1391        dx = .00001
1392        dadx = np.zeros(12)
1393        for i in range(12):
1394            ia = i/3
1395            ix = i%3
1396            Oatoms[ia][ix+1] -= dx
1397            a0 = calcTorsion(Oatoms,SyOps,Amat)
1398            Oatoms[ia][ix+1] += 2*dx
1399            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1400            Oatoms[ia][ix+1] -= dx           
1401        covMatrix = covData['covMatrix']
1402        varyList = covData['varyList']
1403        TorVcov = getVCov(names,varyList,covMatrix)
1404        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
1405        if sig < 0.01:
1406            sig = -0.01
1407   
1408    return Tors,sig
1409       
1410def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
1411    '''default doc string
1412   
1413    :param type name: description
1414   
1415    :returns: type name: description
1416   
1417    '''
1418
1419    def calcDist(Atoms,SyOps,Amat):
1420        XYZ = []
1421        for i,atom in enumerate(Atoms):
1422            Inv,M,T,C,U = SyOps[i]
1423            XYZ.append(np.array(atom[1:4]))
1424            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1425            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1426        V1 = XYZ[1]-XYZ[0]
1427        return np.sqrt(np.sum(V1**2))
1428       
1429    def calcAngle(Atoms,SyOps,Amat):
1430        XYZ = []
1431        for i,atom in enumerate(Atoms):
1432            Inv,M,T,C,U = SyOps[i]
1433            XYZ.append(np.array(atom[1:4]))
1434            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1435            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1436        V1 = XYZ[1]-XYZ[0]
1437        V1 /= np.sqrt(np.sum(V1**2))
1438        V2 = XYZ[1]-XYZ[2]
1439        V2 /= np.sqrt(np.sum(V2**2))
1440        V3 = V2-V1
1441        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1442        return acosd(cang)
1443
1444    def calcTorsion(Atoms,SyOps,Amat):
1445       
1446        XYZ = []
1447        for i,atom in enumerate(Atoms):
1448            Inv,M,T,C,U = SyOps[i]
1449            XYZ.append(np.array(atom[1:4]))
1450            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1451            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1452        V1 = XYZ[1]-XYZ[0]
1453        V2 = XYZ[2]-XYZ[1]
1454        V3 = XYZ[3]-XYZ[2]
1455        V1 /= np.sqrt(np.sum(V1**2))
1456        V2 /= np.sqrt(np.sum(V2**2))
1457        V3 /= np.sqrt(np.sum(V3**2))
1458        M = np.array([V1,V2,V3])
1459        D = nl.det(M)
1460        Ang = 1.0
1461        P12 = np.dot(V1,V2)
1462        P13 = np.dot(V1,V3)
1463        P23 = np.dot(V2,V3)
1464        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1465        return Tors
1466           
1467    Inv = []
1468    SyOps = []
1469    names = []
1470    for i,atom in enumerate(Oatoms):
1471        names += atom[-1]
1472        Op,unit = Atoms[i][-1]
1473        inv = Op/abs(Op)
1474        m,t = SGData['SGOps'][abs(Op)%100-1]
1475        c = SGData['SGCen'][abs(Op)/100]
1476        SyOps.append([inv,m,t,c,unit])
1477    M = len(Oatoms)
1478    if M == 2:
1479        Val = calcDist(Oatoms,SyOps,Amat)
1480    elif M == 3:
1481        Val = calcAngle(Oatoms,SyOps,Amat)
1482    else:
1483        Val = calcTorsion(Oatoms,SyOps,Amat)
1484   
1485    sigVals = [-0.001,-0.01,-0.01]
1486    sig = sigVals[M-3]
1487    if 'covMatrix' in covData:
1488        parmNames = []
1489        dx = .00001
1490        N = M*3
1491        dadx = np.zeros(N)
1492        for i in range(N):
1493            ia = i/3
1494            ix = i%3
1495            Oatoms[ia][ix+1] += dx
1496            if M == 2:
1497                a0 = calcDist(Oatoms,SyOps,Amat)
1498            elif M == 3:
1499                a0 = calcAngle(Oatoms,SyOps,Amat)
1500            else:
1501                a0 = calcTorsion(Oatoms,SyOps,Amat)
1502            Oatoms[ia][ix+1] -= 2*dx
1503            if M == 2:
1504                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1505            elif M == 3:
1506                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1507            else:
1508                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1509        covMatrix = covData['covMatrix']
1510        varyList = covData['varyList']
1511        Vcov = getVCov(names,varyList,covMatrix)
1512        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
1513        if sig < sigVals[M-3]:
1514            sig = sigVals[M-3]
1515   
1516    return Val,sig
1517       
1518def ValEsd(value,esd=0,nTZ=False):
1519    '''Format a floating point number with a given level of precision or
1520    with in crystallographic format with a "esd", as value(esd). If esd is
1521    negative the number is formatted with the level of significant figures
1522    appropriate if abs(esd) were the esd, but the esd is not included.
1523    if the esd is zero, approximately 6 significant figures are printed.
1524    nTZ=True causes "extra" zeros to be removed after the decimal place.
1525    for example:
1526
1527      * "1.235(3)" for value=1.2346 & esd=0.003
1528      * "1.235(3)e4" for value=12346. & esd=30
1529      * "1.235(3)e6" for value=0.12346e7 & esd=3000
1530      * "1.235" for value=1.2346 & esd=-0.003
1531      * "1.240" for value=1.2395 & esd=-0.003
1532      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
1533      * "1.23460" for value=1.2346 & esd=0.0
1534
1535    :param float value: number to be formatted
1536    :param float esd: uncertainty or if esd < 0, specifies level of
1537      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
1538    :param bool nTZ: True to remove trailing zeros (default is False)
1539    :returns: value(esd) or value as a string
1540
1541    '''
1542    # Note: this routine is Python 3 compatible -- I think
1543    cutoff = 3.16228    #=(sqrt(10); same as old GSAS   was 1.95
1544    if math.isnan(value): # invalid value, bail out
1545        return '?'
1546    if math.isnan(esd): # invalid esd, treat as zero
1547        esd = 0
1548        esdoff = 5
1549    elif esd != 0:
1550        # transform the esd to a one or two digit integer
1551        l = math.log10(abs(esd)) % 1.
1552        if l < math.log10(cutoff): l+= 1.       
1553        intesd = int(round(10**l)) # esd as integer
1554        # determine the number of digits offset for the esd
1555        esdoff = int(round(math.log10(intesd*1./abs(esd))))
1556    else:
1557        esdoff = 5
1558    valoff = 0
1559    if abs(value) < abs(esdoff): # value is effectively zero
1560        pass
1561    elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
1562        # where the digit offset is to the left of the decimal place or where too many
1563        # digits are needed
1564        if abs(value) > 1:
1565            valoff = int(math.log10(abs(value)))
1566        elif abs(value) > 0:
1567            valoff = int(math.log10(abs(value))-0.9999999)
1568        else:
1569            valoff = 0
1570    if esd != 0:
1571        if valoff+esdoff < 0:
1572            valoff = esdoff = 0
1573        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
1574    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
1575        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
1576    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
1577        if abs(value) > 0:           
1578            extra = -math.log10(abs(value))
1579        else:
1580            extra = 0
1581        if extra > 0: extra += 1
1582        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
1583    if esd > 0:
1584        out += ("({:d})").format(intesd)  # add the esd
1585    elif nTZ and '.' in out:
1586        out = out.rstrip('0')  # strip zeros to right of decimal
1587        out = out.rstrip('.')  # and decimal place when not needed
1588    if valoff != 0:
1589        out += ("e{:d}").format(valoff) # add an exponent, when needed
1590    return out
1591   
1592################################################################################
1593##### Texture fitting stuff
1594################################################################################
1595
1596def FitTexture(General,Gangls,refData,keyList):
1597    import pytexture as ptx
1598    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
1599   
1600    def printSpHarm(textureData,SHtextureSig):
1601        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1602        names = ['omega','chi','phi']
1603        namstr = '  names :'
1604        ptstr =  '  values:'
1605        sigstr = '  esds  :'
1606        for name in names:
1607            namstr += '%12s'%('Sample '+name)
1608            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1609            if 'Sample '+name in SHtextureSig:
1610                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1611            else:
1612                sigstr += 12*' '
1613        print namstr
1614        print ptstr
1615        print sigstr
1616        print '\n Texture coefficients:'
1617        SHcoeff = textureData['SH Coeff'][1]
1618        SHkeys = SHcoeff.keys()
1619        nCoeff = len(SHcoeff)
1620        nBlock = nCoeff/10+1
1621        iBeg = 0
1622        iFin = min(iBeg+10,nCoeff)
1623        for block in range(nBlock):
1624            namstr = '  names :'
1625            ptstr =  '  values:'
1626            sigstr = '  esds  :'
1627            for name in SHkeys[iBeg:iFin]:
1628                if 'C' in name:
1629                    namstr += '%12s'%(name)
1630                    ptstr += '%12.3f'%(SHcoeff[name])
1631                    if name in SHtextureSig:
1632                        sigstr += '%12.3f'%(SHtextureSig[name])
1633                    else:
1634                        sigstr += 12*' '
1635            print namstr
1636            print ptstr
1637            print sigstr
1638            iBeg += 10
1639            iFin = min(iBeg+10,nCoeff)
1640           
1641    def Dict2Values(parmdict, varylist):
1642        '''Use before call to leastsq to setup list of values for the parameters
1643        in parmdict, as selected by key in varylist'''
1644        return [parmdict[key] for key in varylist] 
1645       
1646    def Values2Dict(parmdict, varylist, values):
1647        ''' Use after call to leastsq to update the parameter dictionary with
1648        values corresponding to keys in varylist'''
1649        parmdict.update(zip(varylist,values))
1650       
1651    def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList):
1652        parmDict.update(zip(varyList,values))
1653        Mat = np.empty(0)
1654        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
1655        for hist in Gangls.keys():
1656            mat = []
1657            for ref in refData[hist]:
1658#                wt = np.sqrt(ref[4])
1659                wt = 1.
1660                ref[6] = 1.
1661                H = ref[:3]
1662                phi,beta = G2lat.CrsAng(H,cell,SGData)
1663                psi,gam,x,x = G2lat.SamAng(ref[3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
1664                for item in parmDict:
1665                    if 'C' in item:
1666                        L,M,N = eval(item.strip('C'))
1667                        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1668                        Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
1669                        Lnorm = G2lat.Lnorm(L)
1670                        ref[6] += parmDict[item]*Lnorm*Kcl*Ksl
1671                mat.append(wt*(ref[5]-ref[6]))
1672            Mat = np.concatenate((Mat,np.array(mat)))
1673        return Mat
1674       
1675    def errSpHarm2(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList):
1676        parmDict.update(zip(varyList,values))
1677        Mat = np.empty(0)
1678        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
1679        for hist in Gangls.keys():
1680            Refs = refData[hist]
1681#            wt = np.sqrt(ref[4])
1682            wt = 1.
1683            Refs[:,6] = 1.
1684            H = Refs[:,:3]
1685            phi,beta = G2lat.CrsAng(H,cell,SGData)
1686            psi,gam,x,x = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
1687            for item in parmDict:
1688                if 'C' in item:
1689                    L,M,N = eval(item.strip('C'))
1690                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1691                    Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
1692                    Lnorm = G2lat.Lnorm(L)
1693                    Refs[:,6] += parmDict[item]*Lnorm*Kcl*Ksl
1694            mat = wt*(Refs[:,5]-Refs[:,6])
1695            Mat = np.concatenate((Mat,mat))
1696        return Mat
1697       
1698    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList):
1699        Mat = np.empty(0)
1700        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
1701        for hist in Gangls.keys():
1702            mat = np.zeros((len(refData[hist]),len(varyList)))
1703            for i,ref in enumerate(refData[hist]):
1704#                wt = np.sqrt(ref[4])
1705                wt = 1.
1706                H = ref[:3]
1707                phi,beta = G2lat.CrsAng(H,cell,SGData)
1708                psi,gam,dPdA,dGdA = G2lat.SamAng(ref[3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
1709                for j,item in enumerate(varyList):
1710                    if 'C' in item:
1711                        L,M,N = eval(item.strip('C'))
1712                        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1713                        Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
1714                        Lnorm = G2lat.Lnorm(L)
1715                        mat[i,j] = -wt*Lnorm*Kcl*Ksl
1716                        for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
1717                            try:
1718                                l = varyList.index(itema)
1719                                mat[i,l] += parmDict[item]*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
1720                            except ValueError:
1721                                pass
1722            if len(Mat):
1723                Mat = np.concatenate((Mat,mat))
1724            else:
1725                Mat = mat
1726        return Mat
1727
1728    print ' Fit texture for '+General['Name']
1729    SGData = General['SGData']
1730    cell = General['Cell'][1:7]
1731    Texture = General['SH Texture']
1732    if not Texture['Order']:
1733        return 'No spherical harmonics coefficients'
1734    varyList = []
1735    parmDict = copy.copy(Texture['SH Coeff'][1])
1736    for item in ['Sample omega','Sample chi','Sample phi']:
1737        parmDict[item] = Texture[item][1]
1738        if Texture[item][0]:
1739            varyList.append(item)
1740    if Texture['SH Coeff'][0]:
1741        varyList += Texture['SH Coeff'][1].keys()
1742    while True:
1743        begin = time.time()
1744        values =  np.array(Dict2Values(parmDict, varyList))
1745        result = so.leastsq(errSpHarm2,values,Dfun=dervSpHarm,full_output=True,
1746            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList))
1747        ncyc = int(result[2]['nfev']/2)
1748        if ncyc:
1749            runtime = time.time()-begin   
1750            chisq = np.sum(result[2]['fvec']**2)
1751            Values2Dict(parmDict, varyList, result[0])
1752            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
1753            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(result[2]['fvec']),' Number of parameters: ',len(varyList)
1754            print 'refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1755            print 'chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)
1756            try:
1757                sig = np.sqrt(np.diag(result[1])*GOF)
1758                if np.any(np.isnan(sig)):
1759                    print '*** Least squares aborted - some invalid esds possible ***'
1760                break                   #refinement succeeded - finish up!
1761            except ValueError:          #result[1] is None on singular matrix
1762                print '**** Refinement failed - singular matrix ****'
1763                return None
1764        else:
1765            break
1766   
1767    for hist in keyList:
1768        print ' Texture corrections for '+hist
1769        for ref in refData[hist]:
1770            print ' %d %d %d %.3f %.3f'%(int(ref[0]),int(ref[1]),int(ref[2]),ref[5],ref[6])
1771    if ncyc:
1772        for parm in parmDict:
1773            if 'C' in parm:
1774                Texture['SH Coeff'][1][parm] = parmDict[parm]
1775            else:
1776                Texture[parm][1] = parmDict[parm] 
1777        sigDict = dict(zip(varyList,sig))
1778        printSpHarm(Texture,sigDict)
1779       
1780    return None
1781   
1782################################################################################
1783##### Fourier & charge flip stuff
1784################################################################################
1785
1786def adjHKLmax(SGData,Hmax):
1787    '''default doc string
1788   
1789    :param type name: description
1790   
1791    :returns: type name: description
1792   
1793    '''
1794    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1795        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
1796        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
1797        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
1798    else:
1799        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
1800        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
1801        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
1802
1803def OmitMap(data,reflDict,pgbar=None):
1804    '''default doc string
1805   
1806    :param type name: description
1807   
1808    :returns: type name: description
1809   
1810    '''
1811    generalData = data['General']
1812    if not generalData['Map']['MapType']:
1813        print '**** ERROR - Fourier map not defined'
1814        return
1815    mapData = generalData['Map']
1816    dmin = mapData['Resolution']
1817    SGData = generalData['SGData']
1818    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1819    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1820    cell = generalData['Cell'][1:8]       
1821    A = G2lat.cell2A(cell[:6])
1822    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1823    adjHKLmax(SGData,Hmax)
1824    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1825    time0 = time.time()
1826    for iref,ref in enumerate(reflDict['RefList']):
1827        if ref[4] >= dmin:
1828            Fosq,Fcsq,ph = ref[8:11]
1829            Uniq = np.inner(ref[:3],SGMT)
1830            Phi = np.inner(ref[:3],SGT)
1831            for i,hkl in enumerate(Uniq):        #uses uniq
1832                hkl = np.asarray(hkl,dtype='i')
1833                dp = 360.*Phi[i]                #and phi
1834                a = cosd(ph+dp)
1835                b = sind(ph+dp)
1836                phasep = complex(a,b)
1837                phasem = complex(a,-b)
1838                Fo = np.sqrt(Fosq)
1839                if '2Fo-Fc' in mapData['MapType']:
1840                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
1841                else:
1842                    F = np.sqrt(Fosq)
1843                h,k,l = hkl+Hmax
1844                Fhkl[h,k,l] = F*phasep
1845                h,k,l = -hkl+Hmax
1846                Fhkl[h,k,l] = F*phasem
1847    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1848    M = np.mgrid[0:4,0:4,0:4]
1849    blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
1850    iBeg = blkIds*rho0.shape/4
1851    iFin = (blkIds+1)*rho0.shape/4
1852    rho_omit = np.zeros_like(rho0)
1853    nBlk = 0
1854    for iB,iF in zip(iBeg,iFin):
1855        rho1 = np.copy(rho0)
1856        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
1857        Fnew = fft.ifftshift(fft.ifftn(rho1))
1858        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
1859        phase = Fnew/np.absolute(Fnew)
1860        OFhkl = np.absolute(Fhkl)*phase
1861        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
1862        rho_omit[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = np.copy(rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]])
1863        nBlk += 1
1864        pgbar.Update(nBlk)
1865    mapData['rho'] = np.real(rho_omit)/cell[6]
1866    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1867    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
1868    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1869    return mapData
1870   
1871def Fourier4DMap(data,reflDict):
1872    '''default doc string
1873   
1874    :param type name: description
1875   
1876    :returns: type name: description
1877   
1878    '''
1879    generalData = data['General']
1880    mapData = generalData['4DmapData']
1881    dmin = mapData['Resolution']
1882    SGData = generalData['SGData']
1883    SSGData = generalData['SSGData']
1884    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1885    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1886    cell = generalData['Cell'][1:8]       
1887    A = G2lat.cell2A(cell[:6])
1888    maxM = generalData['SuperVec'][2]
1889    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
1890    adjHKLmax(SGData,Hmax)
1891    Hmax = np.asarray(Hmax,dtype='i')+1
1892    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1893    time0 = time.time()
1894    for iref,ref in enumerate(reflDict['RefList']):
1895        if ref[5] > dmin:
1896            Fosq,Fcsq,ph = ref[9:12]
1897            Uniq = np.inner(ref[:4],SSGMT)
1898            Phi = np.inner(ref[:4],SSGT)
1899            for i,hkl in enumerate(Uniq):        #uses uniq
1900                hkl = np.asarray(hkl,dtype='i')
1901                dp = 360.*Phi[i]                #and phi
1902                a = cosd(ph+dp)
1903                b = sind(ph+dp)
1904                phasep = complex(a,b)
1905                phasem = complex(a,-b)
1906                if 'Fobs' in mapData['MapType']:
1907                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1908                    h,k,l,m = hkl+Hmax
1909                    Fhkl[h,k,l,m] = F*phasep
1910                    h,k,l,m = -hkl+Hmax
1911                    Fhkl[h,k,l,m] = F*phasem
1912                elif 'delt-F' in mapData['MapType']:
1913                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1914                    h,k,l,m = hkl+Hmax
1915                    Fhkl[h,k,l,m] = dF*phasep
1916                    h,k,l,m = -hkl+Hmax
1917                    Fhkl[h,k,l,m] = dF*phasem
1918    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1919    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1920    mapData['Type'] = reflDict['Type']
1921    mapData['rho'] = np.real(rho)
1922    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1923    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
1924    return mapData
1925
1926def FourierMap(data,reflDict):
1927    '''default doc string
1928   
1929    :param type name: description
1930   
1931    :returns: type name: description
1932   
1933    '''
1934    generalData = data['General']
1935    mapData = generalData['Map']
1936    dmin = mapData['Resolution']
1937    SGData = generalData['SGData']
1938    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1939    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1940    cell = generalData['Cell'][1:8]       
1941    A = G2lat.cell2A(cell[:6])
1942    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1943    adjHKLmax(SGData,Hmax)
1944    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1945#    Fhkl[0,0,0] = generalData['F000X']
1946    time0 = time.time()
1947    for iref,ref in enumerate(reflDict['RefList']):
1948        if ref[4] > dmin:
1949            Fosq,Fcsq,ph = ref[8:11]
1950            Uniq = np.inner(ref[:3],SGMT)
1951            Phi = np.inner(ref[:3],SGT)
1952            for i,hkl in enumerate(Uniq):        #uses uniq
1953                hkl = np.asarray(hkl,dtype='i')
1954                dp = 360.*Phi[i]                #and phi
1955                a = cosd(ph)
1956                b = sind(ph)
1957                phasep = complex(a,b)+dp
1958                phasem = complex(a,-b)-dp
1959                if 'Fobs' in mapData['MapType']:
1960                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1961                    h,k,l = hkl+Hmax
1962                    Fhkl[h,k,l] = F*phasep
1963                    h,k,l = -hkl+Hmax
1964                    Fhkl[h,k,l] = F*phasem
1965                elif 'Fcalc' in mapData['MapType']:
1966                    F = np.sqrt(Fcsq)
1967                    h,k,l = hkl+Hmax
1968                    Fhkl[h,k,l] = F*phasep
1969                    h,k,l = -hkl+Hmax
1970                    Fhkl[h,k,l] = F*phasem
1971                elif 'delt-F' in mapData['MapType']:
1972                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1973                    h,k,l = hkl+Hmax
1974                    Fhkl[h,k,l] = dF*phasep
1975                    h,k,l = -hkl+Hmax
1976                    Fhkl[h,k,l] = dF*phasem
1977                elif '2*Fo-Fc' in mapData['MapType']:
1978                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1979                    h,k,l = hkl+Hmax
1980                    Fhkl[h,k,l] = F*phasep
1981                    h,k,l = -hkl+Hmax
1982                    Fhkl[h,k,l] = F*phasem
1983                elif 'Patterson' in mapData['MapType']:
1984                    h,k,l = hkl+Hmax
1985                    Fhkl[h,k,l] = complex(Fosq,0.)
1986                    h,k,l = -hkl+Hmax
1987                    Fhkl[h,k,l] = complex(Fosq,0.)
1988    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1989    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1990    mapData['Type'] = reflDict['Type']
1991    mapData['rho'] = np.real(rho)
1992    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1993    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
1994    return mapData
1995   
1996# map printing for testing purposes
1997def printRho(SGLaue,rho,rhoMax):                         
1998    '''default doc string
1999   
2000    :param type name: description
2001   
2002    :returns: type name: description
2003   
2004    '''
2005    dim = len(rho.shape)
2006    if dim == 2:
2007        ix,jy = rho.shape
2008        for j in range(jy):
2009            line = ''
2010            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
2011                line += (jy-j)*'  '
2012            for i in range(ix):
2013                r = int(100*rho[i,j]/rhoMax)
2014                line += '%4d'%(r)
2015            print line+'\n'
2016    else:
2017        ix,jy,kz = rho.shape
2018        for k in range(kz):
2019            print 'k = ',k
2020            for j in range(jy):
2021                line = ''
2022                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
2023                    line += (jy-j)*'  '
2024                for i in range(ix):
2025                    r = int(100*rho[i,j,k]/rhoMax)
2026                    line += '%4d'%(r)
2027                print line+'\n'
2028## keep this
2029               
2030def findOffset(SGData,A,Fhkl):   
2031    '''default doc string
2032   
2033    :param type name: description
2034   
2035    :returns: type name: description
2036   
2037    '''
2038    if SGData['SpGrp'] == 'P 1':
2039        return [0,0,0]   
2040    hklShape = Fhkl.shape
2041    hklHalf = np.array(hklShape)/2
2042    sortHKL = np.argsort(Fhkl.flatten())
2043    Fdict = {}
2044    for hkl in sortHKL:
2045        HKL = np.unravel_index(hkl,hklShape)
2046        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
2047        if F == 0.:
2048            break
2049        Fdict['%.6f'%(np.absolute(F))] = hkl
2050    Flist = np.flipud(np.sort(Fdict.keys()))
2051    F = str(1.e6)
2052    i = 0
2053    DH = []
2054    Dphi = []
2055    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2056    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2057    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
2058    for F in Flist:
2059        hkl = np.unravel_index(Fdict[F],hklShape)
2060        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
2061            continue
2062        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
2063        Uniq = np.array(Uniq,dtype='i')
2064        Phi = np.array(Phi)
2065        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
2066        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
2067        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
2068        ang0 = np.angle(Fh0,deg=True)/360.
2069        for H,phi in zip(Uniq,Phi)[1:]:
2070            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
2071            dH = H-hkl
2072            dang = ang-ang0
2073            DH.append(dH)
2074            Dphi.append((dang+.5) % 1.0)
2075        if i > 20 or len(DH) > 30:
2076            break
2077        i += 1
2078    DH = np.array(DH)
2079    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
2080    Dphi = np.array(Dphi)
2081    steps = np.array(hklShape)
2082    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
2083    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
2084    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
2085    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
2086    hist,bins = np.histogram(Mmap,bins=1000)
2087#    for i,item in enumerate(hist[:10]):
2088#        print item,bins[i]
2089    chisq = np.min(Mmap)
2090    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
2091    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
2092#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
2093    return DX
2094   
2095def ChargeFlip(data,reflDict,pgbar):
2096    '''default doc string
2097   
2098    :param type name: description
2099   
2100    :returns: type name: description
2101   
2102    '''
2103    generalData = data['General']
2104    mapData = generalData['Map']
2105    flipData = generalData['Flip']
2106    FFtable = {}
2107    if 'None' not in flipData['Norm element']:
2108        normElem = flipData['Norm element'].upper()
2109        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
2110        for ff in FFs:
2111            if ff['Symbol'] == normElem:
2112                FFtable.update(ff)
2113    dmin = flipData['Resolution']
2114    SGData = generalData['SGData']
2115    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2116    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2117    cell = generalData['Cell'][1:8]       
2118    A = G2lat.cell2A(cell[:6])
2119    Vol = cell[6]
2120    im = 0
2121    if generalData['Type'] in ['modulated','magnetic',]:
2122        im = 1
2123    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2124    adjHKLmax(SGData,Hmax)
2125    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
2126    time0 = time.time()
2127    for iref,ref in enumerate(reflDict['RefList']):
2128        dsp = ref[4+im]
2129        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
2130            continue
2131        if dsp > dmin:
2132            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
2133            if FFtable:
2134                SQ = 0.25/dsp**2
2135                ff *= G2el.ScatFac(FFtable,SQ)[0]
2136            if ref[8+im] > 0.:         #use only +ve Fobs**2
2137                E = np.sqrt(ref[8+im])/ff
2138            else:
2139                E = 0.
2140            ph = ref[10]
2141            ph = rn.uniform(0.,360.)
2142            Uniq = np.inner(ref[:3],SGMT)
2143            Phi = np.inner(ref[:3],SGT)
2144            for i,hkl in enumerate(Uniq):        #uses uniq
2145                hkl = np.asarray(hkl,dtype='i')
2146                dp = 360.*Phi[i]                #and phi
2147                a = cosd(ph+dp)
2148                b = sind(ph+dp)
2149                phasep = complex(a,b)
2150                phasem = complex(a,-b)
2151                h,k,l = hkl+Hmax
2152                Ehkl[h,k,l] = E*phasep
2153                h,k,l = -hkl+Hmax       #Friedel pair refl.
2154                Ehkl[h,k,l] = E*phasem
2155#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
2156    CEhkl = copy.copy(Ehkl)
2157    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
2158    Emask = ma.getmask(MEhkl)
2159    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
2160    Ncyc = 0
2161    old = np.seterr(all='raise')
2162    while True:       
2163        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
2164        CEsig = np.std(CErho)
2165        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
2166        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
2167        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
2168        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
2169        phase = CFhkl/np.absolute(CFhkl)
2170        CEhkl = np.absolute(Ehkl)*phase
2171        Ncyc += 1
2172        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
2173        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
2174        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
2175        if Rcf < 5.:
2176            break
2177        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
2178        if not GoOn or Ncyc > 10000:
2179            break
2180    np.seterr(**old)
2181    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
2182    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
2183    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
2184    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
2185       
2186    mapData['Rcf'] = Rcf
2187    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2188    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2189    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2190    mapData['Type'] = reflDict['Type']
2191    return mapData
2192   
2193def findSSOffset(SGData,SSGData,A,Fhklm):   
2194    '''default doc string
2195   
2196    :param type name: description
2197   
2198    :returns: type name: description
2199   
2200    '''
2201    if SGData['SpGrp'] == 'P 1':
2202        return [0,0,0,0]   
2203    hklmShape = Fhklm.shape
2204    hklmHalf = np.array(hklmShape)/2
2205    sortHKLM = np.argsort(Fhklm.flatten())
2206    Fdict = {}
2207    for hklm in sortHKLM:
2208        HKLM = np.unravel_index(hklm,hklmShape)
2209        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
2210        if F == 0.:
2211            break
2212        Fdict['%.6f'%(np.absolute(F))] = hklm
2213    Flist = np.flipud(np.sort(Fdict.keys()))
2214    F = str(1.e6)
2215    i = 0
2216    DH = []
2217    Dphi = []
2218    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2219    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2220    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
2221    for F in Flist:
2222        hklm = np.unravel_index(Fdict[F],hklmShape)
2223        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
2224            continue
2225        Uniq = np.inner(hklm-hklmHalf,SSGMT)
2226        Phi = np.inner(hklm-hklmHalf,SSGT)
2227        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
2228        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
2229        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
2230        ang0 = np.angle(Fh0,deg=True)/360.
2231        for H,phi in zip(Uniq,Phi)[1:]:
2232            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
2233            dH = H-hklm
2234            dang = ang-ang0
2235            DH.append(dH)
2236            Dphi.append((dang+.5) % 1.0)
2237        if i > 20 or len(DH) > 30:
2238            break
2239        i += 1
2240    DH = np.array(DH)
2241    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
2242    Dphi = np.array(Dphi)
2243    steps = np.array(hklmShape)
2244    X,Y,Z,T = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2],0:1:1./steps[3]]
2245    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
2246    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
2247    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
2248    hist,bins = np.histogram(Mmap,bins=1000)
2249#    for i,item in enumerate(hist[:10]):
2250#        print item,bins[i]
2251    chisq = np.min(Mmap)
2252    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
2253    print ' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
2254#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
2255    return DX
2256   
2257def SSChargeFlip(data,reflDict,pgbar):
2258    '''default doc string
2259   
2260    :param type name: description
2261   
2262    :returns: type name: description
2263   
2264    '''
2265    generalData = data['General']
2266    mapData = generalData['Map']
2267    map4DData = {}
2268    flipData = generalData['Flip']
2269    FFtable = {}
2270    if 'None' not in flipData['Norm element']:
2271        normElem = flipData['Norm element'].upper()
2272        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
2273        for ff in FFs:
2274            if ff['Symbol'] == normElem:
2275                FFtable.update(ff)
2276    dmin = flipData['Resolution']
2277    SGData = generalData['SGData']
2278    SSGData = generalData['SSGData']
2279    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2280    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2281    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2282    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2283    cell = generalData['Cell'][1:8]       
2284    A = G2lat.cell2A(cell[:6])
2285    Vol = cell[6]
2286    maxM = generalData['SuperVec'][2]
2287    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
2288    adjHKLmax(SGData,Hmax)
2289    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
2290    time0 = time.time()
2291    for iref,ref in enumerate(reflDict['RefList']):
2292        dsp = ref[5]
2293        if dsp > dmin:
2294            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
2295            if FFtable:
2296                SQ = 0.25/dsp**2
2297                ff *= G2el.ScatFac(FFtable,SQ)[0]
2298            if ref[9] > 0.:         #use only +ve Fobs**2
2299                E = np.sqrt(ref[9])/ff
2300            else:
2301                E = 0.
2302            ph = ref[11]
2303            ph = rn.uniform(0.,360.)
2304            Uniq = np.inner(ref[:4],SSGMT)
2305            Phi = np.inner(ref[:4],SSGT)
2306            for i,hklm in enumerate(Uniq):        #uses uniq
2307                hklm = np.asarray(hklm,dtype='i')
2308                dp = 360.*Phi[i]                #and phi
2309                a = cosd(ph+dp)
2310                b = sind(ph+dp)
2311                phasep = complex(a,b)
2312                phasem = complex(a,-b)
2313                h,k,l,m = hklm+Hmax
2314                Ehkl[h,k,l,m] = E*phasep
2315                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
2316                Ehkl[h,k,l,m] = E*phasem
2317#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
2318    CEhkl = copy.copy(Ehkl)
2319    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
2320    Emask = ma.getmask(MEhkl)
2321    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
2322    Ncyc = 0
2323    old = np.seterr(all='raise')
2324    while True:       
2325        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
2326        CEsig = np.std(CErho)
2327        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
2328        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
2329        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
2330        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
2331        phase = CFhkl/np.absolute(CFhkl)
2332        CEhkl = np.absolute(Ehkl)*phase
2333        Ncyc += 1
2334        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
2335        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
2336        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
2337        if Rcf < 5.:
2338            break
2339        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
2340        if not GoOn or Ncyc > 10000:
2341            break
2342    np.seterr(**old)
2343    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
2344    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))
2345    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))
2346    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
2347    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
2348
2349    mapData['Rcf'] = Rcf
2350    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2351    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2352    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2353    mapData['Type'] = reflDict['Type']
2354
2355    map4DData['Rcf'] = Rcf
2356    map4DData['rho'] = np.real(np.roll(np.roll(np.roll(np.roll(SSrho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2),roll[3],axis=3))
2357    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
2358    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
2359    map4DData['Type'] = reflDict['Type']
2360    return mapData,map4DData
2361   
2362def SearchMap(generalData,drawingData,Neg=False):
2363    '''Does a search of a density map for peaks meeting the criterion of peak
2364    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
2365    mapData is data['General']['mapData']; the map is also in mapData.
2366
2367    :param generalData: the phase data structure; includes the map
2368    :param drawingData: the drawing data structure
2369    :param Neg:  if True then search for negative peaks (i.e. H-atoms & neutron data)
2370
2371    :returns: (peaks,mags,dzeros) where
2372
2373        * peaks : ndarray
2374          x,y,z positions of the peaks found in the map
2375        * mags : ndarray
2376          the magnitudes of the peaks
2377        * dzeros : ndarray
2378          the distance of the peaks from  the unit cell origin
2379        * dcent : ndarray
2380          the distance of the peaks from  the unit cell center
2381
2382    '''       
2383    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2384   
2385    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
2386   
2387#    def noDuplicate(xyz,peaks,Amat):
2388#        XYZ = np.inner(Amat,xyz)
2389#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2390#            print ' Peak',xyz,' <0.5A from another peak'
2391#            return False
2392#        return True
2393#                           
2394    def fixSpecialPos(xyz,SGData,Amat):
2395        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
2396        X = []
2397        xyzs = [equiv[0] for equiv in equivs]
2398        for x in xyzs:
2399            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
2400                X.append(x)
2401        if len(X) > 1:
2402            return np.average(X,axis=0)
2403        else:
2404            return xyz
2405       
2406    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
2407        Mag,x0,y0,z0,sig = parms
2408        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
2409#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
2410        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
2411       
2412    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
2413        Mag,x0,y0,z0,sig = parms
2414        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2415        return M
2416       
2417    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
2418        Mag,x0,y0,z0,sig = parms
2419        dMdv = np.zeros(([5,]+list(rX.shape)))
2420        delt = .01
2421        for i in range(5):
2422            parms[i] -= delt
2423            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2424            parms[i] += 2.*delt
2425            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2426            parms[i] -= delt
2427            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
2428        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2429        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
2430        dMdv = np.reshape(dMdv,(5,rX.size))
2431        Hess = np.inner(dMdv,dMdv)
2432       
2433        return Vec,Hess
2434       
2435    phaseName = generalData['Name']
2436    SGData = generalData['SGData']
2437    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2438    peaks = []
2439    mags = []
2440    dzeros = []
2441    dcent = []
2442    try:
2443        mapData = generalData['Map']
2444        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
2445        if Neg:
2446            rho = -copy.copy(mapData['rho'])     #flip +/-
2447        else:
2448            rho = copy.copy(mapData['rho'])     #don't mess up original
2449        mapHalf = np.array(rho.shape)/2
2450        res = mapData['Resolution']
2451        incre = np.array(rho.shape,dtype=np.float)
2452        step = max(1.0,1./res)+1
2453        steps = np.array(3*[step,])
2454    except KeyError:
2455        print '**** ERROR - Fourier map not defined'
2456        return peaks,mags
2457    rhoMask = ma.array(rho,mask=(rho<contLevel))
2458    indices = (-1,0,1)
2459    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2460    for roll in rolls:
2461        if np.any(roll):
2462            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
2463    indx = np.transpose(rhoMask.nonzero())
2464    peaks = indx/incre
2465    mags = rhoMask[rhoMask.nonzero()]
2466    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
2467        rho = rollMap(rho,ind)
2468        rMM = mapHalf-steps
2469        rMP = mapHalf+steps+1
2470        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2471        peakInt = np.sum(rhoPeak)*res**3
2472        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2473        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
2474        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
2475            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
2476        x1 = result[0]
2477        if not np.any(x1 < 0):
2478            mag = x1[0]
2479            peak = (np.array(x1[1:4])-ind)/incre
2480        peak = fixSpecialPos(peak,SGData,Amat)
2481        rho = rollMap(rho,-ind)
2482    cent = np.ones(3)*.5     
2483    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
2484    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
2485    if Neg:     #want negative magnitudes for negative peaks
2486        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
2487    else:
2488        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
2489   
2490def sortArray(data,pos,reverse=False):
2491    '''data is a list of items
2492    sort by pos in list; reverse if True
2493    '''
2494    T = []
2495    for i,M in enumerate(data):
2496        try:
2497            T.append((M[pos],i))
2498        except IndexError:
2499            return data
2500    D = dict(zip(T,data))
2501    T.sort()
2502    if reverse:
2503        T.reverse()
2504    X = []
2505    for key in T:
2506        X.append(D[key])
2507    return X
2508
2509def PeaksEquiv(data,Ind):
2510    '''Find the equivalent map peaks for those selected. Works on the
2511    contents of data['Map Peaks'].
2512
2513    :param data: the phase data structure
2514    :param list Ind: list of selected peak indices
2515    :returns: augmented list of peaks including those related by symmetry to the
2516      ones in Ind
2517
2518    '''       
2519    def Duplicate(xyz,peaks,Amat):
2520        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2521            return True
2522        return False
2523                           
2524    generalData = data['General']
2525    cell = generalData['Cell'][1:7]
2526    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2527    A = G2lat.cell2A(cell)
2528    SGData = generalData['SGData']
2529    mapPeaks = data['Map Peaks']
2530    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
2531    Indx = {}
2532    for ind in Ind:
2533        xyz = np.array(mapPeaks[ind][1:4])
2534        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
2535        for jnd,xyz in enumerate(XYZ):       
2536            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
2537    Ind = []
2538    for ind in Indx:
2539        if Indx[ind]:
2540            Ind.append(ind)
2541    return Ind
2542               
2543def PeaksUnique(data,Ind):
2544    '''Finds the symmetry unique set of peaks from those selected. Works on the
2545    contents of data['Map Peaks'].
2546
2547    :param data: the phase data structure
2548    :param list Ind: list of selected peak indices
2549    :returns: the list of symmetry unique peaks from among those given in Ind
2550
2551    '''       
2552#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
2553
2554    def noDuplicate(xyz,peaks,Amat):
2555        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2556            return False
2557        return True
2558                           
2559    generalData = data['General']
2560    cell = generalData['Cell'][1:7]
2561    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2562    A = G2lat.cell2A(cell)
2563    SGData = generalData['SGData']
2564    mapPeaks = data['Map Peaks']
2565    Indx = {}
2566    XYZ = {}
2567    for ind in Ind:
2568        XYZ[ind] = np.array(mapPeaks[ind][1:4])
2569        Indx[ind] = True
2570    for ind in Ind:
2571        if Indx[ind]:
2572            xyz = XYZ[ind]
2573            for jnd in Ind:
2574                if ind != jnd and Indx[jnd]:                       
2575                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
2576                    xyzs = np.array([equiv[0] for equiv in Equiv])
2577                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
2578    Ind = []
2579    for ind in Indx:
2580        if Indx[ind]:
2581            Ind.append(ind)
2582    return Ind
2583   
2584################################################################################
2585##### single peak fitting profile fxn stuff
2586################################################################################
2587
2588def getCWsig(ins,pos):
2589    '''get CW peak profile sigma
2590   
2591    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
2592      as values only
2593    :param float pos: 2-theta of peak
2594    :returns: float getCWsig: peak sigma
2595   
2596    '''
2597    tp = tand(pos/2.0)
2598    return ins['U']*tp**2+ins['V']*tp+ins['W']
2599   
2600def getCWsigDeriv(pos):
2601    '''get derivatives of CW peak profile sigma wrt U,V, & W
2602   
2603    :param float pos: 2-theta of peak
2604   
2605    :returns: list getCWsigDeriv: d(sig)/dU, d(sig)/dV & d(sig)/dW
2606   
2607    '''
2608    tp = tand(pos/2.0)
2609    return tp**2,tp,1.0
2610   
2611def getCWgam(ins,pos):
2612    '''get CW peak profile gamma
2613   
2614    :param dict ins: instrument parameters with at least 'X' & 'Y'
2615      as values only
2616    :param float pos: 2-theta of peak
2617    :returns: float getCWgam: peak gamma
2618   
2619    '''
2620    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
2621   
2622def getCWgamDeriv(pos):
2623    '''get derivatives of CW peak profile gamma wrt X & Y
2624   
2625    :param float pos: 2-theta of peak
2626   
2627    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
2628   
2629    '''
2630    return 1./cosd(pos/2.0),tand(pos/2.0)
2631   
2632def getTOFsig(ins,dsp):
2633    '''get TOF peak profile sigma
2634   
2635    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
2636      as values only
2637    :param float dsp: d-spacing of peak
2638   
2639    :returns: float getTOFsig: peak sigma
2640   
2641    '''
2642    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
2643   
2644def getTOFsigDeriv(dsp):
2645    '''get derivatives of TOF peak profile gamma wrt sig-0, sig-1, & sig-q
2646   
2647    :param float dsp: d-spacing of peak
2648   
2649    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
2650   
2651    '''
2652    return 1.0,dsp**2,dsp**4,1./dsp**2
2653   
2654def getTOFgamma(ins,dsp):
2655    '''get TOF peak profile gamma
2656   
2657    :param dict ins: instrument parameters with at least 'X' & 'Y'
2658      as values only
2659    :param float dsp: d-spacing of peak
2660   
2661    :returns: float getTOFgamma: peak gamma
2662   
2663    '''
2664    return ins['X']*dsp+ins['Y']*dsp**2
2665   
2666def getTOFgammaDeriv(dsp):
2667    '''get derivatives of TOF peak profile gamma wrt X & Y
2668   
2669    :param float dsp: d-spacing of peak
2670   
2671    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
2672   
2673    '''
2674    return dsp,dsp**2
2675   
2676def getTOFbeta(ins,dsp):
2677    '''get TOF peak profile beta
2678   
2679    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
2680      as values only
2681    :param float dsp: d-spacing of peak
2682   
2683    :returns: float getTOFbeta: peak beat
2684   
2685    '''
2686    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
2687   
2688def getTOFbetaDeriv(dsp):
2689    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
2690   
2691    :param float dsp: d-spacing of peak
2692   
2693    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
2694   
2695    '''
2696    return 1.0,1./dsp**4,1./dsp**2
2697   
2698def getTOFalpha(ins,dsp):
2699    '''get TOF peak profile alpha
2700   
2701    :param dict ins: instrument parameters with at least 'alpha'
2702      as values only
2703    :param float dsp: d-spacing of peak
2704   
2705    :returns: flaot getTOFalpha: peak alpha
2706   
2707    '''
2708    return ins['alpha']/dsp
2709   
2710def getTOFalphaDeriv(dsp):
2711    '''get derivatives of TOF peak profile beta wrt alpha
2712   
2713    :param float dsp: d-spacing of peak
2714   
2715    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
2716   
2717    '''
2718    return 1./dsp
2719   
2720def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
2721    '''set starting peak parameters for single peak fits from plot selection or auto selection
2722   
2723    :param dict Parms: instrument parameters dictionary
2724    :param dict Parms2: table lookup for TOF profile coefficients
2725    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
2726    :param float mag: peak top magnitude from pick
2727    :param bool ifQ: True if pos in Q
2728    :param bool useFit: True if use fitted CW Parms values (not defaults)
2729   
2730    :returns: list XY: peak list entry:
2731        for CW: [pos,0,mag,1,sig,0,gam,0]
2732        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2733        NB: mag refinement set by default, all others off
2734   
2735    '''
2736    ind = 0
2737    if useFit:
2738        ind = 1
2739    ins = {}
2740    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
2741        for x in ['U','V','W','X','Y']:
2742            ins[x] = Parms[x][ind]
2743        if ifQ:                              #qplot - convert back to 2-theta
2744            pos = 2.0*asind(pos*wave/(4*math.pi))
2745        sig = getCWsig(ins,pos)
2746        gam = getCWgam(ins,pos)           
2747        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
2748    else:
2749        if ifQ:
2750            dsp = 2.*np.pi/pos
2751            pos = Parms['difC']*dsp
2752        else:
2753            dsp = pos/Parms['difC'][1]
2754        if 'Pdabc' in Parms2:
2755            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
2756                ins[x] = Parms[x][ind]
2757            Pdabc = Parms2['Pdabc'].T
2758            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2759            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2760        else:
2761            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
2762                ins[x] = Parms[x][ind]
2763            alp = getTOFalpha(ins,dsp)
2764            bet = getTOFbeta(ins,dsp)
2765        sig = getTOFsig(ins,dsp)
2766        gam = getTOFgamma(ins,dsp)
2767        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2768    return XY
2769   
2770################################################################################
2771##### MC/SA stuff
2772################################################################################
2773
2774#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2775# Original Author: Travis Oliphant 2002
2776# Bug-fixes in 2006 by Tim Leslie
2777
2778
2779import numpy
2780from numpy import asarray, tan, exp, ones, squeeze, sign, \
2781     all, log, sqrt, pi, shape, array, minimum, where
2782from numpy import random
2783
2784#__all__ = ['anneal']
2785
2786_double_min = numpy.finfo(float).min
2787_double_max = numpy.finfo(float).max
2788class base_schedule(object):
2789    def __init__(self):
2790        self.dwell = 20
2791        self.learn_rate = 0.5
2792        self.lower = -10
2793        self.upper = 10
2794        self.Ninit = 50
2795        self.accepted = 0
2796        self.tests = 0
2797        self.feval = 0
2798        self.k = 0
2799        self.T = None
2800
2801    def init(self, **options):
2802        self.__dict__.update(options)
2803        self.lower = asarray(self.lower)
2804        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2805        self.upper = asarray(self.upper)
2806        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2807        self.k = 0
2808        self.accepted = 0
2809        self.feval = 0
2810        self.tests = 0
2811
2812    def getstart_temp(self, best_state):
2813        """ Find a matching starting temperature and starting parameters vector
2814        i.e. find x0 such that func(x0) = T0.
2815
2816        Parameters
2817        ----------
2818        best_state : _state
2819            A _state object to store the function value and x0 found.
2820
2821        returns
2822        -------
2823        x0 : array
2824            The starting parameters vector.
2825        """
2826
2827        assert(not self.dims is None)
2828        lrange = self.lower
2829        urange = self.upper
2830        fmax = _double_min
2831        fmin = _double_max
2832        for _ in range(self.Ninit):
2833            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2834            fval = self.func(x0, *self.args)
2835            self.feval += 1
2836            if fval > fmax:
2837                fmax = fval
2838            if fval < fmin:
2839                fmin = fval
2840                best_state.cost = fval
2841                best_state.x = array(x0)
2842
2843        self.T0 = (fmax-fmin)*1.5
2844        return best_state.x
2845       
2846    def set_range(self,x0,frac):
2847        delrange = frac*(self.upper-self.lower)
2848        self.upper = x0+delrange
2849        self.lower = x0-delrange
2850
2851    def accept_test(self, dE):
2852        T = self.T
2853        self.tests += 1
2854        if dE < 0:
2855            self.accepted += 1
2856            return 1
2857        p = exp(-dE*1.0/self.boltzmann/T)
2858        if (p > random.uniform(0.0, 1.0)):
2859            self.accepted += 1
2860            return 1
2861        return 0
2862
2863    def update_guess(self, x0):
2864        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2865
2866    def update_temp(self, x0):
2867        pass
2868
2869
2870#  A schedule due to Lester Ingber modified to use bounds - OK
2871class fast_sa(base_schedule):
2872    def init(self, **options):
2873        self.__dict__.update(options)
2874        if self.m is None:
2875            self.m = 1.0
2876        if self.n is None:
2877            self.n = 1.0
2878        self.c = self.m * exp(-self.n * self.quench)
2879
2880    def update_guess(self, x0):
2881        x0 = asarray(x0)
2882        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2883        T = self.T
2884        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
2885        xnew = xc*(self.upper - self.lower)+self.lower
2886        return xnew
2887#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2888#        xc = y*(self.upper - self.lower)
2889#        xnew = x0 + xc
2890#        return xnew
2891
2892    def update_temp(self):
2893        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2894        self.k += 1
2895        return
2896
2897class cauchy_sa(base_schedule):     #modified to use bounds - not good
2898    def update_guess(self, x0):
2899        x0 = asarray(x0)
2900        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2901        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2902        xnew = xc*(self.upper - self.lower)+self.lower
2903        return xnew
2904#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2905#        xc = self.learn_rate * self.T * tan(numbers)
2906#        xnew = x0 + xc
2907#        return xnew
2908
2909    def update_temp(self):
2910        self.T = self.T0/(1+self.k)
2911        self.k += 1
2912        return
2913
2914class boltzmann_sa(base_schedule):
2915#    def update_guess(self, x0):
2916#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2917#        x0 = asarray(x0)
2918#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2919#
2920#        xnew = x0 + xc*std*self.learn_rate
2921#        return xnew
2922
2923    def update_temp(self):
2924        self.k += 1
2925        self.T = self.T0 / log(self.k+1.0)
2926        return
2927
2928class log_sa(base_schedule):        #OK
2929
2930    def init(self,**options):
2931        self.__dict__.update(options)
2932       
2933    def update_guess(self,x0):     #same as default
2934        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2935       
2936    def update_temp(self):
2937        self.k += 1
2938        self.T = self.T0*self.slope**self.k
2939       
2940class _state(object):
2941    def __init__(self):
2942        self.x = None
2943        self.cost = None
2944
2945# TODO:
2946#     allow for general annealing temperature profile
2947#     in that case use update given by alpha and omega and
2948#     variation of all previous updates and temperature?
2949
2950# Simulated annealing
2951
2952def anneal(func, x0, args=(), schedule='fast', full_output=0,
2953           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2954           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2955           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
2956           ranRange=0.10,autoRan=False,dlg=None):
2957    """Minimize a function using simulated annealing.
2958
2959    Schedule is a schedule class implementing the annealing schedule.
2960    Available ones are 'fast', 'cauchy', 'boltzmann'
2961
2962    :param callable func: f(x, \*args)
2963        Function to be optimized.
2964    :param ndarray x0:
2965        Initial guess.
2966    :param tuple args:
2967        Extra parameters to `func`.
2968    :param base_schedule schedule:
2969        Annealing schedule to use (a class).
2970    :param bool full_output:
2971        Whether to return optional outputs.
2972    :param float T0:
2973        Initial Temperature (estimated as 1.2 times the largest
2974        cost-function deviation over random points in the range).
2975    :param float Tf:
2976        Final goal temperature.
2977    :param int maxeval:
2978        Maximum function evaluations.
2979    :param int maxaccept:
2980        Maximum changes to accept.
2981    :param int maxiter:
2982        Maximum cooling iterations.
2983    :param float learn_rate:
2984        Scale constant for adjusting guesses.
2985    :param float boltzmann:
2986        Boltzmann constant in acceptance test
2987        (increase for less stringent test at each temperature).
2988    :param float feps:
2989        Stopping relative error tolerance for the function value in
2990        last four coolings.
2991    :param float quench,m,n:
2992        Parameters to alter fast_sa schedule.
2993    :param float/ndarray lower,upper:
2994        Lower and upper bounds on `x`.
2995    :param int dwell:
2996        The number of times to search the space at each temperature.
2997    :param float slope:
2998        Parameter for log schedule
2999    :param bool ranStart=False:
3000        True for set 10% of ranges about x
3001
3002    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
3003
3004     * xmin (ndarray): Point giving smallest value found.
3005     * Jmin (float): Minimum value of function found.
3006     * T (float): Final temperature.
3007     * feval (int): Number of function evaluations.
3008     * iters (int): Number of cooling iterations.
3009     * accept (int): Number of tests accepted.
3010     * retval (int): Flag indicating stopping condition:
3011
3012              *  0: Points no longer changing
3013              *  1: Cooled to final temperature
3014              *  2: Maximum function evaluations
3015              *  3: Maximum cooling iterations reached
3016              *  4: Maximum accepted query locations reached
3017              *  5: Final point not the minimum amongst encountered points
3018
3019    *Notes*:
3020    Simulated annealing is a random algorithm which uses no derivative
3021    information from the function being optimized. In practice it has
3022    been more useful in discrete optimization than continuous
3023    optimization, as there are usually better algorithms for continuous
3024    optimization problems.
3025
3026    Some experimentation by trying the difference temperature
3027    schedules and altering their parameters is likely required to
3028    obtain good performance.
3029
3030    The randomness in the algorithm comes from random sampling in numpy.
3031    To obtain the same results you can call numpy.random.seed with the
3032    same seed immediately before calling scipy.optimize.anneal.
3033
3034    We give a brief description of how the three temperature schedules
3035    generate new points and vary their temperature. Temperatures are
3036    only updated with iterations in the outer loop. The inner loop is
3037    over xrange(dwell), and new points are generated for
3038    every iteration in the inner loop. (Though whether the proposed
3039    new points are accepted is probabilistic.)
3040
3041    For readability, let d denote the dimension of the inputs to func.
3042    Also, let x_old denote the previous state, and k denote the
3043    iteration number of the outer loop. All other variables not
3044    defined below are input variables to scipy.optimize.anneal itself.
3045
3046    In the 'fast' schedule the updates are ::
3047
3048        u ~ Uniform(0, 1, size=d)
3049        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
3050        xc = y * (upper - lower)
3051        x_new = x_old + xc
3052
3053        c = n * exp(-n * quench)
3054        T_new = T0 * exp(-c * k**quench)
3055
3056
3057    In the 'cauchy' schedule the updates are ::
3058
3059        u ~ Uniform(-pi/2, pi/2, size=d)
3060        xc = learn_rate * T * tan(u)
3061        x_new = x_old + xc
3062
3063        T_new = T0 / (1+k)
3064
3065    In the 'boltzmann' schedule the updates are ::
3066
3067        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
3068        y ~ Normal(0, std, size=d)
3069        x_new = x_old + learn_rate * y
3070
3071        T_new = T0 / log(1+k)
3072
3073    """
3074    x0 = asarray(x0)
3075    lower = asarray(lower)
3076    upper = asarray(upper)
3077
3078    schedule = eval(schedule+'_sa()')
3079    #   initialize the schedule
3080    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
3081                  learn_rate=learn_rate, lower=lower, upper=upper,
3082                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
3083
3084    current_state, last_state, best_state = _state(), _state(), _state()
3085    if ranStart:
3086        schedule.set_range(x0,ranRange)
3087    if T0 is None:
3088        x0 = schedule.getstart_temp(best_state)
3089    else:
3090        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
3091        best_state.x = None
3092        best_state.cost = numpy.Inf
3093
3094    last_state.x = asarray(x0).copy()
3095    fval = func(x0,*args)
3096    schedule.feval += 1
3097    last_state.cost = fval
3098    if last_state.cost < best_state.cost:
3099        best_state.cost = fval
3100        best_state.x = asarray(x0).copy()
3101    schedule.T = schedule.T0
3102    fqueue = [100, 300, 500, 700]
3103    iters = 1
3104    keepGoing = True
3105    bestn = 0
3106    while keepGoing:
3107        retval = 0
3108        for n in xrange(dwell):
3109            current_state.x = schedule.update_guess(last_state.x)
3110            current_state.cost = func(current_state.x,*args)
3111            schedule.feval += 1
3112
3113            dE = current_state.cost - last_state.cost
3114            if schedule.accept_test(dE):
3115                last_state.x = current_state.x.copy()
3116                last_state.cost = current_state.cost
3117                if last_state.cost < best_state.cost:
3118                    best_state.x = last_state.x.copy()
3119                    best_state.cost = last_state.cost
3120                    bestn = n
3121                    if best_state.cost < 1.0 and autoRan:
3122                        schedule.set_range(x0,best_state.cost/2.)                       
3123        if dlg:
3124            GoOn = dlg.Update(min(100.,best_state.cost*100),
3125                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
3126                    'Best trial:',bestn,  \
3127                    'MC/SA Residual:',best_state.cost*100,'%', \
3128                    ))[0]
3129            if not GoOn:
3130                best_state.x = last_state.x.copy()
3131                best_state.cost = last_state.cost
3132                retval = 5
3133        schedule.update_temp()
3134        iters += 1
3135        # Stopping conditions
3136        # 0) last saved values of f from each cooling step
3137        #     are all very similar (effectively cooled)
3138        # 1) Tf is set and we are below it
3139        # 2) maxeval is set and we are past it
3140        # 3) maxiter is set and we are past it
3141        # 4) maxaccept is set and we are past it
3142        # 5) user canceled run via progress bar
3143
3144        fqueue.append(squeeze(last_state.cost))
3145        fqueue.pop(0)
3146        af = asarray(fqueue)*1.0
3147        if retval == 5:
3148            print ' User terminated run; incomplete MC/SA'
3149            keepGoing = False
3150            break
3151        if all(abs((af-af[0])/af[0]) < feps):
3152            retval = 0
3153            if abs(af[-1]-best_state.cost) > feps*10:
3154                retval = 5
3155#                print "Warning: Cooled to %f at %s but this is not" \
3156#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
3157#                      + " the smallest point found."
3158            break
3159        if (Tf is not None) and (schedule.T < Tf):
3160            retval = 1
3161            break
3162        if (maxeval is not None) and (schedule.feval > maxeval):
3163            retval = 2
3164            break
3165        if (iters > maxiter):
3166            print "Warning: Maximum number of iterations exceeded."
3167            retval = 3
3168            break
3169        if (maxaccept is not None) and (schedule.accepted > maxaccept):
3170            retval = 4
3171            break
3172
3173    if full_output:
3174        return best_state.x, best_state.cost, schedule.T, \
3175               schedule.feval, iters, schedule.accepted, retval
3176    else:
3177        return best_state.x, retval
3178
3179def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
3180    outlist = []
3181    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
3182    for n in range(iCyc):
3183        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
3184        outlist.append(result[0])
3185        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
3186    out_q.put(outlist)
3187
3188def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
3189    import multiprocessing as mp
3190   
3191    nprocs = mp.cpu_count()
3192    out_q = mp.Queue()
3193    procs = []
3194    iCyc = np.zeros(nprocs)
3195    for i in range(nCyc):
3196        iCyc[i%nprocs] += 1
3197    for i in range(nprocs):
3198        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
3199        procs.append(p)
3200        p.start()
3201    resultlist = []
3202    for i in range(nprocs):
3203        resultlist += out_q.get()
3204    for p in procs:
3205        p.join()
3206    return resultlist
3207
3208def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
3209    '''default doc string
3210   
3211    :param type name: description
3212   
3213    :returns: type name: description
3214   
3215    '''
3216   
3217    global tsum
3218    tsum = 0.
3219   
3220    def getMDparms(item,pfx,parmDict,varyList):
3221        parmDict[pfx+'MDaxis'] = item['axis']
3222        parmDict[pfx+'MDval'] = item['Coef'][0]
3223        if item['Coef'][1]:
3224            varyList += [pfx+'MDval',]
3225            limits = item['Coef'][2]
3226            lower.append(limits[0])
3227            upper.append(limits[1])
3228                       
3229    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
3230        parmDict[pfx+'Atype'] = item['atType']
3231        aTypes |= set([item['atType'],]) 
3232        pstr = ['Ax','Ay','Az']
3233        XYZ = [0,0,0]
3234        for i in range(3):
3235            name = pfx+pstr[i]
3236            parmDict[name] = item['Pos'][0][i]
3237            XYZ[i] = parmDict[name]
3238            if item['Pos'][1][i]:
3239                varyList += [name,]
3240                limits = item['Pos'][2][i]
3241                lower.append(limits[0])
3242                upper.append(limits[1])
3243        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
3244           
3245    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
3246        parmDict[mfx+'MolCent'] = item['MolCent']
3247        parmDict[mfx+'RBId'] = item['RBId']
3248        pstr = ['Px','Py','Pz']
3249        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
3250        for i in range(3):
3251            name = mfx+pstr[i]
3252            parmDict[name] = item['Pos'][0][i]
3253            if item['Pos'][1][i]:
3254                varyList += [name,]
3255                limits = item['Pos'][2][i]
3256                lower.append(limits[0])
3257                upper.append(limits[1])
3258        AV = item['Ori'][0]
3259        A = AV[0]
3260        V = AV[1:]
3261        for i in range(4):
3262            name = mfx+ostr[i]
3263            if i:
3264                parmDict[name] = V[i-1]
3265            else:
3266                parmDict[name] = A
3267            if item['Ovar'] == 'AV':
3268                varyList += [name,]
3269                limits = item['Ori'][2][i]
3270                lower.append(limits[0])
3271                upper.append(limits[1])
3272            elif item['Ovar'] == 'A' and not i:
3273                varyList += [name,]
3274                limits = item['Ori'][2][i]
3275                lower.append(limits[0])
3276                upper.append(limits[1])
3277        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
3278            for i in range(len(item['Tor'][0])):
3279                name = mfx+'Tor'+str(i)
3280                parmDict[name] = item['Tor'][0][i]
3281                if item['Tor'][1][i]:
3282                    varyList += [name,]
3283                    limits = item['Tor'][2][i]
3284                    lower.append(limits[0])
3285                    upper.append(limits[1])
3286        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
3287        aTypes |= set(atypes)
3288        atNo += len(atypes)
3289        return atNo
3290               
3291    def GetAtomM(Xdata,SGData):
3292        Mdata = []
3293        for xyz in Xdata:
3294            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
3295        return np.array(Mdata)
3296       
3297    def GetAtomT(RBdata,parmDict):
3298        'Needs a doc string'
3299        atNo = parmDict['atNo']
3300        nfixAt = parmDict['nfixAt']
3301        Tdata = atNo*[' ',]
3302        for iatm in range(nfixAt):
3303            parm = ':'+str(iatm)+':Atype'
3304            if parm in parmDict:
3305                Tdata[iatm] = aTypes.index(parmDict[parm])
3306        iatm = nfixAt
3307        for iObj in range(parmDict['nObj']):
3308            pfx = str(iObj)+':'
3309            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3310                if parmDict[pfx+'Type'] == 'Vector':
3311                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3312                    nAtm = len(RBRes['rbVect'][0])
3313                else:       #Residue
3314                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3315                    nAtm = len(RBRes['rbXYZ'])
3316                for i in range(nAtm):
3317                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
3318                    iatm += 1
3319            elif parmDict[pfx+'Type'] == 'Atom':
3320                atNo = parmDict[pfx+'atNo']
3321                parm = pfx+'Atype'              #remove extra ':'
3322                if parm in parmDict:
3323                    Tdata[atNo] = aTypes.index(parmDict[parm])
3324                iatm += 1
3325            else:
3326                continue        #skips March Dollase
3327        return Tdata
3328       
3329    def GetAtomX(RBdata,parmDict):
3330        'Needs a doc string'
3331        Bmat = parmDict['Bmat']
3332        atNo = parmDict['atNo']
3333        nfixAt = parmDict['nfixAt']
3334        Xdata = np.zeros((3,atNo))
3335        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
3336        for iatm in range(nfixAt):
3337            for key in keys:
3338                parm = ':'+str(iatm)+key
3339                if parm in parmDict:
3340                    keys[key][iatm] = parmDict[parm]
3341        iatm = nfixAt
3342        for iObj in range(parmDict['nObj']):
3343            pfx = str(iObj)+':'
3344            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3345                if parmDict[pfx+'Type'] == 'Vector':
3346                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3347                    vecs = RBRes['rbVect']
3348                    mags = RBRes['VectMag']
3349                    Cart = np.zeros_like(vecs[0])
3350                    for vec,mag in zip(vecs,mags):
3351                        Cart += vec*mag
3352                elif parmDict[pfx+'Type'] == 'Residue':
3353                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3354                    Cart = np.array(RBRes['rbXYZ'])
3355                    for itor,seq in enumerate(RBRes['rbSeq']):
3356                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
3357                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
3358                if parmDict[pfx+'MolCent'][1]:
3359                    Cart -= parmDict[pfx+'MolCent'][0]
3360                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3361                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
3362                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
3363                iatm += len(Cart)
3364            elif parmDict[pfx+'Type'] == 'Atom':
3365                atNo = parmDict[pfx+'atNo']
3366                for key in keys:
3367                    parm = pfx+key[1:]              #remove extra ':'
3368                    if parm in parmDict:
3369                        keys[key][atNo] = parmDict[parm]
3370                iatm += 1
3371            else:
3372                continue        #skips March Dollase
3373        return Xdata.T
3374       
3375    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
3376        allX = np.inner(Xdata,SGM)+SGT
3377        allT = np.repeat(Tdata,allX.shape[1])
3378        allM = np.repeat(Mdata,allX.shape[1])
3379        allX = np.reshape(allX,(-1,3))
3380        return allT,allM,allX
3381
3382    def getAllX(Xdata,SGM,SGT):
3383        allX = np.inner(Xdata,SGM)+SGT
3384        allX = np.reshape(allX,(-1,3))
3385        return allX
3386       
3387    def normQuaternions(RBdata,parmDict,varyList,values):
3388        for iObj in range(parmDict['nObj']):
3389            pfx = str(iObj)+':'
3390            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3391                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3392                A,V = Q2AVdeg(Qori)
3393                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
3394                    if i:
3395                        parmDict[pfx+name] = V[i-1]
3396                    else:
3397                        parmDict[pfx+name] = A
3398       
3399    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
3400        ''' Compute structure factors for all h,k,l for phase
3401        input:
3402            refList: [ref] where each ref = h,k,l,m,d,...
3403            rcov:   array[nref,nref] covariance terms between Fo^2 values
3404            ifInv:  bool True if centrosymmetric
3405            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
3406            RBdata: [dict] rigid body dictionary
3407            varyList: [list] names of varied parameters in MC/SA (not used here)           
3408            ParmDict: [dict] problem parameters
3409        puts result F^2 in each ref[5] in refList
3410        returns:
3411            delt-F*rcov*delt-F/sum(Fo^2)^2
3412        '''       
3413        global tsum
3414        t0 = time.time()
3415        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
3416        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
3417        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
3418        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
3419        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
3420        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
3421        refList[5] = Aterm
3422        if not ifInv:
3423            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
3424        if len(cosTable):        #apply MD correction
3425            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
3426        sumFcsq = np.sum(refList[5])
3427        scale = parmDict['sumFosq']/sumFcsq
3428        refList[5] *= scale
3429        refList[6] = refList[4]-refList[5]
3430        M = np.inner(refList[6],np.inner(rcov,refList[6]))
3431        tsum += (time.time()-t0)
3432        return M/np.sum(refList[4]**2)
3433
3434    sq8ln2 = np.sqrt(8*np.log(2))
3435    sq2pi = np.sqrt(2*np.pi)
3436    sq4pi = np.sqrt(4*np.pi)
3437    generalData = data['General']
3438    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3439    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
3440    SGData = generalData['SGData']
3441    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
3442    SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
3443    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
3444    fixAtoms = data['Atoms']                       #if any
3445    cx,ct,cs = generalData['AtomPtrs'][:3]
3446    aTypes = set([])
3447    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
3448    varyList = []
3449    atNo = 0
3450    for atm in fixAtoms:
3451        pfx = ':'+str(atNo)+':'
3452        parmDict[pfx+'Atype'] = atm[ct]
3453        aTypes |= set([atm[ct],])
3454        pstr = ['Ax','Ay','Az']
3455        parmDict[pfx+'Amul'] = atm[cs+1]
3456        for i in range(3):
3457            name = pfx+pstr[i]
3458            parmDict[name] = atm[cx+i]
3459        atNo += 1
3460    parmDict['nfixAt'] = len(fixAtoms)       
3461    MCSA = generalData['MCSA controls']
3462    reflName = MCSA['Data source']
3463    phaseName = generalData['Name']
3464    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
3465    upper = []
3466    lower = []
3467    MDvec = np.zeros(3)
3468    for i,item in enumerate(MCSAObjs):
3469        mfx = str(i)+':'
3470        parmDict[mfx+'Type'] = item['Type']
3471        if item['Type'] == 'MD':
3472            getMDparms(item,mfx,parmDict,varyList)
3473            MDvec = np.array(item['axis'])
3474        elif item['Type'] == 'Atom':
3475            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
3476            parmDict[mfx+'atNo'] = atNo
3477            atNo += 1
3478        elif item['Type'] in ['Residue','Vector']:
3479            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
3480    parmDict['atNo'] = atNo                 #total no. of atoms
3481    parmDict['nObj'] = len(MCSAObjs)
3482    aTypes = list(aTypes)
3483    Tdata = GetAtomT(RBdata,parmDict)
3484    Xdata = GetAtomX(RBdata,parmDict)
3485    Mdata = GetAtomM(Xdata,SGData)
3486    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
3487    FFtables = G2el.GetFFtable(aTypes)
3488    refs = []
3489    allFF = []
3490    cosTable = []
3491    sumFosq = 0
3492    if 'PWDR' in reflName:
3493        for ref in reflData:
3494            h,k,l,m,d,pos,sig,gam,f = ref[:9]
3495            if d >= MCSA['dmin']:
3496                sig = G2pwd.getgamFW(sig,gam)/sq8ln2        #--> sig from FWHM
3497                SQ = 0.25/d**2
3498                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3499                refs.append([h,k,l,m,f*m,pos,sig])
3500                sumFosq += f*m
3501                Heqv = np.inner(np.array([h,k,l]),SGMT)
3502                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3503        nRef = len(refs)
3504        cosTable = np.array(cosTable)**2
3505        rcov = np.zeros((nRef,nRef))
3506        for iref,refI in enumerate(refs):
3507            rcov[iref][iref] = 1./(sq4pi*refI[6])
3508            for jref,refJ in enumerate(refs[:iref]):
3509                t1 = refI[6]**2+refJ[6]**2
3510                t2 = (refJ[5]-refI[5])**2/(2.*t1)
3511                if t2 > 10.:
3512                    rcov[iref][jref] = 0.
3513                else:
3514                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
3515        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3516        Rdiag = np.sqrt(np.diag(rcov))
3517        Rnorm = np.outer(Rdiag,Rdiag)
3518        rcov /= Rnorm
3519    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
3520        vNames = []
3521        pfx = str(data['pId'])+'::PWLref:'
3522        for iref,refI in enumerate(reflData):           #Pawley reflection set
3523            h,k,l,m,d,v,f,s = refI
3524            if d >= MCSA['dmin'] and v:       #skip unrefined ones
3525                vNames.append(pfx+str(iref))
3526                SQ = 0.25/d**2
3527                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3528                refs.append([h,k,l,m,f*m,iref,0.])
3529                sumFosq += f*m
3530                Heqv = np.inner(np.array([h,k,l]),SGMT)
3531                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3532        cosTable = np.array(cosTable)**2
3533        nRef = len(refs)
3534#        if generalData['doPawley'] and (covData['freshCOV'] or  MCSA['newDmin']):
3535        if covData['freshCOV'] or  MCSA['newDmin']:
3536            vList = covData['varyList']
3537            covMatrix = covData['covMatrix']
3538            rcov = getVCov(vNames,vList,covMatrix)
3539            rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3540            Rdiag = np.sqrt(np.diag(rcov))
3541            Rdiag = np.where(Rdiag,Rdiag,1.0)
3542            Rnorm = np.outer(Rdiag,Rdiag)
3543            rcov /= Rnorm
3544            MCSA['rcov'] = rcov
3545            covData['freshCOV'] = False
3546            MCSA['newDmin'] = False
3547        else:
3548            rcov = MCSA['rcov']
3549    elif 'HKLF' in reflName:
3550        for ref in reflData:
3551            [h,k,l,m,d],f = ref[:5],ref[6]
3552            if d >= MCSA['dmin']:
3553                SQ = 0.25/d**2
3554                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3555                refs.append([h,k,l,m,f*m,0.,0.])
3556                sumFosq += f*m
3557        nRef = len(refs)
3558        rcov = np.identity(len(refs))
3559    allFF = np.array(allFF).T
3560    refs = np.array(refs).T
3561    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
3562    print ' Number of parameters varied: %d'%(len(varyList))
3563    parmDict['sumFosq'] = sumFosq
3564    x0 = [parmDict[val] for val in varyList]
3565    ifInv = SGData['SGInv']
3566    # consider replacing anneal with scipy.optimize.basinhopping
3567    results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
3568        schedule=MCSA['Algorithm'], full_output=True,
3569        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
3570        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
3571        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
3572        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
3573        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
3574    M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
3575#    for ref in refs.T:
3576#        print ' %4d %4d %4d %10.3f %10.3f %10.3f'%(int(ref[0]),int(ref[1]),int(ref[2]),ref[4],ref[5],ref[6])
3577#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
3578    Result = [False,False,results[1],results[2],]+list(results[0])
3579    Result.append(varyList)
3580    return Result,tsum
3581
3582       
3583################################################################################
3584##### Quaternion stuff
3585################################################################################
3586
3587def prodQQ(QA,QB):
3588    ''' Grassman quaternion product
3589        QA,QB quaternions; q=r+ai+bj+ck
3590    '''
3591    D = np.zeros(4)
3592    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
3593    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
3594    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
3595    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
3596   
3597#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
3598#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
3599   
3600    return D
3601   
3602def normQ(QA):
3603    ''' get length of quaternion & normalize it
3604        q=r+ai+bj+ck
3605    '''
3606    n = np.sqrt(np.sum(np.array(QA)**2))
3607    return QA/n
3608   
3609def invQ(Q):
3610    '''
3611        get inverse of quaternion
3612        q=r+ai+bj+ck; q* = r-ai-bj-ck
3613    '''
3614    return Q*np.array([1,-1,-1,-1])
3615   
3616def prodQVQ(Q,V):
3617    """
3618    compute the quaternion vector rotation qvq-1 = v'
3619    q=r+ai+bj+ck
3620    """
3621    T2 = Q[0]*Q[1]
3622    T3 = Q[0]*Q[2]
3623    T4 = Q[0]*Q[3]
3624    T5 = -Q[1]*Q[1]
3625    T6 = Q[1]*Q[2]
3626    T7 = Q[1]*Q[3]
3627    T8 = -Q[2]*Q[2]
3628    T9 = Q[2]*Q[3]
3629    T10 = -Q[3]*Q[3]
3630    M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
3631    VP = 2.*np.inner(V,M)
3632    return VP+V
3633   
3634def Q2Mat(Q):
3635    ''' make rotation matrix from quaternion
3636        q=r+ai+bj+ck
3637    '''
3638    QN = normQ(Q)
3639    aa = QN[0]**2
3640    ab = QN[0]*QN[1]
3641    ac = QN[0]*QN[2]
3642    ad = QN[0]*QN[3]
3643    bb = QN[1]**2
3644    bc = QN[1]*QN[2]
3645    bd = QN[1]*QN[3]
3646    cc = QN[2]**2
3647    cd = QN[2]*QN[3]
3648    dd = QN[3]**2
3649    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
3650        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
3651        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
3652    return np.array(M)
3653   
3654def AV2Q(A,V):
3655    ''' convert angle (radians) & vector to quaternion
3656        q=r+ai+bj+ck
3657    '''
3658    Q = np.zeros(4)
3659    d = nl.norm(np.array(V))
3660    if d:
3661        V /= d
3662        if not A:       #==0.
3663            A = 2.*np.pi
3664        p = A/2.
3665        Q[0] = np.cos(p)
3666        Q[1:4] = V*np.sin(p)
3667    else:
3668        Q[3] = 1.
3669    return Q
3670   
3671def AVdeg2Q(A,V):
3672    ''' convert angle (degrees) & vector to quaternion
3673        q=r+ai+bj+ck
3674    '''
3675    Q = np.zeros(4)
3676    d = nl.norm(np.array(V))
3677    if not A:       #== 0.!
3678        A = 360.
3679    if d:
3680        V /= d
3681        p = A/2.
3682        Q[0] = cosd(p)
3683        Q[1:4] = V*sind(p)
3684    else:
3685        Q[3] = 1.
3686    return Q
3687   
3688def Q2AVdeg(Q):
3689    ''' convert quaternion to angle (degrees 0-360) & normalized vector
3690        q=r+ai+bj+ck
3691    '''
3692    A = 2.*acosd(Q[0])
3693    V = np.array(Q[1:])
3694    V /= sind(A/2.)
3695    return A,V
3696   
3697def Q2AV(Q):
3698    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
3699        q=r+ai+bj+ck
3700    '''
3701    A = 2.*np.arccos(Q[0])
3702    V = np.array(Q[1:])
3703    V /= np.sin(A/2.)
3704    return A,V
3705   
3706def randomQ(r0,r1,r2,r3):
3707    ''' create random quaternion from 4 random numbers in range (-1,1)
3708    '''
3709    sum = 0
3710    Q = np.array(4)
3711    Q[0] = r0
3712    sum += Q[0]**2
3713    Q[1] = np.sqrt(1.-sum)*r1
3714    sum += Q[1]**2
3715    Q[2] = np.sqrt(1.-sum)*r2
3716    sum += Q[2]**2
3717    Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
3718    return Q
3719   
3720def randomAVdeg(r0,r1,r2,r3):
3721    ''' create random angle (deg),vector from 4 random number in range (-1,1)
3722    '''
3723    return Q2AVdeg(randomQ(r0,r1,r2,r3))
3724   
3725def makeQuat(A,B,C):
3726    ''' Make quaternion from rotation of A vector to B vector about C axis
3727
3728        :param np.array A,B,C: Cartesian 3-vectors
3729        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
3730    '''
3731
3732    V1 = np.cross(A,C)
3733    V2 = np.cross(B,C)
3734    if nl.norm(V1)*nl.norm(V2):
3735        V1 /= nl.norm(V1)
3736        V2 /= nl.norm(V2)
3737        V3 = np.cross(V1,V2)
3738    else:
3739        V3 = np.zeros(3)
3740    Q = np.array([0.,0.,0.,1.])
3741    D = 0.
3742    if nl.norm(V3):
3743        V3 /= nl.norm(V3)
3744        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
3745        D = np.arccos(D1)/2.0
3746        V1 = C-V3
3747        V2 = C+V3
3748        DM = nl.norm(V1)
3749        DP = nl.norm(V2)
3750        S = np.sin(D)
3751        Q[0] = np.cos(D)
3752        Q[1:] = V3*S
3753        D *= 2.
3754        if DM > DP:
3755            D *= -1.
3756    return Q,D
3757   
3758def annealtests():
3759    from numpy import cos
3760    # minimum expected at ~-0.195
3761    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
3762    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
3763    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
3764    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
3765
3766    # minimum expected at ~[-0.195, -0.1]
3767    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
3768    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='cauchy')
3769    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='fast')
3770    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='boltzmann')
3771
3772
3773if __name__ == '__main__':
3774    annealtests()
Note: See TracBrowser for help on using the repository browser.