source: trunk/GSASIImath.py @ 2522

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

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