source: trunk/GSASIImath.py @ 4176

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

best mag ss sf calc for DyMn6Ge6

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