source: trunk/GSASIImath.py @ 2313

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

remove unused pylab reference

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