source: trunk/GSASIImath.py @ 2764

Last change on this file since 2764 was 2764, checked in by toby, 6 years ago

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