source: trunk/GSASIImath_new.py @ 4693

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

updates to math_new after testing; improve output

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