source: trunk/GSASIImath.py @ 2314

Last change on this file since 2314 was 2314, checked in by vondreele, 6 years ago

change Levenberg/Marquardt? lambda max to 103 from 105 (it was excessive)
fix image loading on tree selection; now any pick of IMG child will reload image

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