source: trunk/GSASIImath.py

Last change on this file was 5320, checked in by toby, 5 weeks ago

more fullrmc work (all but swaps done?); misc docs fixes

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