source: trunk/GSASIImath.py @ 2166

Last change on this file since 2166 was 2166, checked in by vondreele, 6 years ago

add GetMeanWave? to G2math to return Ka-bar

more work on Stacking faults - now runs DIFFaX with data from GSAS-II

result as a powder pattern returned to ycalc.

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