source: trunk/GSASIImath_new.py @ 4780

Last change on this file since 4780 was 4780, checked in by toby, 11 months ago

separate hard & SVD singularities, new code in GSASIImath_new.py not GSASIImath.py

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