source: trunk/GSASIImath.py @ 2448

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

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