source: trunk/GSASIImath.py @ 4421

Last change on this file since 4421 was 4421, checked in by toby, 19 months ago

add R.B. extractor; misc doc updates; plt bug after strip Hs; grid bug: uncompleted edits; mac horz. line fix

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