source: trunk/GSASIImath.py @ 1990

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

Structure factor derivs for SS continued...

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