source: trunk/GSASIImath.py @ 5015

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

changes to comments

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