source: trunk/GSASIImath.py @ 4837

Last change on this file since 4837 was 4830, checked in by vondreele, 8 months ago

show more places in Absorb output
improve image absorption calcs for cylinders now seems to work OK
reorganize Image controls GUI & add option for vertical samples
hide printing of LS cycles if Rrint=False in HessianLSQ
implement (& not use) Klein-Nishina PDF correction
Fix bug in if importer

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