source: trunk/GSASIImath.py @ 4840

Last change on this file since 4840 was 4840, checked in by toby, 8 months ago

LeBail? issues: bug on zero cycles fixed; new LeBail? fit menu item; HessianLSQ: remove extra func eval on 0 cycles; fix noprint options; misc cleanups

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