source: trunk/GSASIImath.py @ 2132

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

more rationalization of FWHM calculations

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