source: trunk/GSASIImath.py @ 4903

Last change on this file since 4903 was 4903, checked in by vondreele, 5 months ago

fix 2 issues in HessianLSQ for SVN LinAlgErrors?
fix a typo in G2plot for SASD plot display
change GetTthAzmDsp? --> GeTthAzmDsp2 in G2plot & G2Image
comment old GetTthAzmDsp? code - fails in recent python/numpy versions; replaced by GetTthAzmDsp2
change GetTthAzmG --> GetTthAzm2 in G2image

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