source: trunk/GSASIImath.py @ 4990

Last change on this file since 4990 was 4990, checked in by toby, 4 months ago

renumber ouches to keep unique; give timing on ouch 4

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