source: trunk/GSASIImath.py @ 4770

Last change on this file since 4770 was 4770, checked in by vondreele, 9 months ago

reinstitute final call to Hess after convergence - required for proper esds

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