source: trunk/GSASIImath.py @ 4213

Last change on this file since 4213 was 4213, checked in by toby, 2 years ago

multiple small changes to allow docs build under 3.x

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