source: trunk/GSASIImath_new.py @ 4683

Last change on this file since 4683 was 4683, checked in by toby, 10 months ago

include trial version of G2math with reworked HessianLSQ

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