source: trunk/GSASIImath.py @ 4789

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

switch to new HessianLSQ; Ouch#4 msg in GUI; misc wx4.1 fixes; docs improvements; switch to new HessianLSQ

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