source: trunk/GSASIImath.py @ 2839

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

fix numpy divide issues inquaternions

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