source: trunk/GSASIImath.py @ 4997

Last change on this file since 4997 was 4997, checked in by vondreele, 3 months ago

move imports to top of G2ctrlGUi code so debugger can find it
fix FillUnitCell? bug when incommensurate mag.
latest revision of incomm. mag. str. factor stuff

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