source: trunk/GSASIImath.py @ 4688

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

correct setting of singular vars in old HessianLSQ

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