source: trunk/GSASIImath_new.py @ 4759

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

wx4.1 & mpl3.3 fixes

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