source: trunk/GSASIImath.py @ 4291

Last change on this file since 4291 was 4291, checked in by vondreele, 3 years ago

new tool in Phase/General? - Compare: compares (P 1 only) structures to ideal octahedra & tetrahedra
gives statistics on bond length, polygon tilts & rotations as well as distortions.

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