source: trunk/GSASIImath.py @ 4519

Last change on this file since 4519 was 4519, checked in by vondreele, 17 months ago

G2fpaGUI - minor cleanup
G2lattice - test on 'T' instead of 'C' so that 'B' PWDR types properly handled
G2math - add getPinkalpha, getPinkbeta % deriv routines for pink beam peak shapes

modify setPeakparams to use them

G2plot - test on 'T' instead of 'C' so that 'B' PWDR types properly handled
G2pwd - add pink beam peak shape function - same as TOF peak shape; scale sig & gam to be in centidegrees

  • add wtFactor to FitPeaks? fix bad PWDR weights for peak fitting so reduced chi2 is more appropriate

G2pwdGUI - test on 'T' instead of 'C' so that 'B' PWDR types properly handled

  • add wtFactor to DoPeakFit? to fix bad PWDR weights for peak fitting so reduced chi2 is more appropriate

G2strMath - latest in incommensurate mag str factors - no real improvement

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