source: trunk/GSASIImath.py @ 1970

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

Add pygauleg to pypowder.for - recompile Win2.7 & Win2.7-64
work on SS structure factor

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