source: trunk/GSASIImath.py @ 4553

Last change on this file since 4553 was 4553, checked in by vondreele, 16 months ago

reimplement new PeaksUnique? code & make more general. Now finds unique peaks closest to x=0, y=0, z=0, origin or center.

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