source: trunk/GSASIImath.py @ 4513

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

fix problem after setting 1-->2 transformation. A block of code neededto be tabbed in once more.
better mag. str. fctr stuff - closer but not right.

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