source: trunk/GSASIImath.py @ 4671

Last change on this file since 4671 was 4671, checked in by vondreele, 12 months ago

add refinable multiplier for a fixed background; multiplier is now normally > 0
usable for peak fitting and Rietved refinement
add the fixed background entry to all background defaults
fix bug in GetDetectorXY fo when cursor outside image - returns [0,0] not None; changes elsewhere to use this
GetTthAzmDsp? now returns explicit list not assumed tuple
put the abs in the nl.qr test for singularities in HessianLSQ
Add 'BF mult' to name list in G2obj
put a try - except TypeError? around setting plot style stuff in PlotPatterns?
Remove the setting of Pattern[0]BackFile? - this was redundant for PWDR
remove picker/pickradius from linescan plot - failed
remove the alternate fixed background definition ('_fixedVary', etc.)
clear the PhaseReaderClass?Drawing? dictionary upon import of phase from a gpx file

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