source: trunk/GSASIImath.py @ 4514

Last change on this file since 4514 was 4514, checked in by vondreele, 17 months ago

remove non gpx entries from recent projects
continue work on magnetism

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