source: trunk/GSASIImath.py @ 4533

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

fix math error in incomm mag str fctr stuff $ expand equations

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