source: trunk/GSASIImath.py @ 4295

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

fix crash on zero cycles (e.g.in a simulation)

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