source: trunk/GSASIImath.py @ 2330

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

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