source: trunk/GSASIImath.py @ 1958

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

add r-values by SS index
add a fade for frac modulation of atoms in drawing
work on modulation structure factor

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