source: trunk/GSASIImath.py @ 2755

Last change on this file since 2755 was 2755, checked in by vondreele, 5 years ago

extend Export HKLs for powder data to include most all columns in GSASII.py
narrow HorizontalLine? display
add a ReloadSubstances? option in case user changed wavelength in SASD & REFD substances
Add a default unit scatter to substances - only in new projects
add modeling of REFD patterns for x-rays (not yet tested - no data) & CW neutrons (tested)
modify XScattDen to give clearly fo, f' & f" components
add NCScattDen for CW neutrons - uses wave to generate real/imaginary components
fix printing of scattering density units - now correct symbols
remove unused imports from GSASIIsasd.py

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