source: trunk/GSASIImath.py @ 1837

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

fix Mass problem on 1st LS if General not visited after atoms added
add (hidden for now) penalty option for texture refinement

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