source: trunk/GSASIImath.py @ 4761

Last change on this file since 4761 was 4761, checked in by toby, 9 months ago

refinements with no variables or no cycles

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