source: trunk/GSASIImath.py @ 4583

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

revise graphics & input for Insert Rigid Body

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