source: trunk/GSASIImath_new.py @ 4690

Last change on this file since 4690 was 4690, checked in by toby, 2 years ago

major updates to trial version of HessianLSQ

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