source: trunk/GSASIImath.py @ 1966

Last change on this file since 1966 was 1966, checked in by vondreele, 7 years ago

new 3D HKL plot commands for selecting view axes & set a zone plot - trouble with rotations tho.
SS structure factor by integration, not Bessel fxns.

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