source: trunk/GSASIImath.py @ 1976

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

fix bug in plotting 3D HKLs
SS str. Fctr. mods - add site fraction & Uij modulation math

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