source: trunk/GSASIImath.py @ 1980

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

refBlk=100 version of SStructureFactor - works

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