source: trunk/GSASIImath.py @ 1839

Last change on this file since 1839 was 1839, checked in by vondreele, 8 years ago

provide graphical display for texture fit, change weighting scheme & suppress negatives

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 130.1 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2015-05-06 21:42:51 +0000 (Wed, 06 May 2015) $
5# $Author: vondreele $
6# $Revision: 1839 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 1839 2015-05-06 21:42:51Z 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: 1839 $")
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,pgbar):
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,pgbar):
1652        parmDict.update(zip(varyList,values))
1653        Mat = np.empty(0)
1654        sumObs = 0
1655        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
1656        for hist in Gangls.keys():
1657            Refs = refData[hist]
1658            Refs[:,5] = np.where(Refs[:,5]>0.,Refs[:,5],0.)
1659            wt = 1./np.sqrt(np.max(Refs[:,5],.25))
1660            sumObs += np.sum(Refs[:,5])
1661            Refs[:,6] = 1.
1662            H = Refs[:,:3]
1663            phi,beta = G2lat.CrsAng(H,cell,SGData)
1664            psi,gam,x,x = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
1665            for item in parmDict:
1666                if 'C' in item:
1667                    L,M,N = eval(item.strip('C'))
1668                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1669                    Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
1670                    Lnorm = G2lat.Lnorm(L)
1671                    Refs[:,6] += parmDict[item]*Lnorm*Kcl*Ksl
1672            mat = wt*(Refs[:,5]-Refs[:,6])
1673            Mat = np.concatenate((Mat,mat))
1674        sumD = np.sum(np.abs(Mat))
1675        R = min(100.,100.*sumD/sumObs)
1676        pgbar.Update(R,newmsg='Residual = %5.2f'%(R))
1677        print ' Residual: %.3f%%'%(R)
1678        return Mat
1679       
1680    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
1681        Mat = np.empty(0)
1682        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
1683        for hist in Gangls.keys():
1684            mat = np.zeros((len(varyList),len(refData[hist])))
1685            Refs = refData[hist]
1686            H = Refs[:,:3]
1687            wt = 1./np.sqrt(np.max(Refs[:,5],.25))
1688            phi,beta = G2lat.CrsAng(H,cell,SGData)
1689            psi,gam,dPdA,dGdA = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
1690            for j,item in enumerate(varyList):
1691                if 'C' in item:
1692                    L,M,N = eval(item.strip('C'))
1693                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1694                    Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
1695                    Lnorm = G2lat.Lnorm(L)
1696                    mat[j] = -wt*Lnorm*Kcl*Ksl
1697                    for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
1698                        try:
1699                            l = varyList.index(itema)
1700                            mat[l] -= parmDict[item]*wt*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
1701                        except ValueError:
1702                            pass
1703            if len(Mat):
1704                Mat = np.concatenate((Mat,mat.T))
1705            else:
1706                Mat = mat.T
1707        print 'deriv'
1708        return Mat
1709
1710    print ' Fit texture for '+General['Name']
1711    SGData = General['SGData']
1712    cell = General['Cell'][1:7]
1713    Texture = General['SH Texture']
1714    if not Texture['Order']:
1715        return 'No spherical harmonics coefficients'
1716    varyList = []
1717    parmDict = copy.copy(Texture['SH Coeff'][1])
1718    for item in ['Sample omega','Sample chi','Sample phi']:
1719        parmDict[item] = Texture[item][1]
1720        if Texture[item][0]:
1721            varyList.append(item)
1722    if Texture['SH Coeff'][0]:
1723        varyList += Texture['SH Coeff'][1].keys()
1724    while True:
1725        begin = time.time()
1726        values =  np.array(Dict2Values(parmDict, varyList))
1727        result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,ftol=1.e-6,
1728            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList,pgbar))
1729        ncyc = int(result[2]['nfev']/2)
1730        if ncyc:
1731            runtime = time.time()-begin   
1732            chisq = np.sum(result[2]['fvec']**2)
1733            Values2Dict(parmDict, varyList, result[0])
1734            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
1735            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(result[2]['fvec']),' Number of parameters: ',len(varyList)
1736            print 'refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
1737            try:
1738                sig = np.sqrt(np.diag(result[1])*GOF)
1739                if np.any(np.isnan(sig)):
1740                    print '*** Least squares aborted - some invalid esds possible ***'
1741                break                   #refinement succeeded - finish up!
1742            except ValueError:          #result[1] is None on singular matrix
1743                print '**** Refinement failed - singular matrix ****'
1744                return None
1745        else:
1746            break
1747   
1748    if ncyc:
1749        for parm in parmDict:
1750            if 'C' in parm:
1751                Texture['SH Coeff'][1][parm] = parmDict[parm]
1752            else:
1753                Texture[parm][1] = parmDict[parm] 
1754        sigDict = dict(zip(varyList,sig))
1755        printSpHarm(Texture,sigDict)
1756       
1757    return None
1758   
1759################################################################################
1760##### Fourier & charge flip stuff
1761################################################################################
1762
1763def adjHKLmax(SGData,Hmax):
1764    '''default doc string
1765   
1766    :param type name: description
1767   
1768    :returns: type name: description
1769   
1770    '''
1771    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1772        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
1773        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
1774        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
1775    else:
1776        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
1777        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
1778        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
1779
1780def OmitMap(data,reflDict,pgbar=None):
1781    '''default doc string
1782   
1783    :param type name: description
1784   
1785    :returns: type name: description
1786   
1787    '''
1788    generalData = data['General']
1789    if not generalData['Map']['MapType']:
1790        print '**** ERROR - Fourier map not defined'
1791        return
1792    mapData = generalData['Map']
1793    dmin = mapData['Resolution']
1794    SGData = generalData['SGData']
1795    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1796    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1797    cell = generalData['Cell'][1:8]       
1798    A = G2lat.cell2A(cell[:6])
1799    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1800    adjHKLmax(SGData,Hmax)
1801    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1802    time0 = time.time()
1803    for iref,ref in enumerate(reflDict['RefList']):
1804        if ref[4] >= dmin:
1805            Fosq,Fcsq,ph = ref[8:11]
1806            Uniq = np.inner(ref[:3],SGMT)
1807            Phi = np.inner(ref[:3],SGT)
1808            for i,hkl in enumerate(Uniq):        #uses uniq
1809                hkl = np.asarray(hkl,dtype='i')
1810                dp = 360.*Phi[i]                #and phi
1811                a = cosd(ph+dp)
1812                b = sind(ph+dp)
1813                phasep = complex(a,b)
1814                phasem = complex(a,-b)
1815                Fo = np.sqrt(Fosq)
1816                if '2Fo-Fc' in mapData['MapType']:
1817                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
1818                else:
1819                    F = np.sqrt(Fosq)
1820                h,k,l = hkl+Hmax
1821                Fhkl[h,k,l] = F*phasep
1822                h,k,l = -hkl+Hmax
1823                Fhkl[h,k,l] = F*phasem
1824    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1825    M = np.mgrid[0:4,0:4,0:4]
1826    blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
1827    iBeg = blkIds*rho0.shape/4
1828    iFin = (blkIds+1)*rho0.shape/4
1829    rho_omit = np.zeros_like(rho0)
1830    nBlk = 0
1831    for iB,iF in zip(iBeg,iFin):
1832        rho1 = np.copy(rho0)
1833        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
1834        Fnew = fft.ifftshift(fft.ifftn(rho1))
1835        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
1836        phase = Fnew/np.absolute(Fnew)
1837        OFhkl = np.absolute(Fhkl)*phase
1838        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
1839        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]])
1840        nBlk += 1
1841        pgbar.Update(nBlk)
1842    mapData['rho'] = np.real(rho_omit)/cell[6]
1843    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1844    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
1845    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1846    return mapData
1847   
1848def Fourier4DMap(data,reflDict):
1849    '''default doc string
1850   
1851    :param type name: description
1852   
1853    :returns: type name: description
1854   
1855    '''
1856    generalData = data['General']
1857    mapData = generalData['4DmapData']
1858    dmin = mapData['Resolution']
1859    SGData = generalData['SGData']
1860    SSGData = generalData['SSGData']
1861    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1862    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1863    cell = generalData['Cell'][1:8]       
1864    A = G2lat.cell2A(cell[:6])
1865    maxM = generalData['SuperVec'][2]
1866    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
1867    adjHKLmax(SGData,Hmax)
1868    Hmax = np.asarray(Hmax,dtype='i')+1
1869    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1870    time0 = time.time()
1871    for iref,ref in enumerate(reflDict['RefList']):
1872        if ref[5] > dmin:
1873            Fosq,Fcsq,ph = ref[9:12]
1874            Uniq = np.inner(ref[:4],SSGMT)
1875            Phi = np.inner(ref[:4],SSGT)
1876            for i,hkl in enumerate(Uniq):        #uses uniq
1877                hkl = np.asarray(hkl,dtype='i')
1878                dp = 360.*Phi[i]                #and phi
1879                a = cosd(ph+dp)
1880                b = sind(ph+dp)
1881                phasep = complex(a,b)
1882                phasem = complex(a,-b)
1883                if 'Fobs' in mapData['MapType']:
1884                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1885                    h,k,l,m = hkl+Hmax
1886                    Fhkl[h,k,l,m] = F*phasep
1887                    h,k,l,m = -hkl+Hmax
1888                    Fhkl[h,k,l,m] = F*phasem
1889                elif 'delt-F' in mapData['MapType']:
1890                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1891                    h,k,l,m = hkl+Hmax
1892                    Fhkl[h,k,l,m] = dF*phasep
1893                    h,k,l,m = -hkl+Hmax
1894                    Fhkl[h,k,l,m] = dF*phasem
1895    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1896    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1897    mapData['Type'] = reflDict['Type']
1898    mapData['rho'] = np.real(rho)
1899    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1900    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
1901    return mapData
1902
1903def FourierMap(data,reflDict):
1904    '''default doc string
1905   
1906    :param type name: description
1907   
1908    :returns: type name: description
1909   
1910    '''
1911    generalData = data['General']
1912    mapData = generalData['Map']
1913    dmin = mapData['Resolution']
1914    SGData = generalData['SGData']
1915    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1916    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1917    cell = generalData['Cell'][1:8]       
1918    A = G2lat.cell2A(cell[:6])
1919    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1920    adjHKLmax(SGData,Hmax)
1921    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1922#    Fhkl[0,0,0] = generalData['F000X']
1923    time0 = time.time()
1924    for iref,ref in enumerate(reflDict['RefList']):
1925        if ref[4] > dmin:
1926            Fosq,Fcsq,ph = ref[8:11]
1927            Uniq = np.inner(ref[:3],SGMT)
1928            Phi = np.inner(ref[:3],SGT)
1929            for i,hkl in enumerate(Uniq):        #uses uniq
1930                hkl = np.asarray(hkl,dtype='i')
1931                dp = 360.*Phi[i]                #and phi
1932                a = cosd(ph)
1933                b = sind(ph)
1934                phasep = complex(a,b)+dp
1935                phasem = complex(a,-b)-dp
1936                if 'Fobs' in mapData['MapType']:
1937                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1938                    h,k,l = hkl+Hmax
1939                    Fhkl[h,k,l] = F*phasep
1940                    h,k,l = -hkl+Hmax
1941                    Fhkl[h,k,l] = F*phasem
1942                elif 'Fcalc' in mapData['MapType']:
1943                    F = np.sqrt(Fcsq)
1944                    h,k,l = hkl+Hmax
1945                    Fhkl[h,k,l] = F*phasep
1946                    h,k,l = -hkl+Hmax
1947                    Fhkl[h,k,l] = F*phasem
1948                elif 'delt-F' in mapData['MapType']:
1949                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1950                    h,k,l = hkl+Hmax
1951                    Fhkl[h,k,l] = dF*phasep
1952                    h,k,l = -hkl+Hmax
1953                    Fhkl[h,k,l] = dF*phasem
1954                elif '2*Fo-Fc' in mapData['MapType']:
1955                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1956                    h,k,l = hkl+Hmax
1957                    Fhkl[h,k,l] = F*phasep
1958                    h,k,l = -hkl+Hmax
1959                    Fhkl[h,k,l] = F*phasem
1960                elif 'Patterson' in mapData['MapType']:
1961                    h,k,l = hkl+Hmax
1962                    Fhkl[h,k,l] = complex(Fosq,0.)
1963                    h,k,l = -hkl+Hmax
1964                    Fhkl[h,k,l] = complex(Fosq,0.)
1965    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1966    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1967    mapData['Type'] = reflDict['Type']
1968    mapData['rho'] = np.real(rho)
1969    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1970    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
1971    return mapData
1972   
1973# map printing for testing purposes
1974def printRho(SGLaue,rho,rhoMax):                         
1975    '''default doc string
1976   
1977    :param type name: description
1978   
1979    :returns: type name: description
1980   
1981    '''
1982    dim = len(rho.shape)
1983    if dim == 2:
1984        ix,jy = rho.shape
1985        for j in range(jy):
1986            line = ''
1987            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1988                line += (jy-j)*'  '
1989            for i in range(ix):
1990                r = int(100*rho[i,j]/rhoMax)
1991                line += '%4d'%(r)
1992            print line+'\n'
1993    else:
1994        ix,jy,kz = rho.shape
1995        for k in range(kz):
1996            print 'k = ',k
1997            for j in range(jy):
1998                line = ''
1999                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
2000                    line += (jy-j)*'  '
2001                for i in range(ix):
2002                    r = int(100*rho[i,j,k]/rhoMax)
2003                    line += '%4d'%(r)
2004                print line+'\n'
2005## keep this
2006               
2007def findOffset(SGData,A,Fhkl):   
2008    '''default doc string
2009   
2010    :param type name: description
2011   
2012    :returns: type name: description
2013   
2014    '''
2015    if SGData['SpGrp'] == 'P 1':
2016        return [0,0,0]   
2017    hklShape = Fhkl.shape
2018    hklHalf = np.array(hklShape)/2
2019    sortHKL = np.argsort(Fhkl.flatten())
2020    Fdict = {}
2021    for hkl in sortHKL:
2022        HKL = np.unravel_index(hkl,hklShape)
2023        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
2024        if F == 0.:
2025            break
2026        Fdict['%.6f'%(np.absolute(F))] = hkl
2027    Flist = np.flipud(np.sort(Fdict.keys()))
2028    F = str(1.e6)
2029    i = 0
2030    DH = []
2031    Dphi = []
2032    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2033    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2034    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
2035    for F in Flist:
2036        hkl = np.unravel_index(Fdict[F],hklShape)
2037        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
2038            continue
2039        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
2040        Uniq = np.array(Uniq,dtype='i')
2041        Phi = np.array(Phi)
2042        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
2043        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
2044        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
2045        ang0 = np.angle(Fh0,deg=True)/360.
2046        for H,phi in zip(Uniq,Phi)[1:]:
2047            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
2048            dH = H-hkl
2049            dang = ang-ang0
2050            DH.append(dH)
2051            Dphi.append((dang+.5) % 1.0)
2052        if i > 20 or len(DH) > 30:
2053            break
2054        i += 1
2055    DH = np.array(DH)
2056    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
2057    Dphi = np.array(Dphi)
2058    steps = np.array(hklShape)
2059    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
2060    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
2061    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
2062    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
2063    hist,bins = np.histogram(Mmap,bins=1000)
2064#    for i,item in enumerate(hist[:10]):
2065#        print item,bins[i]
2066    chisq = np.min(Mmap)
2067    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
2068    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
2069#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
2070    return DX
2071   
2072def ChargeFlip(data,reflDict,pgbar):
2073    '''default doc string
2074   
2075    :param type name: description
2076   
2077    :returns: type name: description
2078   
2079    '''
2080    generalData = data['General']
2081    mapData = generalData['Map']
2082    flipData = generalData['Flip']
2083    FFtable = {}
2084    if 'None' not in flipData['Norm element']:
2085        normElem = flipData['Norm element'].upper()
2086        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
2087        for ff in FFs:
2088            if ff['Symbol'] == normElem:
2089                FFtable.update(ff)
2090    dmin = flipData['Resolution']
2091    SGData = generalData['SGData']
2092    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2093    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2094    cell = generalData['Cell'][1:8]       
2095    A = G2lat.cell2A(cell[:6])
2096    Vol = cell[6]
2097    im = 0
2098    if generalData['Type'] in ['modulated','magnetic',]:
2099        im = 1
2100    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2101    adjHKLmax(SGData,Hmax)
2102    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
2103    time0 = time.time()
2104    for iref,ref in enumerate(reflDict['RefList']):
2105        dsp = ref[4+im]
2106        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
2107            continue
2108        if dsp > dmin:
2109            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
2110            if FFtable:
2111                SQ = 0.25/dsp**2
2112                ff *= G2el.ScatFac(FFtable,SQ)[0]
2113            if ref[8+im] > 0.:         #use only +ve Fobs**2
2114                E = np.sqrt(ref[8+im])/ff
2115            else:
2116                E = 0.
2117            ph = ref[10]
2118            ph = rn.uniform(0.,360.)
2119            Uniq = np.inner(ref[:3],SGMT)
2120            Phi = np.inner(ref[:3],SGT)
2121            for i,hkl in enumerate(Uniq):        #uses uniq
2122                hkl = np.asarray(hkl,dtype='i')
2123                dp = 360.*Phi[i]                #and phi
2124                a = cosd(ph+dp)
2125                b = sind(ph+dp)
2126                phasep = complex(a,b)
2127                phasem = complex(a,-b)
2128                h,k,l = hkl+Hmax
2129                Ehkl[h,k,l] = E*phasep
2130                h,k,l = -hkl+Hmax       #Friedel pair refl.
2131                Ehkl[h,k,l] = E*phasem
2132#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
2133    CEhkl = copy.copy(Ehkl)
2134    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
2135    Emask = ma.getmask(MEhkl)
2136    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
2137    Ncyc = 0
2138    old = np.seterr(all='raise')
2139    while True:       
2140        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
2141        CEsig = np.std(CErho)
2142        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
2143        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
2144        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
2145        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
2146        phase = CFhkl/np.absolute(CFhkl)
2147        CEhkl = np.absolute(Ehkl)*phase
2148        Ncyc += 1
2149        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
2150        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
2151        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
2152        if Rcf < 5.:
2153            break
2154        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
2155        if not GoOn or Ncyc > 10000:
2156            break
2157    np.seterr(**old)
2158    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
2159    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
2160    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
2161    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
2162       
2163    mapData['Rcf'] = Rcf
2164    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2165    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2166    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2167    mapData['Type'] = reflDict['Type']
2168    return mapData
2169   
2170def findSSOffset(SGData,SSGData,A,Fhklm):   
2171    '''default doc string
2172   
2173    :param type name: description
2174   
2175    :returns: type name: description
2176   
2177    '''
2178    if SGData['SpGrp'] == 'P 1':
2179        return [0,0,0,0]   
2180    hklmShape = Fhklm.shape
2181    hklmHalf = np.array(hklmShape)/2
2182    sortHKLM = np.argsort(Fhklm.flatten())
2183    Fdict = {}
2184    for hklm in sortHKLM:
2185        HKLM = np.unravel_index(hklm,hklmShape)
2186        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
2187        if F == 0.:
2188            break
2189        Fdict['%.6f'%(np.absolute(F))] = hklm
2190    Flist = np.flipud(np.sort(Fdict.keys()))
2191    F = str(1.e6)
2192    i = 0
2193    DH = []
2194    Dphi = []
2195    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2196    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2197    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
2198    for F in Flist:
2199        hklm = np.unravel_index(Fdict[F],hklmShape)
2200        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
2201            continue
2202        Uniq = np.inner(hklm-hklmHalf,SSGMT)
2203        Phi = np.inner(hklm-hklmHalf,SSGT)
2204        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
2205        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
2206        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
2207        ang0 = np.angle(Fh0,deg=True)/360.
2208        for H,phi in zip(Uniq,Phi)[1:]:
2209            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
2210            dH = H-hklm
2211            dang = ang-ang0
2212            DH.append(dH)
2213            Dphi.append((dang+.5) % 1.0)
2214        if i > 20 or len(DH) > 30:
2215            break
2216        i += 1
2217    DH = np.array(DH)
2218    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
2219    Dphi = np.array(Dphi)
2220    steps = np.array(hklmShape)
2221    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]]
2222    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
2223    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
2224    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
2225    hist,bins = np.histogram(Mmap,bins=1000)
2226#    for i,item in enumerate(hist[:10]):
2227#        print item,bins[i]
2228    chisq = np.min(Mmap)
2229    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
2230    print ' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
2231#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
2232    return DX
2233   
2234def SSChargeFlip(data,reflDict,pgbar):
2235    '''default doc string
2236   
2237    :param type name: description
2238   
2239    :returns: type name: description
2240   
2241    '''
2242    generalData = data['General']
2243    mapData = generalData['Map']
2244    map4DData = {}
2245    flipData = generalData['Flip']
2246    FFtable = {}
2247    if 'None' not in flipData['Norm element']:
2248        normElem = flipData['Norm element'].upper()
2249        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
2250        for ff in FFs:
2251            if ff['Symbol'] == normElem:
2252                FFtable.update(ff)
2253    dmin = flipData['Resolution']
2254    SGData = generalData['SGData']
2255    SSGData = generalData['SSGData']
2256    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2257    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2258    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2259    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2260    cell = generalData['Cell'][1:8]       
2261    A = G2lat.cell2A(cell[:6])
2262    Vol = cell[6]
2263    maxM = generalData['SuperVec'][2]
2264    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
2265    adjHKLmax(SGData,Hmax)
2266    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
2267    time0 = time.time()
2268    for iref,ref in enumerate(reflDict['RefList']):
2269        dsp = ref[5]
2270        if dsp > dmin:
2271            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
2272            if FFtable:
2273                SQ = 0.25/dsp**2
2274                ff *= G2el.ScatFac(FFtable,SQ)[0]
2275            if ref[9] > 0.:         #use only +ve Fobs**2
2276                E = np.sqrt(ref[9])/ff
2277            else:
2278                E = 0.
2279            ph = ref[11]
2280            ph = rn.uniform(0.,360.)
2281            Uniq = np.inner(ref[:4],SSGMT)
2282            Phi = np.inner(ref[:4],SSGT)
2283            for i,hklm in enumerate(Uniq):        #uses uniq
2284                hklm = np.asarray(hklm,dtype='i')
2285                dp = 360.*Phi[i]                #and phi
2286                a = cosd(ph+dp)
2287                b = sind(ph+dp)
2288                phasep = complex(a,b)
2289                phasem = complex(a,-b)
2290                h,k,l,m = hklm+Hmax
2291                Ehkl[h,k,l,m] = E*phasep
2292                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
2293                Ehkl[h,k,l,m] = E*phasem
2294#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
2295    CEhkl = copy.copy(Ehkl)
2296    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
2297    Emask = ma.getmask(MEhkl)
2298    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
2299    Ncyc = 0
2300    old = np.seterr(all='raise')
2301    while True:       
2302        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
2303        CEsig = np.std(CErho)
2304        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
2305        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
2306        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
2307        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
2308        phase = CFhkl/np.absolute(CFhkl)
2309        CEhkl = np.absolute(Ehkl)*phase
2310        Ncyc += 1
2311        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
2312        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
2313        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
2314        if Rcf < 5.:
2315            break
2316        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
2317        if not GoOn or Ncyc > 10000:
2318            break
2319    np.seterr(**old)
2320    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
2321    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))
2322    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))
2323    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
2324    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
2325
2326    mapData['Rcf'] = Rcf
2327    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2328    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2329    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2330    mapData['Type'] = reflDict['Type']
2331
2332    map4DData['Rcf'] = Rcf
2333    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))
2334    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
2335    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
2336    map4DData['Type'] = reflDict['Type']
2337    return mapData,map4DData
2338   
2339def SearchMap(generalData,drawingData,Neg=False):
2340    '''Does a search of a density map for peaks meeting the criterion of peak
2341    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
2342    mapData is data['General']['mapData']; the map is also in mapData.
2343
2344    :param generalData: the phase data structure; includes the map
2345    :param drawingData: the drawing data structure
2346    :param Neg:  if True then search for negative peaks (i.e. H-atoms & neutron data)
2347
2348    :returns: (peaks,mags,dzeros) where
2349
2350        * peaks : ndarray
2351          x,y,z positions of the peaks found in the map
2352        * mags : ndarray
2353          the magnitudes of the peaks
2354        * dzeros : ndarray
2355          the distance of the peaks from  the unit cell origin
2356        * dcent : ndarray
2357          the distance of the peaks from  the unit cell center
2358
2359    '''       
2360    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2361   
2362    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
2363   
2364#    def noDuplicate(xyz,peaks,Amat):
2365#        XYZ = np.inner(Amat,xyz)
2366#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2367#            print ' Peak',xyz,' <0.5A from another peak'
2368#            return False
2369#        return True
2370#                           
2371    def fixSpecialPos(xyz,SGData,Amat):
2372        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
2373        X = []
2374        xyzs = [equiv[0] for equiv in equivs]
2375        for x in xyzs:
2376            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
2377                X.append(x)
2378        if len(X) > 1:
2379            return np.average(X,axis=0)
2380        else:
2381            return xyz
2382       
2383    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
2384        Mag,x0,y0,z0,sig = parms
2385        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
2386#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
2387        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
2388       
2389    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
2390        Mag,x0,y0,z0,sig = parms
2391        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2392        return M
2393       
2394    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
2395        Mag,x0,y0,z0,sig = parms
2396        dMdv = np.zeros(([5,]+list(rX.shape)))
2397        delt = .01
2398        for i in range(5):
2399            parms[i] -= delt
2400            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2401            parms[i] += 2.*delt
2402            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2403            parms[i] -= delt
2404            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
2405        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2406        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
2407        dMdv = np.reshape(dMdv,(5,rX.size))
2408        Hess = np.inner(dMdv,dMdv)
2409       
2410        return Vec,Hess
2411       
2412    phaseName = generalData['Name']
2413    SGData = generalData['SGData']
2414    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2415    peaks = []
2416    mags = []
2417    dzeros = []
2418    dcent = []
2419    try:
2420        mapData = generalData['Map']
2421        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
2422        if Neg:
2423            rho = -copy.copy(mapData['rho'])     #flip +/-
2424        else:
2425            rho = copy.copy(mapData['rho'])     #don't mess up original
2426        mapHalf = np.array(rho.shape)/2
2427        res = mapData['Resolution']
2428        incre = np.array(rho.shape,dtype=np.float)
2429        step = max(1.0,1./res)+1
2430        steps = np.array(3*[step,])
2431    except KeyError:
2432        print '**** ERROR - Fourier map not defined'
2433        return peaks,mags
2434    rhoMask = ma.array(rho,mask=(rho<contLevel))
2435    indices = (-1,0,1)
2436    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2437    for roll in rolls:
2438        if np.any(roll):
2439            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
2440    indx = np.transpose(rhoMask.nonzero())
2441    peaks = indx/incre
2442    mags = rhoMask[rhoMask.nonzero()]
2443    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
2444        rho = rollMap(rho,ind)
2445        rMM = mapHalf-steps
2446        rMP = mapHalf+steps+1
2447        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2448        peakInt = np.sum(rhoPeak)*res**3
2449        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2450        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
2451        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
2452            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
2453        x1 = result[0]
2454        if not np.any(x1 < 0):
2455            mag = x1[0]
2456            peak = (np.array(x1[1:4])-ind)/incre
2457        peak = fixSpecialPos(peak,SGData,Amat)
2458        rho = rollMap(rho,-ind)
2459    cent = np.ones(3)*.5     
2460    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
2461    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
2462    if Neg:     #want negative magnitudes for negative peaks
2463        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
2464    else:
2465        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
2466   
2467def sortArray(data,pos,reverse=False):
2468    '''data is a list of items
2469    sort by pos in list; reverse if True
2470    '''
2471    T = []
2472    for i,M in enumerate(data):
2473        try:
2474            T.append((M[pos],i))
2475        except IndexError:
2476            return data
2477    D = dict(zip(T,data))
2478    T.sort()
2479    if reverse:
2480        T.reverse()
2481    X = []
2482    for key in T:
2483        X.append(D[key])
2484    return X
2485
2486def PeaksEquiv(data,Ind):
2487    '''Find the equivalent map peaks for those selected. Works on the
2488    contents of data['Map Peaks'].
2489
2490    :param data: the phase data structure
2491    :param list Ind: list of selected peak indices
2492    :returns: augmented list of peaks including those related by symmetry to the
2493      ones in Ind
2494
2495    '''       
2496    def Duplicate(xyz,peaks,Amat):
2497        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2498            return True
2499        return False
2500                           
2501    generalData = data['General']
2502    cell = generalData['Cell'][1:7]
2503    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2504    A = G2lat.cell2A(cell)
2505    SGData = generalData['SGData']
2506    mapPeaks = data['Map Peaks']
2507    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
2508    Indx = {}
2509    for ind in Ind:
2510        xyz = np.array(mapPeaks[ind][1:4])
2511        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
2512        for jnd,xyz in enumerate(XYZ):       
2513            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
2514    Ind = []
2515    for ind in Indx:
2516        if Indx[ind]:
2517            Ind.append(ind)
2518    return Ind
2519               
2520def PeaksUnique(data,Ind):
2521    '''Finds the symmetry unique set of peaks from those selected. Works on the
2522    contents of data['Map Peaks'].
2523
2524    :param data: the phase data structure
2525    :param list Ind: list of selected peak indices
2526    :returns: the list of symmetry unique peaks from among those given in Ind
2527
2528    '''       
2529#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
2530
2531    def noDuplicate(xyz,peaks,Amat):
2532        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2533            return False
2534        return True
2535                           
2536    generalData = data['General']
2537    cell = generalData['Cell'][1:7]
2538    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2539    A = G2lat.cell2A(cell)
2540    SGData = generalData['SGData']
2541    mapPeaks = data['Map Peaks']
2542    Indx = {}
2543    XYZ = {}
2544    for ind in Ind:
2545        XYZ[ind] = np.array(mapPeaks[ind][1:4])
2546        Indx[ind] = True
2547    for ind in Ind:
2548        if Indx[ind]:
2549            xyz = XYZ[ind]
2550            for jnd in Ind:
2551                if ind != jnd and Indx[jnd]:                       
2552                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
2553                    xyzs = np.array([equiv[0] for equiv in Equiv])
2554                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
2555    Ind = []
2556    for ind in Indx:
2557        if Indx[ind]:
2558            Ind.append(ind)
2559    return Ind
2560   
2561################################################################################
2562##### single peak fitting profile fxn stuff
2563################################################################################
2564
2565def getCWsig(ins,pos):
2566    '''get CW peak profile sigma
2567   
2568    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
2569      as values only
2570    :param float pos: 2-theta of peak
2571    :returns: float getCWsig: peak sigma
2572   
2573    '''
2574    tp = tand(pos/2.0)
2575    return ins['U']*tp**2+ins['V']*tp+ins['W']
2576   
2577def getCWsigDeriv(pos):
2578    '''get derivatives of CW peak profile sigma wrt U,V, & W
2579   
2580    :param float pos: 2-theta of peak
2581   
2582    :returns: list getCWsigDeriv: d(sig)/dU, d(sig)/dV & d(sig)/dW
2583   
2584    '''
2585    tp = tand(pos/2.0)
2586    return tp**2,tp,1.0
2587   
2588def getCWgam(ins,pos):
2589    '''get CW peak profile gamma
2590   
2591    :param dict ins: instrument parameters with at least 'X' & 'Y'
2592      as values only
2593    :param float pos: 2-theta of peak
2594    :returns: float getCWgam: peak gamma
2595   
2596    '''
2597    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
2598   
2599def getCWgamDeriv(pos):
2600    '''get derivatives of CW peak profile gamma wrt X & Y
2601   
2602    :param float pos: 2-theta of peak
2603   
2604    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
2605   
2606    '''
2607    return 1./cosd(pos/2.0),tand(pos/2.0)
2608   
2609def getTOFsig(ins,dsp):
2610    '''get TOF peak profile sigma
2611   
2612    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
2613      as values only
2614    :param float dsp: d-spacing of peak
2615   
2616    :returns: float getTOFsig: peak sigma
2617   
2618    '''
2619    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
2620   
2621def getTOFsigDeriv(dsp):
2622    '''get derivatives of TOF peak profile gamma wrt sig-0, sig-1, & sig-q
2623   
2624    :param float dsp: d-spacing of peak
2625   
2626    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
2627   
2628    '''
2629    return 1.0,dsp**2,dsp**4,1./dsp**2
2630   
2631def getTOFgamma(ins,dsp):
2632    '''get TOF peak profile gamma
2633   
2634    :param dict ins: instrument parameters with at least 'X' & 'Y'
2635      as values only
2636    :param float dsp: d-spacing of peak
2637   
2638    :returns: float getTOFgamma: peak gamma
2639   
2640    '''
2641    return ins['X']*dsp+ins['Y']*dsp**2
2642   
2643def getTOFgammaDeriv(dsp):
2644    '''get derivatives of TOF peak profile gamma wrt X & Y
2645   
2646    :param float dsp: d-spacing of peak
2647   
2648    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
2649   
2650    '''
2651    return dsp,dsp**2
2652   
2653def getTOFbeta(ins,dsp):
2654    '''get TOF peak profile beta
2655   
2656    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
2657      as values only
2658    :param float dsp: d-spacing of peak
2659   
2660    :returns: float getTOFbeta: peak beat
2661   
2662    '''
2663    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
2664   
2665def getTOFbetaDeriv(dsp):
2666    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
2667   
2668    :param float dsp: d-spacing of peak
2669   
2670    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
2671   
2672    '''
2673    return 1.0,1./dsp**4,1./dsp**2
2674   
2675def getTOFalpha(ins,dsp):
2676    '''get TOF peak profile alpha
2677   
2678    :param dict ins: instrument parameters with at least 'alpha'
2679      as values only
2680    :param float dsp: d-spacing of peak
2681   
2682    :returns: flaot getTOFalpha: peak alpha
2683   
2684    '''
2685    return ins['alpha']/dsp
2686   
2687def getTOFalphaDeriv(dsp):
2688    '''get derivatives of TOF peak profile beta wrt alpha
2689   
2690    :param float dsp: d-spacing of peak
2691   
2692    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
2693   
2694    '''
2695    return 1./dsp
2696   
2697def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
2698    '''set starting peak parameters for single peak fits from plot selection or auto selection
2699   
2700    :param dict Parms: instrument parameters dictionary
2701    :param dict Parms2: table lookup for TOF profile coefficients
2702    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
2703    :param float mag: peak top magnitude from pick
2704    :param bool ifQ: True if pos in Q
2705    :param bool useFit: True if use fitted CW Parms values (not defaults)
2706   
2707    :returns: list XY: peak list entry:
2708        for CW: [pos,0,mag,1,sig,0,gam,0]
2709        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2710        NB: mag refinement set by default, all others off
2711   
2712    '''
2713    ind = 0
2714    if useFit:
2715        ind = 1
2716    ins = {}
2717    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
2718        for x in ['U','V','W','X','Y']:
2719            ins[x] = Parms[x][ind]
2720        if ifQ:                              #qplot - convert back to 2-theta
2721            pos = 2.0*asind(pos*wave/(4*math.pi))
2722        sig = getCWsig(ins,pos)
2723        gam = getCWgam(ins,pos)           
2724        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
2725    else:
2726        if ifQ:
2727            dsp = 2.*np.pi/pos
2728            pos = Parms['difC']*dsp
2729        else:
2730            dsp = pos/Parms['difC'][1]
2731        if 'Pdabc' in Parms2:
2732            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
2733                ins[x] = Parms[x][ind]
2734            Pdabc = Parms2['Pdabc'].T
2735            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2736            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2737        else:
2738            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
2739                ins[x] = Parms[x][ind]
2740            alp = getTOFalpha(ins,dsp)
2741            bet = getTOFbeta(ins,dsp)
2742        sig = getTOFsig(ins,dsp)
2743        gam = getTOFgamma(ins,dsp)
2744        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2745    return XY
2746   
2747################################################################################
2748##### MC/SA stuff
2749################################################################################
2750
2751#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2752# Original Author: Travis Oliphant 2002
2753# Bug-fixes in 2006 by Tim Leslie
2754
2755
2756import numpy
2757from numpy import asarray, tan, exp, ones, squeeze, sign, \
2758     all, log, sqrt, pi, shape, array, minimum, where
2759from numpy import random
2760
2761#__all__ = ['anneal']
2762
2763_double_min = numpy.finfo(float).min
2764_double_max = numpy.finfo(float).max
2765class base_schedule(object):
2766    def __init__(self):
2767        self.dwell = 20
2768        self.learn_rate = 0.5
2769        self.lower = -10
2770        self.upper = 10
2771        self.Ninit = 50
2772        self.accepted = 0
2773        self.tests = 0
2774        self.feval = 0
2775        self.k = 0
2776        self.T = None
2777
2778    def init(self, **options):
2779        self.__dict__.update(options)
2780        self.lower = asarray(self.lower)
2781        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2782        self.upper = asarray(self.upper)
2783        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2784        self.k = 0
2785        self.accepted = 0
2786        self.feval = 0
2787        self.tests = 0
2788
2789    def getstart_temp(self, best_state):
2790        """ Find a matching starting temperature and starting parameters vector
2791        i.e. find x0 such that func(x0) = T0.
2792
2793        Parameters
2794        ----------
2795        best_state : _state
2796            A _state object to store the function value and x0 found.
2797
2798        returns
2799        -------
2800        x0 : array
2801            The starting parameters vector.
2802        """
2803
2804        assert(not self.dims is None)
2805        lrange = self.lower
2806        urange = self.upper
2807        fmax = _double_min
2808        fmin = _double_max
2809        for _ in range(self.Ninit):
2810            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2811            fval = self.func(x0, *self.args)
2812            self.feval += 1
2813            if fval > fmax:
2814                fmax = fval
2815            if fval < fmin:
2816                fmin = fval
2817                best_state.cost = fval
2818                best_state.x = array(x0)
2819
2820        self.T0 = (fmax-fmin)*1.5
2821        return best_state.x
2822       
2823    def set_range(self,x0,frac):
2824        delrange = frac*(self.upper-self.lower)
2825        self.upper = x0+delrange
2826        self.lower = x0-delrange
2827
2828    def accept_test(self, dE):
2829        T = self.T
2830        self.tests += 1
2831        if dE < 0:
2832            self.accepted += 1
2833            return 1
2834        p = exp(-dE*1.0/self.boltzmann/T)
2835        if (p > random.uniform(0.0, 1.0)):
2836            self.accepted += 1
2837            return 1
2838        return 0
2839
2840    def update_guess(self, x0):
2841        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2842
2843    def update_temp(self, x0):
2844        pass
2845
2846
2847#  A schedule due to Lester Ingber modified to use bounds - OK
2848class fast_sa(base_schedule):
2849    def init(self, **options):
2850        self.__dict__.update(options)
2851        if self.m is None:
2852            self.m = 1.0
2853        if self.n is None:
2854            self.n = 1.0
2855        self.c = self.m * exp(-self.n * self.quench)
2856
2857    def update_guess(self, x0):
2858        x0 = asarray(x0)
2859        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2860        T = self.T
2861        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
2862        xnew = xc*(self.upper - self.lower)+self.lower
2863        return xnew
2864#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2865#        xc = y*(self.upper - self.lower)
2866#        xnew = x0 + xc
2867#        return xnew
2868
2869    def update_temp(self):
2870        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2871        self.k += 1
2872        return
2873
2874class cauchy_sa(base_schedule):     #modified to use bounds - not good
2875    def update_guess(self, x0):
2876        x0 = asarray(x0)
2877        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2878        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2879        xnew = xc*(self.upper - self.lower)+self.lower
2880        return xnew
2881#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2882#        xc = self.learn_rate * self.T * tan(numbers)
2883#        xnew = x0 + xc
2884#        return xnew
2885
2886    def update_temp(self):
2887        self.T = self.T0/(1+self.k)
2888        self.k += 1
2889        return
2890
2891class boltzmann_sa(base_schedule):
2892#    def update_guess(self, x0):
2893#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2894#        x0 = asarray(x0)
2895#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2896#
2897#        xnew = x0 + xc*std*self.learn_rate
2898#        return xnew
2899
2900    def update_temp(self):
2901        self.k += 1
2902        self.T = self.T0 / log(self.k+1.0)
2903        return
2904
2905class log_sa(base_schedule):        #OK
2906
2907    def init(self,**options):
2908        self.__dict__.update(options)
2909       
2910    def update_guess(self,x0):     #same as default
2911        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2912       
2913    def update_temp(self):
2914        self.k += 1
2915        self.T = self.T0*self.slope**self.k
2916       
2917class _state(object):
2918    def __init__(self):
2919        self.x = None
2920        self.cost = None
2921
2922# TODO:
2923#     allow for general annealing temperature profile
2924#     in that case use update given by alpha and omega and
2925#     variation of all previous updates and temperature?
2926
2927# Simulated annealing
2928
2929def anneal(func, x0, args=(), schedule='fast', full_output=0,
2930           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2931           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2932           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
2933           ranRange=0.10,autoRan=False,dlg=None):
2934    """Minimize a function using simulated annealing.
2935
2936    Schedule is a schedule class implementing the annealing schedule.
2937    Available ones are 'fast', 'cauchy', 'boltzmann'
2938
2939    :param callable func: f(x, \*args)
2940        Function to be optimized.
2941    :param ndarray x0:
2942        Initial guess.
2943    :param tuple args:
2944        Extra parameters to `func`.
2945    :param base_schedule schedule:
2946        Annealing schedule to use (a class).
2947    :param bool full_output:
2948        Whether to return optional outputs.
2949    :param float T0:
2950        Initial Temperature (estimated as 1.2 times the largest
2951        cost-function deviation over random points in the range).
2952    :param float Tf:
2953        Final goal temperature.
2954    :param int maxeval:
2955        Maximum function evaluations.
2956    :param int maxaccept:
2957        Maximum changes to accept.
2958    :param int maxiter:
2959        Maximum cooling iterations.
2960    :param float learn_rate:
2961        Scale constant for adjusting guesses.
2962    :param float boltzmann:
2963        Boltzmann constant in acceptance test
2964        (increase for less stringent test at each temperature).
2965    :param float feps:
2966        Stopping relative error tolerance for the function value in
2967        last four coolings.
2968    :param float quench,m,n:
2969        Parameters to alter fast_sa schedule.
2970    :param float/ndarray lower,upper:
2971        Lower and upper bounds on `x`.
2972    :param int dwell:
2973        The number of times to search the space at each temperature.
2974    :param float slope:
2975        Parameter for log schedule
2976    :param bool ranStart=False:
2977        True for set 10% of ranges about x
2978
2979    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2980
2981     * xmin (ndarray): Point giving smallest value found.
2982     * Jmin (float): Minimum value of function found.
2983     * T (float): Final temperature.
2984     * feval (int): Number of function evaluations.
2985     * iters (int): Number of cooling iterations.
2986     * accept (int): Number of tests accepted.
2987     * retval (int): Flag indicating stopping condition:
2988
2989              *  0: Points no longer changing
2990              *  1: Cooled to final temperature
2991              *  2: Maximum function evaluations
2992              *  3: Maximum cooling iterations reached
2993              *  4: Maximum accepted query locations reached
2994              *  5: Final point not the minimum amongst encountered points
2995
2996    *Notes*:
2997    Simulated annealing is a random algorithm which uses no derivative
2998    information from the function being optimized. In practice it has
2999    been more useful in discrete optimization than continuous
3000    optimization, as there are usually better algorithms for continuous
3001    optimization problems.
3002
3003    Some experimentation by trying the difference temperature
3004    schedules and altering their parameters is likely required to
3005    obtain good performance.
3006
3007    The randomness in the algorithm comes from random sampling in numpy.
3008    To obtain the same results you can call numpy.random.seed with the
3009    same seed immediately before calling scipy.optimize.anneal.
3010
3011    We give a brief description of how the three temperature schedules
3012    generate new points and vary their temperature. Temperatures are
3013    only updated with iterations in the outer loop. The inner loop is
3014    over xrange(dwell), and new points are generated for
3015    every iteration in the inner loop. (Though whether the proposed
3016    new points are accepted is probabilistic.)
3017
3018    For readability, let d denote the dimension of the inputs to func.
3019    Also, let x_old denote the previous state, and k denote the
3020    iteration number of the outer loop. All other variables not
3021    defined below are input variables to scipy.optimize.anneal itself.
3022
3023    In the 'fast' schedule the updates are ::
3024
3025        u ~ Uniform(0, 1, size=d)
3026        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
3027        xc = y * (upper - lower)
3028        x_new = x_old + xc
3029
3030        c = n * exp(-n * quench)
3031        T_new = T0 * exp(-c * k**quench)
3032
3033
3034    In the 'cauchy' schedule the updates are ::
3035
3036        u ~ Uniform(-pi/2, pi/2, size=d)
3037        xc = learn_rate * T * tan(u)
3038        x_new = x_old + xc
3039
3040        T_new = T0 / (1+k)
3041
3042    In the 'boltzmann' schedule the updates are ::
3043
3044        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
3045        y ~ Normal(0, std, size=d)
3046        x_new = x_old + learn_rate * y
3047
3048        T_new = T0 / log(1+k)
3049
3050    """
3051    x0 = asarray(x0)
3052    lower = asarray(lower)
3053    upper = asarray(upper)
3054
3055    schedule = eval(schedule+'_sa()')
3056    #   initialize the schedule
3057    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
3058                  learn_rate=learn_rate, lower=lower, upper=upper,
3059                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
3060
3061    current_state, last_state, best_state = _state(), _state(), _state()
3062    if ranStart:
3063        schedule.set_range(x0,ranRange)
3064    if T0 is None:
3065        x0 = schedule.getstart_temp(best_state)
3066    else:
3067        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
3068        best_state.x = None
3069        best_state.cost = numpy.Inf
3070
3071    last_state.x = asarray(x0).copy()
3072    fval = func(x0,*args)
3073    schedule.feval += 1
3074    last_state.cost = fval
3075    if last_state.cost < best_state.cost:
3076        best_state.cost = fval
3077        best_state.x = asarray(x0).copy()
3078    schedule.T = schedule.T0
3079    fqueue = [100, 300, 500, 700]
3080    iters = 1
3081    keepGoing = True
3082    bestn = 0
3083    while keepGoing:
3084        retval = 0
3085        for n in xrange(dwell):
3086            current_state.x = schedule.update_guess(last_state.x)
3087            current_state.cost = func(current_state.x,*args)
3088            schedule.feval += 1
3089
3090            dE = current_state.cost - last_state.cost
3091            if schedule.accept_test(dE):
3092                last_state.x = current_state.x.copy()
3093                last_state.cost = current_state.cost
3094                if last_state.cost < best_state.cost:
3095                    best_state.x = last_state.x.copy()
3096                    best_state.cost = last_state.cost
3097                    bestn = n
3098                    if best_state.cost < 1.0 and autoRan:
3099                        schedule.set_range(x0,best_state.cost/2.)                       
3100        if dlg:
3101            GoOn = dlg.Update(min(100.,best_state.cost*100),
3102                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
3103                    'Best trial:',bestn,  \
3104                    'MC/SA Residual:',best_state.cost*100,'%', \
3105                    ))[0]
3106            if not GoOn:
3107                best_state.x = last_state.x.copy()
3108                best_state.cost = last_state.cost
3109                retval = 5
3110        schedule.update_temp()
3111        iters += 1
3112        # Stopping conditions
3113        # 0) last saved values of f from each cooling step
3114        #     are all very similar (effectively cooled)
3115        # 1) Tf is set and we are below it
3116        # 2) maxeval is set and we are past it
3117        # 3) maxiter is set and we are past it
3118        # 4) maxaccept is set and we are past it
3119        # 5) user canceled run via progress bar
3120
3121        fqueue.append(squeeze(last_state.cost))
3122        fqueue.pop(0)
3123        af = asarray(fqueue)*1.0
3124        if retval == 5:
3125            print ' User terminated run; incomplete MC/SA'
3126            keepGoing = False
3127            break
3128        if all(abs((af-af[0])/af[0]) < feps):
3129            retval = 0
3130            if abs(af[-1]-best_state.cost) > feps*10:
3131                retval = 5
3132#                print "Warning: Cooled to %f at %s but this is not" \
3133#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
3134#                      + " the smallest point found."
3135            break
3136        if (Tf is not None) and (schedule.T < Tf):
3137            retval = 1
3138            break
3139        if (maxeval is not None) and (schedule.feval > maxeval):
3140            retval = 2
3141            break
3142        if (iters > maxiter):
3143            print "Warning: Maximum number of iterations exceeded."
3144            retval = 3
3145            break
3146        if (maxaccept is not None) and (schedule.accepted > maxaccept):
3147            retval = 4
3148            break
3149
3150    if full_output:
3151        return best_state.x, best_state.cost, schedule.T, \
3152               schedule.feval, iters, schedule.accepted, retval
3153    else:
3154        return best_state.x, retval
3155
3156def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
3157    outlist = []
3158    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
3159    for n in range(iCyc):
3160        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
3161        outlist.append(result[0])
3162        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
3163    out_q.put(outlist)
3164
3165def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
3166    import multiprocessing as mp
3167   
3168    nprocs = mp.cpu_count()
3169    out_q = mp.Queue()
3170    procs = []
3171    iCyc = np.zeros(nprocs)
3172    for i in range(nCyc):
3173        iCyc[i%nprocs] += 1
3174    for i in range(nprocs):
3175        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
3176        procs.append(p)
3177        p.start()
3178    resultlist = []
3179    for i in range(nprocs):
3180        resultlist += out_q.get()
3181    for p in procs:
3182        p.join()
3183    return resultlist
3184
3185def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
3186    '''default doc string
3187   
3188    :param type name: description
3189   
3190    :returns: type name: description
3191   
3192    '''
3193   
3194    global tsum
3195    tsum = 0.
3196   
3197    def getMDparms(item,pfx,parmDict,varyList):
3198        parmDict[pfx+'MDaxis'] = item['axis']
3199        parmDict[pfx+'MDval'] = item['Coef'][0]
3200        if item['Coef'][1]:
3201            varyList += [pfx+'MDval',]
3202            limits = item['Coef'][2]
3203            lower.append(limits[0])
3204            upper.append(limits[1])
3205                       
3206    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
3207        parmDict[pfx+'Atype'] = item['atType']
3208        aTypes |= set([item['atType'],]) 
3209        pstr = ['Ax','Ay','Az']
3210        XYZ = [0,0,0]
3211        for i in range(3):
3212            name = pfx+pstr[i]
3213            parmDict[name] = item['Pos'][0][i]
3214            XYZ[i] = parmDict[name]
3215            if item['Pos'][1][i]:
3216                varyList += [name,]
3217                limits = item['Pos'][2][i]
3218                lower.append(limits[0])
3219                upper.append(limits[1])
3220        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
3221           
3222    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
3223        parmDict[mfx+'MolCent'] = item['MolCent']
3224        parmDict[mfx+'RBId'] = item['RBId']
3225        pstr = ['Px','Py','Pz']
3226        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
3227        for i in range(3):
3228            name = mfx+pstr[i]
3229            parmDict[name] = item['Pos'][0][i]
3230            if item['Pos'][1][i]:
3231                varyList += [name,]
3232                limits = item['Pos'][2][i]
3233                lower.append(limits[0])
3234                upper.append(limits[1])
3235        AV = item['Ori'][0]
3236        A = AV[0]
3237        V = AV[1:]
3238        for i in range(4):
3239            name = mfx+ostr[i]
3240            if i:
3241                parmDict[name] = V[i-1]
3242            else:
3243                parmDict[name] = A
3244            if item['Ovar'] == 'AV':
3245                varyList += [name,]
3246                limits = item['Ori'][2][i]
3247                lower.append(limits[0])
3248                upper.append(limits[1])
3249            elif item['Ovar'] == 'A' and not i:
3250                varyList += [name,]
3251                limits = item['Ori'][2][i]
3252                lower.append(limits[0])
3253                upper.append(limits[1])
3254        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
3255            for i in range(len(item['Tor'][0])):
3256                name = mfx+'Tor'+str(i)
3257                parmDict[name] = item['Tor'][0][i]
3258                if item['Tor'][1][i]:
3259                    varyList += [name,]
3260                    limits = item['Tor'][2][i]
3261                    lower.append(limits[0])
3262                    upper.append(limits[1])
3263        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
3264        aTypes |= set(atypes)
3265        atNo += len(atypes)
3266        return atNo
3267               
3268    def GetAtomM(Xdata,SGData):
3269        Mdata = []
3270        for xyz in Xdata:
3271            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
3272        return np.array(Mdata)
3273       
3274    def GetAtomT(RBdata,parmDict):
3275        'Needs a doc string'
3276        atNo = parmDict['atNo']
3277        nfixAt = parmDict['nfixAt']
3278        Tdata = atNo*[' ',]
3279        for iatm in range(nfixAt):
3280            parm = ':'+str(iatm)+':Atype'
3281            if parm in parmDict:
3282                Tdata[iatm] = aTypes.index(parmDict[parm])
3283        iatm = nfixAt
3284        for iObj in range(parmDict['nObj']):
3285            pfx = str(iObj)+':'
3286            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3287                if parmDict[pfx+'Type'] == 'Vector':
3288                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3289                    nAtm = len(RBRes['rbVect'][0])
3290                else:       #Residue
3291                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3292                    nAtm = len(RBRes['rbXYZ'])
3293                for i in range(nAtm):
3294                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
3295                    iatm += 1
3296            elif parmDict[pfx+'Type'] == 'Atom':
3297                atNo = parmDict[pfx+'atNo']
3298                parm = pfx+'Atype'              #remove extra ':'
3299                if parm in parmDict:
3300                    Tdata[atNo] = aTypes.index(parmDict[parm])
3301                iatm += 1
3302            else:
3303                continue        #skips March Dollase
3304        return Tdata
3305       
3306    def GetAtomX(RBdata,parmDict):
3307        'Needs a doc string'
3308        Bmat = parmDict['Bmat']
3309        atNo = parmDict['atNo']
3310        nfixAt = parmDict['nfixAt']
3311        Xdata = np.zeros((3,atNo))
3312        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
3313        for iatm in range(nfixAt):
3314            for key in keys:
3315                parm = ':'+str(iatm)+key
3316                if parm in parmDict:
3317                    keys[key][iatm] = parmDict[parm]
3318        iatm = nfixAt
3319        for iObj in range(parmDict['nObj']):
3320            pfx = str(iObj)+':'
3321            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3322                if parmDict[pfx+'Type'] == 'Vector':
3323                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3324                    vecs = RBRes['rbVect']
3325                    mags = RBRes['VectMag']
3326                    Cart = np.zeros_like(vecs[0])
3327                    for vec,mag in zip(vecs,mags):
3328                        Cart += vec*mag
3329                elif parmDict[pfx+'Type'] == 'Residue':
3330                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3331                    Cart = np.array(RBRes['rbXYZ'])
3332                    for itor,seq in enumerate(RBRes['rbSeq']):
3333                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
3334                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
3335                if parmDict[pfx+'MolCent'][1]:
3336                    Cart -= parmDict[pfx+'MolCent'][0]
3337                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3338                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
3339                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
3340                iatm += len(Cart)
3341            elif parmDict[pfx+'Type'] == 'Atom':
3342                atNo = parmDict[pfx+'atNo']
3343                for key in keys:
3344                    parm = pfx+key[1:]              #remove extra ':'
3345                    if parm in parmDict:
3346                        keys[key][atNo] = parmDict[parm]
3347                iatm += 1
3348            else:
3349                continue        #skips March Dollase
3350        return Xdata.T
3351       
3352    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
3353        allX = np.inner(Xdata,SGM)+SGT
3354        allT = np.repeat(Tdata,allX.shape[1])
3355        allM = np.repeat(Mdata,allX.shape[1])
3356        allX = np.reshape(allX,(-1,3))
3357        return allT,allM,allX
3358
3359    def getAllX(Xdata,SGM,SGT):
3360        allX = np.inner(Xdata,SGM)+SGT
3361        allX = np.reshape(allX,(-1,3))
3362        return allX
3363       
3364    def normQuaternions(RBdata,parmDict,varyList,values):
3365        for iObj in range(parmDict['nObj']):
3366            pfx = str(iObj)+':'
3367            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3368                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3369                A,V = Q2AVdeg(Qori)
3370                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
3371                    if i:
3372                        parmDict[pfx+name] = V[i-1]
3373                    else:
3374                        parmDict[pfx+name] = A
3375       
3376    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
3377        ''' Compute structure factors for all h,k,l for phase
3378        input:
3379            refList: [ref] where each ref = h,k,l,m,d,...
3380            rcov:   array[nref,nref] covariance terms between Fo^2 values
3381            ifInv:  bool True if centrosymmetric
3382            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
3383            RBdata: [dict] rigid body dictionary
3384            varyList: [list] names of varied parameters in MC/SA (not used here)           
3385            ParmDict: [dict] problem parameters
3386        puts result F^2 in each ref[5] in refList
3387        returns:
3388            delt-F*rcov*delt-F/sum(Fo^2)^2
3389        '''       
3390        global tsum
3391        t0 = time.time()
3392        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
3393        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
3394        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
3395        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
3396        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
3397        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
3398        refList[5] = Aterm
3399        if not ifInv:
3400            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
3401        if len(cosTable):        #apply MD correction
3402            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
3403        sumFcsq = np.sum(refList[5])
3404        scale = parmDict['sumFosq']/sumFcsq
3405        refList[5] *= scale
3406        refList[6] = refList[4]-refList[5]
3407        M = np.inner(refList[6],np.inner(rcov,refList[6]))
3408        tsum += (time.time()-t0)
3409        return M/np.sum(refList[4]**2)
3410
3411    sq8ln2 = np.sqrt(8*np.log(2))
3412    sq2pi = np.sqrt(2*np.pi)
3413    sq4pi = np.sqrt(4*np.pi)
3414    generalData = data['General']
3415    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3416    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
3417    SGData = generalData['SGData']
3418    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
3419    SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
3420    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
3421    fixAtoms = data['Atoms']                       #if any
3422    cx,ct,cs = generalData['AtomPtrs'][:3]
3423    aTypes = set([])
3424    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
3425    varyList = []
3426    atNo = 0
3427    for atm in fixAtoms:
3428        pfx = ':'+str(atNo)+':'
3429        parmDict[pfx+'Atype'] = atm[ct]
3430        aTypes |= set([atm[ct],])
3431        pstr = ['Ax','Ay','Az']
3432        parmDict[pfx+'Amul'] = atm[cs+1]
3433        for i in range(3):
3434            name = pfx+pstr[i]
3435            parmDict[name] = atm[cx+i]
3436        atNo += 1
3437    parmDict['nfixAt'] = len(fixAtoms)       
3438    MCSA = generalData['MCSA controls']
3439    reflName = MCSA['Data source']
3440    phaseName = generalData['Name']
3441    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
3442    upper = []
3443    lower = []
3444    MDvec = np.zeros(3)
3445    for i,item in enumerate(MCSAObjs):
3446        mfx = str(i)+':'
3447        parmDict[mfx+'Type'] = item['Type']
3448        if item['Type'] == 'MD':
3449            getMDparms(item,mfx,parmDict,varyList)
3450            MDvec = np.array(item['axis'])
3451        elif item['Type'] == 'Atom':
3452            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
3453            parmDict[mfx+'atNo'] = atNo
3454            atNo += 1
3455        elif item['Type'] in ['Residue','Vector']:
3456            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
3457    parmDict['atNo'] = atNo                 #total no. of atoms
3458    parmDict['nObj'] = len(MCSAObjs)
3459    aTypes = list(aTypes)
3460    Tdata = GetAtomT(RBdata,parmDict)
3461    Xdata = GetAtomX(RBdata,parmDict)
3462    Mdata = GetAtomM(Xdata,SGData)
3463    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
3464    FFtables = G2el.GetFFtable(aTypes)
3465    refs = []
3466    allFF = []
3467    cosTable = []
3468    sumFosq = 0
3469    if 'PWDR' in reflName:
3470        for ref in reflData:
3471            h,k,l,m,d,pos,sig,gam,f = ref[:9]
3472            if d >= MCSA['dmin']:
3473                sig = G2pwd.getgamFW(sig,gam)/sq8ln2        #--> sig from FWHM
3474                SQ = 0.25/d**2
3475                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3476                refs.append([h,k,l,m,f*m,pos,sig])
3477                sumFosq += f*m
3478                Heqv = np.inner(np.array([h,k,l]),SGMT)
3479                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3480        nRef = len(refs)
3481        cosTable = np.array(cosTable)**2
3482        rcov = np.zeros((nRef,nRef))
3483        for iref,refI in enumerate(refs):
3484            rcov[iref][iref] = 1./(sq4pi*refI[6])
3485            for jref,refJ in enumerate(refs[:iref]):
3486                t1 = refI[6]**2+refJ[6]**2
3487                t2 = (refJ[5]-refI[5])**2/(2.*t1)
3488                if t2 > 10.:
3489                    rcov[iref][jref] = 0.
3490                else:
3491                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
3492        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3493        Rdiag = np.sqrt(np.diag(rcov))
3494        Rnorm = np.outer(Rdiag,Rdiag)
3495        rcov /= Rnorm
3496    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
3497        vNames = []
3498        pfx = str(data['pId'])+'::PWLref:'
3499        for iref,refI in enumerate(reflData):           #Pawley reflection set
3500            h,k,l,m,d,v,f,s = refI
3501            if d >= MCSA['dmin'] and v:       #skip unrefined ones
3502                vNames.append(pfx+str(iref))
3503                SQ = 0.25/d**2
3504                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3505                refs.append([h,k,l,m,f*m,iref,0.])
3506                sumFosq += f*m
3507                Heqv = np.inner(np.array([h,k,l]),SGMT)
3508                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3509        cosTable = np.array(cosTable)**2
3510        nRef = len(refs)
3511#        if generalData['doPawley'] and (covData['freshCOV'] or  MCSA['newDmin']):
3512        if covData['freshCOV'] or  MCSA['newDmin']:
3513            vList = covData['varyList']
3514            covMatrix = covData['covMatrix']
3515            rcov = getVCov(vNames,vList,covMatrix)
3516            rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3517            Rdiag = np.sqrt(np.diag(rcov))
3518            Rdiag = np.where(Rdiag,Rdiag,1.0)
3519            Rnorm = np.outer(Rdiag,Rdiag)
3520            rcov /= Rnorm
3521            MCSA['rcov'] = rcov
3522            covData['freshCOV'] = False
3523            MCSA['newDmin'] = False
3524        else:
3525            rcov = MCSA['rcov']
3526    elif 'HKLF' in reflName:
3527        for ref in reflData:
3528            [h,k,l,m,d],f = ref[:5],ref[6]
3529            if d >= MCSA['dmin']:
3530                SQ = 0.25/d**2
3531                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3532                refs.append([h,k,l,m,f*m,0.,0.])
3533                sumFosq += f*m
3534        nRef = len(refs)
3535        rcov = np.identity(len(refs))
3536    allFF = np.array(allFF).T
3537    refs = np.array(refs).T
3538    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
3539    print ' Number of parameters varied: %d'%(len(varyList))
3540    parmDict['sumFosq'] = sumFosq
3541    x0 = [parmDict[val] for val in varyList]
3542    ifInv = SGData['SGInv']
3543    # consider replacing anneal with scipy.optimize.basinhopping
3544    results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
3545        schedule=MCSA['Algorithm'], full_output=True,
3546        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
3547        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
3548        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
3549        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
3550        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
3551    M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
3552#    for ref in refs.T:
3553#        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])
3554#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
3555    Result = [False,False,results[1],results[2],]+list(results[0])
3556    Result.append(varyList)
3557    return Result,tsum
3558
3559       
3560################################################################################
3561##### Quaternion stuff
3562################################################################################
3563
3564def prodQQ(QA,QB):
3565    ''' Grassman quaternion product
3566        QA,QB quaternions; q=r+ai+bj+ck
3567    '''
3568    D = np.zeros(4)
3569    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
3570    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
3571    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
3572    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
3573   
3574#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
3575#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
3576   
3577    return D
3578   
3579def normQ(QA):
3580    ''' get length of quaternion & normalize it
3581        q=r+ai+bj+ck
3582    '''
3583    n = np.sqrt(np.sum(np.array(QA)**2))
3584    return QA/n
3585   
3586def invQ(Q):
3587    '''
3588        get inverse of quaternion
3589        q=r+ai+bj+ck; q* = r-ai-bj-ck
3590    '''
3591    return Q*np.array([1,-1,-1,-1])
3592   
3593def prodQVQ(Q,V):
3594    """
3595    compute the quaternion vector rotation qvq-1 = v'
3596    q=r+ai+bj+ck
3597    """
3598    T2 = Q[0]*Q[1]
3599    T3 = Q[0]*Q[2]
3600    T4 = Q[0]*Q[3]
3601    T5 = -Q[1]*Q[1]
3602    T6 = Q[1]*Q[2]
3603    T7 = Q[1]*Q[3]
3604    T8 = -Q[2]*Q[2]
3605    T9 = Q[2]*Q[3]
3606    T10 = -Q[3]*Q[3]
3607    M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
3608    VP = 2.*np.inner(V,M)
3609    return VP+V
3610   
3611def Q2Mat(Q):
3612    ''' make rotation matrix from quaternion
3613        q=r+ai+bj+ck
3614    '''
3615    QN = normQ(Q)
3616    aa = QN[0]**2
3617    ab = QN[0]*QN[1]
3618    ac = QN[0]*QN[2]
3619    ad = QN[0]*QN[3]
3620    bb = QN[1]**2
3621    bc = QN[1]*QN[2]
3622    bd = QN[1]*QN[3]
3623    cc = QN[2]**2
3624    cd = QN[2]*QN[3]
3625    dd = QN[3]**2
3626    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
3627        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
3628        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
3629    return np.array(M)
3630   
3631def AV2Q(A,V):
3632    ''' convert angle (radians) & vector to quaternion
3633        q=r+ai+bj+ck
3634    '''
3635    Q = np.zeros(4)
3636    d = nl.norm(np.array(V))
3637    if d:
3638        V /= d
3639        if not A:       #==0.
3640            A = 2.*np.pi
3641        p = A/2.
3642        Q[0] = np.cos(p)
3643        Q[1:4] = V*np.sin(p)
3644    else:
3645        Q[3] = 1.
3646    return Q
3647   
3648def AVdeg2Q(A,V):
3649    ''' convert angle (degrees) & vector to quaternion
3650        q=r+ai+bj+ck
3651    '''
3652    Q = np.zeros(4)
3653    d = nl.norm(np.array(V))
3654    if not A:       #== 0.!
3655        A = 360.
3656    if d:
3657        V /= d
3658        p = A/2.
3659        Q[0] = cosd(p)
3660        Q[1:4] = V*sind(p)
3661    else:
3662        Q[3] = 1.
3663    return Q
3664   
3665def Q2AVdeg(Q):
3666    ''' convert quaternion to angle (degrees 0-360) & normalized vector
3667        q=r+ai+bj+ck
3668    '''
3669    A = 2.*acosd(Q[0])
3670    V = np.array(Q[1:])
3671    V /= sind(A/2.)
3672    return A,V
3673   
3674def Q2AV(Q):
3675    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
3676        q=r+ai+bj+ck
3677    '''
3678    A = 2.*np.arccos(Q[0])
3679    V = np.array(Q[1:])
3680    V /= np.sin(A/2.)
3681    return A,V
3682   
3683def randomQ(r0,r1,r2,r3):
3684    ''' create random quaternion from 4 random numbers in range (-1,1)
3685    '''
3686    sum = 0
3687    Q = np.array(4)
3688    Q[0] = r0
3689    sum += Q[0]**2
3690    Q[1] = np.sqrt(1.-sum)*r1
3691    sum += Q[1]**2
3692    Q[2] = np.sqrt(1.-sum)*r2
3693    sum += Q[2]**2
3694    Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
3695    return Q
3696   
3697def randomAVdeg(r0,r1,r2,r3):
3698    ''' create random angle (deg),vector from 4 random number in range (-1,1)
3699    '''
3700    return Q2AVdeg(randomQ(r0,r1,r2,r3))
3701   
3702def makeQuat(A,B,C):
3703    ''' Make quaternion from rotation of A vector to B vector about C axis
3704
3705        :param np.array A,B,C: Cartesian 3-vectors
3706        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
3707    '''
3708
3709    V1 = np.cross(A,C)
3710    V2 = np.cross(B,C)
3711    if nl.norm(V1)*nl.norm(V2):
3712        V1 /= nl.norm(V1)
3713        V2 /= nl.norm(V2)
3714        V3 = np.cross(V1,V2)
3715    else:
3716        V3 = np.zeros(3)
3717    Q = np.array([0.,0.,0.,1.])
3718    D = 0.
3719    if nl.norm(V3):
3720        V3 /= nl.norm(V3)
3721        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
3722        D = np.arccos(D1)/2.0
3723        V1 = C-V3
3724        V2 = C+V3
3725        DM = nl.norm(V1)
3726        DP = nl.norm(V2)
3727        S = np.sin(D)
3728        Q[0] = np.cos(D)
3729        Q[1:] = V3*S
3730        D *= 2.
3731        if DM > DP:
3732            D *= -1.
3733    return Q,D
3734   
3735def annealtests():
3736    from numpy import cos
3737    # minimum expected at ~-0.195
3738    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
3739    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
3740    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
3741    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
3742
3743    # minimum expected at ~[-0.195, -0.1]
3744    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
3745    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')
3746    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')
3747    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')
3748
3749
3750if __name__ == '__main__':
3751    annealtests()
Note: See TracBrowser for help on using the repository browser.