source: trunk/GSASIImath.py @ 1992

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

fix SS position derivative error

  • 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 20:35:11 +0000 (Wed, 07 Oct 2015) $
5# $Author: vondreele $
6# $Revision: 1992 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 1992 2015-10-07 20:35:11Z 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: 1992 $")
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]),np.concatenate((dGdAf,dGdBf),-1), \
1065        np.concatenate((dGdAx,dGdBx),-1),np.concatenate((dGdAu,dGdBu),-1)      #ops X atoms
1066   
1067def posFourier(tau,psin,pcos,smul):
1068    A = np.array([ps[:,np.newaxis]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)])*smul
1069    B = np.array([pc[:,np.newaxis]*np.cos(2*np.pi*(i+1)*tau) for i,pc in enumerate(pcos)])
1070    return np.sum(A,axis=0)+np.sum(B,axis=0)
1071   
1072def posSawtooth(tau,Toff,slopes):
1073    Tau = (tau-Toff)%1.
1074    A = slopes[:,np.newaxis]*Tau
1075    return A
1076
1077def posZigZag(tau,Toff,slopes):
1078    Tau = (tau-Toff)%1.
1079    A = np.where(Tau <= 0.5,slopes[:,np.newaxis]*Tau,slopes[:,np.newaxis]*(1.-Tau))
1080    return A
1081   
1082def fracCrenel(tau,Toff,Twid):
1083    Tau = (tau-Toff)%1.
1084    A = np.where(Tau<Twid,1.,0.)
1085    return A
1086   
1087def fracFourier(tau,fsin,fcos):
1088    A = np.array([fs[:,np.newaxis]*np.sin(2.*np.pi*(i+1)*tau) for i,fs in enumerate(fsin)])
1089    B = np.array([fc[:,np.newaxis]*np.cos(2.*np.pi*(i+1)*tau) for i,fc in enumerate(fcos)])
1090    return np.sum(A,axis=0)+np.sum(B,axis=0)
1091   
1092def ApplyModulation(data,tau):
1093    '''Applies modulation to drawing atom positions & Uijs for given tau
1094    '''
1095    generalData = data['General']
1096    SGData = generalData['SGData']
1097    SSGData = generalData['SSGData']
1098    cx,ct,cs,cia = generalData['AtomPtrs']
1099    drawingData = data['Drawing']
1100    dcx,dct,dcs,dci = drawingData['atomPtrs']
1101    atoms = data['Atoms']
1102    drawAtoms = drawingData['Atoms']
1103    Fade = np.zeros(len(drawAtoms))
1104    for atom in atoms:   
1105        atxyz = G2spc.MoveToUnitCell(np.array(atom[cx:cx+3]))[0]
1106        atuij = np.array(atom[cia+2:cia+8])
1107        waveType = atom[-1]['SS1']['waveType']
1108        Sfrac = atom[-1]['SS1']['Sfrac']
1109        Spos = atom[-1]['SS1']['Spos']
1110        Sadp = atom[-1]['SS1']['Sadp']
1111        indx = FindAtomIndexByIDs(drawAtoms,dci,[atom[cia+8],],True)
1112        for ind in indx:
1113            drawatom = drawAtoms[ind]
1114            opr = drawatom[dcs-1]
1115            sop,ssop,icent = G2spc.OpsfromStringOps(opr,SGData,SSGData)
1116            sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(tau,sop,ssop,atxyz)
1117            smul = sdet         # n-fold vs m operator on wave
1118            tauT *= icent       #invert wave on -1
1119            wave = np.zeros(3)
1120            uwave = np.zeros(6)
1121            #how do I handle Sfrac? - fade the atoms?
1122            if len(Sfrac):
1123                scof = []
1124                ccof = []
1125                for i,sfrac in enumerate(Sfrac):
1126                    if not i and 'Crenel' in waveType:
1127                        Fade[ind] += fracCrenel(tauT,sfrac[0][0],sfrac[0][1])
1128                    else:
1129                        scof.append(sfrac[0][0])
1130                        ccof.append(sfrac[0][1])
1131                    if len(scof):
1132                        Fade[ind] += np.sum(fracFourier(tauT,scof,ccof))                           
1133            if len(Spos):
1134                scof = []
1135                ccof = []
1136                for i,spos in enumerate(Spos):
1137                    if waveType in ['Sawtooth','ZigZag'] and not i:
1138                        Toff = spos[0][0]
1139                        slopes = np.array(spos[0][1:])
1140                        if waveType == 'Sawtooth':
1141                            wave = posSawtooth(tauT,Toff,slopes)
1142                        elif waveType == 'ZigZag':
1143                            wave = posZigZag(tauT,Toff,slopes)
1144                    else:
1145                        scof.append(spos[0][:3])
1146                        ccof.append(spos[0][3:])
1147                wave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof),smul),axis=1)
1148            if len(Sadp):
1149                scof = []
1150                ccof = []
1151                for i,sadp in enumerate(Sadp):
1152                    scof.append(sadp[0][:6])
1153                    ccof.append(sadp[0][6:])
1154                uwave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof),smul),axis=1)
1155            if atom[cia] == 'A':                   
1156                X,U = G2spc.ApplyStringOps(opr,SGData,atxyz+wave,atuij+uwave)
1157                drawatom[dcx:dcx+3] = X
1158                drawatom[dci-6:dci] = U
1159            else:
1160                X = G2spc.ApplyStringOps(opr,SGData,atxyz+wave)
1161                drawatom[dcx:dcx+3] = X
1162    return drawAtoms,Fade
1163   
1164# gauleg.py Gauss Legendre numerical quadrature, x and w computation
1165# integrate from a to b using n evaluations of the function f(x) 
1166# usage: from gauleg import gaulegf         
1167#        x,w = gaulegf( a, b, n)                               
1168#        area = 0.0                                           
1169#        for i in range(1,n+1):          #  yes, 1..n                   
1170#          area += w[i]*f(x[i])                                   
1171
1172import math
1173def gaulegf(a, b, n):
1174  x = range(n+1) # x[0] unused
1175  w = range(n+1) # w[0] unused
1176  eps = 3.0E-14
1177  m = (n+1)/2
1178  xm = 0.5*(b+a)
1179  xl = 0.5*(b-a)
1180  for i in range(1,m+1):
1181    z = math.cos(3.141592654*(i-0.25)/(n+0.5))
1182    while True:
1183      p1 = 1.0
1184      p2 = 0.0
1185      for j in range(1,n+1):
1186        p3 = p2
1187        p2 = p1
1188        p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
1189
1190      pp = n*(z*p1-p2)/(z*z-1.0)
1191      z1 = z
1192      z = z1 - p1/pp
1193      if abs(z-z1) <= eps:
1194            break
1195
1196    x[i] = xm - xl*z
1197    x[n+1-i] = xm + xl*z
1198    w[i] = 2.0*xl/((1.0-z*z)*pp*pp)
1199    w[n+1-i] = w[i]
1200  return np.array(x), np.array(w)
1201# end gaulegf
1202   
1203   
1204def BessJn(nmax,x):
1205    ''' compute Bessel function J(n,x) from scipy routine & recurrance relation
1206    returns sequence of J(n,x) for n in range [-nmax...0...nmax]
1207   
1208    :param  integer nmax: maximul order for Jn(x)
1209    :param  float x: argument for Jn(x)
1210   
1211    :returns numpy array: [J(-nmax,x)...J(0,x)...J(nmax,x)]
1212   
1213    '''
1214    import scipy.special as sp
1215    bessJn = np.zeros(2*nmax+1)
1216    bessJn[nmax] = sp.j0(x)
1217    bessJn[nmax+1] = sp.j1(x)
1218    bessJn[nmax-1] = -bessJn[nmax+1]
1219    for i in range(2,nmax+1):
1220        bessJn[i+nmax] = 2*(i-1)*bessJn[nmax+i-1]/x-bessJn[nmax+i-2]
1221        bessJn[nmax-i] = bessJn[i+nmax]*(-1)**i
1222    return bessJn
1223   
1224def BessIn(nmax,x):
1225    ''' compute modified Bessel function I(n,x) from scipy routines & recurrance relation
1226    returns sequence of I(n,x) for n in range [-nmax...0...nmax]
1227   
1228    :param  integer nmax: maximul order for In(x)
1229    :param  float x: argument for In(x)
1230   
1231    :returns numpy array: [I(-nmax,x)...I(0,x)...I(nmax,x)]
1232   
1233    '''
1234    import scipy.special as sp
1235    bessIn = np.zeros(2*nmax+1)
1236    bessIn[nmax] = sp.i0(x)
1237    bessIn[nmax+1] = sp.i1(x)
1238    bessIn[nmax-1] = bessIn[nmax+1]
1239    for i in range(2,nmax+1):
1240        bessIn[i+nmax] = bessIn[nmax+i-2]-2*(i-1)*bessIn[nmax+i-1]/x
1241        bessIn[nmax-i] = bessIn[i+nmax]
1242    return bessIn
1243       
1244   
1245################################################################################
1246##### distance, angle, planes, torsion stuff
1247################################################################################
1248
1249def getSyXYZ(XYZ,ops,SGData):
1250    '''default doc string
1251   
1252    :param type name: description
1253   
1254    :returns: type name: description
1255   
1256    '''
1257    XYZout = np.zeros_like(XYZ)
1258    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
1259        if op == '1':
1260            XYZout[i] = xyz
1261        else:
1262            oprs = op.split('+')
1263            unit = [0,0,0]
1264            if len(oprs)>1:
1265                unit = np.array(list(eval(oprs[1])))
1266            syop =int(oprs[0])
1267            inv = syop/abs(syop)
1268            syop *= inv
1269            cent = syop/100
1270            syop %= 100
1271            syop -= 1
1272            M,T = SGData['SGOps'][syop]
1273            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
1274    return XYZout
1275   
1276def getRestDist(XYZ,Amat):
1277    '''default doc string
1278   
1279    :param type name: description
1280   
1281    :returns: type name: description
1282   
1283    '''
1284    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
1285   
1286def getRestDeriv(Func,XYZ,Amat,ops,SGData):
1287    '''default doc string
1288   
1289    :param type name: description
1290   
1291    :returns: type name: description
1292   
1293    '''
1294    deriv = np.zeros((len(XYZ),3))
1295    dx = 0.00001
1296    for j,xyz in enumerate(XYZ):
1297        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1298            XYZ[j] -= x
1299            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1300            XYZ[j] += 2*x
1301            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1302            XYZ[j] -= x
1303            deriv[j][i] = (d1-d2)/(2*dx)
1304    return deriv.flatten()
1305
1306def getRestAngle(XYZ,Amat):
1307    '''default doc string
1308   
1309    :param type name: description
1310   
1311    :returns: type name: description
1312   
1313    '''
1314   
1315    def calcVec(Ox,Tx,Amat):
1316        return np.inner(Amat,(Tx-Ox))
1317
1318    VecA = calcVec(XYZ[1],XYZ[0],Amat)
1319    VecA /= np.sqrt(np.sum(VecA**2))
1320    VecB = calcVec(XYZ[1],XYZ[2],Amat)
1321    VecB /= np.sqrt(np.sum(VecB**2))
1322    edge = VecB-VecA
1323    edge = np.sum(edge**2)
1324    angle = (2.-edge)/2.
1325    angle = max(angle,-1.)
1326    return acosd(angle)
1327   
1328def getRestPlane(XYZ,Amat):
1329    '''default doc string
1330   
1331    :param type name: description
1332   
1333    :returns: type name: description
1334   
1335    '''
1336    sumXYZ = np.zeros(3)
1337    for xyz in XYZ:
1338        sumXYZ += xyz
1339    sumXYZ /= len(XYZ)
1340    XYZ = np.array(XYZ)-sumXYZ
1341    XYZ = np.inner(Amat,XYZ).T
1342    Zmat = np.zeros((3,3))
1343    for i,xyz in enumerate(XYZ):
1344        Zmat += np.outer(xyz.T,xyz)
1345    Evec,Emat = nl.eig(Zmat)
1346    Evec = np.sqrt(Evec)/(len(XYZ)-3)
1347    Order = np.argsort(Evec)
1348    return Evec[Order[0]]
1349   
1350def getRestChiral(XYZ,Amat):   
1351    '''default doc string
1352   
1353    :param type name: description
1354   
1355    :returns: type name: description
1356   
1357    '''
1358    VecA = np.empty((3,3))   
1359    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
1360    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
1361    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
1362    return nl.det(VecA)
1363   
1364def getRestTorsion(XYZ,Amat):
1365    '''default doc string
1366   
1367    :param type name: description
1368   
1369    :returns: type name: description
1370   
1371    '''
1372    VecA = np.empty((3,3))
1373    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
1374    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
1375    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
1376    D = nl.det(VecA)
1377    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
1378    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
1379    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
1380    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
1381    Ang = 1.0
1382    if abs(P12) < 1.0 and abs(P23) < 1.0:
1383        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
1384    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
1385    return TOR
1386   
1387def calcTorsionEnergy(TOR,Coeff=[]):
1388    '''default doc string
1389   
1390    :param type name: description
1391   
1392    :returns: type name: description
1393   
1394    '''
1395    sum = 0.
1396    if len(Coeff):
1397        cof = np.reshape(Coeff,(3,3)).T
1398        delt = TOR-cof[1]
1399        delt = np.where(delt<-180.,delt+360.,delt)
1400        delt = np.where(delt>180.,delt-360.,delt)
1401        term = -cof[2]*delt**2
1402        val = cof[0]*np.exp(term/1000.0)
1403        pMax = cof[0][np.argmin(val)]
1404        Eval = np.sum(val)
1405        sum = Eval-pMax
1406    return sum,Eval
1407
1408def getTorsionDeriv(XYZ,Amat,Coeff):
1409    '''default doc string
1410   
1411    :param type name: description
1412   
1413    :returns: type name: description
1414   
1415    '''
1416    deriv = np.zeros((len(XYZ),3))
1417    dx = 0.00001
1418    for j,xyz in enumerate(XYZ):
1419        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1420            XYZ[j] -= x
1421            tor = getRestTorsion(XYZ,Amat)
1422            p1,d1 = calcTorsionEnergy(tor,Coeff)
1423            XYZ[j] += 2*x
1424            tor = getRestTorsion(XYZ,Amat)
1425            p2,d2 = calcTorsionEnergy(tor,Coeff)           
1426            XYZ[j] -= x
1427            deriv[j][i] = (p2-p1)/(2*dx)
1428    return deriv.flatten()
1429
1430def getRestRama(XYZ,Amat):
1431    '''Computes a pair of torsion angles in a 5 atom string
1432   
1433    :param nparray XYZ: crystallographic coordinates of 5 atoms
1434    :param nparray Amat: crystal to cartesian transformation matrix
1435   
1436    :returns: list (phi,psi) two torsion angles in degrees
1437   
1438    '''
1439    phi = getRestTorsion(XYZ[:5],Amat)
1440    psi = getRestTorsion(XYZ[1:],Amat)
1441    return phi,psi
1442   
1443def calcRamaEnergy(phi,psi,Coeff=[]):
1444    '''Computes pseudo potential energy from a pair of torsion angles and a
1445    numerical description of the potential energy surface. Used to create
1446    penalty function in LS refinement:     
1447    :math:`Eval(\\phi,\\psi) = C[0]*exp(-V/1000)`
1448
1449    where :math:`V = -C[3] * (\\phi-C[1])^2 - C[4]*(\\psi-C[2])^2 - 2*(\\phi-C[1])*(\\psi-C[2])`
1450   
1451    :param float phi: first torsion angle (:math:`\\phi`)
1452    :param float psi: second torsion angle (:math:`\\psi`)
1453    :param list Coeff: pseudo potential coefficients
1454   
1455    :returns: list (sum,Eval): pseudo-potential difference from minimum & value;
1456      sum is used for penalty function.
1457   
1458    '''
1459    sum = 0.
1460    Eval = 0.
1461    if len(Coeff):
1462        cof = Coeff.T
1463        dPhi = phi-cof[1]
1464        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
1465        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
1466        dPsi = psi-cof[2]
1467        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
1468        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
1469        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
1470        val = cof[0]*np.exp(val/1000.)
1471        pMax = cof[0][np.argmin(val)]
1472        Eval = np.sum(val)
1473        sum = Eval-pMax
1474    return sum,Eval
1475
1476def getRamaDeriv(XYZ,Amat,Coeff):
1477    '''Computes numerical derivatives of torsion angle pair pseudo potential
1478    with respect of crystallographic atom coordinates of the 5 atom sequence
1479   
1480    :param nparray XYZ: crystallographic coordinates of 5 atoms
1481    :param nparray Amat: crystal to cartesian transformation matrix
1482    :param list Coeff: pseudo potential coefficients
1483   
1484    :returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
1485     crystallographic xyz coordinates.
1486   
1487    '''
1488    deriv = np.zeros((len(XYZ),3))
1489    dx = 0.00001
1490    for j,xyz in enumerate(XYZ):
1491        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1492            XYZ[j] -= x
1493            phi,psi = getRestRama(XYZ,Amat)
1494            p1,d1 = calcRamaEnergy(phi,psi,Coeff)
1495            XYZ[j] += 2*x
1496            phi,psi = getRestRama(XYZ,Amat)
1497            p2,d2 = calcRamaEnergy(phi,psi,Coeff)
1498            XYZ[j] -= x
1499            deriv[j][i] = (p2-p1)/(2*dx)
1500    return deriv.flatten()
1501
1502def getRestPolefig(ODFln,SamSym,Grid):
1503    '''default doc string
1504   
1505    :param type name: description
1506   
1507    :returns: type name: description
1508   
1509    '''
1510    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
1511    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
1512    R = np.where(R <= 1.,2.*atand(R),0.0)
1513    Z = np.zeros_like(R)
1514    Z = G2lat.polfcal(ODFln,SamSym,R,P)
1515    Z = np.reshape(Z,(Grid,Grid))
1516    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
1517
1518def getRestPolefigDerv(HKL,Grid,SHCoeff):
1519    '''default doc string
1520   
1521    :param type name: description
1522   
1523    :returns: type name: description
1524   
1525    '''
1526    pass
1527       
1528def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
1529    '''default doc string
1530   
1531    :param type name: description
1532   
1533    :returns: type name: description
1534   
1535    '''
1536    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
1537        TxT = inv*(np.inner(M,Tx)+T)+C+U
1538        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
1539       
1540    inv = Top/abs(Top)
1541    cent = abs(Top)/100
1542    op = abs(Top)%100-1
1543    M,T = SGData['SGOps'][op]
1544    C = SGData['SGCen'][cent]
1545    dx = .00001
1546    deriv = np.zeros(6)
1547    for i in [0,1,2]:
1548        Oxyz[i] -= dx
1549        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1550        Oxyz[i] += 2*dx
1551        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1552        Oxyz[i] -= dx
1553        Txyz[i] -= dx
1554        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1555        Txyz[i] += 2*dx
1556        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1557        Txyz[i] -= dx
1558    return deriv
1559   
1560def getAngSig(VA,VB,Amat,SGData,covData={}):
1561    '''default doc string
1562   
1563    :param type name: description
1564   
1565    :returns: type name: description
1566   
1567    '''
1568    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
1569        TxT = inv*(np.inner(M,Tx)+T)+C+U
1570        return np.inner(Amat,(TxT-Ox))
1571       
1572    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
1573        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
1574        VecA /= np.sqrt(np.sum(VecA**2))
1575        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
1576        VecB /= np.sqrt(np.sum(VecB**2))
1577        edge = VecB-VecA
1578        edge = np.sum(edge**2)
1579        angle = (2.-edge)/2.
1580        angle = max(angle,-1.)
1581        return acosd(angle)
1582       
1583    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
1584    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
1585    invA = invB = 1
1586    invA = TopA/abs(TopA)
1587    invB = TopB/abs(TopB)
1588    centA = abs(TopA)/100
1589    centB = abs(TopB)/100
1590    opA = abs(TopA)%100-1
1591    opB = abs(TopB)%100-1
1592    MA,TA = SGData['SGOps'][opA]
1593    MB,TB = SGData['SGOps'][opB]
1594    CA = SGData['SGCen'][centA]
1595    CB = SGData['SGCen'][centB]
1596    if 'covMatrix' in covData:
1597        covMatrix = covData['covMatrix']
1598        varyList = covData['varyList']
1599        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
1600        dx = .00001
1601        dadx = np.zeros(9)
1602        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1603        for i in [0,1,2]:
1604            OxA[i] -= dx
1605            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1606            OxA[i] += 2*dx
1607            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1608            OxA[i] -= dx
1609           
1610            TxA[i] -= dx
1611            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1612            TxA[i] += 2*dx
1613            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1614            TxA[i] -= dx
1615           
1616            TxB[i] -= dx
1617            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1618            TxB[i] += 2*dx
1619            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1620            TxB[i] -= dx
1621           
1622        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1623        if sigAng < 0.01:
1624            sigAng = 0.0
1625        return Ang,sigAng
1626    else:
1627        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
1628
1629def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
1630    '''default doc string
1631   
1632    :param type name: description
1633   
1634    :returns: type name: description
1635   
1636    '''
1637    def calcDist(Atoms,SyOps,Amat):
1638        XYZ = []
1639        for i,atom in enumerate(Atoms):
1640            Inv,M,T,C,U = SyOps[i]
1641            XYZ.append(np.array(atom[1:4]))
1642            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1643            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1644        V1 = XYZ[1]-XYZ[0]
1645        return np.sqrt(np.sum(V1**2))
1646       
1647    Inv = []
1648    SyOps = []
1649    names = []
1650    for i,atom in enumerate(Oatoms):
1651        names += atom[-1]
1652        Op,unit = Atoms[i][-1]
1653        inv = Op/abs(Op)
1654        m,t = SGData['SGOps'][abs(Op)%100-1]
1655        c = SGData['SGCen'][abs(Op)/100]
1656        SyOps.append([inv,m,t,c,unit])
1657    Dist = calcDist(Oatoms,SyOps,Amat)
1658   
1659    sig = -0.001
1660    if 'covMatrix' in covData:
1661        parmNames = []
1662        dx = .00001
1663        dadx = np.zeros(6)
1664        for i in range(6):
1665            ia = i/3
1666            ix = i%3
1667            Oatoms[ia][ix+1] += dx
1668            a0 = calcDist(Oatoms,SyOps,Amat)
1669            Oatoms[ia][ix+1] -= 2*dx
1670            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1671        covMatrix = covData['covMatrix']
1672        varyList = covData['varyList']
1673        DistVcov = getVCov(names,varyList,covMatrix)
1674        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
1675        if sig < 0.001:
1676            sig = -0.001
1677   
1678    return Dist,sig
1679
1680def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
1681    '''default doc string
1682   
1683    :param type name: description
1684   
1685    :returns: type name: description
1686   
1687    '''
1688
1689    def calcAngle(Atoms,SyOps,Amat):
1690        XYZ = []
1691        for i,atom in enumerate(Atoms):
1692            Inv,M,T,C,U = SyOps[i]
1693            XYZ.append(np.array(atom[1:4]))
1694            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1695            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1696        V1 = XYZ[1]-XYZ[0]
1697        V1 /= np.sqrt(np.sum(V1**2))
1698        V2 = XYZ[1]-XYZ[2]
1699        V2 /= np.sqrt(np.sum(V2**2))
1700        V3 = V2-V1
1701        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1702        return acosd(cang)
1703
1704    Inv = []
1705    SyOps = []
1706    names = []
1707    for i,atom in enumerate(Oatoms):
1708        names += atom[-1]
1709        Op,unit = Atoms[i][-1]
1710        inv = Op/abs(Op)
1711        m,t = SGData['SGOps'][abs(Op)%100-1]
1712        c = SGData['SGCen'][abs(Op)/100]
1713        SyOps.append([inv,m,t,c,unit])
1714    Angle = calcAngle(Oatoms,SyOps,Amat)
1715   
1716    sig = -0.01
1717    if 'covMatrix' in covData:
1718        parmNames = []
1719        dx = .00001
1720        dadx = np.zeros(9)
1721        for i in range(9):
1722            ia = i/3
1723            ix = i%3
1724            Oatoms[ia][ix+1] += dx
1725            a0 = calcAngle(Oatoms,SyOps,Amat)
1726            Oatoms[ia][ix+1] -= 2*dx
1727            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1728        covMatrix = covData['covMatrix']
1729        varyList = covData['varyList']
1730        AngVcov = getVCov(names,varyList,covMatrix)
1731        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1732        if sig < 0.01:
1733            sig = -0.01
1734   
1735    return Angle,sig
1736
1737def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
1738    '''default doc string
1739   
1740    :param type name: description
1741   
1742    :returns: type name: description
1743   
1744    '''
1745
1746    def calcTorsion(Atoms,SyOps,Amat):
1747       
1748        XYZ = []
1749        for i,atom in enumerate(Atoms):
1750            Inv,M,T,C,U = SyOps[i]
1751            XYZ.append(np.array(atom[1:4]))
1752            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1753            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1754        V1 = XYZ[1]-XYZ[0]
1755        V2 = XYZ[2]-XYZ[1]
1756        V3 = XYZ[3]-XYZ[2]
1757        V1 /= np.sqrt(np.sum(V1**2))
1758        V2 /= np.sqrt(np.sum(V2**2))
1759        V3 /= np.sqrt(np.sum(V3**2))
1760        M = np.array([V1,V2,V3])
1761        D = nl.det(M)
1762        Ang = 1.0
1763        P12 = np.dot(V1,V2)
1764        P13 = np.dot(V1,V3)
1765        P23 = np.dot(V2,V3)
1766        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1767        return Tors
1768           
1769    Inv = []
1770    SyOps = []
1771    names = []
1772    for i,atom in enumerate(Oatoms):
1773        names += atom[-1]
1774        Op,unit = Atoms[i][-1]
1775        inv = Op/abs(Op)
1776        m,t = SGData['SGOps'][abs(Op)%100-1]
1777        c = SGData['SGCen'][abs(Op)/100]
1778        SyOps.append([inv,m,t,c,unit])
1779    Tors = calcTorsion(Oatoms,SyOps,Amat)
1780   
1781    sig = -0.01
1782    if 'covMatrix' in covData:
1783        parmNames = []
1784        dx = .00001
1785        dadx = np.zeros(12)
1786        for i in range(12):
1787            ia = i/3
1788            ix = i%3
1789            Oatoms[ia][ix+1] -= dx
1790            a0 = calcTorsion(Oatoms,SyOps,Amat)
1791            Oatoms[ia][ix+1] += 2*dx
1792            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1793            Oatoms[ia][ix+1] -= dx           
1794        covMatrix = covData['covMatrix']
1795        varyList = covData['varyList']
1796        TorVcov = getVCov(names,varyList,covMatrix)
1797        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
1798        if sig < 0.01:
1799            sig = -0.01
1800   
1801    return Tors,sig
1802       
1803def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
1804    '''default doc string
1805   
1806    :param type name: description
1807   
1808    :returns: type name: description
1809   
1810    '''
1811
1812    def calcDist(Atoms,SyOps,Amat):
1813        XYZ = []
1814        for i,atom in enumerate(Atoms):
1815            Inv,M,T,C,U = SyOps[i]
1816            XYZ.append(np.array(atom[1:4]))
1817            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1818            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1819        V1 = XYZ[1]-XYZ[0]
1820        return np.sqrt(np.sum(V1**2))
1821       
1822    def calcAngle(Atoms,SyOps,Amat):
1823        XYZ = []
1824        for i,atom in enumerate(Atoms):
1825            Inv,M,T,C,U = SyOps[i]
1826            XYZ.append(np.array(atom[1:4]))
1827            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1828            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1829        V1 = XYZ[1]-XYZ[0]
1830        V1 /= np.sqrt(np.sum(V1**2))
1831        V2 = XYZ[1]-XYZ[2]
1832        V2 /= np.sqrt(np.sum(V2**2))
1833        V3 = V2-V1
1834        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1835        return acosd(cang)
1836
1837    def calcTorsion(Atoms,SyOps,Amat):
1838       
1839        XYZ = []
1840        for i,atom in enumerate(Atoms):
1841            Inv,M,T,C,U = SyOps[i]
1842            XYZ.append(np.array(atom[1:4]))
1843            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1844            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1845        V1 = XYZ[1]-XYZ[0]
1846        V2 = XYZ[2]-XYZ[1]
1847        V3 = XYZ[3]-XYZ[2]
1848        V1 /= np.sqrt(np.sum(V1**2))
1849        V2 /= np.sqrt(np.sum(V2**2))
1850        V3 /= np.sqrt(np.sum(V3**2))
1851        M = np.array([V1,V2,V3])
1852        D = nl.det(M)
1853        Ang = 1.0
1854        P12 = np.dot(V1,V2)
1855        P13 = np.dot(V1,V3)
1856        P23 = np.dot(V2,V3)
1857        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1858        return Tors
1859           
1860    Inv = []
1861    SyOps = []
1862    names = []
1863    for i,atom in enumerate(Oatoms):
1864        names += atom[-1]
1865        Op,unit = Atoms[i][-1]
1866        inv = Op/abs(Op)
1867        m,t = SGData['SGOps'][abs(Op)%100-1]
1868        c = SGData['SGCen'][abs(Op)/100]
1869        SyOps.append([inv,m,t,c,unit])
1870    M = len(Oatoms)
1871    if M == 2:
1872        Val = calcDist(Oatoms,SyOps,Amat)
1873    elif M == 3:
1874        Val = calcAngle(Oatoms,SyOps,Amat)
1875    else:
1876        Val = calcTorsion(Oatoms,SyOps,Amat)
1877   
1878    sigVals = [-0.001,-0.01,-0.01]
1879    sig = sigVals[M-3]
1880    if 'covMatrix' in covData:
1881        parmNames = []
1882        dx = .00001
1883        N = M*3
1884        dadx = np.zeros(N)
1885        for i in range(N):
1886            ia = i/3
1887            ix = i%3
1888            Oatoms[ia][ix+1] += dx
1889            if M == 2:
1890                a0 = calcDist(Oatoms,SyOps,Amat)
1891            elif M == 3:
1892                a0 = calcAngle(Oatoms,SyOps,Amat)
1893            else:
1894                a0 = calcTorsion(Oatoms,SyOps,Amat)
1895            Oatoms[ia][ix+1] -= 2*dx
1896            if M == 2:
1897                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1898            elif M == 3:
1899                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1900            else:
1901                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1902        covMatrix = covData['covMatrix']
1903        varyList = covData['varyList']
1904        Vcov = getVCov(names,varyList,covMatrix)
1905        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
1906        if sig < sigVals[M-3]:
1907            sig = sigVals[M-3]
1908   
1909    return Val,sig
1910       
1911def ValEsd(value,esd=0,nTZ=False):
1912    '''Format a floating point number with a given level of precision or
1913    with in crystallographic format with a "esd", as value(esd). If esd is
1914    negative the number is formatted with the level of significant figures
1915    appropriate if abs(esd) were the esd, but the esd is not included.
1916    if the esd is zero, approximately 6 significant figures are printed.
1917    nTZ=True causes "extra" zeros to be removed after the decimal place.
1918    for example:
1919
1920      * "1.235(3)" for value=1.2346 & esd=0.003
1921      * "1.235(3)e4" for value=12346. & esd=30
1922      * "1.235(3)e6" for value=0.12346e7 & esd=3000
1923      * "1.235" for value=1.2346 & esd=-0.003
1924      * "1.240" for value=1.2395 & esd=-0.003
1925      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
1926      * "1.23460" for value=1.2346 & esd=0.0
1927
1928    :param float value: number to be formatted
1929    :param float esd: uncertainty or if esd < 0, specifies level of
1930      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
1931    :param bool nTZ: True to remove trailing zeros (default is False)
1932    :returns: value(esd) or value as a string
1933
1934    '''
1935    # Note: this routine is Python 3 compatible -- I think
1936    cutoff = 3.16228    #=(sqrt(10); same as old GSAS   was 1.95
1937    if math.isnan(value): # invalid value, bail out
1938        return '?'
1939    if math.isnan(esd): # invalid esd, treat as zero
1940        esd = 0
1941        esdoff = 5
1942    elif esd != 0:
1943        # transform the esd to a one or two digit integer
1944        l = math.log10(abs(esd)) % 1.
1945        if l < math.log10(cutoff): l+= 1.       
1946        intesd = int(round(10**l)) # esd as integer
1947        # determine the number of digits offset for the esd
1948        esdoff = int(round(math.log10(intesd*1./abs(esd))))
1949    else:
1950        esdoff = 5
1951    valoff = 0
1952    if abs(value) < abs(esdoff): # value is effectively zero
1953        pass
1954    elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
1955        # where the digit offset is to the left of the decimal place or where too many
1956        # digits are needed
1957        if abs(value) > 1:
1958            valoff = int(math.log10(abs(value)))
1959        elif abs(value) > 0:
1960            valoff = int(math.log10(abs(value))-0.9999999)
1961        else:
1962            valoff = 0
1963    if esd != 0:
1964        if valoff+esdoff < 0:
1965            valoff = esdoff = 0
1966        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
1967    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
1968        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
1969    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
1970        if abs(value) > 0:           
1971            extra = -math.log10(abs(value))
1972        else:
1973            extra = 0
1974        if extra > 0: extra += 1
1975        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
1976    if esd > 0:
1977        out += ("({:d})").format(intesd)  # add the esd
1978    elif nTZ and '.' in out:
1979        out = out.rstrip('0')  # strip zeros to right of decimal
1980        out = out.rstrip('.')  # and decimal place when not needed
1981    if valoff != 0:
1982        out += ("e{:d}").format(valoff) # add an exponent, when needed
1983    return out
1984   
1985################################################################################
1986##### Texture fitting stuff
1987################################################################################
1988
1989def FitTexture(General,Gangls,refData,keyList,pgbar):
1990    import pytexture as ptx
1991    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
1992   
1993    def printSpHarm(textureData,SHtextureSig):
1994        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1995        names = ['omega','chi','phi']
1996        namstr = '  names :'
1997        ptstr =  '  values:'
1998        sigstr = '  esds  :'
1999        for name in names:
2000            namstr += '%12s'%('Sample '+name)
2001            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2002            if 'Sample '+name in SHtextureSig:
2003                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2004            else:
2005                sigstr += 12*' '
2006        print namstr
2007        print ptstr
2008        print sigstr
2009        print '\n Texture coefficients:'
2010        SHcoeff = textureData['SH Coeff'][1]
2011        SHkeys = SHcoeff.keys()
2012        nCoeff = len(SHcoeff)
2013        nBlock = nCoeff/10+1
2014        iBeg = 0
2015        iFin = min(iBeg+10,nCoeff)
2016        for block in range(nBlock):
2017            namstr = '  names :'
2018            ptstr =  '  values:'
2019            sigstr = '  esds  :'
2020            for name in SHkeys[iBeg:iFin]:
2021                if 'C' in name:
2022                    namstr += '%12s'%(name)
2023                    ptstr += '%12.3f'%(SHcoeff[name])
2024                    if name in SHtextureSig:
2025                        sigstr += '%12.3f'%(SHtextureSig[name])
2026                    else:
2027                        sigstr += 12*' '
2028            print namstr
2029            print ptstr
2030            print sigstr
2031            iBeg += 10
2032            iFin = min(iBeg+10,nCoeff)
2033           
2034    def Dict2Values(parmdict, varylist):
2035        '''Use before call to leastsq to setup list of values for the parameters
2036        in parmdict, as selected by key in varylist'''
2037        return [parmdict[key] for key in varylist] 
2038       
2039    def Values2Dict(parmdict, varylist, values):
2040        ''' Use after call to leastsq to update the parameter dictionary with
2041        values corresponding to keys in varylist'''
2042        parmdict.update(zip(varylist,values))
2043       
2044    def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2045        parmDict.update(zip(varyList,values))
2046        Mat = np.empty(0)
2047        sumObs = 0
2048        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
2049        for hist in Gangls.keys():
2050            Refs = refData[hist]
2051            Refs[:,5] = np.where(Refs[:,5]>0.,Refs[:,5],0.)
2052            wt = 1./np.sqrt(np.max(Refs[:,4],.25))
2053#            wt = 1./np.max(Refs[:,4],.25)
2054            sumObs += np.sum(wt*Refs[:,5])
2055            Refs[:,6] = 1.
2056            H = Refs[:,:3]
2057            phi,beta = G2lat.CrsAng(H,cell,SGData)
2058            psi,gam,x,x = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2059            for item in parmDict:
2060                if 'C' in item:
2061                    L,M,N = eval(item.strip('C'))
2062                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2063                    Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
2064                    Lnorm = G2lat.Lnorm(L)
2065                    Refs[:,6] += parmDict[item]*Lnorm*Kcl*Ksl
2066            mat = wt*(Refs[:,5]-Refs[:,6])
2067            Mat = np.concatenate((Mat,mat))
2068        sumD = np.sum(np.abs(Mat))
2069        R = min(100.,100.*sumD/sumObs)
2070        pgbar.Update(R,newmsg='Residual = %5.2f'%(R))
2071        print ' Residual: %.3f%%'%(R)
2072        return Mat
2073       
2074    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2075        Mat = np.empty(0)
2076        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
2077        for hist in Gangls.keys():
2078            mat = np.zeros((len(varyList),len(refData[hist])))
2079            Refs = refData[hist]
2080            H = Refs[:,:3]
2081            wt = 1./np.sqrt(np.max(Refs[:,4],.25))
2082#            wt = 1./np.max(Refs[:,4],.25)
2083            phi,beta = G2lat.CrsAng(H,cell,SGData)
2084            psi,gam,dPdA,dGdA = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2085            for j,item in enumerate(varyList):
2086                if 'C' in item:
2087                    L,M,N = eval(item.strip('C'))
2088                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2089                    Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
2090                    Lnorm = G2lat.Lnorm(L)
2091                    mat[j] = -wt*Lnorm*Kcl*Ksl
2092                    for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
2093                        try:
2094                            l = varyList.index(itema)
2095                            mat[l] -= parmDict[item]*wt*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
2096                        except ValueError:
2097                            pass
2098            if len(Mat):
2099                Mat = np.concatenate((Mat,mat.T))
2100            else:
2101                Mat = mat.T
2102        print 'deriv'
2103        return Mat
2104
2105    print ' Fit texture for '+General['Name']
2106    SGData = General['SGData']
2107    cell = General['Cell'][1:7]
2108    Texture = General['SH Texture']
2109    if not Texture['Order']:
2110        return 'No spherical harmonics coefficients'
2111    varyList = []
2112    parmDict = copy.copy(Texture['SH Coeff'][1])
2113    for item in ['Sample omega','Sample chi','Sample phi']:
2114        parmDict[item] = Texture[item][1]
2115        if Texture[item][0]:
2116            varyList.append(item)
2117    if Texture['SH Coeff'][0]:
2118        varyList += Texture['SH Coeff'][1].keys()
2119    while True:
2120        begin = time.time()
2121        values =  np.array(Dict2Values(parmDict, varyList))
2122        result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,ftol=1.e-6,
2123            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList,pgbar))
2124        ncyc = int(result[2]['nfev']/2)
2125        if ncyc:
2126            runtime = time.time()-begin   
2127            chisq = np.sum(result[2]['fvec']**2)
2128            Values2Dict(parmDict, varyList, result[0])
2129            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
2130            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(result[2]['fvec']),' Number of parameters: ',len(varyList)
2131            print 'refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2132            try:
2133                sig = np.sqrt(np.diag(result[1])*GOF)
2134                if np.any(np.isnan(sig)):
2135                    print '*** Least squares aborted - some invalid esds possible ***'
2136                break                   #refinement succeeded - finish up!
2137            except ValueError:          #result[1] is None on singular matrix
2138                print '**** Refinement failed - singular matrix ****'
2139                return None
2140        else:
2141            break
2142   
2143    if ncyc:
2144        for parm in parmDict:
2145            if 'C' in parm:
2146                Texture['SH Coeff'][1][parm] = parmDict[parm]
2147            else:
2148                Texture[parm][1] = parmDict[parm] 
2149        sigDict = dict(zip(varyList,sig))
2150        printSpHarm(Texture,sigDict)
2151       
2152    return None
2153   
2154################################################################################
2155##### Fourier & charge flip stuff
2156################################################################################
2157
2158def adjHKLmax(SGData,Hmax):
2159    '''default doc string
2160   
2161    :param type name: description
2162   
2163    :returns: type name: description
2164   
2165    '''
2166    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
2167        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
2168        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
2169        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
2170    else:
2171        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
2172        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
2173        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
2174
2175def OmitMap(data,reflDict,pgbar=None):
2176    '''default doc string
2177   
2178    :param type name: description
2179   
2180    :returns: type name: description
2181   
2182    '''
2183    generalData = data['General']
2184    if not generalData['Map']['MapType']:
2185        print '**** ERROR - Fourier map not defined'
2186        return
2187    mapData = generalData['Map']
2188    dmin = mapData['Resolution']
2189    SGData = generalData['SGData']
2190    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2191    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2192    cell = generalData['Cell'][1:8]       
2193    A = G2lat.cell2A(cell[:6])
2194    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2195    adjHKLmax(SGData,Hmax)
2196    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2197    time0 = time.time()
2198    for iref,ref in enumerate(reflDict['RefList']):
2199        if ref[4] >= dmin:
2200            Fosq,Fcsq,ph = ref[8:11]
2201            Uniq = np.inner(ref[:3],SGMT)
2202            Phi = np.inner(ref[:3],SGT)
2203            for i,hkl in enumerate(Uniq):        #uses uniq
2204                hkl = np.asarray(hkl,dtype='i')
2205                dp = 360.*Phi[i]                #and phi
2206                a = cosd(ph+dp)
2207                b = sind(ph+dp)
2208                phasep = complex(a,b)
2209                phasem = complex(a,-b)
2210                Fo = np.sqrt(Fosq)
2211                if '2Fo-Fc' in mapData['MapType']:
2212                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
2213                else:
2214                    F = np.sqrt(Fosq)
2215                h,k,l = hkl+Hmax
2216                Fhkl[h,k,l] = F*phasep
2217                h,k,l = -hkl+Hmax
2218                Fhkl[h,k,l] = F*phasem
2219    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
2220    M = np.mgrid[0:4,0:4,0:4]
2221    blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
2222    iBeg = blkIds*rho0.shape/4
2223    iFin = (blkIds+1)*rho0.shape/4
2224    rho_omit = np.zeros_like(rho0)
2225    nBlk = 0
2226    for iB,iF in zip(iBeg,iFin):
2227        rho1 = np.copy(rho0)
2228        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
2229        Fnew = fft.ifftshift(fft.ifftn(rho1))
2230        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
2231        phase = Fnew/np.absolute(Fnew)
2232        OFhkl = np.absolute(Fhkl)*phase
2233        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
2234        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]])
2235        nBlk += 1
2236        pgbar.Update(nBlk)
2237    mapData['rho'] = np.real(rho_omit)/cell[6]
2238    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2239    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2240    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
2241    return mapData
2242   
2243def Fourier4DMap(data,reflDict):
2244    '''default doc string
2245   
2246    :param type name: description
2247   
2248    :returns: type name: description
2249   
2250    '''
2251    generalData = data['General']
2252    map4DData = generalData['4DmapData']
2253    mapData = generalData['Map']
2254    dmin = mapData['Resolution']
2255    SGData = generalData['SGData']
2256    SSGData = generalData['SSGData']
2257    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2258    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2259    cell = generalData['Cell'][1:8]       
2260    A = G2lat.cell2A(cell[:6])
2261    maxM = generalData['SuperVec'][2]
2262    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
2263    adjHKLmax(SGData,Hmax)
2264    Hmax = np.asarray(Hmax,dtype='i')+1
2265    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2266    time0 = time.time()
2267    for iref,ref in enumerate(reflDict['RefList']):
2268        if ref[5] > dmin:
2269            Fosq,Fcsq,ph = ref[9:12]
2270            Fosq = np.where(Fosq>0.,Fosq,0.)    #can't use Fo^2 < 0
2271            Uniq = np.inner(ref[:4],SSGMT)
2272            Phi = np.inner(ref[:4],SSGT)
2273            for i,hkl in enumerate(Uniq):        #uses uniq
2274                hkl = np.asarray(hkl,dtype='i')
2275                dp = 360.*Phi[i]                #and phi
2276                a = cosd(ph+dp)
2277                b = sind(ph+dp)
2278                phasep = complex(a,b)
2279                phasem = complex(a,-b)
2280                if 'Fobs' in mapData['MapType']:
2281                    F = np.sqrt(Fosq)
2282                    h,k,l,m = hkl+Hmax
2283                    Fhkl[h,k,l,m] = F*phasep
2284                    h,k,l,m = -hkl+Hmax
2285                    Fhkl[h,k,l,m] = F*phasem
2286                elif 'Fcalc' in mapData['MapType']:
2287                    F = np.sqrt(Fcsq)
2288                    h,k,l,m = hkl+Hmax
2289                    Fhkl[h,k,l,m] = F*phasep
2290                    h,k,l,m = -hkl+Hmax
2291                    Fhkl[h,k,l,m] = F*phasem                   
2292                elif 'delt-F' in mapData['MapType']:
2293                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
2294                    h,k,l,m = hkl+Hmax
2295                    Fhkl[h,k,l,m] = dF*phasep
2296                    h,k,l,m = -hkl+Hmax
2297                    Fhkl[h,k,l,m] = dF*phasem
2298    SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6]                #4D map
2299    rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6]    #3D map
2300    map4DData['rho'] = np.real(SSrho)
2301    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
2302    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
2303    map4DData['Type'] = reflDict['Type']
2304    mapData['Type'] = reflDict['Type']
2305    mapData['rho'] = np.real(rho)
2306    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2307    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2308    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
2309
2310def FourierMap(data,reflDict):
2311    '''default doc string
2312   
2313    :param type name: description
2314   
2315    :returns: type name: description
2316   
2317    '''
2318    generalData = data['General']
2319    mapData = generalData['Map']
2320    dmin = mapData['Resolution']
2321    SGData = generalData['SGData']
2322    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2323    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2324    cell = generalData['Cell'][1:8]       
2325    A = G2lat.cell2A(cell[:6])
2326    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2327    adjHKLmax(SGData,Hmax)
2328    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2329#    Fhkl[0,0,0] = generalData['F000X']
2330    time0 = time.time()
2331    for iref,ref in enumerate(reflDict['RefList']):
2332        if ref[4] > dmin:
2333            Fosq,Fcsq,ph = ref[8:11]
2334            Uniq = np.inner(ref[:3],SGMT)
2335            Phi = np.inner(ref[:3],SGT)
2336            for i,hkl in enumerate(Uniq):        #uses uniq
2337                hkl = np.asarray(hkl,dtype='i')
2338                dp = 360.*Phi[i]                #and phi
2339                a = cosd(ph+dp)
2340                b = sind(ph+dp)
2341                phasep = complex(a,b)
2342                phasem = complex(a,-b)
2343                if 'Fobs' in mapData['MapType']:
2344                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
2345                    h,k,l = hkl+Hmax
2346                    Fhkl[h,k,l] = F*phasep
2347                    h,k,l = -hkl+Hmax
2348                    Fhkl[h,k,l] = F*phasem
2349                elif 'Fcalc' in mapData['MapType']:
2350                    F = np.sqrt(Fcsq)
2351                    h,k,l = hkl+Hmax
2352                    Fhkl[h,k,l] = F*phasep
2353                    h,k,l = -hkl+Hmax
2354                    Fhkl[h,k,l] = F*phasem
2355                elif 'delt-F' in mapData['MapType']:
2356                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
2357                    h,k,l = hkl+Hmax
2358                    Fhkl[h,k,l] = dF*phasep
2359                    h,k,l = -hkl+Hmax
2360                    Fhkl[h,k,l] = dF*phasem
2361                elif '2*Fo-Fc' in mapData['MapType']:
2362                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
2363                    h,k,l = hkl+Hmax
2364                    Fhkl[h,k,l] = F*phasep
2365                    h,k,l = -hkl+Hmax
2366                    Fhkl[h,k,l] = F*phasem
2367                elif 'Patterson' in mapData['MapType']:
2368                    h,k,l = hkl+Hmax
2369                    Fhkl[h,k,l] = complex(Fosq,0.)
2370                    h,k,l = -hkl+Hmax
2371                    Fhkl[h,k,l] = complex(Fosq,0.)
2372    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
2373    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
2374    mapData['Type'] = reflDict['Type']
2375    mapData['rho'] = np.real(rho)
2376    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2377    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2378   
2379# map printing for testing purposes
2380def printRho(SGLaue,rho,rhoMax):                         
2381    '''default doc string
2382   
2383    :param type name: description
2384   
2385    :returns: type name: description
2386   
2387    '''
2388    dim = len(rho.shape)
2389    if dim == 2:
2390        ix,jy = rho.shape
2391        for j in range(jy):
2392            line = ''
2393            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
2394                line += (jy-j)*'  '
2395            for i in range(ix):
2396                r = int(100*rho[i,j]/rhoMax)
2397                line += '%4d'%(r)
2398            print line+'\n'
2399    else:
2400        ix,jy,kz = rho.shape
2401        for k in range(kz):
2402            print 'k = ',k
2403            for j in range(jy):
2404                line = ''
2405                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
2406                    line += (jy-j)*'  '
2407                for i in range(ix):
2408                    r = int(100*rho[i,j,k]/rhoMax)
2409                    line += '%4d'%(r)
2410                print line+'\n'
2411## keep this
2412               
2413def findOffset(SGData,A,Fhkl):   
2414    '''default doc string
2415   
2416    :param type name: description
2417   
2418    :returns: type name: description
2419   
2420    '''
2421    if SGData['SpGrp'] == 'P 1':
2422        return [0,0,0]   
2423    hklShape = Fhkl.shape
2424    hklHalf = np.array(hklShape)/2
2425    sortHKL = np.argsort(Fhkl.flatten())
2426    Fdict = {}
2427    for hkl in sortHKL:
2428        HKL = np.unravel_index(hkl,hklShape)
2429        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
2430        if F == 0.:
2431            break
2432        Fdict['%.6f'%(np.absolute(F))] = hkl
2433    Flist = np.flipud(np.sort(Fdict.keys()))
2434    F = str(1.e6)
2435    i = 0
2436    DH = []
2437    Dphi = []
2438    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2439    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2440    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
2441    for F in Flist:
2442        hkl = np.unravel_index(Fdict[F],hklShape)
2443        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
2444            continue
2445        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
2446        Uniq = np.array(Uniq,dtype='i')
2447        Phi = np.array(Phi)
2448        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
2449        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
2450        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
2451        ang0 = np.angle(Fh0,deg=True)/360.
2452        for H,phi in zip(Uniq,Phi)[1:]:
2453            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
2454            dH = H-hkl
2455            dang = ang-ang0
2456            DH.append(dH)
2457            Dphi.append((dang+.5) % 1.0)
2458        if i > 20 or len(DH) > 30:
2459            break
2460        i += 1
2461    DH = np.array(DH)
2462    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
2463    Dphi = np.array(Dphi)
2464    steps = np.array(hklShape)
2465    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
2466    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
2467    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
2468    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
2469    hist,bins = np.histogram(Mmap,bins=1000)
2470#    for i,item in enumerate(hist[:10]):
2471#        print item,bins[i]
2472    chisq = np.min(Mmap)
2473    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
2474    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
2475#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
2476    return DX
2477   
2478def ChargeFlip(data,reflDict,pgbar):
2479    '''default doc string
2480   
2481    :param type name: description
2482   
2483    :returns: type name: description
2484   
2485    '''
2486    generalData = data['General']
2487    mapData = generalData['Map']
2488    flipData = generalData['Flip']
2489    FFtable = {}
2490    if 'None' not in flipData['Norm element']:
2491        normElem = flipData['Norm element'].upper()
2492        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
2493        for ff in FFs:
2494            if ff['Symbol'] == normElem:
2495                FFtable.update(ff)
2496    dmin = flipData['Resolution']
2497    SGData = generalData['SGData']
2498    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2499    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2500    cell = generalData['Cell'][1:8]       
2501    A = G2lat.cell2A(cell[:6])
2502    Vol = cell[6]
2503    im = 0
2504    if generalData['Type'] in ['modulated','magnetic',]:
2505        im = 1
2506    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2507    adjHKLmax(SGData,Hmax)
2508    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
2509    time0 = time.time()
2510    for iref,ref in enumerate(reflDict['RefList']):
2511        dsp = ref[4+im]
2512        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
2513            continue
2514        if dsp > dmin:
2515            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
2516            if FFtable:
2517                SQ = 0.25/dsp**2
2518                ff *= G2el.ScatFac(FFtable,SQ)[0]
2519            if ref[8+im] > 0.:         #use only +ve Fobs**2
2520                E = np.sqrt(ref[8+im])/ff
2521            else:
2522                E = 0.
2523            ph = ref[10]
2524            ph = rn.uniform(0.,360.)
2525            Uniq = np.inner(ref[:3],SGMT)
2526            Phi = np.inner(ref[:3],SGT)
2527            for i,hkl in enumerate(Uniq):        #uses uniq
2528                hkl = np.asarray(hkl,dtype='i')
2529                dp = 360.*Phi[i]                #and phi
2530                a = cosd(ph+dp)
2531                b = sind(ph+dp)
2532                phasep = complex(a,b)
2533                phasem = complex(a,-b)
2534                h,k,l = hkl+Hmax
2535                Ehkl[h,k,l] = E*phasep
2536                h,k,l = -hkl+Hmax       #Friedel pair refl.
2537                Ehkl[h,k,l] = E*phasem
2538#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
2539    CEhkl = copy.copy(Ehkl)
2540    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
2541    Emask = ma.getmask(MEhkl)
2542    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
2543    Ncyc = 0
2544    old = np.seterr(all='raise')
2545    while True:       
2546        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
2547        CEsig = np.std(CErho)
2548        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
2549        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
2550        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
2551        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
2552        phase = CFhkl/np.absolute(CFhkl)
2553        CEhkl = np.absolute(Ehkl)*phase
2554        Ncyc += 1
2555        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
2556        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
2557        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
2558        if Rcf < 5.:
2559            break
2560        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
2561        if not GoOn or Ncyc > 10000:
2562            break
2563    np.seterr(**old)
2564    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
2565    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
2566    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
2567    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
2568       
2569    mapData['Rcf'] = Rcf
2570    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2571    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2572    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2573    mapData['Type'] = reflDict['Type']
2574    return mapData
2575   
2576def findSSOffset(SGData,SSGData,A,Fhklm):   
2577    '''default doc string
2578   
2579    :param type name: description
2580   
2581    :returns: type name: description
2582   
2583    '''
2584    if SGData['SpGrp'] == 'P 1':
2585        return [0,0,0,0]   
2586    hklmShape = Fhklm.shape
2587    hklmHalf = np.array(hklmShape)/2
2588    sortHKLM = np.argsort(Fhklm.flatten())
2589    Fdict = {}
2590    for hklm in sortHKLM:
2591        HKLM = np.unravel_index(hklm,hklmShape)
2592        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
2593        if F == 0.:
2594            break
2595        Fdict['%.6f'%(np.absolute(F))] = hklm
2596    Flist = np.flipud(np.sort(Fdict.keys()))
2597    F = str(1.e6)
2598    i = 0
2599    DH = []
2600    Dphi = []
2601    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2602    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2603    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
2604    for F in Flist:
2605        hklm = np.unravel_index(Fdict[F],hklmShape)
2606        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
2607            continue
2608        Uniq = np.inner(hklm-hklmHalf,SSGMT)
2609        Phi = np.inner(hklm-hklmHalf,SSGT)
2610        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
2611        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
2612        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
2613        ang0 = np.angle(Fh0,deg=True)/360.
2614        for H,phi in zip(Uniq,Phi)[1:]:
2615            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
2616            dH = H-hklm
2617            dang = ang-ang0
2618            DH.append(dH)
2619            Dphi.append((dang+.5) % 1.0)
2620        if i > 20 or len(DH) > 30:
2621            break
2622        i += 1
2623    DH = np.array(DH)
2624    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
2625    Dphi = np.array(Dphi)
2626    steps = np.array(hklmShape)
2627    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]]
2628    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
2629    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
2630    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
2631    hist,bins = np.histogram(Mmap,bins=1000)
2632#    for i,item in enumerate(hist[:10]):
2633#        print item,bins[i]
2634    chisq = np.min(Mmap)
2635    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
2636    print ' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
2637#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
2638    return DX
2639   
2640def SSChargeFlip(data,reflDict,pgbar):
2641    '''default doc string
2642   
2643    :param type name: description
2644   
2645    :returns: type name: description
2646   
2647    '''
2648    generalData = data['General']
2649    mapData = generalData['Map']
2650    map4DData = {}
2651    flipData = generalData['Flip']
2652    FFtable = {}
2653    if 'None' not in flipData['Norm element']:
2654        normElem = flipData['Norm element'].upper()
2655        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
2656        for ff in FFs:
2657            if ff['Symbol'] == normElem:
2658                FFtable.update(ff)
2659    dmin = flipData['Resolution']
2660    SGData = generalData['SGData']
2661    SSGData = generalData['SSGData']
2662    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2663    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2664    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2665    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2666    cell = generalData['Cell'][1:8]       
2667    A = G2lat.cell2A(cell[:6])
2668    Vol = cell[6]
2669    maxM = generalData['SuperVec'][2]
2670    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
2671    adjHKLmax(SGData,Hmax)
2672    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
2673    time0 = time.time()
2674    for iref,ref in enumerate(reflDict['RefList']):
2675        dsp = ref[5]
2676        if dsp > dmin:
2677            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
2678            if FFtable:
2679                SQ = 0.25/dsp**2
2680                ff *= G2el.ScatFac(FFtable,SQ)[0]
2681            if ref[9] > 0.:         #use only +ve Fobs**2
2682                E = np.sqrt(ref[9])/ff
2683            else:
2684                E = 0.
2685            ph = ref[11]
2686            ph = rn.uniform(0.,360.)
2687            Uniq = np.inner(ref[:4],SSGMT)
2688            Phi = np.inner(ref[:4],SSGT)
2689            for i,hklm in enumerate(Uniq):        #uses uniq
2690                hklm = np.asarray(hklm,dtype='i')
2691                dp = 360.*Phi[i]                #and phi
2692                a = cosd(ph+dp)
2693                b = sind(ph+dp)
2694                phasep = complex(a,b)
2695                phasem = complex(a,-b)
2696                h,k,l,m = hklm+Hmax
2697                Ehkl[h,k,l,m] = E*phasep
2698                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
2699                Ehkl[h,k,l,m] = E*phasem
2700#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
2701    CEhkl = copy.copy(Ehkl)
2702    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
2703    Emask = ma.getmask(MEhkl)
2704    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
2705    Ncyc = 0
2706    old = np.seterr(all='raise')
2707    while True:       
2708        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
2709        CEsig = np.std(CErho)
2710        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
2711        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
2712        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
2713        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
2714        phase = CFhkl/np.absolute(CFhkl)
2715        CEhkl = np.absolute(Ehkl)*phase
2716        Ncyc += 1
2717        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
2718        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
2719        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
2720        if Rcf < 5.:
2721            break
2722        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
2723        if not GoOn or Ncyc > 10000:
2724            break
2725    np.seterr(**old)
2726    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
2727    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))
2728    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))
2729    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
2730    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
2731
2732    mapData['Rcf'] = Rcf
2733    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2734    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2735    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2736    mapData['Type'] = reflDict['Type']
2737
2738    map4DData['Rcf'] = Rcf
2739    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))
2740    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
2741    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
2742    map4DData['Type'] = reflDict['Type']
2743    return mapData,map4DData
2744   
2745def getRho(xyz,mapData):
2746    ''' get scattering density at a point by 8-point interpolation
2747    param xyz: coordinate to be probed
2748    param: mapData: dict of map data
2749   
2750    :returns: density at xyz
2751    '''
2752    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2753    if not len(mapData):
2754        return 0.0
2755    rho = copy.copy(mapData['rho'])     #don't mess up original
2756    if not len(rho):
2757        return 0.0
2758    mapShape = np.array(rho.shape)
2759    mapStep = 1./mapShape
2760    X = np.array(xyz)%1.    #get into unit cell
2761    I = np.array(X*mapShape,dtype='int')
2762    D = X-I*mapStep         #position inside map cell
2763    D12 = D[0]*D[1]
2764    D13 = D[0]*D[2]
2765    D23 = D[1]*D[2]
2766    D123 = np.prod(D)
2767    Rho = rollMap(rho,-I)    #shifts map so point is in corner
2768    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]+  \
2769        Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \
2770        Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \
2771        Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123)   
2772    return R
2773       
2774def SearchMap(generalData,drawingData,Neg=False):
2775    '''Does a search of a density map for peaks meeting the criterion of peak
2776    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
2777    mapData is data['General']['mapData']; the map is also in mapData.
2778
2779    :param generalData: the phase data structure; includes the map
2780    :param drawingData: the drawing data structure
2781    :param Neg:  if True then search for negative peaks (i.e. H-atoms & neutron data)
2782
2783    :returns: (peaks,mags,dzeros) where
2784
2785        * peaks : ndarray
2786          x,y,z positions of the peaks found in the map
2787        * mags : ndarray
2788          the magnitudes of the peaks
2789        * dzeros : ndarray
2790          the distance of the peaks from  the unit cell origin
2791        * dcent : ndarray
2792          the distance of the peaks from  the unit cell center
2793
2794    '''       
2795    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2796   
2797    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
2798   
2799#    def noDuplicate(xyz,peaks,Amat):
2800#        XYZ = np.inner(Amat,xyz)
2801#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2802#            print ' Peak',xyz,' <0.5A from another peak'
2803#            return False
2804#        return True
2805#                           
2806    def fixSpecialPos(xyz,SGData,Amat):
2807        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
2808        X = []
2809        xyzs = [equiv[0] for equiv in equivs]
2810        for x in xyzs:
2811            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
2812                X.append(x)
2813        if len(X) > 1:
2814            return np.average(X,axis=0)
2815        else:
2816            return xyz
2817       
2818    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
2819        Mag,x0,y0,z0,sig = parms
2820        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
2821#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
2822        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
2823       
2824    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
2825        Mag,x0,y0,z0,sig = parms
2826        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2827        return M
2828       
2829    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
2830        Mag,x0,y0,z0,sig = parms
2831        dMdv = np.zeros(([5,]+list(rX.shape)))
2832        delt = .01
2833        for i in range(5):
2834            parms[i] -= delt
2835            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2836            parms[i] += 2.*delt
2837            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2838            parms[i] -= delt
2839            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
2840        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2841        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
2842        dMdv = np.reshape(dMdv,(5,rX.size))
2843        Hess = np.inner(dMdv,dMdv)
2844       
2845        return Vec,Hess
2846       
2847    phaseName = generalData['Name']
2848    SGData = generalData['SGData']
2849    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2850    peaks = []
2851    mags = []
2852    dzeros = []
2853    dcent = []
2854    try:
2855        mapData = generalData['Map']
2856        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
2857        if Neg:
2858            rho = -copy.copy(mapData['rho'])     #flip +/-
2859        else:
2860            rho = copy.copy(mapData['rho'])     #don't mess up original
2861        mapHalf = np.array(rho.shape)/2
2862        res = mapData['Resolution']
2863        incre = np.array(rho.shape,dtype=np.float)
2864        step = max(1.0,1./res)+1
2865        steps = np.array(3*[step,])
2866    except KeyError:
2867        print '**** ERROR - Fourier map not defined'
2868        return peaks,mags
2869    rhoMask = ma.array(rho,mask=(rho<contLevel))
2870    indices = (-1,0,1)
2871    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2872    for roll in rolls:
2873        if np.any(roll):
2874            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
2875    indx = np.transpose(rhoMask.nonzero())
2876    peaks = indx/incre
2877    mags = rhoMask[rhoMask.nonzero()]
2878    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
2879        rho = rollMap(rho,ind)
2880        rMM = mapHalf-steps
2881        rMP = mapHalf+steps+1
2882        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2883        peakInt = np.sum(rhoPeak)*res**3
2884        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2885        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
2886        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
2887            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
2888        x1 = result[0]
2889        if not np.any(x1 < 0):
2890            mag = x1[0]
2891            peak = (np.array(x1[1:4])-ind)/incre
2892        peak = fixSpecialPos(peak,SGData,Amat)
2893        rho = rollMap(rho,-ind)
2894    cent = np.ones(3)*.5     
2895    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
2896    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
2897    if Neg:     #want negative magnitudes for negative peaks
2898        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
2899    else:
2900        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
2901   
2902def sortArray(data,pos,reverse=False):
2903    '''data is a list of items
2904    sort by pos in list; reverse if True
2905    '''
2906    T = []
2907    for i,M in enumerate(data):
2908        try:
2909            T.append((M[pos],i))
2910        except IndexError:
2911            return data
2912    D = dict(zip(T,data))
2913    T.sort()
2914    if reverse:
2915        T.reverse()
2916    X = []
2917    for key in T:
2918        X.append(D[key])
2919    return X
2920
2921def PeaksEquiv(data,Ind):
2922    '''Find the equivalent map peaks for those selected. Works on the
2923    contents of data['Map Peaks'].
2924
2925    :param data: the phase data structure
2926    :param list Ind: list of selected peak indices
2927    :returns: augmented list of peaks including those related by symmetry to the
2928      ones in Ind
2929
2930    '''       
2931    def Duplicate(xyz,peaks,Amat):
2932        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2933            return True
2934        return False
2935                           
2936    generalData = data['General']
2937    cell = generalData['Cell'][1:7]
2938    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2939    A = G2lat.cell2A(cell)
2940    SGData = generalData['SGData']
2941    mapPeaks = data['Map Peaks']
2942    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
2943    Indx = {}
2944    for ind in Ind:
2945        xyz = np.array(mapPeaks[ind][1:4])
2946        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
2947        for jnd,xyz in enumerate(XYZ):       
2948            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
2949    Ind = []
2950    for ind in Indx:
2951        if Indx[ind]:
2952            Ind.append(ind)
2953    return Ind
2954               
2955def PeaksUnique(data,Ind):
2956    '''Finds the symmetry unique set of peaks from those selected. Works on the
2957    contents of data['Map Peaks'].
2958
2959    :param data: the phase data structure
2960    :param list Ind: list of selected peak indices
2961    :returns: the list of symmetry unique peaks from among those given in Ind
2962
2963    '''       
2964#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
2965
2966    def noDuplicate(xyz,peaks,Amat):
2967        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2968            return False
2969        return True
2970                           
2971    generalData = data['General']
2972    cell = generalData['Cell'][1:7]
2973    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2974    A = G2lat.cell2A(cell)
2975    SGData = generalData['SGData']
2976    mapPeaks = data['Map Peaks']
2977    Indx = {}
2978    XYZ = {}
2979    for ind in Ind:
2980        XYZ[ind] = np.array(mapPeaks[ind][1:4])
2981        Indx[ind] = True
2982    for ind in Ind:
2983        if Indx[ind]:
2984            xyz = XYZ[ind]
2985            for jnd in Ind:
2986                if ind != jnd and Indx[jnd]:                       
2987                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
2988                    xyzs = np.array([equiv[0] for equiv in Equiv])
2989                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
2990    Ind = []
2991    for ind in Indx:
2992        if Indx[ind]:
2993            Ind.append(ind)
2994    return Ind
2995   
2996################################################################################
2997##### single peak fitting profile fxn stuff
2998################################################################################
2999
3000def getCWsig(ins,pos):
3001    '''get CW peak profile sigma
3002   
3003    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
3004      as values only
3005    :param float pos: 2-theta of peak
3006    :returns: float getCWsig: peak sigma
3007   
3008    '''
3009    tp = tand(pos/2.0)
3010    return ins['U']*tp**2+ins['V']*tp+ins['W']
3011   
3012def getCWsigDeriv(pos):
3013    '''get derivatives of CW peak profile sigma wrt U,V, & W
3014   
3015    :param float pos: 2-theta of peak
3016   
3017    :returns: list getCWsigDeriv: d(sig)/dU, d(sig)/dV & d(sig)/dW
3018   
3019    '''
3020    tp = tand(pos/2.0)
3021    return tp**2,tp,1.0
3022   
3023def getCWgam(ins,pos):
3024    '''get CW peak profile gamma
3025   
3026    :param dict ins: instrument parameters with at least 'X' & 'Y'
3027      as values only
3028    :param float pos: 2-theta of peak
3029    :returns: float getCWgam: peak gamma
3030   
3031    '''
3032    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
3033   
3034def getCWgamDeriv(pos):
3035    '''get derivatives of CW peak profile gamma wrt X & Y
3036   
3037    :param float pos: 2-theta of peak
3038   
3039    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
3040   
3041    '''
3042    return 1./cosd(pos/2.0),tand(pos/2.0)
3043   
3044def getTOFsig(ins,dsp):
3045    '''get TOF peak profile sigma
3046   
3047    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
3048      as values only
3049    :param float dsp: d-spacing of peak
3050   
3051    :returns: float getTOFsig: peak sigma
3052   
3053    '''
3054    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
3055   
3056def getTOFsigDeriv(dsp):
3057    '''get derivatives of TOF peak profile gamma wrt sig-0, sig-1, & sig-q
3058   
3059    :param float dsp: d-spacing of peak
3060   
3061    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
3062   
3063    '''
3064    return 1.0,dsp**2,dsp**4,1./dsp**2
3065   
3066def getTOFgamma(ins,dsp):
3067    '''get TOF peak profile gamma
3068   
3069    :param dict ins: instrument parameters with at least 'X' & 'Y'
3070      as values only
3071    :param float dsp: d-spacing of peak
3072   
3073    :returns: float getTOFgamma: peak gamma
3074   
3075    '''
3076    return ins['X']*dsp+ins['Y']*dsp**2
3077   
3078def getTOFgammaDeriv(dsp):
3079    '''get derivatives of TOF peak profile gamma wrt X & Y
3080   
3081    :param float dsp: d-spacing of peak
3082   
3083    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
3084   
3085    '''
3086    return dsp,dsp**2
3087   
3088def getTOFbeta(ins,dsp):
3089    '''get TOF peak profile beta
3090   
3091    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
3092      as values only
3093    :param float dsp: d-spacing of peak
3094   
3095    :returns: float getTOFbeta: peak beat
3096   
3097    '''
3098    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
3099   
3100def getTOFbetaDeriv(dsp):
3101    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
3102   
3103    :param float dsp: d-spacing of peak
3104   
3105    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
3106   
3107    '''
3108    return 1.0,1./dsp**4,1./dsp**2
3109   
3110def getTOFalpha(ins,dsp):
3111    '''get TOF peak profile alpha
3112   
3113    :param dict ins: instrument parameters with at least 'alpha'
3114      as values only
3115    :param float dsp: d-spacing of peak
3116   
3117    :returns: flaot getTOFalpha: peak alpha
3118   
3119    '''
3120    return ins['alpha']/dsp
3121   
3122def getTOFalphaDeriv(dsp):
3123    '''get derivatives of TOF peak profile beta wrt alpha
3124   
3125    :param float dsp: d-spacing of peak
3126   
3127    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
3128   
3129    '''
3130    return 1./dsp
3131   
3132def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
3133    '''set starting peak parameters for single peak fits from plot selection or auto selection
3134   
3135    :param dict Parms: instrument parameters dictionary
3136    :param dict Parms2: table lookup for TOF profile coefficients
3137    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
3138    :param float mag: peak top magnitude from pick
3139    :param bool ifQ: True if pos in Q
3140    :param bool useFit: True if use fitted CW Parms values (not defaults)
3141   
3142    :returns: list XY: peak list entry:
3143        for CW: [pos,0,mag,1,sig,0,gam,0]
3144        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3145        NB: mag refinement set by default, all others off
3146   
3147    '''
3148    ind = 0
3149    if useFit:
3150        ind = 1
3151    ins = {}
3152    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
3153        for x in ['U','V','W','X','Y']:
3154            ins[x] = Parms[x][ind]
3155        if ifQ:                              #qplot - convert back to 2-theta
3156            pos = 2.0*asind(pos*wave/(4*math.pi))
3157        sig = getCWsig(ins,pos)
3158        gam = getCWgam(ins,pos)           
3159        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
3160    else:
3161        if ifQ:
3162            dsp = 2.*np.pi/pos
3163            pos = Parms['difC']*dsp
3164        else:
3165            dsp = pos/Parms['difC'][1]
3166        if 'Pdabc' in Parms2:
3167            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
3168                ins[x] = Parms[x][ind]
3169            Pdabc = Parms2['Pdabc'].T
3170            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
3171            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
3172        else:
3173            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
3174                ins[x] = Parms[x][ind]
3175            alp = getTOFalpha(ins,dsp)
3176            bet = getTOFbeta(ins,dsp)
3177        sig = getTOFsig(ins,dsp)
3178        gam = getTOFgamma(ins,dsp)
3179        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3180    return XY
3181   
3182################################################################################
3183##### MC/SA stuff
3184################################################################################
3185
3186#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
3187# Original Author: Travis Oliphant 2002
3188# Bug-fixes in 2006 by Tim Leslie
3189
3190
3191import numpy
3192from numpy import asarray, tan, exp, ones, squeeze, sign, \
3193     all, log, sqrt, pi, shape, array, minimum, where
3194from numpy import random
3195
3196#__all__ = ['anneal']
3197
3198_double_min = numpy.finfo(float).min
3199_double_max = numpy.finfo(float).max
3200class base_schedule(object):
3201    def __init__(self):
3202        self.dwell = 20
3203        self.learn_rate = 0.5
3204        self.lower = -10
3205        self.upper = 10
3206        self.Ninit = 50
3207        self.accepted = 0
3208        self.tests = 0
3209        self.feval = 0
3210        self.k = 0
3211        self.T = None
3212
3213    def init(self, **options):
3214        self.__dict__.update(options)
3215        self.lower = asarray(self.lower)
3216        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
3217        self.upper = asarray(self.upper)
3218        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
3219        self.k = 0
3220        self.accepted = 0
3221        self.feval = 0
3222        self.tests = 0
3223
3224    def getstart_temp(self, best_state):
3225        """ Find a matching starting temperature and starting parameters vector
3226        i.e. find x0 such that func(x0) = T0.
3227
3228        Parameters
3229        ----------
3230        best_state : _state
3231            A _state object to store the function value and x0 found.
3232
3233        returns
3234        -------
3235        x0 : array
3236            The starting parameters vector.
3237        """
3238
3239        assert(not self.dims is None)
3240        lrange = self.lower
3241        urange = self.upper
3242        fmax = _double_min
3243        fmin = _double_max
3244        for _ in range(self.Ninit):
3245            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
3246            fval = self.func(x0, *self.args)
3247            self.feval += 1
3248            if fval > fmax:
3249                fmax = fval
3250            if fval < fmin:
3251                fmin = fval
3252                best_state.cost = fval
3253                best_state.x = array(x0)
3254
3255        self.T0 = (fmax-fmin)*1.5
3256        return best_state.x
3257       
3258    def set_range(self,x0,frac):
3259        delrange = frac*(self.upper-self.lower)
3260        self.upper = x0+delrange
3261        self.lower = x0-delrange
3262
3263    def accept_test(self, dE):
3264        T = self.T
3265        self.tests += 1
3266        if dE < 0:
3267            self.accepted += 1
3268            return 1
3269        p = exp(-dE*1.0/self.boltzmann/T)
3270        if (p > random.uniform(0.0, 1.0)):
3271            self.accepted += 1
3272            return 1
3273        return 0
3274
3275    def update_guess(self, x0):
3276        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
3277
3278    def update_temp(self, x0):
3279        pass
3280
3281
3282#  A schedule due to Lester Ingber modified to use bounds - OK
3283class fast_sa(base_schedule):
3284    def init(self, **options):
3285        self.__dict__.update(options)
3286        if self.m is None:
3287            self.m = 1.0
3288        if self.n is None:
3289            self.n = 1.0
3290        self.c = self.m * exp(-self.n * self.quench)
3291
3292    def update_guess(self, x0):
3293        x0 = asarray(x0)
3294        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
3295        T = self.T
3296        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
3297        xnew = xc*(self.upper - self.lower)+self.lower
3298        return xnew
3299#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
3300#        xc = y*(self.upper - self.lower)
3301#        xnew = x0 + xc
3302#        return xnew
3303
3304    def update_temp(self):
3305        self.T = self.T0*exp(-self.c * self.k**(self.quench))
3306        self.k += 1
3307        return
3308
3309class cauchy_sa(base_schedule):     #modified to use bounds - not good
3310    def update_guess(self, x0):
3311        x0 = asarray(x0)
3312        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
3313        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
3314        xnew = xc*(self.upper - self.lower)+self.lower
3315        return xnew
3316#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
3317#        xc = self.learn_rate * self.T * tan(numbers)
3318#        xnew = x0 + xc
3319#        return xnew
3320
3321    def update_temp(self):
3322        self.T = self.T0/(1+self.k)
3323        self.k += 1
3324        return
3325
3326class boltzmann_sa(base_schedule):
3327#    def update_guess(self, x0):
3328#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
3329#        x0 = asarray(x0)
3330#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
3331#
3332#        xnew = x0 + xc*std*self.learn_rate
3333#        return xnew
3334
3335    def update_temp(self):
3336        self.k += 1
3337        self.T = self.T0 / log(self.k+1.0)
3338        return
3339
3340class log_sa(base_schedule):        #OK
3341
3342    def init(self,**options):
3343        self.__dict__.update(options)
3344       
3345    def update_guess(self,x0):     #same as default
3346        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
3347       
3348    def update_temp(self):
3349        self.k += 1
3350        self.T = self.T0*self.slope**self.k
3351       
3352class _state(object):
3353    def __init__(self):
3354        self.x = None
3355        self.cost = None
3356
3357# TODO:
3358#     allow for general annealing temperature profile
3359#     in that case use update given by alpha and omega and
3360#     variation of all previous updates and temperature?
3361
3362# Simulated annealing
3363
3364def anneal(func, x0, args=(), schedule='fast', full_output=0,
3365           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
3366           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
3367           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
3368           ranRange=0.10,autoRan=False,dlg=None):
3369    """Minimize a function using simulated annealing.
3370
3371    Schedule is a schedule class implementing the annealing schedule.
3372    Available ones are 'fast', 'cauchy', 'boltzmann'
3373
3374    :param callable func: f(x, \*args)
3375        Function to be optimized.
3376    :param ndarray x0:
3377        Initial guess.
3378    :param tuple args:
3379        Extra parameters to `func`.
3380    :param base_schedule schedule:
3381        Annealing schedule to use (a class).
3382    :param bool full_output:
3383        Whether to return optional outputs.
3384    :param float T0:
3385        Initial Temperature (estimated as 1.2 times the largest
3386        cost-function deviation over random points in the range).
3387    :param float Tf:
3388        Final goal temperature.
3389    :param int maxeval:
3390        Maximum function evaluations.
3391    :param int maxaccept:
3392        Maximum changes to accept.
3393    :param int maxiter:
3394        Maximum cooling iterations.
3395    :param float learn_rate:
3396        Scale constant for adjusting guesses.
3397    :param float boltzmann:
3398        Boltzmann constant in acceptance test
3399        (increase for less stringent test at each temperature).
3400    :param float feps:
3401        Stopping relative error tolerance for the function value in
3402        last four coolings.
3403    :param float quench,m,n:
3404        Parameters to alter fast_sa schedule.
3405    :param float/ndarray lower,upper:
3406        Lower and upper bounds on `x`.
3407    :param int dwell:
3408        The number of times to search the space at each temperature.
3409    :param float slope:
3410        Parameter for log schedule
3411    :param bool ranStart=False:
3412        True for set 10% of ranges about x
3413
3414    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
3415
3416     * xmin (ndarray): Point giving smallest value found.
3417     * Jmin (float): Minimum value of function found.
3418     * T (float): Final temperature.
3419     * feval (int): Number of function evaluations.
3420     * iters (int): Number of cooling iterations.
3421     * accept (int): Number of tests accepted.
3422     * retval (int): Flag indicating stopping condition:
3423
3424              *  0: Points no longer changing
3425              *  1: Cooled to final temperature
3426              *  2: Maximum function evaluations
3427              *  3: Maximum cooling iterations reached
3428              *  4: Maximum accepted query locations reached
3429              *  5: Final point not the minimum amongst encountered points
3430
3431    *Notes*:
3432    Simulated annealing is a random algorithm which uses no derivative
3433    information from the function being optimized. In practice it has
3434    been more useful in discrete optimization than continuous
3435    optimization, as there are usually better algorithms for continuous
3436    optimization problems.
3437
3438    Some experimentation by trying the difference temperature
3439    schedules and altering their parameters is likely required to
3440    obtain good performance.
3441
3442    The randomness in the algorithm comes from random sampling in numpy.
3443    To obtain the same results you can call numpy.random.seed with the
3444    same seed immediately before calling scipy.optimize.anneal.
3445
3446    We give a brief description of how the three temperature schedules
3447    generate new points and vary their temperature. Temperatures are
3448    only updated with iterations in the outer loop. The inner loop is
3449    over xrange(dwell), and new points are generated for
3450    every iteration in the inner loop. (Though whether the proposed
3451    new points are accepted is probabilistic.)
3452
3453    For readability, let d denote the dimension of the inputs to func.
3454    Also, let x_old denote the previous state, and k denote the
3455    iteration number of the outer loop. All other variables not
3456    defined below are input variables to scipy.optimize.anneal itself.
3457
3458    In the 'fast' schedule the updates are ::
3459
3460        u ~ Uniform(0, 1, size=d)
3461        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
3462        xc = y * (upper - lower)
3463        x_new = x_old + xc
3464
3465        c = n * exp(-n * quench)
3466        T_new = T0 * exp(-c * k**quench)
3467
3468
3469    In the 'cauchy' schedule the updates are ::
3470
3471        u ~ Uniform(-pi/2, pi/2, size=d)
3472        xc = learn_rate * T * tan(u)
3473        x_new = x_old + xc
3474
3475        T_new = T0 / (1+k)
3476
3477    In the 'boltzmann' schedule the updates are ::
3478
3479        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
3480        y ~ Normal(0, std, size=d)
3481        x_new = x_old + learn_rate * y
3482
3483        T_new = T0 / log(1+k)
3484
3485    """
3486    x0 = asarray(x0)
3487    lower = asarray(lower)
3488    upper = asarray(upper)
3489
3490    schedule = eval(schedule+'_sa()')
3491    #   initialize the schedule
3492    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
3493                  learn_rate=learn_rate, lower=lower, upper=upper,
3494                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
3495
3496    current_state, last_state, best_state = _state(), _state(), _state()
3497    if ranStart:
3498        schedule.set_range(x0,ranRange)
3499    if T0 is None:
3500        x0 = schedule.getstart_temp(best_state)
3501    else:
3502        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
3503        best_state.x = None
3504        best_state.cost = numpy.Inf
3505
3506    last_state.x = asarray(x0).copy()
3507    fval = func(x0,*args)
3508    schedule.feval += 1
3509    last_state.cost = fval
3510    if last_state.cost < best_state.cost:
3511        best_state.cost = fval
3512        best_state.x = asarray(x0).copy()
3513    schedule.T = schedule.T0
3514    fqueue = [100, 300, 500, 700]
3515    iters = 1
3516    keepGoing = True
3517    bestn = 0
3518    while keepGoing:
3519        retval = 0
3520        for n in xrange(dwell):
3521            current_state.x = schedule.update_guess(last_state.x)
3522            current_state.cost = func(current_state.x,*args)
3523            schedule.feval += 1
3524
3525            dE = current_state.cost - last_state.cost
3526            if schedule.accept_test(dE):
3527                last_state.x = current_state.x.copy()
3528                last_state.cost = current_state.cost
3529                if last_state.cost < best_state.cost:
3530                    best_state.x = last_state.x.copy()
3531                    best_state.cost = last_state.cost
3532                    bestn = n
3533                    if best_state.cost < 1.0 and autoRan:
3534                        schedule.set_range(x0,best_state.cost/2.)                       
3535        if dlg:
3536            GoOn = dlg.Update(min(100.,best_state.cost*100),
3537                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
3538                    'Best trial:',bestn,  \
3539                    'MC/SA Residual:',best_state.cost*100,'%', \
3540                    ))[0]
3541            if not GoOn:
3542                best_state.x = last_state.x.copy()
3543                best_state.cost = last_state.cost
3544                retval = 5
3545        schedule.update_temp()
3546        iters += 1
3547        # Stopping conditions
3548        # 0) last saved values of f from each cooling step
3549        #     are all very similar (effectively cooled)
3550        # 1) Tf is set and we are below it
3551        # 2) maxeval is set and we are past it
3552        # 3) maxiter is set and we are past it
3553        # 4) maxaccept is set and we are past it
3554        # 5) user canceled run via progress bar
3555
3556        fqueue.append(squeeze(last_state.cost))
3557        fqueue.pop(0)
3558        af = asarray(fqueue)*1.0
3559        if retval == 5:
3560            print ' User terminated run; incomplete MC/SA'
3561            keepGoing = False
3562            break
3563        if all(abs((af-af[0])/af[0]) < feps):
3564            retval = 0
3565            if abs(af[-1]-best_state.cost) > feps*10:
3566                retval = 5
3567#                print "Warning: Cooled to %f at %s but this is not" \
3568#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
3569#                      + " the smallest point found."
3570            break
3571        if (Tf is not None) and (schedule.T < Tf):
3572            retval = 1
3573            break
3574        if (maxeval is not None) and (schedule.feval > maxeval):
3575            retval = 2
3576            break
3577        if (iters > maxiter):
3578            print "Warning: Maximum number of iterations exceeded."
3579            retval = 3
3580            break
3581        if (maxaccept is not None) and (schedule.accepted > maxaccept):
3582            retval = 4
3583            break
3584
3585    if full_output:
3586        return best_state.x, best_state.cost, schedule.T, \
3587               schedule.feval, iters, schedule.accepted, retval
3588    else:
3589        return best_state.x, retval
3590
3591def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
3592    outlist = []
3593    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
3594    for n in range(iCyc):
3595        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
3596        outlist.append(result[0])
3597        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
3598    out_q.put(outlist)
3599
3600def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
3601    import multiprocessing as mp
3602   
3603    nprocs = mp.cpu_count()
3604    out_q = mp.Queue()
3605    procs = []
3606    iCyc = np.zeros(nprocs)
3607    for i in range(nCyc):
3608        iCyc[i%nprocs] += 1
3609    for i in range(nprocs):
3610        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
3611        procs.append(p)
3612        p.start()
3613    resultlist = []
3614    for i in range(nprocs):
3615        resultlist += out_q.get()
3616    for p in procs:
3617        p.join()
3618    return resultlist
3619
3620def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
3621    '''default doc string
3622   
3623    :param type name: description
3624   
3625    :returns: type name: description
3626   
3627    '''
3628   
3629    global tsum
3630    tsum = 0.
3631   
3632    def getMDparms(item,pfx,parmDict,varyList):
3633        parmDict[pfx+'MDaxis'] = item['axis']
3634        parmDict[pfx+'MDval'] = item['Coef'][0]
3635        if item['Coef'][1]:
3636            varyList += [pfx+'MDval',]
3637            limits = item['Coef'][2]
3638            lower.append(limits[0])
3639            upper.append(limits[1])
3640                       
3641    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
3642        parmDict[pfx+'Atype'] = item['atType']
3643        aTypes |= set([item['atType'],]) 
3644        pstr = ['Ax','Ay','Az']
3645        XYZ = [0,0,0]
3646        for i in range(3):
3647            name = pfx+pstr[i]
3648            parmDict[name] = item['Pos'][0][i]
3649            XYZ[i] = parmDict[name]
3650            if item['Pos'][1][i]:
3651                varyList += [name,]
3652                limits = item['Pos'][2][i]
3653                lower.append(limits[0])
3654                upper.append(limits[1])
3655        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
3656           
3657    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
3658        parmDict[mfx+'MolCent'] = item['MolCent']
3659        parmDict[mfx+'RBId'] = item['RBId']
3660        pstr = ['Px','Py','Pz']
3661        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
3662        for i in range(3):
3663            name = mfx+pstr[i]
3664            parmDict[name] = item['Pos'][0][i]
3665            if item['Pos'][1][i]:
3666                varyList += [name,]
3667                limits = item['Pos'][2][i]
3668                lower.append(limits[0])
3669                upper.append(limits[1])
3670        AV = item['Ori'][0]
3671        A = AV[0]
3672        V = AV[1:]
3673        for i in range(4):
3674            name = mfx+ostr[i]
3675            if i:
3676                parmDict[name] = V[i-1]
3677            else:
3678                parmDict[name] = A
3679            if item['Ovar'] == 'AV':
3680                varyList += [name,]
3681                limits = item['Ori'][2][i]
3682                lower.append(limits[0])
3683                upper.append(limits[1])
3684            elif item['Ovar'] == 'A' and not i:
3685                varyList += [name,]
3686                limits = item['Ori'][2][i]
3687                lower.append(limits[0])
3688                upper.append(limits[1])
3689        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
3690            for i in range(len(item['Tor'][0])):
3691                name = mfx+'Tor'+str(i)
3692                parmDict[name] = item['Tor'][0][i]
3693                if item['Tor'][1][i]:
3694                    varyList += [name,]
3695                    limits = item['Tor'][2][i]
3696                    lower.append(limits[0])
3697                    upper.append(limits[1])
3698        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
3699        aTypes |= set(atypes)
3700        atNo += len(atypes)
3701        return atNo
3702               
3703    def GetAtomM(Xdata,SGData):
3704        Mdata = []
3705        for xyz in Xdata:
3706            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
3707        return np.array(Mdata)
3708       
3709    def GetAtomT(RBdata,parmDict):
3710        'Needs a doc string'
3711        atNo = parmDict['atNo']
3712        nfixAt = parmDict['nfixAt']
3713        Tdata = atNo*[' ',]
3714        for iatm in range(nfixAt):
3715            parm = ':'+str(iatm)+':Atype'
3716            if parm in parmDict:
3717                Tdata[iatm] = aTypes.index(parmDict[parm])
3718        iatm = nfixAt
3719        for iObj in range(parmDict['nObj']):
3720            pfx = str(iObj)+':'
3721            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3722                if parmDict[pfx+'Type'] == 'Vector':
3723                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3724                    nAtm = len(RBRes['rbVect'][0])
3725                else:       #Residue
3726                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3727                    nAtm = len(RBRes['rbXYZ'])
3728                for i in range(nAtm):
3729                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
3730                    iatm += 1
3731            elif parmDict[pfx+'Type'] == 'Atom':
3732                atNo = parmDict[pfx+'atNo']
3733                parm = pfx+'Atype'              #remove extra ':'
3734                if parm in parmDict:
3735                    Tdata[atNo] = aTypes.index(parmDict[parm])
3736                iatm += 1
3737            else:
3738                continue        #skips March Dollase
3739        return Tdata
3740       
3741    def GetAtomX(RBdata,parmDict):
3742        'Needs a doc string'
3743        Bmat = parmDict['Bmat']
3744        atNo = parmDict['atNo']
3745        nfixAt = parmDict['nfixAt']
3746        Xdata = np.zeros((3,atNo))
3747        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
3748        for iatm in range(nfixAt):
3749            for key in keys:
3750                parm = ':'+str(iatm)+key
3751                if parm in parmDict:
3752                    keys[key][iatm] = parmDict[parm]
3753        iatm = nfixAt
3754        for iObj in range(parmDict['nObj']):
3755            pfx = str(iObj)+':'
3756            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3757                if parmDict[pfx+'Type'] == 'Vector':
3758                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3759                    vecs = RBRes['rbVect']
3760                    mags = RBRes['VectMag']
3761                    Cart = np.zeros_like(vecs[0])
3762                    for vec,mag in zip(vecs,mags):
3763                        Cart += vec*mag
3764                elif parmDict[pfx+'Type'] == 'Residue':
3765                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3766                    Cart = np.array(RBRes['rbXYZ'])
3767                    for itor,seq in enumerate(RBRes['rbSeq']):
3768                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
3769                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
3770                if parmDict[pfx+'MolCent'][1]:
3771                    Cart -= parmDict[pfx+'MolCent'][0]
3772                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3773                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
3774                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
3775                iatm += len(Cart)
3776            elif parmDict[pfx+'Type'] == 'Atom':
3777                atNo = parmDict[pfx+'atNo']
3778                for key in keys:
3779                    parm = pfx+key[1:]              #remove extra ':'
3780                    if parm in parmDict:
3781                        keys[key][atNo] = parmDict[parm]
3782                iatm += 1
3783            else:
3784                continue        #skips March Dollase
3785        return Xdata.T
3786       
3787    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
3788        allX = np.inner(Xdata,SGM)+SGT
3789        allT = np.repeat(Tdata,allX.shape[1])
3790        allM = np.repeat(Mdata,allX.shape[1])
3791        allX = np.reshape(allX,(-1,3))
3792        return allT,allM,allX
3793
3794    def getAllX(Xdata,SGM,SGT):
3795        allX = np.inner(Xdata,SGM)+SGT
3796        allX = np.reshape(allX,(-1,3))
3797        return allX
3798       
3799    def normQuaternions(RBdata,parmDict,varyList,values):
3800        for iObj in range(parmDict['nObj']):
3801            pfx = str(iObj)+':'
3802            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3803                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3804                A,V = Q2AVdeg(Qori)
3805                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
3806                    if i:
3807                        parmDict[pfx+name] = V[i-1]
3808                    else:
3809                        parmDict[pfx+name] = A
3810       
3811    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
3812        ''' Compute structure factors for all h,k,l for phase
3813        input:
3814            refList: [ref] where each ref = h,k,l,m,d,...
3815            rcov:   array[nref,nref] covariance terms between Fo^2 values
3816            ifInv:  bool True if centrosymmetric
3817            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
3818            RBdata: [dict] rigid body dictionary
3819            varyList: [list] names of varied parameters in MC/SA (not used here)           
3820            ParmDict: [dict] problem parameters
3821        puts result F^2 in each ref[5] in refList
3822        returns:
3823            delt-F*rcov*delt-F/sum(Fo^2)^2
3824        '''       
3825        global tsum
3826        t0 = time.time()
3827        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
3828        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
3829        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
3830        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
3831        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
3832        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
3833        refList[5] = Aterm
3834        if not ifInv:
3835            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
3836        if len(cosTable):        #apply MD correction
3837            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
3838        sumFcsq = np.sum(refList[5])
3839        scale = parmDict['sumFosq']/sumFcsq
3840        refList[5] *= scale
3841        refList[6] = refList[4]-refList[5]
3842        M = np.inner(refList[6],np.inner(rcov,refList[6]))
3843        tsum += (time.time()-t0)
3844        return M/np.sum(refList[4]**2)
3845
3846    sq8ln2 = np.sqrt(8*np.log(2))
3847    sq2pi = np.sqrt(2*np.pi)
3848    sq4pi = np.sqrt(4*np.pi)
3849    generalData = data['General']
3850    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3851    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
3852    SGData = generalData['SGData']
3853    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
3854    SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
3855    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
3856    fixAtoms = data['Atoms']                       #if any
3857    cx,ct,cs = generalData['AtomPtrs'][:3]
3858    aTypes = set([])
3859    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
3860    varyList = []
3861    atNo = 0
3862    for atm in fixAtoms:
3863        pfx = ':'+str(atNo)+':'
3864        parmDict[pfx+'Atype'] = atm[ct]
3865        aTypes |= set([atm[ct],])
3866        pstr = ['Ax','Ay','Az']
3867        parmDict[pfx+'Amul'] = atm[cs+1]
3868        for i in range(3):
3869            name = pfx+pstr[i]
3870            parmDict[name] = atm[cx+i]
3871        atNo += 1
3872    parmDict['nfixAt'] = len(fixAtoms)       
3873    MCSA = generalData['MCSA controls']
3874    reflName = MCSA['Data source']
3875    phaseName = generalData['Name']
3876    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
3877    upper = []
3878    lower = []
3879    MDvec = np.zeros(3)
3880    for i,item in enumerate(MCSAObjs):
3881        mfx = str(i)+':'
3882        parmDict[mfx+'Type'] = item['Type']
3883        if item['Type'] == 'MD':
3884            getMDparms(item,mfx,parmDict,varyList)
3885            MDvec = np.array(item['axis'])
3886        elif item['Type'] == 'Atom':
3887            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
3888            parmDict[mfx+'atNo'] = atNo
3889            atNo += 1
3890        elif item['Type'] in ['Residue','Vector']:
3891            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
3892    parmDict['atNo'] = atNo                 #total no. of atoms
3893    parmDict['nObj'] = len(MCSAObjs)
3894    aTypes = list(aTypes)
3895    Tdata = GetAtomT(RBdata,parmDict)
3896    Xdata = GetAtomX(RBdata,parmDict)
3897    Mdata = GetAtomM(Xdata,SGData)
3898    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
3899    FFtables = G2el.GetFFtable(aTypes)
3900    refs = []
3901    allFF = []
3902    cosTable = []
3903    sumFosq = 0
3904    if 'PWDR' in reflName:
3905        for ref in reflData:
3906            h,k,l,m,d,pos,sig,gam,f = ref[:9]
3907            if d >= MCSA['dmin']:
3908                sig = G2pwd.getgamFW(sig,gam)/sq8ln2        #--> sig from FWHM
3909                SQ = 0.25/d**2
3910                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3911                refs.append([h,k,l,m,f*m,pos,sig])
3912                sumFosq += f*m
3913                Heqv = np.inner(np.array([h,k,l]),SGMT)
3914                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3915        nRef = len(refs)
3916        cosTable = np.array(cosTable)**2
3917        rcov = np.zeros((nRef,nRef))
3918        for iref,refI in enumerate(refs):
3919            rcov[iref][iref] = 1./(sq4pi*refI[6])
3920            for jref,refJ in enumerate(refs[:iref]):
3921                t1 = refI[6]**2+refJ[6]**2
3922                t2 = (refJ[5]-refI[5])**2/(2.*t1)
3923                if t2 > 10.:
3924                    rcov[iref][jref] = 0.
3925                else:
3926                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
3927        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3928        Rdiag = np.sqrt(np.diag(rcov))
3929        Rnorm = np.outer(Rdiag,Rdiag)
3930        rcov /= Rnorm
3931    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
3932        vNames = []
3933        pfx = str(data['pId'])+'::PWLref:'
3934        for iref,refI in enumerate(reflData):           #Pawley reflection set
3935            h,k,l,m,d,v,f,s = refI
3936            if d >= MCSA['dmin'] and v:       #skip unrefined ones
3937                vNames.append(pfx+str(iref))
3938                SQ = 0.25/d**2
3939                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3940                refs.append([h,k,l,m,f*m,iref,0.])
3941                sumFosq += f*m
3942                Heqv = np.inner(np.array([h,k,l]),SGMT)
3943                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3944        cosTable = np.array(cosTable)**2
3945        nRef = len(refs)
3946#        if generalData['doPawley'] and (covData['freshCOV'] or  MCSA['newDmin']):
3947        if covData['freshCOV'] or  MCSA['newDmin']:
3948            vList = covData['varyList']
3949            covMatrix = covData['covMatrix']
3950            rcov = getVCov(vNames,vList,covMatrix)
3951            rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3952            Rdiag = np.sqrt(np.diag(rcov))
3953            Rdiag = np.where(Rdiag,Rdiag,1.0)
3954            Rnorm = np.outer(Rdiag,Rdiag)
3955            rcov /= Rnorm
3956            MCSA['rcov'] = rcov
3957            covData['freshCOV'] = False
3958            MCSA['newDmin'] = False
3959        else:
3960            rcov = MCSA['rcov']
3961    elif 'HKLF' in reflName:
3962        for ref in reflData:
3963            [h,k,l,m,d],f = ref[:5],ref[6]
3964            if d >= MCSA['dmin']:
3965                SQ = 0.25/d**2
3966                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3967                refs.append([h,k,l,m,f*m,0.,0.])
3968                sumFosq += f*m
3969        nRef = len(refs)
3970        rcov = np.identity(len(refs))
3971    allFF = np.array(allFF).T
3972    refs = np.array(refs).T
3973    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
3974    print ' Number of parameters varied: %d'%(len(varyList))
3975    parmDict['sumFosq'] = sumFosq
3976    x0 = [parmDict[val] for val in varyList]
3977    ifInv = SGData['SGInv']
3978    # consider replacing anneal with scipy.optimize.basinhopping
3979    results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
3980        schedule=MCSA['Algorithm'], full_output=True,
3981        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
3982        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
3983        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
3984        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
3985        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
3986    M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
3987#    for ref in refs.T:
3988#        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])
3989#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
3990    Result = [False,False,results[1],results[2],]+list(results[0])
3991    Result.append(varyList)
3992    return Result,tsum
3993
3994       
3995################################################################################
3996##### Quaternion stuff
3997################################################################################
3998
3999def prodQQ(QA,QB):
4000    ''' Grassman quaternion product
4001        QA,QB quaternions; q=r+ai+bj+ck
4002    '''
4003    D = np.zeros(4)
4004    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
4005    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
4006    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
4007    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
4008   
4009#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
4010#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
4011   
4012    return D
4013   
4014def normQ(QA):
4015    ''' get length of quaternion & normalize it
4016        q=r+ai+bj+ck
4017    '''
4018    n = np.sqrt(np.sum(np.array(QA)**2))
4019    return QA/n
4020   
4021def invQ(Q):
4022    '''
4023        get inverse of quaternion
4024        q=r+ai+bj+ck; q* = r-ai-bj-ck
4025    '''
4026    return Q*np.array([1,-1,-1,-1])
4027   
4028def prodQVQ(Q,V):
4029    """
4030    compute the quaternion vector rotation qvq-1 = v'
4031    q=r+ai+bj+ck
4032    """
4033    T2 = Q[0]*Q[1]
4034    T3 = Q[0]*Q[2]
4035    T4 = Q[0]*Q[3]
4036    T5 = -Q[1]*Q[1]
4037    T6 = Q[1]*Q[2]
4038    T7 = Q[1]*Q[3]
4039    T8 = -Q[2]*Q[2]
4040    T9 = Q[2]*Q[3]
4041    T10 = -Q[3]*Q[3]
4042    M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
4043    VP = 2.*np.inner(V,M)
4044    return VP+V
4045   
4046def Q2Mat(Q):
4047    ''' make rotation matrix from quaternion
4048        q=r+ai+bj+ck
4049    '''
4050    QN = normQ(Q)
4051    aa = QN[0]**2
4052    ab = QN[0]*QN[1]
4053    ac = QN[0]*QN[2]
4054    ad = QN[0]*QN[3]
4055    bb = QN[1]**2
4056    bc = QN[1]*QN[2]
4057    bd = QN[1]*QN[3]
4058    cc = QN[2]**2
4059    cd = QN[2]*QN[3]
4060    dd = QN[3]**2
4061    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
4062        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
4063        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
4064    return np.array(M)
4065   
4066def AV2Q(A,V):
4067    ''' convert angle (radians) & vector to quaternion
4068        q=r+ai+bj+ck
4069    '''
4070    Q = np.zeros(4)
4071    d = nl.norm(np.array(V))
4072    if d:
4073        V /= d
4074        if not A:       #==0.
4075            A = 2.*np.pi
4076        p = A/2.
4077        Q[0] = np.cos(p)
4078        Q[1:4] = V*np.sin(p)
4079    else:
4080        Q[3] = 1.
4081    return Q
4082   
4083def AVdeg2Q(A,V):
4084    ''' convert angle (degrees) & vector to quaternion
4085        q=r+ai+bj+ck
4086    '''
4087    Q = np.zeros(4)
4088    d = nl.norm(np.array(V))
4089    if not A:       #== 0.!
4090        A = 360.
4091    if d:
4092        V /= d
4093        p = A/2.
4094        Q[0] = cosd(p)
4095        Q[1:4] = V*sind(p)
4096    else:
4097        Q[3] = 1.
4098    return Q
4099   
4100def Q2AVdeg(Q):
4101    ''' convert quaternion to angle (degrees 0-360) & normalized vector
4102        q=r+ai+bj+ck
4103    '''
4104    A = 2.*acosd(Q[0])
4105    V = np.array(Q[1:])
4106    V /= sind(A/2.)
4107    return A,V
4108   
4109def Q2AV(Q):
4110    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
4111        q=r+ai+bj+ck
4112    '''
4113    A = 2.*np.arccos(Q[0])
4114    V = np.array(Q[1:])
4115    V /= np.sin(A/2.)
4116    return A,V
4117   
4118def randomQ(r0,r1,r2,r3):
4119    ''' create random quaternion from 4 random numbers in range (-1,1)
4120    '''
4121    sum = 0
4122    Q = np.array(4)
4123    Q[0] = r0
4124    sum += Q[0]**2
4125    Q[1] = np.sqrt(1.-sum)*r1
4126    sum += Q[1]**2
4127    Q[2] = np.sqrt(1.-sum)*r2
4128    sum += Q[2]**2
4129    Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
4130    return Q
4131   
4132def randomAVdeg(r0,r1,r2,r3):
4133    ''' create random angle (deg),vector from 4 random number in range (-1,1)
4134    '''
4135    return Q2AVdeg(randomQ(r0,r1,r2,r3))
4136   
4137def makeQuat(A,B,C):
4138    ''' Make quaternion from rotation of A vector to B vector about C axis
4139
4140        :param np.array A,B,C: Cartesian 3-vectors
4141        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
4142    '''
4143
4144    V1 = np.cross(A,C)
4145    V2 = np.cross(B,C)
4146    if nl.norm(V1)*nl.norm(V2):
4147        V1 /= nl.norm(V1)
4148        V2 /= nl.norm(V2)
4149        V3 = np.cross(V1,V2)
4150    else:
4151        V3 = np.zeros(3)
4152    Q = np.array([0.,0.,0.,1.])
4153    D = 0.
4154    if nl.norm(V3):
4155        V3 /= nl.norm(V3)
4156        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
4157        D = np.arccos(D1)/2.0
4158        V1 = C-V3
4159        V2 = C+V3
4160        DM = nl.norm(V1)
4161        DP = nl.norm(V2)
4162        S = np.sin(D)
4163        Q[0] = np.cos(D)
4164        Q[1:] = V3*S
4165        D *= 2.
4166        if DM > DP:
4167            D *= -1.
4168    return Q,D
4169   
4170def annealtests():
4171    from numpy import cos
4172    # minimum expected at ~-0.195
4173    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
4174    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
4175    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
4176    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
4177
4178    # minimum expected at ~[-0.195, -0.1]
4179    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
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='cauchy')
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='fast')
4182    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')
4183
4184
4185if __name__ == '__main__':
4186    annealtests()
Note: See TracBrowser for help on using the repository browser.