source: trunk/GSASIImath.py @ 4760

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

fix zero cycle bug

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