source: trunk/GSASIImath.py @ 5013

Last change on this file since 5013 was 5013, checked in by vondreele, 4 months ago

remove invalid plot drag message - just console clutter
current best effort for incomm mag str fctrs

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