source: trunk/GSASIImath.py @ 4182

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

mag. str. fctr. changes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 191.3 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2019-10-18 20:01:27 +0000 (Fri, 18 Oct 2019) $
5# $Author: vondreele $
6# $Revision: 4182 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 4182 2019-10-18 20:01:27Z 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: 4182 $")
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    if SGData['SGGray']:
1381        SGMT = np.vstack((SGMT,SGMT))
1382        Sinv =np.vstack((Sinv,Sinv))
1383        SGT = np.vstack((SGT,SGT+.5))%1.
1384    detSM = nl.det(SGMT)
1385    mst = Sinv[:,3,:3]
1386    epsinv = Sinv[:,3,3]
1387    phi = np.inner(XYZ,modQ).T
1388    TA = np.sum(mst[nxs,:,:]*(XYZ-SGT[:,:3][nxs,:,:]),axis=-1).T
1389    tauT =  TA[nxs,:,:] + epsinv[nxs,:,nxs]*(glTau[:,nxs,nxs]-SGT[:,3][nxs,:,nxs]+phi[nxs,:,:])
1390    modind = np.arange(nWaves)+1.
1391    phase = (modind[:,nxs,nxs]*tauT)     #Nops,Natm,Nwave
1392    psin = np.sin(twopi*phase)
1393    pcos = np.cos(twopi*phase)
1394    MmodA = np.sum(Am[nxs,nxs,:,:,:]*pcos[:,:,:,nxs,nxs],axis=3)    #cos term
1395    MmodB = np.sum(Bm[nxs,nxs,:,:,:]*psin[:,:,:,nxs,nxs],axis=3)    #sin term
1396#    if SGData['SGGray']:
1397#        MmodA = -np.sum(SGMT[nxs,:,nxs,:,:]*MmodA[:,:,:,nxs,:],axis=-1)*detSM[nxs,:,nxs,nxs]
1398#        MmodB = -np.sum(SGMT[nxs,:,nxs,:,:]*MmodB[:,:,:,nxs,:],axis=-1)*detSM[nxs,:,nxs,nxs]
1399#    else:
1400    MmodA = np.sum(SGMT[nxs,:,nxs,:,:]*MmodA[:,:,:,nxs,:],axis=-1)*SGData['MagMom'][nxs,:,nxs,nxs]
1401    MmodB = np.sum(SGMT[nxs,:,nxs,:,:]*MmodB[:,:,:,nxs,:],axis=-1)*SGData['MagMom'][nxs,:,nxs,nxs]
1402    return MmodB,MmodA    #Ntau,Nops,Natm,,Mxyz; cos & sin parts; sum matches drawn atom moments
1403       
1404def Modulation(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
1405    '''
1406    H: array nRefBlk x ops X hklt
1407    HP: array nRefBlk x ops X hklt proj to hkl
1408    nWaves: list number of waves for frac, pos, uij & mag
1409    Fmod: array 2 x atoms x waves    (sin,cos terms)
1410    Xmod: array atoms X 3 X ngl
1411    Umod: array atoms x 3x3 x ngl
1412    glTau,glWt: arrays Gauss-Lorentzian pos & wts
1413    '''
1414   
1415    if nWaves[2]:       #uij (adp) waves
1416        if len(HP.shape) > 2:
1417            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1418        else:
1419            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1420    else:
1421        HbH = 1.0
1422    HdotX = np.inner(HP,Xmod)                   #refBlk x ops x atoms X ngl
1423    if len(H.shape) > 2:
1424        D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
1425        HdotXD = twopi*(HdotX+D[:,:,nxs,:])
1426    else:
1427        D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
1428        HdotXD = twopi*(HdotX+D[:,nxs,:])
1429    cosHA = np.sum(Fmod*HbH*np.cos(HdotXD)*glWt,axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
1430    sinHA = np.sum(Fmod*HbH*np.sin(HdotXD)*glWt,axis=-1)       #imag part; ditto
1431    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
1432   
1433#def MagModulation(H,HP,nWaves,Fmod,Xmod,Umod,Mmod,glTau,glWt):
1434#    '''
1435#    H: array nRefBlk x ops X hklt
1436#    HP: array nRefBlk x ops X hklt proj to hkl
1437#    nWaves: list number of waves for frac, pos, uij & mag
1438#    Fmod: array 2 x atoms x waves    (sin,cos terms)
1439#    Xmod: array atoms X 3 X ngl
1440#    Umod: array atoms x 3x3 x ngl
1441#    Mmod: array atoms x 3x3 x ngl
1442#    glTau,glWt: arrays Gauss-Lorentzian pos & wts
1443#    '''
1444#   
1445#    if nWaves[2]:       #uij (adp) waves
1446#        if len(HP.shape) > 2:
1447#            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1448#        else:
1449#            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1450#    else:
1451#        HbH = 1.0
1452#    HdotX = np.inner(HP,Xmod)                   #refBlk x ops x atoms X ngl
1453#    if len(H.shape) > 2:
1454#        D = H[:,:,3:]*glTau[nxs,nxs,:]              #m*e*tau; refBlk x ops X ngl
1455#        HdotXD = twopi*(HdotX+D[:,:,nxs,:])
1456#    else:
1457#        D = H[:,3:]*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
1458#        HdotXD = twopi*(HdotX+D[:,nxs,:])
1459#    M = np.swapaxes(Mmod,1,2)
1460#    cosHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.cos(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
1461#    sinHA = np.sum(M[nxs,nxs,:,:,:]*(Fmod*HbH*np.sin(HdotXD)[:,:,:,nxs,:]*glWt),axis=-1)       #imag part; ditto
1462#    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
1463#
1464def ModulationTw(H,HP,nWaves,Fmod,Xmod,Umod,glTau,glWt):
1465    '''
1466    H: array nRefBlk x tw x ops X hklt
1467    HP: array nRefBlk x tw x ops X hklt proj to hkl
1468    Fmod: array 2 x atoms x waves    (sin,cos terms)
1469    Xmod: array atoms X ngl X 3
1470    Umod: array atoms x ngl x 3x3
1471    glTau,glWt: arrays Gauss-Lorentzian pos & wts
1472    '''
1473   
1474    if nWaves[2]:
1475        if len(HP.shape) > 3:   #Blocks of reflections
1476            HbH = np.exp(-np.sum(HP[:,:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1477        else:   #single reflections
1478            HbH = np.exp(-np.sum(HP[:,nxs,nxs,:]*np.inner(HP,Umod),axis=-1)) # refBlk x ops x atoms x ngl add Overhauser corr.?
1479    else:
1480        HbH = 1.0
1481    HdotX = np.inner(HP,Xmod)                   #refBlk x tw x ops x atoms X ngl
1482    if len(H.shape) > 3:
1483        D = glTau*H[:,:,:,3:,nxs]              #m*e*tau; refBlk x tw x ops X ngl
1484        HdotXD = twopi*(HdotX+D[:,:,:,nxs,:])
1485    else:
1486        D = H*glTau[nxs,:]              #m*e*tau; refBlk x ops X ngl
1487        HdotXD = twopi*(HdotX+D[:,nxs,:])
1488    cosHA = np.sum(Fmod*HbH*np.cos(HdotXD)*glWt,axis=-1)       #real part; refBlk X ops x atoms; sum for G-L integration
1489    sinHA = np.sum(Fmod*HbH*np.sin(HdotXD)*glWt,axis=-1)       #imag part; ditto
1490    return np.array([cosHA,sinHA])             # 2 x refBlk x SGops x atoms
1491   
1492def makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast):
1493    '''
1494    Only for Fourier waves for fraction, position & adp (probably not used for magnetism)
1495    FSSdata: array 2 x atoms x waves    (sin,cos terms)
1496    XSSdata: array 2x3 x atoms X waves (sin,cos terms)
1497    USSdata: array 2x6 x atoms X waves (sin,cos terms)
1498    Mast: array orthogonalization matrix for Uij
1499    '''
1500    glTau,glWt = pwd.pygauleg(0.,1.,ngl)         #get Gauss-Legendre intervals & weights
1501    waveShapes = [FSSdata.T.shape,XSSdata.T.shape,USSdata.T.shape]
1502    Af = np.array(FSSdata[0]).T    #sin frac mods x waves x atoms
1503    Bf = np.array(FSSdata[1]).T    #cos frac mods...
1504    Ax = np.array(XSSdata[:3]).T   #...cos pos mods x waves x atoms
1505    Bx = np.array(XSSdata[3:]).T   #...cos pos mods
1506    Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T   #atoms x waves x sin Uij mods
1507    Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T   #...cos Uij mods
1508    nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]] 
1509    StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))    #atoms x waves x 3 x ngl
1510    CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl))
1511    ZtauXt = np.zeros((Ax.shape[0],2,3,ngl))               #atoms x Tminmax x 3 x ngl
1512    ZtauXx = np.zeros((Ax.shape[0],3,ngl))               #atoms x XYZmax x ngl
1513    for iatm in range(Ax.shape[0]):
1514        nx = 0
1515        if 'ZigZag' in waveTypes[iatm]:
1516            nx = 1
1517        elif 'Block' in waveTypes[iatm]:
1518            nx = 1
1519        tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau  #Xwaves x ngl
1520        if nx:   
1521            StauX[iatm][nx:] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X ngl
1522            CtauX[iatm][nx:] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
1523        else:
1524            StauX[iatm] = np.ones_like(Ax)[iatm,:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:]   #atoms X waves X 3(xyz) X ngl
1525            CtauX[iatm] = np.ones_like(Bx)[iatm,:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:]   #ditto
1526    if nWaves[0]:
1527        tauF = np.arange(1.,nWaves[0]+1)[:,nxs]*glTau  #Fwaves x ngl
1528        StauF = np.ones_like(Af)[:,:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf
1529        CtauF = np.ones_like(Bf)[:,:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf
1530    else:
1531        StauF = 1.0
1532        CtauF = 1.0
1533    if nWaves[2]:
1534        tauU = np.arange(1.,nWaves[2]+1)[:,nxs]*glTau     #Uwaves x ngl
1535        StauU = np.ones_like(Au)[:,:,:,:,nxs]*np.sin(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmodA/dAu
1536        CtauU = np.ones_like(Bu)[:,:,:,:,nxs]*np.cos(twopi*tauU)[nxs,:,nxs,nxs,:]   #also dUmodB/dBu
1537        UmodA = Au[:,:,:,:,nxs]*StauU #atoms x waves x 3x3 x ngl
1538        UmodB = Bu[:,:,:,:,nxs]*CtauU #ditto
1539#derivs need to be ops x atoms x waves x 6uij; ops x atoms x waves x ngl x 6uij before sum
1540        StauU = np.rollaxis(np.rollaxis(np.swapaxes(StauU,2,4),-1),-1)
1541        CtauU = np.rollaxis(np.rollaxis(np.swapaxes(CtauU,2,4),-1),-1)
1542    else:
1543        StauU = 1.0
1544        CtauU = 1.0
1545        UmodA = 0.
1546        UmodB = 0.
1547    return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB
1548   
1549def ModulationDerv(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt):
1550    '''
1551    Compute Fourier modulation derivatives
1552    H: array ops X hklt proj to hkl
1553    HP: array ops X hklt proj to hkl
1554    Hij: array 2pi^2[a*^2h^2 b*^2k^2 c*^2l^2 a*b*hk a*c*hl b*c*kl] of projected hklm to hkl space
1555    '''
1556   
1557    Mf = [H.shape[0],]+list(waveShapes[0])    #=[ops,atoms,waves,2] (sin+cos frac mods)
1558    dGdMfC = np.zeros(Mf)
1559    dGdMfS = np.zeros(Mf)
1560    Mx = [H.shape[0],]+list(waveShapes[1])   #=[ops,atoms,waves,6] (sin+cos pos mods)
1561    dGdMxC = np.zeros(Mx)
1562    dGdMxS = np.zeros(Mx)
1563    Mu = [H.shape[0],]+list(waveShapes[2])    #=[ops,atoms,waves,12] (sin+cos Uij mods)
1564    dGdMuC = np.zeros(Mu)
1565    dGdMuS = np.zeros(Mu)
1566   
1567    D = twopi*H[:,3][:,nxs]*glTau[nxs,:]              #m*e*tau; ops X ngl
1568    HdotX = twopi*np.inner(HP,Xmod)        #ops x atoms X ngl
1569    HdotXD = HdotX+D[:,nxs,:]
1570    if nWaves[2]:
1571        Umod = np.swapaxes((UmodAB),2,4)      #atoms x waves x ngl x 3x3 (symmetric so I can do this!)
1572        HuH = np.sum(HP[:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #ops x atoms x waves x ngl
1573        HuH = np.sum(HP[:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1)    #ops x atoms x waves x ngl
1574        HbH = np.exp(-np.sum(HuH,axis=-2)) # ops x atoms x ngl; sum waves - OK vs Modulation version
1575#        part1 = -np.exp(-HuH)*Fmod[nxs,:,nxs,:]    #ops x atoms x waves x ngl
1576        part1 = -np.exp(-HuH)*Fmod    #ops x atoms x waves x ngl
1577        dUdAu = Hij[:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[0]),0,4)[nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6sinUij
1578        dUdBu = Hij[:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[1]),0,4)[nxs,:,:,:,:]    #ops x atoms x waves x ngl x 6cosUij
1579        dGdMuCa = np.sum(part1[:,:,:,:,nxs]*dUdAu*np.cos(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
1580        dGdMuCb = np.sum(part1[:,:,:,:,nxs]*dUdBu*np.cos(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1581        dGdMuC = np.concatenate((dGdMuCa,dGdMuCb),axis=-1)   #ops x atoms x waves x 12uij
1582        dGdMuSa = np.sum(part1[:,:,:,:,nxs]*dUdAu*np.sin(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   #ops x atoms x waves x 6uij; G-L sum
1583        dGdMuSb = np.sum(part1[:,:,:,:,nxs]*dUdBu*np.sin(HdotXD)[:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1584        dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1)   #ops x atoms x waves x 12uij
1585    else:
1586        HbH = np.ones_like(HdotX)
1587    dHdXA = twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[0],-1,-2)[nxs,:,:,:,:]    #ops x atoms x sine waves x ngl x xyz
1588    dHdXB = twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[1],-1,-2)[nxs,:,:,:,:]    #ditto - cos waves
1589# ops x atoms x waves x 2xyz - real part - good
1590#    dGdMxCa = -np.sum((Fmod[nxs,:,:]*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1591#    dGdMxCb = -np.sum((Fmod[nxs,:,:]*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1592    dGdMxCa = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1593    dGdMxCb = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)
1594    dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1)
1595# ops x atoms x waves x 2xyz - imag part - good
1596#    dGdMxSa = np.sum((Fmod[nxs,:,:]*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
1597#    dGdMxSb = np.sum((Fmod[nxs,:,:]*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
1598    dGdMxSa = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
1599    dGdMxSb = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2)   
1600    dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1)
1601    return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS]
1602       
1603def posFourier(tau,psin,pcos):
1604    A = np.array([ps[:,nxs]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)])
1605    B = np.array([pc[:,nxs]*np.cos(2*np.pi*(i+1)*tau) for i,pc in enumerate(pcos)])
1606    return np.sum(A,axis=0)+np.sum(B,axis=0)
1607   
1608def posZigZag(T,Tmm,Xmax):
1609    DT = Tmm[1]-Tmm[0]
1610    Su = 2.*Xmax/DT
1611    Sd = 2.*Xmax/(1.-DT)
1612    A = np.array([np.where(  0.< (t-Tmm[0])%1. <= DT, -Xmax+Su*((t-Tmm[0])%1.), Xmax-Sd*((t-Tmm[1])%1.)) for t in T])
1613    return A
1614   
1615#def posZigZagDerv(T,Tmm,Xmax):
1616#    DT = Tmm[1]-Tmm[0]
1617#    Su = 2.*Xmax/DT
1618#    Sd = 2.*Xmax/(1.-DT)
1619#    dAdT = np.zeros((2,3,len(T)))
1620#    dAdT[0] = np.array([np.where(Tmm[0] < t <= Tmm[1],Su*(t-Tmm[0]-1)/DT,-Sd*(t-Tmm[1])/(1.-DT)) for t in T]).T
1621#    dAdT[1] = np.array([np.where(Tmm[0] < t <= Tmm[1],-Su*(t-Tmm[0])/DT,Sd*(t-Tmm[1])/(1.-DT)) for t in T]).T
1622#    dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-1.+2.*(t-Tmm[0])/DT,1.-2.*(t-Tmm[1])%1./DT) for t in T])
1623#    return dAdT,dAdX
1624
1625def posBlock(T,Tmm,Xmax):
1626    A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-Xmax,Xmax) for t in T])
1627    return A
1628   
1629#def posBlockDerv(T,Tmm,Xmax):
1630#    dAdT = np.zeros((2,3,len(T)))
1631#    ind = np.searchsorted(T,Tmm)
1632#    dAdT[0,:,ind[0]] = -Xmax/len(T)
1633#    dAdT[1,:,ind[1]] = Xmax/len(T)
1634#    dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t <= Tmm[1],-1.,1.) for t in T])  #OK
1635#    return dAdT,dAdX
1636   
1637def fracCrenel(tau,Toff,Twid):
1638    Tau = (tau-Toff)%1.
1639    A = np.where(Tau<Twid,1.,0.)
1640    return A
1641   
1642def fracFourier(tau,fsin,fcos):
1643    if len(fsin) == 1:
1644        A = np.array([fsin[0]*np.sin(2.*np.pi*tau)])
1645        B = np.array([fcos[0]*np.cos(2.*np.pi*tau)])
1646    else:
1647        A = np.array([fs[:,nxs]*np.sin(2.*np.pi*(i+1)*tau) for i,fs in enumerate(fsin)])
1648        B = np.array([fc[:,nxs]*np.cos(2.*np.pi*(i+1)*tau) for i,fc in enumerate(fcos)])
1649    return np.sum(A,axis=0)+np.sum(B,axis=0)
1650
1651def ApplyModulation(data,tau):
1652    '''Applies modulation to drawing atom positions & Uijs for given tau
1653    '''
1654    generalData = data['General']
1655    cell = generalData['Cell'][1:7]
1656    G,g = G2lat.cell2Gmat(cell)
1657    SGData = generalData['SGData']
1658    SSGData = generalData['SSGData']
1659    cx,ct,cs,cia = generalData['AtomPtrs']
1660    drawingData = data['Drawing']
1661    modul = generalData['SuperVec'][0]
1662    dcx,dct,dcs,dci = drawingData['atomPtrs']
1663    atoms = data['Atoms']
1664    drawAtoms = drawingData['Atoms']
1665    Fade = np.ones(len(drawAtoms))
1666    for atom in atoms:
1667        atxyz = np.array(atom[cx:cx+3])
1668        atuij = np.array(atom[cia+2:cia+8])
1669        Sfrac = atom[-1]['SS1']['Sfrac']
1670        Spos = atom[-1]['SS1']['Spos']
1671        Sadp = atom[-1]['SS1']['Sadp']
1672        if generalData['Type'] == 'magnetic':
1673            Smag = atom[-1]['SS1']['Smag']
1674            atmom = np.array(atom[cx+4:cx+7])
1675        indx = FindAtomIndexByIDs(drawAtoms,dci,[atom[cia+8],],True)
1676        for ind in indx:
1677            drawatom = drawAtoms[ind]
1678            opr = drawatom[dcs-1]
1679            sop,ssop,icent,cent,unit = G2spc.OpsfromStringOps(opr,SGData,SSGData)
1680            drxyz = (np.inner(sop[0],atxyz)+sop[1]+cent)*icent+np.array(unit)
1681            tauT = G2spc.getTauT(tau,sop,ssop,drxyz,modul)[-1]
1682            tauT *= icent       #invert wave on -1
1683#            print(tau,tauT,opr,G2spc.MT2text(sop).replace(' ',''),G2spc.SSMT2text(ssop).replace(' ',''))
1684            wave = np.zeros(3)
1685            uwave = np.zeros(6)
1686            mom = np.zeros(3)
1687            if len(Sfrac):
1688                scof = []
1689                ccof = []
1690                waveType = Sfrac[0]
1691                for i,sfrac in enumerate(Sfrac[1:]):
1692                    if not i and 'Crenel' in waveType:
1693                        Fade[ind] += fracCrenel(tauT,sfrac[0][0],sfrac[0][1])
1694                    else:
1695                        scof.append(sfrac[0][0])
1696                        ccof.append(sfrac[0][1])
1697                    if len(scof):
1698                        Fade[ind] += np.sum(fracFourier(tauT,scof,ccof))                           
1699            if len(Spos):
1700                scof = []
1701                ccof = []
1702                waveType = Spos[0]
1703                for i,spos in enumerate(Spos[1:]):
1704                    if waveType in ['ZigZag','Block'] and not i:
1705                        Tminmax = spos[0][:2]
1706                        XYZmax = np.array(spos[0][2:5])
1707                        if waveType == 'Block':
1708                            wave = np.array(posBlock([tauT,],Tminmax,XYZmax))[0]
1709                        elif waveType == 'ZigZag':
1710                            wave = np.array(posZigZag([tauT,],Tminmax,XYZmax))[0]
1711                    else:
1712                        scof.append(spos[0][:3])
1713                        ccof.append(spos[0][3:])
1714                if len(scof):
1715                    wave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
1716            if generalData['Type'] == 'magnetic' and len(Smag):
1717                scof = []
1718                ccof = []
1719                waveType = Smag[0]
1720                for i,spos in enumerate(Smag[1:]):
1721                    scof.append(spos[0][:3])
1722                    ccof.append(spos[0][3:])
1723                if len(scof):
1724                    mom += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1)
1725            if len(Sadp):
1726                scof = []
1727                ccof = []
1728                waveType = Sadp[0]
1729                for i,sadp in enumerate(Sadp[1:]):
1730                    scof.append(sadp[0][:6])
1731                    ccof.append(sadp[0][6:])
1732                ures = posFourier(tauT,np.array(scof),np.array(ccof))
1733                if np.any(ures):
1734                    uwave += np.sum(ures,axis=1)
1735            if atom[cia] == 'A':                   
1736                X,U = G2spc.ApplyStringOps(opr,SGData,atxyz+wave,atuij+uwave)
1737                drawatom[dcx:dcx+3] = X
1738                drawatom[dci-6:dci] = U
1739            else:
1740                X = G2spc.ApplyStringOps(opr,SGData,atxyz+wave)
1741                drawatom[dcx:dcx+3] = X
1742            if generalData['Type'] == 'magnetic':
1743                M = G2spc.ApplyStringOpsMom(opr,SGData,atmom+mom)
1744                drawatom[dcx+3:dcx+6] = M
1745    return drawAtoms,Fade
1746   
1747# gauleg.py Gauss Legendre numerical quadrature, x and w computation
1748# integrate from a to b using n evaluations of the function f(x) 
1749# usage: from gauleg import gaulegf         
1750#        x,w = gaulegf( a, b, n)                               
1751#        area = 0.0                                           
1752#        for i in range(1,n+1):          #  yes, 1..n                   
1753#          area += w[i]*f(x[i])                                   
1754
1755def gaulegf(a, b, n):
1756    x = range(n+1) # x[0] unused
1757    w = range(n+1) # w[0] unused
1758    eps = 3.0E-14
1759    m = (n+1)/2
1760    xm = 0.5*(b+a)
1761    xl = 0.5*(b-a)
1762    for i in range(1,m+1):
1763        z = math.cos(3.141592654*(i-0.25)/(n+0.5))
1764        while True:
1765            p1 = 1.0
1766            p2 = 0.0
1767            for j in range(1,n+1):
1768                p3 = p2
1769                p2 = p1
1770                p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
1771       
1772            pp = n*(z*p1-p2)/(z*z-1.0)
1773            z1 = z
1774            z = z1 - p1/pp
1775            if abs(z-z1) <= eps:
1776                break
1777
1778        x[i] = xm - xl*z
1779        x[n+1-i] = xm + xl*z
1780        w[i] = 2.0*xl/((1.0-z*z)*pp*pp)
1781        w[n+1-i] = w[i]
1782    return np.array(x), np.array(w)
1783# end gaulegf
1784   
1785   
1786def BessJn(nmax,x):
1787    ''' compute Bessel function J(n,x) from scipy routine & recurrance relation
1788    returns sequence of J(n,x) for n in range [-nmax...0...nmax]
1789   
1790    :param  integer nmax: maximul order for Jn(x)
1791    :param  float x: argument for Jn(x)
1792   
1793    :returns numpy array: [J(-nmax,x)...J(0,x)...J(nmax,x)]
1794   
1795    '''
1796    import scipy.special as sp
1797    bessJn = np.zeros(2*nmax+1)
1798    bessJn[nmax] = sp.j0(x)
1799    bessJn[nmax+1] = sp.j1(x)
1800    bessJn[nmax-1] = -bessJn[nmax+1]
1801    for i in range(2,nmax+1):
1802        bessJn[i+nmax] = 2*(i-1)*bessJn[nmax+i-1]/x-bessJn[nmax+i-2]
1803        bessJn[nmax-i] = bessJn[i+nmax]*(-1)**i
1804    return bessJn
1805   
1806def BessIn(nmax,x):
1807    ''' compute modified Bessel function I(n,x) from scipy routines & recurrance relation
1808    returns sequence of I(n,x) for n in range [-nmax...0...nmax]
1809   
1810    :param  integer nmax: maximul order for In(x)
1811    :param  float x: argument for In(x)
1812   
1813    :returns numpy array: [I(-nmax,x)...I(0,x)...I(nmax,x)]
1814   
1815    '''
1816    import scipy.special as sp
1817    bessIn = np.zeros(2*nmax+1)
1818    bessIn[nmax] = sp.i0(x)
1819    bessIn[nmax+1] = sp.i1(x)
1820    bessIn[nmax-1] = bessIn[nmax+1]
1821    for i in range(2,nmax+1):
1822        bessIn[i+nmax] = bessIn[nmax+i-2]-2*(i-1)*bessIn[nmax+i-1]/x
1823        bessIn[nmax-i] = bessIn[i+nmax]
1824    return bessIn
1825       
1826   
1827################################################################################
1828##### distance, angle, planes, torsion stuff
1829################################################################################
1830
1831def CalcDist(distance_dict, distance_atoms, parmDict):
1832    if not len(parmDict):
1833        return 0.
1834    pId = distance_dict['pId']
1835    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1836    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1837    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[0])] for x in ['x','y','z']]
1838    Txyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[1])] for x in ['x','y','z']]
1839    inv = 1
1840    symNo = distance_dict['symNo']
1841    if symNo < 0:
1842        inv = -1
1843        symNo *= -1
1844    cen = symNo//100
1845    op = symNo%100-1
1846    M,T = distance_dict['SGData']['SGOps'][op]
1847    D = T*inv+distance_dict['SGData']['SGCen'][cen]
1848    D += distance_dict['cellNo']
1849    Txyz = np.inner(M*inv,Txyz)+D
1850    dist = np.sqrt(np.sum(np.inner(Amat,(Txyz-Oxyz))**2))
1851#    GSASIIpath.IPyBreak()
1852    return dist   
1853   
1854def CalcDistDeriv(distance_dict, distance_atoms, parmDict):
1855    if not len(parmDict):
1856        return None
1857    pId = distance_dict['pId']
1858    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1859    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1860    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[0])] for x in ['x','y','z']]
1861    Txyz = [parmDict['%s::A%s:%d'%(pId,x,distance_atoms[1])] for x in ['x','y','z']]
1862    symNo = distance_dict['symNo']
1863    Tunit = distance_dict['cellNo']
1864    SGData = distance_dict['SGData']   
1865    deriv = getDistDerv(Oxyz,Txyz,Amat,Tunit,symNo,SGData)
1866    return deriv
1867   
1868def CalcAngle(angle_dict, angle_atoms, parmDict):
1869    if not len(parmDict):
1870        return 0.
1871    pId = angle_dict['pId']
1872    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1873    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1874    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[0])] for x in ['x','y','z']]
1875    Axyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][0])] for x in ['x','y','z']]
1876    Bxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][1])] for x in ['x','y','z']]
1877    ABxyz = [Axyz,Bxyz]
1878    symNo = angle_dict['symNo']
1879    vec = np.zeros((2,3))
1880    for i in range(2):
1881        inv = 1
1882        if symNo[i] < 0:
1883            inv = -1
1884        cen = inv*symNo[i]//100
1885        op = inv*symNo[i]%100-1
1886        M,T = angle_dict['SGData']['SGOps'][op]
1887        D = T*inv+angle_dict['SGData']['SGCen'][cen]
1888        D += angle_dict['cellNo'][i]
1889        ABxyz[i] = np.inner(M*inv,ABxyz[i])+D
1890        vec[i] = np.inner(Amat,(ABxyz[i]-Oxyz))
1891        dist = np.sqrt(np.sum(vec[i]**2))
1892        if not dist:
1893            return 0.
1894        vec[i] /= dist
1895    angle = acosd(np.sum(vec[0]*vec[1]))
1896#    GSASIIpath.IPyBreak()
1897    return angle
1898
1899def CalcAngleDeriv(angle_dict, angle_atoms, parmDict):
1900    if not len(parmDict):
1901        return None
1902    pId = angle_dict['pId']
1903    A = [parmDict['%s::A%d'%(pId,i)] for i in range(6)]
1904    Amat = G2lat.cell2AB(G2lat.A2cell(A))[0]
1905    Oxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[0])] for x in ['x','y','z']]
1906    Axyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][0])] for x in ['x','y','z']]
1907    Bxyz = [parmDict['%s::A%s:%d'%(pId,x,angle_atoms[1][1])] for x in ['x','y','z']]
1908    symNo = angle_dict['symNo']
1909    Tunit = angle_dict['cellNo']
1910    SGData = angle_dict['SGData']   
1911    deriv = getAngleDerv(Oxyz,Axyz,Bxyz,Amat,Tunit,symNo,SGData)
1912    return deriv
1913
1914def getSyXYZ(XYZ,ops,SGData):
1915    '''default doc
1916   
1917    :param type name: description
1918   
1919    :returns: type name: description
1920   
1921    '''
1922    XYZout = np.zeros_like(XYZ)
1923    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
1924        if op == '1':
1925            XYZout[i] = xyz
1926        else:
1927            oprs = op.split('+')
1928            unit = [0,0,0]
1929            if len(oprs)>1:
1930                unit = np.array(list(eval(oprs[1])))
1931            syop =int(oprs[0])
1932            inv = syop//abs(syop)
1933            syop *= inv
1934            cent = syop//100
1935            syop %= 100
1936            syop -= 1
1937            M,T = SGData['SGOps'][syop]
1938            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
1939    return XYZout
1940   
1941def getRestDist(XYZ,Amat):
1942    '''default doc string
1943   
1944    :param type name: description
1945   
1946    :returns: type name: description
1947   
1948    '''
1949    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
1950   
1951def getRestDeriv(Func,XYZ,Amat,ops,SGData):
1952    '''default doc string
1953   
1954    :param type name: description
1955   
1956    :returns: type name: description
1957   
1958    '''
1959    deriv = np.zeros((len(XYZ),3))
1960    dx = 0.00001
1961    for j,xyz in enumerate(XYZ):
1962        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
1963            XYZ[j] -= x
1964            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1965            XYZ[j] += 2*x
1966            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
1967            XYZ[j] -= x
1968            deriv[j][i] = (d1-d2)/(2*dx)
1969    return deriv.flatten()
1970
1971def getRestAngle(XYZ,Amat):
1972    '''default doc string
1973   
1974    :param type name: description
1975   
1976    :returns: type name: description
1977   
1978    '''
1979   
1980    def calcVec(Ox,Tx,Amat):
1981        return np.inner(Amat,(Tx-Ox))
1982
1983    VecA = calcVec(XYZ[1],XYZ[0],Amat)
1984    VecA /= np.sqrt(np.sum(VecA**2))
1985    VecB = calcVec(XYZ[1],XYZ[2],Amat)
1986    VecB /= np.sqrt(np.sum(VecB**2))
1987    edge = VecB-VecA
1988    edge = np.sum(edge**2)
1989    angle = (2.-edge)/2.
1990    angle = max(angle,-1.)
1991    return acosd(angle)
1992   
1993def getRestPlane(XYZ,Amat):
1994    '''default doc string
1995   
1996    :param type name: description
1997   
1998    :returns: type name: description
1999   
2000    '''
2001    sumXYZ = np.zeros(3)
2002    for xyz in XYZ:
2003        sumXYZ += xyz
2004    sumXYZ /= len(XYZ)
2005    XYZ = np.array(XYZ)-sumXYZ
2006    XYZ = np.inner(Amat,XYZ).T
2007    Zmat = np.zeros((3,3))
2008    for i,xyz in enumerate(XYZ):
2009        Zmat += np.outer(xyz.T,xyz)
2010    Evec,Emat = nl.eig(Zmat)
2011    Evec = np.sqrt(Evec)/(len(XYZ)-3)
2012    Order = np.argsort(Evec)
2013    return Evec[Order[0]]
2014   
2015def getRestChiral(XYZ,Amat):   
2016    '''default doc string
2017   
2018    :param type name: description
2019   
2020    :returns: type name: description
2021   
2022    '''
2023    VecA = np.empty((3,3))   
2024    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
2025    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
2026    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
2027    return nl.det(VecA)
2028   
2029def getRestTorsion(XYZ,Amat):
2030    '''default doc string
2031   
2032    :param type name: description
2033   
2034    :returns: type name: description
2035   
2036    '''
2037    VecA = np.empty((3,3))
2038    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
2039    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
2040    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
2041    D = nl.det(VecA)
2042    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
2043    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
2044    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
2045    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
2046    Ang = 1.0
2047    if abs(P12) < 1.0 and abs(P23) < 1.0:
2048        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
2049    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
2050    return TOR
2051   
2052def calcTorsionEnergy(TOR,Coeff=[]):
2053    '''default doc string
2054   
2055    :param type name: description
2056   
2057    :returns: type name: description
2058   
2059    '''
2060    sum = 0.
2061    if len(Coeff):
2062        cof = np.reshape(Coeff,(3,3)).T
2063        delt = TOR-cof[1]
2064        delt = np.where(delt<-180.,delt+360.,delt)
2065        delt = np.where(delt>180.,delt-360.,delt)
2066        term = -cof[2]*delt**2
2067        val = cof[0]*np.exp(term/1000.0)
2068        pMax = cof[0][np.argmin(val)]
2069        Eval = np.sum(val)
2070        sum = Eval-pMax
2071    return sum,Eval
2072
2073def getTorsionDeriv(XYZ,Amat,Coeff):
2074    '''default doc string
2075   
2076    :param type name: description
2077   
2078    :returns: type name: description
2079   
2080    '''
2081    deriv = np.zeros((len(XYZ),3))
2082    dx = 0.00001
2083    for j,xyz in enumerate(XYZ):
2084        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
2085            XYZ[j] -= x
2086            tor = getRestTorsion(XYZ,Amat)
2087            p1,d1 = calcTorsionEnergy(tor,Coeff)
2088            XYZ[j] += 2*x
2089            tor = getRestTorsion(XYZ,Amat)
2090            p2,d2 = calcTorsionEnergy(tor,Coeff)           
2091            XYZ[j] -= x
2092            deriv[j][i] = (p2-p1)/(2*dx)
2093    return deriv.flatten()
2094
2095def getRestRama(XYZ,Amat):
2096    '''Computes a pair of torsion angles in a 5 atom string
2097   
2098    :param nparray XYZ: crystallographic coordinates of 5 atoms
2099    :param nparray Amat: crystal to cartesian transformation matrix
2100   
2101    :returns: list (phi,psi) two torsion angles in degrees
2102   
2103    '''
2104    phi = getRestTorsion(XYZ[:5],Amat)
2105    psi = getRestTorsion(XYZ[1:],Amat)
2106    return phi,psi
2107   
2108def calcRamaEnergy(phi,psi,Coeff=[]):
2109    '''Computes pseudo potential energy from a pair of torsion angles and a
2110    numerical description of the potential energy surface. Used to create
2111    penalty function in LS refinement:     
2112    :math:`Eval(\\phi,\\psi) = C[0]*exp(-V/1000)`
2113
2114    where :math:`V = -C[3] * (\\phi-C[1])^2 - C[4]*(\\psi-C[2])^2 - 2*(\\phi-C[1])*(\\psi-C[2])`
2115   
2116    :param float phi: first torsion angle (:math:`\\phi`)
2117    :param float psi: second torsion angle (:math:`\\psi`)
2118    :param list Coeff: pseudo potential coefficients
2119   
2120    :returns: list (sum,Eval): pseudo-potential difference from minimum & value;
2121      sum is used for penalty function.
2122   
2123    '''
2124    sum = 0.
2125    Eval = 0.
2126    if len(Coeff):
2127        cof = Coeff.T
2128        dPhi = phi-cof[1]
2129        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
2130        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
2131        dPsi = psi-cof[2]
2132        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
2133        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
2134        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
2135        val = cof[0]*np.exp(val/1000.)
2136        pMax = cof[0][np.argmin(val)]
2137        Eval = np.sum(val)
2138        sum = Eval-pMax
2139    return sum,Eval
2140
2141def getRamaDeriv(XYZ,Amat,Coeff):
2142    '''Computes numerical derivatives of torsion angle pair pseudo potential
2143    with respect of crystallographic atom coordinates of the 5 atom sequence
2144   
2145    :param nparray XYZ: crystallographic coordinates of 5 atoms
2146    :param nparray Amat: crystal to cartesian transformation matrix
2147    :param list Coeff: pseudo potential coefficients
2148   
2149    :returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
2150     crystallographic xyz coordinates.
2151   
2152    '''
2153    deriv = np.zeros((len(XYZ),3))
2154    dx = 0.00001
2155    for j,xyz in enumerate(XYZ):
2156        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
2157            XYZ[j] -= x
2158            phi,psi = getRestRama(XYZ,Amat)
2159            p1,d1 = calcRamaEnergy(phi,psi,Coeff)
2160            XYZ[j] += 2*x
2161            phi,psi = getRestRama(XYZ,Amat)
2162            p2,d2 = calcRamaEnergy(phi,psi,Coeff)
2163            XYZ[j] -= x
2164            deriv[j][i] = (p2-p1)/(2*dx)
2165    return deriv.flatten()
2166
2167def getRestPolefig(ODFln,SamSym,Grid):
2168    '''default doc string
2169   
2170    :param type name: description
2171   
2172    :returns: type name: description
2173   
2174    '''
2175    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
2176    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
2177    R = np.where(R <= 1.,2.*atand(R),0.0)
2178    Z = np.zeros_like(R)
2179    Z = G2lat.polfcal(ODFln,SamSym,R,P)
2180    Z = np.reshape(Z,(Grid,Grid))
2181    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
2182
2183def getRestPolefigDerv(HKL,Grid,SHCoeff):
2184    '''default doc string
2185   
2186    :param type name: description
2187   
2188    :returns: type name: description
2189   
2190    '''
2191    pass
2192       
2193def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
2194    '''default doc string
2195   
2196    :param type name: description
2197   
2198    :returns: type name: description
2199   
2200    '''
2201    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
2202        TxT = inv*(np.inner(M,Tx)+T)+C+U
2203        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
2204       
2205    inv = Top/abs(Top)
2206    cent = abs(Top)//100
2207    op = abs(Top)%100-1
2208    M,T = SGData['SGOps'][op]
2209    C = SGData['SGCen'][cent]
2210    dx = .00001
2211    deriv = np.zeros(6)
2212    for i in [0,1,2]:
2213        Oxyz[i] -= dx
2214        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
2215        Oxyz[i] += 2*dx
2216        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
2217        Oxyz[i] -= dx
2218        Txyz[i] -= dx
2219        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
2220        Txyz[i] += 2*dx
2221        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
2222        Txyz[i] -= dx
2223    return deriv
2224   
2225def getAngleDerv(Oxyz,Axyz,Bxyz,Amat,Tunit,symNo,SGData):
2226   
2227    def calcAngle(Oxyz,ABxyz,Amat,Tunit,symNo,SGData):
2228        vec = np.zeros((2,3))
2229        for i in range(2):
2230            inv = 1
2231            if symNo[i] < 0:
2232                inv = -1
2233            cen = inv*symNo[i]//100
2234            op = inv*symNo[i]%100-1
2235            M,T = SGData['SGOps'][op]
2236            D = T*inv+SGData['SGCen'][cen]
2237            D += Tunit[i]
2238            ABxyz[i] = np.inner(M*inv,ABxyz[i])+D
2239            vec[i] = np.inner(Amat,(ABxyz[i]-Oxyz))
2240            dist = np.sqrt(np.sum(vec[i]**2))
2241            if not dist:
2242                return 0.
2243            vec[i] /= dist
2244        angle = acosd(np.sum(vec[0]*vec[1]))
2245    #    GSASIIpath.IPyBreak()
2246        return angle
2247       
2248    dx = .00001
2249    deriv = np.zeros(9)
2250    for i in [0,1,2]:
2251        Oxyz[i] -= dx
2252        a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2253        Oxyz[i] += 2*dx
2254        deriv[i] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2255        Oxyz[i] -= dx
2256        Axyz[i] -= dx
2257        a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2258        Axyz[i] += 2*dx
2259        deriv[i+3] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2260        Axyz[i] -= dx
2261        Bxyz[i] -= dx
2262        a0 = calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)
2263        Bxyz[i] += 2*dx
2264        deriv[i+6] = (calcAngle(Oxyz,[Axyz,Bxyz],Amat,Tunit,symNo,SGData)-a0)/(2.*dx)
2265        Bxyz[i] -= dx
2266    return deriv
2267   
2268def getAngSig(VA,VB,Amat,SGData,covData={}):
2269    '''default doc string
2270   
2271    :param type name: description
2272   
2273    :returns: type name: description
2274   
2275    '''
2276    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
2277        TxT = inv*(np.inner(M,Tx)+T)+C+U
2278        return np.inner(Amat,(TxT-Ox))
2279       
2280    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
2281        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
2282        VecA /= np.sqrt(np.sum(VecA**2))
2283        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
2284        VecB /= np.sqrt(np.sum(VecB**2))
2285        edge = VecB-VecA
2286        edge = np.sum(edge**2)
2287        angle = (2.-edge)/2.
2288        angle = max(angle,-1.)
2289        return acosd(angle)
2290       
2291    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
2292    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
2293    invA = invB = 1
2294    invA = TopA//abs(TopA)
2295    invB = TopB//abs(TopB)
2296    centA = abs(TopA)//100
2297    centB = abs(TopB)//100
2298    opA = abs(TopA)%100-1
2299    opB = abs(TopB)%100-1
2300    MA,TA = SGData['SGOps'][opA]
2301    MB,TB = SGData['SGOps'][opB]
2302    CA = SGData['SGCen'][centA]
2303    CB = SGData['SGCen'][centB]
2304    if 'covMatrix' in covData:
2305        covMatrix = covData['covMatrix']
2306        varyList = covData['varyList']
2307        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
2308        dx = .00001
2309        dadx = np.zeros(9)
2310        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2311        for i in [0,1,2]:
2312            OxA[i] -= dx
2313            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2314            OxA[i] += 2*dx
2315            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2316            OxA[i] -= dx
2317           
2318            TxA[i] -= dx
2319            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2320            TxA[i] += 2*dx
2321            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2322            TxA[i] -= dx
2323           
2324            TxB[i] -= dx
2325            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
2326            TxB[i] += 2*dx
2327            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
2328            TxB[i] -= dx
2329           
2330        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
2331        if sigAng < 0.01:
2332            sigAng = 0.0
2333        return Ang,sigAng
2334    else:
2335        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
2336
2337def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
2338    '''default doc string
2339   
2340    :param type name: description
2341   
2342    :returns: type name: description
2343   
2344    '''
2345    def calcDist(Atoms,SyOps,Amat):
2346        XYZ = []
2347        for i,atom in enumerate(Atoms):
2348            Inv,M,T,C,U = SyOps[i]
2349            XYZ.append(np.array(atom[1:4]))
2350            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2351            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2352        V1 = XYZ[1]-XYZ[0]
2353        return np.sqrt(np.sum(V1**2))
2354       
2355    SyOps = []
2356    names = []
2357    for i,atom in enumerate(Oatoms):
2358        names += atom[-1]
2359        Op,unit = Atoms[i][-1]
2360        inv = Op//abs(Op)
2361        m,t = SGData['SGOps'][abs(Op)%100-1]
2362        c = SGData['SGCen'][abs(Op)//100]
2363        SyOps.append([inv,m,t,c,unit])
2364    Dist = calcDist(Oatoms,SyOps,Amat)
2365   
2366    sig = -0.001
2367    if 'covMatrix' in covData:
2368        dx = .00001
2369        dadx = np.zeros(6)
2370        for i in range(6):
2371            ia = i//3
2372            ix = i%3
2373            Oatoms[ia][ix+1] += dx
2374            a0 = calcDist(Oatoms,SyOps,Amat)
2375            Oatoms[ia][ix+1] -= 2*dx
2376            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2377        covMatrix = covData['covMatrix']
2378        varyList = covData['varyList']
2379        DistVcov = getVCov(names,varyList,covMatrix)
2380        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
2381        if sig < 0.001:
2382            sig = -0.001
2383   
2384    return Dist,sig
2385
2386def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
2387    '''default doc string
2388   
2389    :param type name: description
2390   
2391    :returns: type name: description
2392   
2393    '''
2394
2395    def calcAngle(Atoms,SyOps,Amat):
2396        XYZ = []
2397        for i,atom in enumerate(Atoms):
2398            Inv,M,T,C,U = SyOps[i]
2399            XYZ.append(np.array(atom[1:4]))
2400            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2401            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2402        V1 = XYZ[1]-XYZ[0]
2403        V1 /= np.sqrt(np.sum(V1**2))
2404        V2 = XYZ[1]-XYZ[2]
2405        V2 /= np.sqrt(np.sum(V2**2))
2406        V3 = V2-V1
2407        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
2408        return acosd(cang)
2409
2410    SyOps = []
2411    names = []
2412    for i,atom in enumerate(Oatoms):
2413        names += atom[-1]
2414        Op,unit = Atoms[i][-1]
2415        inv = Op//abs(Op)
2416        m,t = SGData['SGOps'][abs(Op)%100-1]
2417        c = SGData['SGCen'][abs(Op)//100]
2418        SyOps.append([inv,m,t,c,unit])
2419    Angle = calcAngle(Oatoms,SyOps,Amat)
2420   
2421    sig = -0.01
2422    if 'covMatrix' in covData:
2423        dx = .00001
2424        dadx = np.zeros(9)
2425        for i in range(9):
2426            ia = i//3
2427            ix = i%3
2428            Oatoms[ia][ix+1] += dx
2429            a0 = calcAngle(Oatoms,SyOps,Amat)
2430            Oatoms[ia][ix+1] -= 2*dx
2431            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2432        covMatrix = covData['covMatrix']
2433        varyList = covData['varyList']
2434        AngVcov = getVCov(names,varyList,covMatrix)
2435        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
2436        if sig < 0.01:
2437            sig = -0.01
2438   
2439    return Angle,sig
2440
2441def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
2442    '''default doc string
2443   
2444    :param type name: description
2445   
2446    :returns: type name: description
2447   
2448    '''
2449
2450    def calcTorsion(Atoms,SyOps,Amat):
2451       
2452        XYZ = []
2453        for i,atom in enumerate(Atoms):
2454            Inv,M,T,C,U = SyOps[i]
2455            XYZ.append(np.array(atom[1:4]))
2456            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2457            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2458        V1 = XYZ[1]-XYZ[0]
2459        V2 = XYZ[2]-XYZ[1]
2460        V3 = XYZ[3]-XYZ[2]
2461        V1 /= np.sqrt(np.sum(V1**2))
2462        V2 /= np.sqrt(np.sum(V2**2))
2463        V3 /= np.sqrt(np.sum(V3**2))
2464        M = np.array([V1,V2,V3])
2465        D = nl.det(M)
2466        P12 = np.dot(V1,V2)
2467        P13 = np.dot(V1,V3)
2468        P23 = np.dot(V2,V3)
2469        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
2470        return Tors
2471           
2472    SyOps = []
2473    names = []
2474    for i,atom in enumerate(Oatoms):
2475        names += atom[-1]
2476        Op,unit = Atoms[i][-1]
2477        inv = Op//abs(Op)
2478        m,t = SGData['SGOps'][abs(Op)%100-1]
2479        c = SGData['SGCen'][abs(Op)//100]
2480        SyOps.append([inv,m,t,c,unit])
2481    Tors = calcTorsion(Oatoms,SyOps,Amat)
2482   
2483    sig = -0.01
2484    if 'covMatrix' in covData:
2485        dx = .00001
2486        dadx = np.zeros(12)
2487        for i in range(12):
2488            ia = i//3
2489            ix = i%3
2490            Oatoms[ia][ix+1] -= dx
2491            a0 = calcTorsion(Oatoms,SyOps,Amat)
2492            Oatoms[ia][ix+1] += 2*dx
2493            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2494            Oatoms[ia][ix+1] -= dx           
2495        covMatrix = covData['covMatrix']
2496        varyList = covData['varyList']
2497        TorVcov = getVCov(names,varyList,covMatrix)
2498        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
2499        if sig < 0.01:
2500            sig = -0.01
2501   
2502    return Tors,sig
2503       
2504def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
2505    '''default doc string
2506   
2507    :param type name: description
2508   
2509    :returns: type name: description
2510   
2511    '''
2512
2513    def calcDist(Atoms,SyOps,Amat):
2514        XYZ = []
2515        for i,atom in enumerate(Atoms):
2516            Inv,M,T,C,U = SyOps[i]
2517            XYZ.append(np.array(atom[1:4]))
2518            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2519            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2520        V1 = XYZ[1]-XYZ[0]
2521        return np.sqrt(np.sum(V1**2))
2522       
2523    def calcAngle(Atoms,SyOps,Amat):
2524        XYZ = []
2525        for i,atom in enumerate(Atoms):
2526            Inv,M,T,C,U = SyOps[i]
2527            XYZ.append(np.array(atom[1:4]))
2528            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2529            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2530        V1 = XYZ[1]-XYZ[0]
2531        V1 /= np.sqrt(np.sum(V1**2))
2532        V2 = XYZ[1]-XYZ[2]
2533        V2 /= np.sqrt(np.sum(V2**2))
2534        V3 = V2-V1
2535        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
2536        return acosd(cang)
2537
2538    def calcTorsion(Atoms,SyOps,Amat):
2539       
2540        XYZ = []
2541        for i,atom in enumerate(Atoms):
2542            Inv,M,T,C,U = SyOps[i]
2543            XYZ.append(np.array(atom[1:4]))
2544            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
2545            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
2546        V1 = XYZ[1]-XYZ[0]
2547        V2 = XYZ[2]-XYZ[1]
2548        V3 = XYZ[3]-XYZ[2]
2549        V1 /= np.sqrt(np.sum(V1**2))
2550        V2 /= np.sqrt(np.sum(V2**2))
2551        V3 /= np.sqrt(np.sum(V3**2))
2552        M = np.array([V1,V2,V3])
2553        D = nl.det(M)
2554        P12 = np.dot(V1,V2)
2555        P13 = np.dot(V1,V3)
2556        P23 = np.dot(V2,V3)
2557        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
2558        return Tors
2559           
2560    SyOps = []
2561    names = []
2562    for i,atom in enumerate(Oatoms):
2563        names += atom[-1]
2564        Op,unit = Atoms[i][-1]
2565        inv = Op//abs(Op)
2566        m,t = SGData['SGOps'][abs(Op)%100-1]
2567        c = SGData['SGCen'][abs(Op)//100]
2568        SyOps.append([inv,m,t,c,unit])
2569    M = len(Oatoms)
2570    if M == 2:
2571        Val = calcDist(Oatoms,SyOps,Amat)
2572    elif M == 3:
2573        Val = calcAngle(Oatoms,SyOps,Amat)
2574    else:
2575        Val = calcTorsion(Oatoms,SyOps,Amat)
2576   
2577    sigVals = [-0.001,-0.01,-0.01]
2578    sig = sigVals[M-3]
2579    if 'covMatrix' in covData:
2580        dx = .00001
2581        N = M*3
2582        dadx = np.zeros(N)
2583        for i in range(N):
2584            ia = i//3
2585            ix = i%3
2586            Oatoms[ia][ix+1] += dx
2587            if M == 2:
2588                a0 = calcDist(Oatoms,SyOps,Amat)
2589            elif M == 3:
2590                a0 = calcAngle(Oatoms,SyOps,Amat)
2591            else:
2592                a0 = calcTorsion(Oatoms,SyOps,Amat)
2593            Oatoms[ia][ix+1] -= 2*dx
2594            if M == 2:
2595                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
2596            elif M == 3:
2597                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
2598            else:
2599                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
2600        covMatrix = covData['covMatrix']
2601        varyList = covData['varyList']
2602        Vcov = getVCov(names,varyList,covMatrix)
2603        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
2604        if sig < sigVals[M-3]:
2605            sig = sigVals[M-3]
2606   
2607    return Val,sig
2608       
2609def ValEsd(value,esd=0,nTZ=False):
2610    '''Format a floating point number with a given level of precision or
2611    with in crystallographic format with a "esd", as value(esd). If esd is
2612    negative the number is formatted with the level of significant figures
2613    appropriate if abs(esd) were the esd, but the esd is not included.
2614    if the esd is zero, approximately 6 significant figures are printed.
2615    nTZ=True causes "extra" zeros to be removed after the decimal place.
2616    for example:
2617
2618      * "1.235(3)" for value=1.2346 & esd=0.003
2619      * "1.235(3)e4" for value=12346. & esd=30
2620      * "1.235(3)e6" for value=0.12346e7 & esd=3000
2621      * "1.235" for value=1.2346 & esd=-0.003
2622      * "1.240" for value=1.2395 & esd=-0.003
2623      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
2624      * "1.23460" for value=1.2346 & esd=0.0
2625
2626    :param float value: number to be formatted
2627    :param float esd: uncertainty or if esd < 0, specifies level of
2628      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
2629    :param bool nTZ: True to remove trailing zeros (default is False)
2630    :returns: value(esd) or value as a string
2631
2632    '''
2633    # Note: this routine is Python 3 compatible -- I think
2634    cutoff = 3.16228    #=(sqrt(10); same as old GSAS   was 1.95
2635    if math.isnan(value): # invalid value, bail out
2636        return '?'
2637    if math.isnan(esd): # invalid esd, treat as zero
2638        esd = 0
2639        esdoff = 5
2640#    if esd < 1.e-5:
2641#        esd = 0
2642#        esdoff = 5
2643    elif esd != 0:
2644        # transform the esd to a one or two digit integer
2645        l = math.log10(abs(esd)) % 1.
2646        if l < math.log10(cutoff): l+= 1.       
2647        intesd = int(round(10**l)) # esd as integer
2648        # determine the number of digits offset for the esd
2649        esdoff = int(round(math.log10(intesd*1./abs(esd))))
2650    else:
2651        esdoff = 5
2652    valoff = 0
2653    if abs(value) < abs(esdoff): # value is effectively zero
2654        pass
2655    elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
2656        # where the digit offset is to the left of the decimal place or where too many
2657        # digits are needed
2658        if abs(value) > 1:
2659            valoff = int(math.log10(abs(value)))
2660        elif abs(value) > 0:
2661            valoff = int(math.log10(abs(value))-0.9999999)
2662        else:
2663            valoff = 0
2664    if esd != 0:
2665        if valoff+esdoff < 0:
2666            valoff = esdoff = 0
2667        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
2668    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
2669        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
2670    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
2671        if abs(value) > 0:           
2672            extra = -math.log10(abs(value))
2673        else:
2674            extra = 0
2675        if extra > 0: extra += 1
2676        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
2677    if esd > 0:
2678        out += ("({:d})").format(intesd)  # add the esd
2679    elif nTZ and '.' in out:
2680        out = out.rstrip('0')  # strip zeros to right of decimal
2681        out = out.rstrip('.')  # and decimal place when not needed
2682    if valoff != 0:
2683        out += ("e{:d}").format(valoff) # add an exponent, when needed
2684    return out
2685   
2686###############################################################################
2687##### Protein validation - "ERRATV2" analysis
2688###############################################################################
2689
2690def validProtein(Phase,old):
2691   
2692    def sumintact(intact):
2693        return {'CC':intact['CC'],'NN':intact['NN'],'OO':intact['OO'],
2694        'CN':(intact['CN']+intact['NC']),'CO':(intact['CO']+intact['OC']),
2695        'NO':(intact['NO']+intact['ON'])}
2696       
2697    resNames = ['ALA','ARG','ASN','ASP','CYS','GLN','GLU','GLY','HIS','ILE',
2698        'LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL','MSE']
2699# data from errat.f
2700    b1_old = np.array([ 
2701        [1154.343,  600.213, 1051.018, 1132.885,  960.738],
2702        [600.213, 1286.818, 1282.042,  957.156,  612.789],
2703        [1051.018, 1282.042, 3519.471,  991.974, 1226.491],
2704        [1132.885,  957.156,  991.974, 1798.672,  820.355],
2705        [960.738,  612.789, 1226.491,  820.355, 2428.966]
2706        ])
2707    avg_old = np.array([ 0.225, 0.281, 0.071, 0.237, 0.044])    #Table 1 3.5A Obsd. Fr. p 1513
2708# data taken from erratv2.ccp
2709    b1 = np.array([
2710          [5040.279078850848200,        3408.805141583649400,   4152.904423767300600,   4236.200004171890200,   5054.781210204625500], 
2711          [3408.805141583648900,        8491.906094010220800,   5958.881777877950300,   1521.387352718486200,   4304.078200827221700], 
2712          [4152.904423767301500,        5958.881777877952100,   7637.167089335050100,   6620.715738223072500,   5287.691183798410700], 
2713          [4236.200004171890200,        1521.387352718486200,   6620.715738223072500,   18368.343774298410000,  4050.797811118806700], 
2714          [5054.781210204625500,        4304.078200827220800,   5287.691183798409800,   4050.797811118806700,   6666.856740479164700]])
2715    avg = np.array([0.192765509919262, 0.195575208778518, 0.275322406824210, 0.059102357035642, 0.233154192767480])
2716    General = Phase['General']
2717    Amat,Bmat = G2lat.cell2AB(General['Cell'][1:7])
2718    cx,ct,cs,cia = General['AtomPtrs']
2719    Atoms = Phase['Atoms']
2720    cartAtoms = []
2721    xyzmin = 999.*np.ones(3)
2722    xyzmax = -999.*np.ones(3)
2723    #select residue atoms,S,Se --> O make cartesian
2724    for atom in Atoms:
2725        if atom[1] in resNames:
2726            cartAtoms.append(atom[:cx+3])
2727            if atom[4].strip() in ['S','Se']:
2728                if not old:
2729                    continue        #S,Se skipped for erratv2?
2730                cartAtoms[-1][3] = 'Os'
2731                cartAtoms[-1][4] = 'O'
2732            cartAtoms[-1][cx:cx+3] = np.inner(Amat,cartAtoms[-1][cx:cx+3])
2733            cartAtoms[-1].append(atom[cia+8])
2734    XYZ = np.array([atom[cx:cx+3] for atom in cartAtoms])
2735    xyzmin = np.array([np.min(XYZ.T[i]) for i in [0,1,2]])
2736    xyzmax = np.array([np.max(XYZ.T[i]) for i in [0,1,2]])
2737    nbox = list(np.array(np.ceil((xyzmax-xyzmin)/4.),dtype=int))+[15,]
2738    Boxes = np.zeros(nbox,dtype=int)
2739    iBox = np.array([np.trunc((XYZ.T[i]-xyzmin[i])/4.) for i in [0,1,2]],dtype=int).T
2740    for ib,box in enumerate(iBox):  #put in a try for too many atoms in box (IndexError)?
2741        try:
2742            Boxes[box[0],box[1],box[2],0] += 1
2743            Boxes[box[0],box[1],box[2],Boxes[box[0],box[1],box[2],0]] = ib
2744        except IndexError:
2745            G2fil.G2Print('Error: too many atoms in box' )
2746            continue
2747    #Box content checks with errat.f $ erratv2.cpp ibox1 arrays
2748    indices = (-1,0,1)
2749    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices]) 
2750    dsmax = 3.75**2
2751    if old:
2752        dsmax = 3.5**2
2753    chains = []
2754    resIntAct = []
2755    chainIntAct = []
2756    res = []
2757    resNames = []
2758    resIDs = {}
2759    resname = []
2760    resID = {}
2761    newChain = True
2762    intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2763    for ia,atom in enumerate(cartAtoms):
2764        jntact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2765        if atom[2] not in chains:   #get chain id & save residue sequence from last chain
2766            chains.append(atom[2])
2767            if len(resIntAct):
2768                resIntAct.append(sumintact(intact))
2769                chainIntAct.append(resIntAct)
2770                resNames += resname
2771                resIDs.update(resID)
2772                res = []
2773                resname = []
2774                resID = {}
2775                resIntAct = []
2776                intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2777                newChain = True
2778        if atom[0] not in res:  #new residue, get residue no.
2779            if res and int(res[-1]) != int(atom[0])-1:  #a gap in chain - not new chain
2780                intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2781                ires = int(res[-1])
2782                for i in range(int(atom[0])-ires-1):
2783                    res.append(str(ires+i+1))
2784                    resname.append('')
2785                    resIntAct.append(sumintact(intact))
2786            res.append(atom[0])
2787            name = '%s-%s%s'%(atom[2],atom[0],atom[1])
2788            resname.append(name)
2789            resID[name] = atom[-1]
2790            if not newChain:
2791                resIntAct.append(sumintact(intact))
2792            intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2793            newChain = False
2794        ibox = iBox[ia]         #box location of atom
2795        tgts = []
2796        for unit in Units:      #assemble list of all possible target atoms
2797            jbox = ibox+unit
2798            if np.all(jbox>=0) and np.all((jbox-nbox[:3])<0):               
2799                tgts += list(Boxes[jbox[0],jbox[1],jbox[2]])
2800        tgts = list(set(tgts))
2801        tgts = [tgt for tgt in tgts if atom[:3] != cartAtoms[tgt][:3]]    #exclude same residue
2802        tgts = [tgt for tgt in tgts if np.sum((XYZ[ia]-XYZ[tgt])**2) < dsmax]
2803        ires = int(atom[0])
2804        if old:
2805            if atom[3].strip() == 'C':
2806                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) in [ires-1,ires+1])]
2807            elif atom[3].strip() == 'N':
2808                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() in ['C','CA'] and int(cartAtoms[tgt][0]) in [ires-1,ires+1])]
2809            elif atom[3].strip() == 'CA':
2810                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) in [ires-1,ires+1])]
2811        else:
2812            tgts = [tgt for tgt in tgts if not int(cartAtoms[tgt][0]) in [ires+1,ires+2,ires+3,ires+4,ires+5,ires+6,ires+7,ires+8]]
2813            if atom[3].strip() == 'C':
2814                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'N' and int(cartAtoms[tgt][0]) == ires+1)]
2815            elif atom[3].strip() == 'N':
2816                tgts = [tgt for tgt in tgts if not (cartAtoms[tgt][3].strip() == 'C' and int(cartAtoms[tgt][0]) == ires-1)]
2817        for tgt in tgts:
2818            dsqt = np.sqrt(np.sum((XYZ[ia]-XYZ[tgt])**2))
2819            mult = 1.0
2820            if dsqt > 3.25 and not old:
2821                mult = 2.*(3.75-dsqt)
2822            intype = atom[4].strip()+cartAtoms[tgt][4].strip()
2823            if 'S' not in intype:
2824                intact[intype] += mult
2825                jntact[intype] += mult
2826#        print ia,atom[0]+atom[1]+atom[3],tgts,jntact['CC'],jntact['CN']+jntact['NC'],jntact['CO']+jntact['OC'],jntact['NN'],jntact['NO']+jntact['ON']
2827    resNames += resname
2828    resIDs.update(resID)
2829    resIntAct.append(sumintact(intact))
2830    chainIntAct.append(resIntAct)
2831    chainProb = []
2832    for ich,chn in enumerate(chains):
2833        IntAct = chainIntAct[ich]
2834        nRes = len(IntAct)
2835        Probs = [0.,0.,0.,0.]   #skip 1st 4 residues in chain
2836        for i in range(4,nRes-4):
2837            if resNames[i]:
2838                mtrx = np.zeros(5)
2839                summ = 0.
2840                for j in range(i-4,i+5):
2841                    summ += np.sum(np.array(list(IntAct[j].values())))
2842                    if old:
2843                        mtrx[0] += IntAct[j]['CC']
2844                        mtrx[1] += IntAct[j]['CO']
2845                        mtrx[2] += IntAct[j]['NN']
2846                        mtrx[3] += IntAct[j]['NO']
2847                        mtrx[4] += IntAct[j]['OO']
2848                    else:
2849                        mtrx[0] += IntAct[j]['CC']
2850                        mtrx[1] += IntAct[j]['CN']
2851                        mtrx[2] += IntAct[j]['CO']
2852                        mtrx[3] += IntAct[j]['NN']
2853                        mtrx[4] += IntAct[j]['NO']
2854                mtrx /= summ
2855    #            print i+1,mtrx*summ
2856                if old:
2857                    mtrx -= avg_old
2858                    prob = np.inner(np.inner(mtrx,b1_old),mtrx)
2859                else:
2860                    mtrx -= avg
2861                    prob = np.inner(np.inner(mtrx,b1),mtrx)
2862            else:       #skip the gaps
2863                prob = 0.0
2864            Probs.append(prob)
2865        Probs += 4*[0.,]        #skip last 4 residues in chain
2866        chainProb += Probs
2867    return resNames,chainProb,resIDs
2868   
2869################################################################################
2870##### Texture fitting stuff
2871################################################################################
2872
2873def FitTexture(General,Gangls,refData,keyList,pgbar):
2874    import pytexture as ptx
2875    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2876   
2877    def printSpHarm(textureData,SHtextureSig):
2878        print ('\n Spherical harmonics texture: Order:' + str(textureData['Order']))
2879        names = ['omega','chi','phi']
2880        namstr = '  names :'
2881        ptstr =  '  values:'
2882        sigstr = '  esds  :'
2883        for name in names:
2884            namstr += '%12s'%('Sample '+name)
2885            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2886            if 'Sample '+name in SHtextureSig:
2887                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2888            else:
2889                sigstr += 12*' '
2890        print (namstr)
2891        print (ptstr)
2892        print (sigstr)
2893        print ('\n Texture coefficients:')
2894        SHcoeff = textureData['SH Coeff'][1]
2895        SHkeys = list(SHcoeff.keys())
2896        nCoeff = len(SHcoeff)
2897        nBlock = nCoeff//10+1
2898        iBeg = 0
2899        iFin = min(iBeg+10,nCoeff)
2900        for block in range(nBlock):
2901            namstr = '  names :'
2902            ptstr =  '  values:'
2903            sigstr = '  esds  :'
2904            for name in SHkeys[iBeg:iFin]:
2905                if 'C' in name:
2906                    namstr += '%12s'%(name)
2907                    ptstr += '%12.3f'%(SHcoeff[name])
2908                    if name in SHtextureSig:
2909                        sigstr += '%12.3f'%(SHtextureSig[name])
2910                    else:
2911                        sigstr += 12*' '
2912            print (namstr)
2913            print (ptstr)
2914            print (sigstr)
2915            iBeg += 10
2916            iFin = min(iBeg+10,nCoeff)
2917           
2918    def Dict2Values(parmdict, varylist):
2919        '''Use before call to leastsq to setup list of values for the parameters
2920        in parmdict, as selected by key in varylist'''
2921        return [parmdict[key] for key in varylist] 
2922       
2923    def Values2Dict(parmdict, varylist, values):
2924        ''' Use after call to leastsq to update the parameter dictionary with
2925        values corresponding to keys in varylist'''
2926        parmdict.update(list(zip(varylist,values)))
2927       
2928    def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2929        parmDict.update(list(zip(varyList,values)))
2930        Mat = np.empty(0)
2931        sumObs = 0
2932        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
2933        for hist in Gangls.keys():
2934            Refs = refData[hist]
2935            Refs[:,5] = np.where(Refs[:,5]>0.,Refs[:,5],0.)
2936            wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2937#            wt = 1./np.max(Refs[:,4],.25)
2938            sumObs += np.sum(wt*Refs[:,5])
2939            Refs[:,6] = 1.
2940            H = Refs[:,:3]
2941            phi,beta = G2lat.CrsAng(H,cell,SGData)
2942            psi,gam,x,x = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2943            for item in parmDict:
2944                if 'C' in item:
2945                    L,M,N = eval(item.strip('C'))
2946                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2947                    Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
2948                    Lnorm = G2lat.Lnorm(L)
2949                    Refs[:,6] += parmDict[item]*Lnorm*Kcl*Ksl
2950            mat = wt*(Refs[:,5]-Refs[:,6])
2951            Mat = np.concatenate((Mat,mat))
2952        sumD = np.sum(np.abs(Mat))
2953        R = min(100.,100.*sumD/sumObs)
2954        pgbar.Update(R,newmsg='Residual = %5.2f'%(R))
2955        print (' Residual: %.3f%%'%(R))
2956        return Mat
2957       
2958    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2959        Mat = np.empty(0)
2960        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
2961        for hist in Gangls.keys():
2962            mat = np.zeros((len(varyList),len(refData[hist])))
2963            Refs = refData[hist]
2964            H = Refs[:,:3]
2965            wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2966#            wt = 1./np.max(Refs[:,4],.25)
2967            phi,beta = G2lat.CrsAng(H,cell,SGData)
2968            psi,gam,dPdA,dGdA = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2969            for j,item in enumerate(varyList):
2970                if 'C' in item:
2971                    L,M,N = eval(item.strip('C'))
2972                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2973                    Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
2974                    Lnorm = G2lat.Lnorm(L)
2975                    mat[j] = -wt*Lnorm*Kcl*Ksl
2976                    for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
2977                        try:
2978                            l = varyList.index(itema)
2979                            mat[l] -= parmDict[item]*wt*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
2980                        except ValueError:
2981                            pass
2982            if len(Mat):
2983                Mat = np.concatenate((Mat,mat.T))
2984            else:
2985                Mat = mat.T
2986        print ('deriv')
2987        return Mat
2988
2989    print (' Fit texture for '+General['Name'])
2990    SGData = General['SGData']
2991    cell = General['Cell'][1:7]
2992    Texture = General['SH Texture']
2993    if not Texture['Order']:
2994        return 'No spherical harmonics coefficients'
2995    varyList = []
2996    parmDict = copy.copy(Texture['SH Coeff'][1])
2997    for item in ['Sample omega','Sample chi','Sample phi']:
2998        parmDict[item] = Texture[item][1]
2999        if Texture[item][0]:
3000            varyList.append(item)
3001    if Texture['SH Coeff'][0]:
3002        varyList += list(Texture['SH Coeff'][1].keys())
3003    while True:
3004        begin = time.time()
3005        values =  np.array(Dict2Values(parmDict, varyList))
3006        result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,ftol=1.e-6,
3007            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList,pgbar))
3008        ncyc = int(result[2]['nfev']//2)
3009        if ncyc:
3010            runtime = time.time()-begin   
3011            chisq = np.sum(result[2]['fvec']**2)
3012            Values2Dict(parmDict, varyList, result[0])
3013            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
3014            G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(result[2]['fvec']),len(varyList)))
3015            G2fil.G2Print ('refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc))
3016            try:
3017                sig = np.sqrt(np.diag(result[1])*GOF)
3018                if np.any(np.isnan(sig)):
3019                    G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***', mode='error')
3020                break                   #refinement succeeded - finish up!
3021            except ValueError:          #result[1] is None on singular matrix
3022                G2fil.G2Print ('**** Refinement failed - singular matrix ****', mode='error')
3023                return None
3024        else:
3025            break
3026   
3027    if ncyc:
3028        for parm in parmDict:
3029            if 'C' in parm:
3030                Texture['SH Coeff'][1][parm] = parmDict[parm]
3031            else:
3032                Texture[parm][1] = parmDict[parm] 
3033        sigDict = dict(zip(varyList,sig))
3034        printSpHarm(Texture,sigDict)
3035       
3036    return None
3037   
3038################################################################################
3039##### Fourier & charge flip stuff
3040################################################################################
3041
3042def adjHKLmax(SGData,Hmax):
3043    '''default doc string
3044   
3045    :param type name: description
3046   
3047    :returns: type name: description
3048   
3049    '''
3050    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
3051        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
3052        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
3053        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
3054    else:
3055        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
3056        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
3057        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
3058
3059def OmitMap(data,reflDict,pgbar=None):
3060    '''default doc string
3061   
3062    :param type name: description
3063   
3064    :returns: type name: description
3065   
3066    '''
3067    generalData = data['General']
3068    if not generalData['Map']['MapType']:
3069        G2fil.G2Print ('**** ERROR - Fourier map not defined')
3070        return
3071    mapData = generalData['Map']
3072    dmin = mapData['GridStep']*2.
3073    SGData = generalData['SGData']
3074    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3075    SGT = np.array([ops[1] for ops in SGData['SGOps']])
3076    cell = generalData['Cell'][1:8]       
3077    A = G2lat.cell2A(cell[:6])
3078    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3079    adjHKLmax(SGData,Hmax)
3080    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3081    time0 = time.time()
3082    for iref,ref in enumerate(reflDict['RefList']):
3083        if ref[4] >= dmin:
3084            Fosq,Fcsq,ph = ref[8:11]
3085            Uniq = np.inner(ref[:3],SGMT)
3086            Phi = np.inner(ref[:3],SGT)
3087            for i,hkl in enumerate(Uniq):        #uses uniq
3088                hkl = np.asarray(hkl,dtype='i')
3089                dp = 360.*Phi[i]                #and phi
3090                a = cosd(ph+dp)
3091                b = sind(ph+dp)
3092                phasep = complex(a,b)
3093                phasem = complex(a,-b)
3094                if '2Fo-Fc' in mapData['MapType']:
3095                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
3096                else:
3097                    F = np.sqrt(Fosq)
3098                h,k,l = hkl+Hmax
3099                Fhkl[h,k,l] = F*phasep
3100                h,k,l = -hkl+Hmax
3101                Fhkl[h,k,l] = F*phasem
3102    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
3103    M = np.mgrid[0:4,0:4,0:4]
3104    blkIds = np.array(list(zip(M[0].flatten(),M[1].flatten(),M[2].flatten())))
3105    iBeg = blkIds*rho0.shape//4
3106    iFin = (blkIds+1)*rho0.shape//4
3107    rho_omit = np.zeros_like(rho0)
3108    nBlk = 0
3109    for iB,iF in zip(iBeg,iFin):
3110        rho1 = np.copy(rho0)
3111        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
3112        Fnew = fft.ifftshift(fft.ifftn(rho1))
3113        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
3114        phase = Fnew/np.absolute(Fnew)
3115        OFhkl = np.absolute(Fhkl)*phase
3116        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
3117        rho_omit[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = np.copy(rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]])
3118        nBlk += 1
3119        pgbar.Update(nBlk)
3120    mapData['rho'] = np.real(rho_omit)/cell[6]
3121    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3122    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3123    G2fil.G2Print ('Omit map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape)))
3124    return mapData
3125   
3126def FourierMap(data,reflDict):
3127    '''default doc string
3128   
3129    :param type name: description
3130   
3131    :returns: type name: description
3132   
3133    '''
3134    generalData = data['General']
3135    mapData = generalData['Map']
3136    dmin = mapData['GridStep']*2.
3137    SGData = generalData['SGData']
3138    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3139    SGT = np.array([ops[1] for ops in SGData['SGOps']])
3140    cell = generalData['Cell'][1:8]       
3141    A = G2lat.cell2A(cell[:6])
3142    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3143    adjHKLmax(SGData,Hmax)
3144    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3145#    Fhkl[0,0,0] = generalData['F000X']
3146    time0 = time.time()
3147    for iref,ref in enumerate(reflDict['RefList']):
3148        if ref[4] > dmin:
3149            Fosq,Fcsq,ph = ref[8:11]
3150            Uniq = np.inner(ref[:3],SGMT)
3151            Phi = np.inner(ref[:3],SGT)
3152            for i,hkl in enumerate(Uniq):        #uses uniq
3153                hkl = np.asarray(hkl,dtype='i')
3154                dp = 360.*Phi[i]                #and phi
3155                a = cosd(ph+dp)
3156                b = sind(ph+dp)
3157                phasep = complex(a,b)
3158                phasem = complex(a,-b)
3159                if 'Fobs' in mapData['MapType']:
3160                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
3161                    h,k,l = hkl+Hmax
3162                    Fhkl[h,k,l] = F*phasep
3163                    h,k,l = -hkl+Hmax
3164                    Fhkl[h,k,l] = F*phasem
3165                elif 'Fcalc' in mapData['MapType']:
3166                    F = np.sqrt(Fcsq)
3167                    h,k,l = hkl+Hmax
3168                    Fhkl[h,k,l] = F*phasep
3169                    h,k,l = -hkl+Hmax
3170                    Fhkl[h,k,l] = F*phasem
3171                elif 'delt-F' in mapData['MapType']:
3172                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
3173                    h,k,l = hkl+Hmax
3174                    Fhkl[h,k,l] = dF*phasep
3175                    h,k,l = -hkl+Hmax
3176                    Fhkl[h,k,l] = dF*phasem
3177                elif '2*Fo-Fc' in mapData['MapType']:
3178                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
3179                    h,k,l = hkl+Hmax
3180                    Fhkl[h,k,l] = F*phasep
3181                    h,k,l = -hkl+Hmax
3182                    Fhkl[h,k,l] = F*phasem
3183                elif 'Patterson' in mapData['MapType']:
3184                    h,k,l = hkl+Hmax
3185                    Fhkl[h,k,l] = complex(Fosq,0.)
3186                    h,k,l = -hkl+Hmax
3187                    Fhkl[h,k,l] = complex(Fosq,0.)
3188    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
3189    G2fil.G2Print ('Fourier map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape)))
3190    mapData['Type'] = reflDict['Type']
3191    mapData['rho'] = np.real(rho)
3192    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3193    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3194   
3195def Fourier4DMap(data,reflDict):
3196    '''default doc string
3197   
3198    :param type name: description
3199   
3200    :returns: type name: description
3201   
3202    '''
3203    generalData = data['General']
3204    map4DData = generalData['4DmapData']
3205    mapData = generalData['Map']
3206    dmin = mapData['GridStep']*2.
3207    SGData = generalData['SGData']
3208    SSGData = generalData['SSGData']
3209    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3210    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3211    cell = generalData['Cell'][1:8]       
3212    A = G2lat.cell2A(cell[:6])
3213    maxM = 4
3214    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
3215    adjHKLmax(SGData,Hmax)
3216    Hmax = np.asarray(Hmax,dtype='i')+1
3217    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3218    time0 = time.time()
3219    for iref,ref in enumerate(reflDict['RefList']):
3220        if ref[5] > dmin:
3221            Fosq,Fcsq,ph = ref[9:12]
3222            Fosq = np.where(Fosq>0.,Fosq,0.)    #can't use Fo^2 < 0
3223            Uniq = np.inner(ref[:4],SSGMT)
3224            Phi = np.inner(ref[:4],SSGT)
3225            for i,hkl in enumerate(Uniq):        #uses uniq
3226                hkl = np.asarray(hkl,dtype='i')
3227                dp = 360.*Phi[i]                #and phi
3228                a = cosd(ph+dp)
3229                b = sind(ph+dp)
3230                phasep = complex(a,b)
3231                phasem = complex(a,-b)
3232                if 'Fobs' in mapData['MapType']:
3233                    F = np.sqrt(Fosq)
3234                    h,k,l,m = hkl+Hmax
3235                    Fhkl[h,k,l,m] = F*phasep
3236                    h,k,l,m = -hkl+Hmax
3237                    Fhkl[h,k,l,m] = F*phasem
3238                elif 'Fcalc' in mapData['MapType']:
3239                    F = np.sqrt(Fcsq)
3240                    h,k,l,m = hkl+Hmax
3241                    Fhkl[h,k,l,m] = F*phasep
3242                    h,k,l,m = -hkl+Hmax
3243                    Fhkl[h,k,l,m] = F*phasem                   
3244                elif 'delt-F' in mapData['MapType']:
3245                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
3246                    h,k,l,m = hkl+Hmax
3247                    Fhkl[h,k,l,m] = dF*phasep
3248                    h,k,l,m = -hkl+Hmax
3249                    Fhkl[h,k,l,m] = dF*phasem
3250    SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6]          #4D map
3251    rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6]    #3D map
3252    map4DData['rho'] = np.real(SSrho)
3253    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3254    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3255    map4DData['Type'] = reflDict['Type']
3256    mapData['Type'] = reflDict['Type']
3257    mapData['rho'] = np.real(rho)
3258    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3259    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3260    G2fil.G2Print ('Fourier map time: %.4f no. elements: %d dimensions: %s'%(time.time()-time0,Fhkl.size,str(Fhkl.shape)))
3261
3262# map printing for testing purposes
3263def printRho(SGLaue,rho,rhoMax):                         
3264    '''default doc string
3265   
3266    :param type name: description
3267   
3268    :returns: type name: description
3269   
3270    '''
3271    dim = len(rho.shape)
3272    if dim == 2:
3273        ix,jy = rho.shape
3274        for j in range(jy):
3275            line = ''
3276            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3277                line += (jy-j)*'  '
3278            for i in range(ix):
3279                r = int(100*rho[i,j]/rhoMax)
3280                line += '%4d'%(r)
3281            print (line+'\n')
3282    else:
3283        ix,jy,kz = rho.shape
3284        for k in range(kz):
3285            print ('k = %d'%k)
3286            for j in range(jy):
3287                line = ''
3288                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3289                    line += (jy-j)*'  '
3290                for i in range(ix):
3291                    r = int(100*rho[i,j,k]/rhoMax)
3292                    line += '%4d'%(r)
3293                print (line+'\n')
3294## keep this
3295               
3296def findOffset(SGData,A,Fhkl):   
3297    '''default doc string
3298   
3299    :param type name: description
3300   
3301    :returns: type name: description
3302   
3303    '''
3304    if SGData['SpGrp'] == 'P 1':
3305        return [0,0,0]   
3306    hklShape = Fhkl.shape
3307    hklHalf = np.array(hklShape)/2
3308    sortHKL = np.argsort(Fhkl.flatten())
3309    Fdict = {}
3310    for hkl in sortHKL:
3311        HKL = np.unravel_index(hkl,hklShape)
3312        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
3313        if F == 0.:
3314            break
3315        Fdict['%.6f'%(np.absolute(F))] = hkl
3316    Flist = np.flipud(np.sort(list(Fdict.keys())))
3317    F = str(1.e6)
3318    i = 0
3319    DH = []
3320    Dphi = []
3321    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3322    for F in Flist:
3323        hkl = np.unravel_index(Fdict[F],hklShape)
3324        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
3325            continue
3326        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
3327        Uniq = np.array(Uniq,dtype='i')
3328        Phi = np.array(Phi)
3329        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
3330        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
3331        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
3332        ang0 = np.angle(Fh0,deg=True)/360.
3333        for H,phi in list(zip(Uniq,Phi))[1:]:
3334            ang = (np.angle(Fhkl[int(H[0]),int(H[1]),int(H[2])],deg=True)/360.-phi)
3335            dH = H-hkl
3336            dang = ang-ang0
3337            DH.append(dH)
3338            Dphi.append((dang+.5) % 1.0)
3339        if i > 20 or len(DH) > 30:
3340            break
3341        i += 1
3342    DH = np.array(DH)
3343    G2fil.G2Print (' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist)))
3344    Dphi = np.array(Dphi)
3345    steps = np.array(hklShape)
3346    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
3347    XYZ = np.array(list(zip(X.flatten(),Y.flatten(),Z.flatten())))
3348    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
3349    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3350    hist,bins = np.histogram(Mmap,bins=1000)
3351#    for i,item in enumerate(hist[:10]):
3352#        print item,bins[i]
3353    chisq = np.min(Mmap)
3354    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3355    G2fil.G2Print (' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2]))
3356#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
3357    return DX
3358   
3359def ChargeFlip(data,reflDict,pgbar):
3360    '''default doc string
3361   
3362    :param type name: description
3363   
3364    :returns: type name: description
3365   
3366    '''
3367    generalData = data['General']
3368    mapData = generalData['Map']
3369    flipData = generalData['Flip']
3370    FFtable = {}
3371    if 'None' not in flipData['Norm element']:
3372        normElem = flipData['Norm element'].upper()
3373        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3374        for ff in FFs:
3375            if ff['Symbol'] == normElem:
3376                FFtable.update(ff)
3377    dmin = flipData['GridStep']*2.
3378    SGData = generalData['SGData']
3379    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3380    SGT = np.array([ops[1] for ops in SGData['SGOps']])
3381    cell = generalData['Cell'][1:8]       
3382    A = G2lat.cell2A(cell[:6])
3383    Vol = cell[6]
3384    im = 0
3385    if generalData['Modulated'] == True:
3386        im = 1
3387    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3388    adjHKLmax(SGData,Hmax)
3389    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3390    time0 = time.time()
3391    for iref,ref in enumerate(reflDict['RefList']):
3392        dsp = ref[4+im]
3393        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
3394            continue
3395        if dsp > dmin:
3396            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3397            if FFtable:
3398                SQ = 0.25/dsp**2
3399                ff *= G2el.ScatFac(FFtable,SQ)[0]
3400            if ref[8+im] > 0.:         #use only +ve Fobs**2
3401                E = np.sqrt(ref[8+im])/ff
3402            else:
3403                E = 0.
3404            ph = ref[10]
3405            ph = rn.uniform(0.,360.)
3406            Uniq = np.inner(ref[:3],SGMT)
3407            Phi = np.inner(ref[:3],SGT)
3408            for i,hkl in enumerate(Uniq):        #uses uniq
3409                hkl = np.asarray(hkl,dtype='i')
3410                dp = 360.*Phi[i]                #and phi
3411                a = cosd(ph+dp)
3412                b = sind(ph+dp)
3413                phasep = complex(a,b)
3414                phasem = complex(a,-b)
3415                h,k,l = hkl+Hmax
3416                Ehkl[h,k,l] = E*phasep
3417                h,k,l = -hkl+Hmax
3418                Ehkl[h,k,l] = E*phasem
3419#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3420    testHKL = np.array(flipData['testHKL'])+Hmax
3421    CEhkl = copy.copy(Ehkl)
3422    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3423    Emask = ma.getmask(MEhkl)
3424    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3425    Ncyc = 0
3426    old = np.seterr(all='raise')
3427    twophases = []
3428    while True:       
3429        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3430        CEsig = np.std(CErho)
3431        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3432        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3433        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3434        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3435        phase = CFhkl/np.absolute(CFhkl)
3436        twophases.append([np.angle(phase[h,k,l]) for h,k,l in testHKL])
3437        CEhkl = np.absolute(Ehkl)*phase
3438        Ncyc += 1
3439        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3440        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3441        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3442        if Rcf < 5.:
3443            break
3444        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3445        if not GoOn or Ncyc > 10000:
3446            break
3447    np.seterr(**old)
3448    G2fil.G2Print (' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size))
3449    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.  #? to get on same scale as e-map
3450    G2fil.G2Print (' No.cycles = %d Residual Rcf =%8.3f%s Map size: %s'%(Ncyc,Rcf,'%',str(CErho.shape)))
3451    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3452       
3453    mapData['Rcf'] = Rcf
3454    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3455    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3456    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3457    mapData['Type'] = reflDict['Type']
3458    return mapData,twophases
3459   
3460def findSSOffset(SGData,SSGData,A,Fhklm):   
3461    '''default doc string
3462   
3463    :param type name: description
3464   
3465    :returns: type name: description
3466   
3467    '''
3468    if SGData['SpGrp'] == 'P 1':
3469        return [0,0,0,0]   
3470    hklmShape = Fhklm.shape
3471    hklmHalf = np.array(hklmShape)/2
3472    sortHKLM = np.argsort(Fhklm.flatten())
3473    Fdict = {}
3474    for hklm in sortHKLM:
3475        HKLM = np.unravel_index(hklm,hklmShape)
3476        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
3477        if F == 0.:
3478            break
3479        Fdict['%.6f'%(np.absolute(F))] = hklm
3480    Flist = np.flipud(np.sort(list(Fdict.keys())))
3481    F = str(1.e6)
3482    i = 0
3483    DH = []
3484    Dphi = []
3485    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3486    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3487    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3488    for F in Flist:
3489        hklm = np.unravel_index(Fdict[F],hklmShape)
3490        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
3491            continue
3492        Uniq = np.inner(hklm-hklmHalf,SSGMT)
3493        Phi = np.inner(hklm-hklmHalf,SSGT)
3494        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
3495        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
3496        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
3497        ang0 = np.angle(Fh0,deg=True)/360.
3498        for H,phi in list(zip(Uniq,Phi))[1:]:
3499            H = np.array(H,dtype=int)
3500            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
3501            dH = H-hklm
3502            dang = ang-ang0
3503            DH.append(dH)
3504            Dphi.append((dang+.5) % 1.0)
3505        if i > 20 or len(DH) > 30:
3506            break
3507        i += 1
3508    DH = np.array(DH)
3509    G2fil.G2Print (' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist)))
3510    Dphi = np.array(Dphi)
3511    steps = np.array(hklmShape)
3512    X,Y,Z,T = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2],0:1:1./steps[3]]
3513    XYZT = np.array(list(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten())))
3514    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
3515    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3516    hist,bins = np.histogram(Mmap,bins=1000)
3517#    for i,item in enumerate(hist[:10]):
3518#        print item,bins[i]
3519    chisq = np.min(Mmap)
3520    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3521    G2fil.G2Print (' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3]))
3522#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
3523    return DX
3524   
3525def SSChargeFlip(data,reflDict,pgbar):
3526    '''default doc string
3527   
3528    :param type name: description
3529   
3530    :returns: type name: description
3531   
3532    '''
3533    generalData = data['General']
3534    mapData = generalData['Map']
3535    map4DData = {}
3536    flipData = generalData['Flip']
3537    FFtable = {}
3538    if 'None' not in flipData['Norm element']:
3539        normElem = flipData['Norm element'].upper()
3540        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3541        for ff in FFs:
3542            if ff['Symbol'] == normElem:
3543                FFtable.update(ff)
3544    dmin = flipData['GridStep']*2.
3545    SGData = generalData['SGData']
3546    SSGData = generalData['SSGData']
3547    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3548    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3549    cell = generalData['Cell'][1:8]       
3550    A = G2lat.cell2A(cell[:6])
3551    Vol = cell[6]
3552    maxM = 4
3553    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
3554    adjHKLmax(SGData,Hmax)
3555    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3556    time0 = time.time()
3557    for iref,ref in enumerate(reflDict['RefList']):
3558        dsp = ref[5]
3559        if dsp > dmin:
3560            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3561            if FFtable:
3562                SQ = 0.25/dsp**2
3563                ff *= G2el.ScatFac(FFtable,SQ)[0]
3564            if ref[9] > 0.:         #use only +ve Fobs**2
3565                E = np.sqrt(ref[9])/ff
3566            else:
3567                E = 0.
3568            ph = ref[11]
3569            ph = rn.uniform(0.,360.)
3570            Uniq = np.inner(ref[:4],SSGMT)
3571            Phi = np.inner(ref[:4],SSGT)
3572            for i,hklm in enumerate(Uniq):        #uses uniq
3573                hklm = np.asarray(hklm,dtype='i')
3574                dp = 360.*Phi[i]                #and phi
3575                a = cosd(ph+dp)
3576                b = sind(ph+dp)
3577                phasep = complex(a,b)
3578                phasem = complex(a,-b)
3579                h,k,l,m = hklm+Hmax
3580                Ehkl[h,k,l,m] = E*phasep
3581                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
3582                Ehkl[h,k,l,m] = E*phasem
3583#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3584    CEhkl = copy.copy(Ehkl)
3585    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3586    Emask = ma.getmask(MEhkl)
3587    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3588    Ncyc = 0
3589    old = np.seterr(all='raise')
3590    while True:       
3591        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3592        CEsig = np.std(CErho)
3593        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3594        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3595        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3596        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3597        phase = CFhkl/np.absolute(CFhkl)
3598        CEhkl = np.absolute(Ehkl)*phase
3599        Ncyc += 1
3600        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3601        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3602        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3603        if Rcf < 5.:
3604            break
3605        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3606        if not GoOn or Ncyc > 10000:
3607            break
3608    np.seterr(**old)
3609    G2fil.G2Print (' Charge flip time: %.4f no. elements: %d'%(time.time()-time0,Ehkl.size))
3610    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))/10.    #? to get on same scale as e-map
3611    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.                  #? ditto
3612    G2fil.G2Print (' No.cycles = %d Residual Rcf =%8.3f%s Map size: %s'%(Ncyc,Rcf,'%',str(CErho.shape)))
3613    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3614
3615    mapData['Rcf'] = Rcf
3616    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3617    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3618    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3619    mapData['Type'] = reflDict['Type']
3620
3621    map4DData['Rcf'] = Rcf
3622    map4DData['rho'] = np.real(np.roll(np.roll(np.roll(np.roll(SSrho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2),roll[3],axis=3))
3623    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3624    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3625    map4DData['Type'] = reflDict['Type']
3626    return mapData,map4DData
3627   
3628def getRho(xyz,mapData):
3629    ''' get scattering density at a point by 8-point interpolation
3630    param xyz: coordinate to be probed
3631    param: mapData: dict of map data
3632   
3633    :returns: density at xyz
3634    '''
3635    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3636    if not len(mapData):
3637        return 0.0
3638    rho = copy.copy(mapData['rho'])     #don't mess up original
3639    if not len(rho):
3640        return 0.0
3641    mapShape = np.array(rho.shape)
3642    mapStep = 1./mapShape
3643    X = np.array(xyz)%1.    #get into unit cell
3644    I = np.array(X*mapShape,dtype='int')
3645    D = X-I*mapStep         #position inside map cell
3646    D12 = D[0]*D[1]
3647    D13 = D[0]*D[2]
3648    D23 = D[1]*D[2]
3649    D123 = np.prod(D)
3650    Rho = rollMap(rho,-I)    #shifts map so point is in corner
3651    R = Rho[0,0,0]*(1.-np.sum(D))+Rho[1,0,0]*D[0]+Rho[0,1,0]*D[1]+Rho[0,0,1]*D[2]+  \
3652        Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \
3653        Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \
3654        Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123)   
3655    return R
3656       
3657def getRhos(XYZ,rho):
3658    ''' get scattering density at an array of point by 8-point interpolation
3659    this is faster than gerRho which is only used for single points. However, getRhos is
3660    replaced by scipy.ndimage.interpolation.map_coordinates which does a better job &  is just as fast.
3661    Thus, getRhos is unused in GSAS-II at this time.
3662    param xyz:  array coordinates to be probed Nx3
3663    param: rho: array copy of map (NB: don't use original!)
3664   
3665    :returns: density at xyz
3666    '''
3667    def getBoxes(rho,I):
3668        Rhos = np.zeros((2,2,2))
3669        Mx,My,Mz = rho.shape
3670        Ix,Iy,Iz = I
3671        Rhos = np.array([[[rho[Ix%Mx,Iy%My,Iz%Mz],rho[Ix%Mx,Iy%My,(Iz+1)%Mz]],
3672                          [rho[Ix%Mx,(Iy+1)%My,Iz%Mz],rho[Ix%Mx,(Iy+1)%My,(Iz+1)%Mz]]],
3673                        [[rho[(Ix+1)%Mx,Iy%My,Iz%Mz],rho[(Ix+1)%Mx,Iy%My,(Iz+1)%Mz]],
3674                          [rho[(Ix+1)%Mx,(Iy+1)%My,Iz%Mz],rho[(Ix+1)%Mx,(Iy+1)%My,(Iz+1)%Mz]]]])
3675        return Rhos
3676       
3677    Blk = 400     #400 doesn't seem to matter
3678    nBlk = len(XYZ)//Blk        #select Blk so this is an exact divide
3679    mapShape = np.array(rho.shape)
3680    mapStep = 1./mapShape
3681    X = XYZ%1.    #get into unit cell
3682    iBeg = 0
3683    R = np.zeros(len(XYZ))
3684#actually a lot faster!
3685    for iblk in range(nBlk):
3686        iFin = iBeg+Blk
3687        Xs = X[iBeg:iFin]
3688        I = np.array(np.rint(Xs*mapShape),dtype='int')
3689        Rhos = np.array([getBoxes(rho,i) for i in I])
3690        Ds = Xs-I*mapStep
3691        RIJs = Rhos[:,0,:2,:2]*(1.-Ds[:,0][:,nxs,nxs])
3692        RIs = RIJs[:,0]*(1.-Ds[:,1][:,nxs])+RIJs[:,1]*Ds[:,1][:,nxs]
3693        R[iBeg:iFin] = RIs[:,0]*(1.-Ds[:,2])+RIs[:,1]*Ds[:,2]
3694        iBeg += Blk
3695    return R
3696       
3697def SearchMap(generalData,drawingData,Neg=False):
3698    '''Does a search of a density map for peaks meeting the criterion of peak
3699    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
3700    mapData is data['General']['mapData']; the map is also in mapData.
3701
3702    :param generalData: the phase data structure; includes the map
3703    :param drawingData: the drawing data structure
3704    :param Neg:  if True then search for negative peaks (i.e. H-atoms & neutron data)
3705
3706    :returns: (peaks,mags,dzeros) where
3707
3708        * peaks : ndarray
3709          x,y,z positions of the peaks found in the map
3710        * mags : ndarray
3711          the magnitudes of the peaks
3712        * dzeros : ndarray
3713          the distance of the peaks from  the unit cell origin
3714        * dcent : ndarray
3715          the distance of the peaks from  the unit cell center
3716
3717    '''       
3718    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3719   
3720    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
3721   
3722    def fixSpecialPos(xyz,SGData,Amat):
3723        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
3724        X = []
3725        xyzs = [equiv[0] for equiv in equivs]
3726        for x in xyzs:
3727            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
3728                X.append(x)
3729        if len(X) > 1:
3730            return np.average(X,axis=0)
3731        else:
3732            return xyz
3733       
3734    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
3735        Mag,x0,y0,z0,sig = parms
3736        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
3737#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
3738        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
3739       
3740    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
3741        Mag,x0,y0,z0,sig = parms
3742        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3743        return M
3744       
3745    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
3746        Mag,x0,y0,z0,sig = parms
3747        dMdv = np.zeros(([5,]+list(rX.shape)))
3748        delt = .01
3749        for i in range(5):
3750            parms[i] -= delt
3751            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3752            parms[i] += 2.*delt
3753            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3754            parms[i] -= delt
3755            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
3756        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3757        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
3758        dMdv = np.reshape(dMdv,(5,rX.size))
3759        Hess = np.inner(dMdv,dMdv)
3760       
3761        return Vec,Hess
3762       
3763    SGData = generalData['SGData']
3764    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3765    peaks = []
3766    mags = []
3767    dzeros = []
3768    dcent = []
3769    try:
3770        mapData = generalData['Map']
3771        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
3772        if Neg:
3773            rho = -copy.copy(mapData['rho'])     #flip +/-
3774        else:
3775            rho = copy.copy(mapData['rho'])     #don't mess up original
3776        mapHalf = np.array(rho.shape)/2
3777        res = mapData['GridStep']*2.
3778        incre = np.array(rho.shape,dtype=np.float)
3779        step = max(1.0,1./res)+1
3780        steps = np.array((3*[step,]),dtype='int32')
3781    except KeyError:
3782        G2fil.G2Print ('**** ERROR - Fourier map not defined')
3783        return peaks,mags
3784    rhoMask = ma.array(rho,mask=(rho<contLevel))
3785    indices = (-1,0,1)
3786    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3787    for roll in rolls:
3788        if np.any(roll):
3789            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
3790    indx = np.transpose(rhoMask.nonzero())
3791    peaks = indx/incre
3792    mags = rhoMask[rhoMask.nonzero()]
3793    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
3794        rho = rollMap(rho,ind)
3795        rMM = mapHalf-steps
3796        rMP = mapHalf+steps+1
3797        rhoPeak = rho[int(rMM[0]):int(rMP[0]),int(rMM[1]):int(rMP[1]),int(rMM[2]):int(rMP[2])]
3798        peakInt = np.sum(rhoPeak)*res**3
3799        rX,rY,rZ = np.mgrid[int(rMM[0]):int(rMP[0]),int(rMM[1]):int(rMP[1]),int(rMM[2]):int(rMP[2])]
3800        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
3801        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
3802            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
3803        x1 = result[0]
3804        if not np.any(x1 < 0):
3805            peak = (np.array(x1[1:4])-ind)/incre
3806        peak = fixSpecialPos(peak,SGData,Amat)
3807        rho = rollMap(rho,-ind)
3808    cent = np.ones(3)*.5     
3809    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
3810    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
3811    if Neg:     #want negative magnitudes for negative peaks
3812        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3813    else:
3814        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3815   
3816def sortArray(data,pos,reverse=False):
3817    '''data is a list of items
3818    sort by pos in list; reverse if True
3819    '''
3820    T = []
3821    for i,M in enumerate(data):
3822        try:
3823            T.append((M[pos],i))
3824        except IndexError:
3825            return data
3826    D = dict(zip(T,data))
3827    T.sort()
3828    if reverse:
3829        T.reverse()
3830    X = []
3831    for key in T:
3832        X.append(D[key])
3833    return X
3834
3835def PeaksEquiv(data,Ind):
3836    '''Find the equivalent map peaks for those selected. Works on the
3837    contents of data['Map Peaks'].
3838
3839    :param data: the phase data structure
3840    :param list Ind: list of selected peak indices
3841    :returns: augmented list of peaks including those related by symmetry to the
3842      ones in Ind
3843
3844    '''       
3845    def Duplicate(xyz,peaks,Amat):
3846        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3847            return True
3848        return False
3849                           
3850    generalData = data['General']
3851    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3852    SGData = generalData['SGData']
3853    mapPeaks = data['Map Peaks']
3854    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
3855    Indx = {}
3856    for ind in Ind:
3857        xyz = np.array(mapPeaks[ind][1:4])
3858        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
3859        for jnd,xyz in enumerate(XYZ):       
3860            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
3861    Ind = []
3862    for ind in Indx:
3863        if Indx[ind]:
3864            Ind.append(ind)
3865    return Ind
3866               
3867def PeaksUnique(data,Ind):
3868    '''Finds the symmetry unique set of peaks from those selected. Works on the
3869    contents of data['Map Peaks'].
3870
3871    :param data: the phase data structure
3872    :param list Ind: list of selected peak indices
3873    :returns: the list of symmetry unique peaks from among those given in Ind
3874
3875    '''       
3876#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
3877
3878    def noDuplicate(xyz,peaks,Amat):
3879        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3880            return False
3881        return True
3882                           
3883    generalData = data['General']
3884    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3885    SGData = generalData['SGData']
3886    mapPeaks = data['Map Peaks']
3887    Indx = {}
3888    XYZ = {}
3889    for ind in Ind:
3890        XYZ[ind] = np.array(mapPeaks[ind][1:4])
3891        Indx[ind] = True
3892    for ind in Ind:
3893        if Indx[ind]:
3894            xyz = XYZ[ind]
3895            for jnd in Ind:
3896                if ind != jnd and Indx[jnd]:                       
3897                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
3898                    xyzs = np.array([equiv[0] for equiv in Equiv])
3899                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
3900    Ind = []
3901    for ind in Indx:
3902        if Indx[ind]:
3903            Ind.append(ind)
3904    return Ind
3905
3906################################################################################
3907##### Dysnomia setup & return stuff
3908################################################################################
3909   
3910
3911   
3912################################################################################
3913##### single peak fitting profile fxn stuff
3914################################################################################
3915
3916def getCWsig(ins,pos):
3917    '''get CW peak profile sigma^2
3918   
3919    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
3920      as values only
3921    :param float pos: 2-theta of peak
3922    :returns: float getCWsig: peak sigma^2
3923   
3924    '''
3925    tp = tand(pos/2.0)
3926    return ins['U']*tp**2+ins['V']*tp+ins['W']
3927   
3928def getCWsigDeriv(pos):
3929    '''get derivatives of CW peak profile sigma^2 wrt U,V, & W
3930   
3931    :param float pos: 2-theta of peak
3932   
3933    :returns: list getCWsigDeriv: d(sig^2)/dU, d(sig)/dV & d(sig)/dW
3934   
3935    '''
3936    tp = tand(pos/2.0)
3937    return tp**2,tp,1.0
3938   
3939def getCWgam(ins,pos):
3940    '''get CW peak profile gamma
3941   
3942    :param dict ins: instrument parameters with at least 'X', 'Y' & 'Z'
3943      as values only
3944    :param float pos: 2-theta of peak
3945    :returns: float getCWgam: peak gamma
3946   
3947    '''
3948    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)+ins.get('Z',0.0)
3949   
3950def getCWgamDeriv(pos):
3951    '''get derivatives of CW peak profile gamma wrt X, Y & Z
3952   
3953    :param float pos: 2-theta of peak
3954   
3955    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
3956   
3957    '''
3958    return 1./cosd(pos/2.0),tand(pos/2.0),1.0
3959   
3960def getTOFsig(ins,dsp):
3961    '''get TOF peak profile sigma^2
3962   
3963    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
3964      as values only
3965    :param float dsp: d-spacing of peak
3966   
3967    :returns: float getTOFsig: peak sigma^2
3968   
3969    '''
3970    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']*dsp
3971   
3972def getTOFsigDeriv(dsp):
3973    '''get derivatives of TOF peak profile sigma^2 wrt sig-0, sig-1, & sig-q
3974   
3975    :param float dsp: d-spacing of peak
3976   
3977    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
3978   
3979    '''
3980    return 1.0,dsp**2,dsp**4,dsp
3981   
3982def getTOFgamma(ins,dsp):
3983    '''get TOF peak profile gamma
3984   
3985    :param dict ins: instrument parameters with at least 'X', 'Y' & 'Z'
3986      as values only
3987    :param float dsp: d-spacing of peak
3988   
3989    :returns: float getTOFgamma: peak gamma
3990   
3991    '''
3992    return ins['Z']+ins['X']*dsp+ins['Y']*dsp**2
3993   
3994def getTOFgammaDeriv(dsp):
3995    '''get derivatives of TOF peak profile gamma wrt X, Y & Z
3996   
3997    :param float dsp: d-spacing of peak
3998   
3999    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
4000   
4001    '''
4002    return dsp,dsp**2,1.0
4003   
4004def getTOFbeta(ins,dsp):
4005    '''get TOF peak profile beta
4006   
4007    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
4008      as values only
4009    :param float dsp: d-spacing of peak
4010   
4011    :returns: float getTOFbeta: peak beat
4012   
4013    '''
4014    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
4015   
4016def getTOFbetaDeriv(dsp):
4017    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
4018   
4019    :param float dsp: d-spacing of peak
4020   
4021    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
4022   
4023    '''
4024    return 1.0,1./dsp**4,1./dsp**2
4025   
4026def getTOFalpha(ins,dsp):
4027    '''get TOF peak profile alpha
4028   
4029    :param dict ins: instrument parameters with at least 'alpha'
4030      as values only
4031    :param float dsp: d-spacing of peak
4032   
4033    :returns: flaot getTOFalpha: peak alpha
4034   
4035    '''
4036    return ins['alpha']/dsp
4037   
4038def getTOFalphaDeriv(dsp):
4039    '''get derivatives of TOF peak profile beta wrt alpha
4040   
4041    :param float dsp: d-spacing of peak
4042   
4043    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
4044   
4045    '''
4046    return 1./dsp
4047   
4048def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
4049    '''set starting peak parameters for single peak fits from plot selection or auto selection
4050   
4051    :param dict Parms: instrument parameters dictionary
4052    :param dict Parms2: table lookup for TOF profile coefficients
4053    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
4054    :param float mag: peak top magnitude from pick
4055    :param bool ifQ: True if pos in Q
4056    :param bool useFit: True if use fitted CW Parms values (not defaults)
4057   
4058    :returns: list XY: peak list entry:
4059        for CW: [pos,0,mag,1,sig,0,gam,0]
4060        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
4061        NB: mag refinement set by default, all others off
4062   
4063    '''
4064    ind = 0
4065    if useFit:
4066        ind = 1
4067    ins = {}
4068    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an else
4069        for x in ['U','V','W','X','Y','Z']:
4070            ins[x] = Parms.get(x,[0.0,0.0])[ind]
4071        if ifQ:                              #qplot - convert back to 2-theta
4072            pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi))
4073        sig = getCWsig(ins,pos)
4074        gam = getCWgam(ins,pos)           
4075        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
4076    else:
4077        if ifQ:
4078            dsp = 2.*np.pi/pos
4079            pos = Parms['difC']*dsp
4080        else:
4081            dsp = pos/Parms['difC'][1]
4082        if 'Pdabc' in Parms2:
4083            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y','Z']:
4084                ins[x] = Parms.get(x,[0.0,0.0])[ind]
4085            Pdabc = Parms2['Pdabc'].T
4086            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
4087            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
4088        else:
4089            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y','Z']:
4090                ins[x] = Parms.get(x,[0.0,0.0])[ind]
4091            alp = getTOFalpha(ins,dsp)
4092            bet = getTOFbeta(ins,dsp)
4093        sig = getTOFsig(ins,dsp)
4094        gam = getTOFgamma(ins,dsp)
4095        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
4096    return XY
4097   
4098################################################################################
4099##### MC/SA stuff
4100################################################################################
4101
4102#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
4103# Original Author: Travis Oliphant 2002
4104# Bug-fixes in 2006 by Tim Leslie
4105
4106
4107import numpy
4108from numpy import asarray, exp, squeeze, sign, \
4109     all, shape, array, where
4110from numpy import random
4111
4112#__all__ = ['anneal']
4113
4114_double_min = numpy.finfo(float).min
4115_double_max = numpy.finfo(float).max
4116class base_schedule(object):
4117    def __init__(self):
4118        self.dwell = 20
4119        self.lower = 0.
4120        self.upper = 1.
4121        self.Ninit = 50
4122        self.accepted = 0
4123        self.tests = 0
4124        self.feval = 0
4125        self.k = 0
4126        self.T = None
4127
4128    def init(self, **options):
4129        self.__dict__.update(options)
4130        self.lower = asarray(self.lower)
4131        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
4132        self.upper = asarray(self.upper)
4133        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
4134        self.k = 0
4135        self.accepted = 0
4136        self.feval = 0
4137        self.tests = 0
4138
4139    def getstart_temp(self, best_state):
4140        """ Find a matching starting temperature and starting parameters vector
4141        i.e. find x0 such that func(x0) = T0.
4142
4143        :param best_state: _state
4144            A _state object to store the function value and x0 found.
4145
4146        :returns: x0 : array
4147            The starting parameters vector.
4148        """
4149
4150        assert(not self.dims is None)
4151        lrange = self.lower
4152        urange = self.upper
4153        fmax = _double_min
4154        fmin = _double_max
4155        for _ in range(self.Ninit):
4156            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
4157            fval = self.func(x0, *self.args)
4158            self.feval += 1
4159            if fval > fmax:
4160                fmax = fval
4161            if fval < fmin:
4162                fmin = fval
4163                best_state.cost = fval
4164                best_state.x = array(x0)
4165
4166        self.T0 = (fmax-fmin)*1.5
4167        return best_state.x
4168       
4169    def set_range(self,x0,frac):
4170        delrange = frac*(self.upper-self.lower)
4171        self.upper = x0+delrange
4172        self.lower = x0-delrange
4173
4174    def accept_test(self, dE):
4175        T = self.T
4176        self.tests += 1
4177        if dE < 0:
4178            self.accepted += 1
4179            return 1
4180        p = exp(-dE*1.0/T)
4181        if (p > random.uniform(0.0, 1.0)):
4182            self.accepted += 1
4183            return 1
4184        return 0
4185
4186    def update_guess(self, x0):
4187        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
4188
4189    def update_temp(self, x0):
4190        pass
4191
4192class fast_sa(base_schedule):
4193    def init(self, **options):
4194        self.__dict__.update(options)
4195
4196    def update_guess(self, x0):
4197        x0 = asarray(x0)
4198        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4199        T = self.T
4200        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
4201        xnew = xc*(self.upper - self.lower)+self.lower
4202        return xnew
4203
4204    def update_temp(self):
4205        self.T = self.T0*exp(-self.c * self.k**(self.quench))
4206        self.k += 1
4207        return
4208
4209class log_sa(base_schedule):        #OK
4210
4211    def init(self,**options):
4212        self.__dict__.update(options)
4213       
4214    def update_guess(self,x0):     #same as default #TODO - is this a reasonable update procedure?
4215        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4216        T = self.T
4217        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
4218        xnew = xc*(self.upper - self.lower)+self.lower
4219        return xnew
4220       
4221    def update_temp(self):
4222        self.k += 1
4223        self.T = self.T0*self.slope**self.k
4224       
4225class _state(object):
4226    def __init__(self):
4227        self.x = None
4228        self.cost = None
4229
4230def makeTsched(data):
4231    if data['Algorithm'] == 'fast':
4232        sched = fast_sa()
4233        sched.quench = data['fast parms'][0]
4234        sched.c = data['fast parms'][1]
4235    elif data['Algorithm'] == 'log':
4236        sched = log_sa()
4237        sched.slope = data['log slope']
4238    sched.T0 = data['Annealing'][0]
4239    if not sched.T0:
4240        sched.T0 = 50.
4241    Tf = data['Annealing'][1]
4242    if not Tf:
4243        Tf = 0.001
4244    Tsched = [sched.T0,]
4245    while Tsched[-1] > Tf:
4246        sched.update_temp()
4247        Tsched.append(sched.T)
4248    return Tsched[1:]
4249   
4250def anneal(func, x0, args=(), schedule='fast', 
4251           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
4252           feps=1e-6, quench=1.0, c=1.0,
4253           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
4254           ranRange=0.10,autoRan=False,dlg=None):
4255    """Minimize a function using simulated annealing.
4256
4257    Schedule is a schedule class implementing the annealing schedule.
4258    Available ones are 'fast', 'cauchy', 'boltzmann'
4259
4260    :param callable func: f(x, \*args)
4261        Function to be optimized.
4262    :param ndarray x0:
4263        Initial guess.
4264    :param tuple args:
4265        Extra parameters to `func`.
4266    :param base_schedule schedule:
4267        Annealing schedule to use (a class).
4268    :param float T0:
4269        Initial Temperature (estimated as 1.2 times the largest
4270        cost-function deviation over random points in the range).
4271    :param float Tf:
4272        Final goal temperature.
4273    :param int maxeval:
4274        Maximum function evaluations.
4275    :param int maxaccept:
4276        Maximum changes to accept.
4277    :param int maxiter:
4278        Maximum cooling iterations.
4279    :param float feps:
4280        Stopping relative error tolerance for the function value in
4281        last four coolings.
4282    :param float quench,c:
4283        Parameters to alter fast_sa schedule.
4284    :param float/ndarray lower,upper:
4285        Lower and upper bounds on `x`.
4286    :param int dwell:
4287        The number of times to search the space at each temperature.
4288    :param float slope:
4289        Parameter for log schedule
4290    :param bool ranStart=False:
4291        True for set 10% of ranges about x
4292
4293    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
4294
4295     * xmin (ndarray): Point giving smallest value found.
4296     * Jmin (float): Minimum value of function found.
4297     * T (float): Final temperature.
4298     * feval (int): Number of function evaluations.
4299     * iters (int): Number of cooling iterations.
4300     * accept (int): Number of tests accepted.
4301     * retval (int): Flag indicating stopping condition:
4302
4303              *  0: Points no longer changing
4304              *  1: Cooled to final temperature
4305              *  2: Maximum function evaluations
4306              *  3: Maximum cooling iterations reached
4307              *  4: Maximum accepted query locations reached
4308              *  5: Final point not the minimum amongst encountered points
4309
4310    *Notes*:
4311    Simulated annealing is a random algorithm which uses no derivative
4312    information from the function being optimized. In practice it has
4313    been more useful in discrete optimization than continuous
4314    optimization, as there are usually better algorithms for continuous
4315    optimization problems.
4316
4317    Some experimentation by trying the difference temperature
4318    schedules and altering their parameters is likely required to
4319    obtain good performance.
4320
4321    The randomness in the algorithm comes from random sampling in numpy.
4322    To obtain the same results you can call numpy.random.seed with the
4323    same seed immediately before calling scipy.optimize.anneal.
4324
4325    We give a brief description of how the three temperature schedules
4326    generate new points and vary their temperature. Temperatures are
4327    only updated with iterations in the outer loop. The inner loop is
4328    over range(dwell), and new points are generated for
4329    every iteration in the inner loop. (Though whether the proposed
4330    new points are accepted is probabilistic.)
4331
4332    For readability, let d denote the dimension of the inputs to func.
4333    Also, let x_old denote the previous state, and k denote the
4334    iteration number of the outer loop. All other variables not
4335    defined below are input variables to scipy.optimize.anneal itself.
4336
4337    In the 'fast' schedule the updates are ::
4338
4339        u ~ Uniform(0, 1, size=d)
4340        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
4341        xc = y * (upper - lower)
4342        x_new = x_old + xc
4343
4344        T_new = T0 * exp(-c * k**quench)
4345
4346    """
4347   
4348    ''' Scipy license:
4349        Copyright (c) 2001, 2002 Enthought, Inc.
4350    All rights reserved.
4351   
4352    Copyright (c) 2003-2016 SciPy Developers.
4353    All rights reserved.
4354   
4355    Redistribution and use in source and binary forms, with or without
4356    modification, are permitted provided that the following conditions are met:
4357   
4358      a. Redistributions of source code must retain the above copyright notice,
4359         this list of conditions and the following disclaimer.
4360      b. Redistributions in binary form must reproduce the above copyright
4361         notice, this list of conditions and the following disclaimer in the
4362         documentation and/or other materials provided with the distribution.
4363      c. Neither the name of Enthought nor the names of the SciPy Developers
4364         may be used to endorse or promote products derived from this software
4365         without specific prior written permission.
4366    '''
4367    x0 = asarray(x0)
4368    lower = asarray(lower)
4369    upper = asarray(upper)
4370
4371    schedule = eval(schedule+'_sa()')
4372    #   initialize the schedule
4373    schedule.init(dims=shape(x0),func=func,args=args,T0=T0,lower=lower, upper=upper,
4374        c=c, quench=quench, dwell=dwell, slope=slope)
4375
4376    current_state, last_state, best_state = _state(), _state(), _state()
4377    if ranStart:
4378        schedule.set_range(x0,ranRange)
4379    if T0 is None:
4380        x0 = schedule.getstart_temp(best_state)
4381    else:
4382        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
4383        best_state.x = None
4384        best_state.cost = numpy.Inf
4385
4386    last_state.x = asarray(x0).copy()
4387    fval = func(x0,*args)
4388    schedule.feval += 1
4389    last_state.cost = fval
4390    if last_state.cost < best_state.cost:
4391        best_state.cost = fval
4392        best_state.x = asarray(x0).copy()
4393    schedule.T = schedule.T0
4394    fqueue = [100, 300, 500, 700]
4395    iters = 1
4396    keepGoing = True
4397    bestn = 0
4398    while keepGoing:
4399        retval = 0
4400        for n in range(dwell):
4401            current_state.x = schedule.update_guess(last_state.x)
4402            current_state.cost = func(current_state.x,*args)
4403            schedule.feval += 1
4404
4405            dE = current_state.cost - last_state.cost
4406            if schedule.accept_test(dE):
4407                last_state.x = current_state.x.copy()
4408                last_state.cost = current_state.cost
4409                if last_state.cost < best_state.cost:
4410                    best_state.x = last_state.x.copy()
4411                    best_state.cost = last_state.cost
4412                    bestn = n
4413                    if best_state.cost < 1.0 and autoRan:
4414                        schedule.set_range(x0,best_state.cost/2.)                       
4415        if dlg:
4416            GoOn = dlg.Update(min(100.,best_state.cost*100),
4417                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
4418                    'Best trial:',bestn,  \
4419                    'MC/SA Residual:',best_state.cost*100,'%', \
4420                    ))[0]
4421            if not GoOn:
4422                best_state.x = last_state.x.copy()
4423                best_state.cost = last_state.cost
4424                retval = 5
4425        schedule.update_temp()
4426        iters += 1
4427        # Stopping conditions
4428        # 0) last saved values of f from each cooling step
4429        #     are all very similar (effectively cooled)
4430        # 1) Tf is set and we are below it
4431        # 2) maxeval is set and we are past it
4432        # 3) maxiter is set and we are past it
4433        # 4) maxaccept is set and we are past it
4434        # 5) user canceled run via progress bar
4435
4436        fqueue.append(squeeze(last_state.cost))
4437        fqueue.pop(0)
4438        af = asarray(fqueue)*1.0
4439        if retval == 5:
4440            G2fil.G2Print ('Error: User terminated run; incomplete MC/SA')
4441            keepGoing = False
4442            break
4443        if all(abs((af-af[0])/af[0]) < feps):
4444            retval = 0
4445            if abs(af[-1]-best_state.cost) > feps*10:
4446                retval = 5
4447                G2fil.G2Print (" Warning: Cooled to %.4f > selected Tmin %.4f in %d steps"%(squeeze(last_state.cost),Tf,iters-1))
4448            break
4449        if (Tf is not None) and (schedule.T < Tf):
4450#            print ' Minimum T reached in %d steps'%(iters-1)
4451            retval = 1
4452            break
4453        if (maxeval is not None) and (schedule.feval > maxeval):
4454            retval = 2
4455            break
4456        if (iters > maxiter):
4457            G2fil.G2Print  (" Warning: Maximum number of iterations exceeded.")
4458            retval = 3
4459            break
4460        if (maxaccept is not None) and (schedule.accepted > maxaccept):
4461            retval = 4
4462            break
4463
4464    return best_state.x, best_state.cost, schedule.T, \
4465           schedule.feval, iters, schedule.accepted, retval
4466
4467def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,nprocess=-1):
4468    outlist = []
4469    timelist = []
4470    nsflist = []
4471    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
4472    for n in range(iCyc):
4473        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None,False)         #mcsa result,time,rcov
4474        outlist.append(result[0])
4475        timelist.append(result[1])
4476        nsflist.append(result[2])
4477        G2fil.G2Print (' MC/SA final fit: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1]))
4478    out_q.put(outlist)
4479    out_t.put(timelist)
4480    out_n.put(nsflist)
4481
4482def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData,nprocs):
4483    import multiprocessing as mp
4484   
4485    out_q = mp.Queue()
4486    out_t = mp.Queue()
4487    out_n = mp.Queue()
4488    procs = []
4489    totsftime = 0.
4490    totnsf = 0
4491    iCyc = np.zeros(nprocs)
4492    for i in range(nCyc):
4493        iCyc[i%nprocs] += 1
4494    for i in range(nprocs):
4495        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,i))
4496        procs.append(p)
4497        p.start()
4498    resultlist = []
4499    for i in range(nprocs):
4500        resultlist += out_q.get()
4501        totsftime += np.sum(out_t.get())
4502        totnsf += np.sum(out_n.get())
4503    for p in procs:
4504        p.join()
4505    return resultlist,totsftime,totnsf
4506
4507def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar,start=True):
4508    '''default doc string
4509   
4510    :param type name: description
4511   
4512    :returns: type name: description
4513    '''
4514   
4515    class RandomDisplacementBounds(object):
4516        """random displacement with bounds"""
4517        def __init__(self, xmin, xmax, stepsize=0.5):
4518            self.xmin = xmin
4519            self.xmax = xmax
4520            self.stepsize = stepsize
4521   
4522        def __call__(self, x):
4523            """take a random step but ensure the new position is within the bounds"""
4524            while True:
4525                # this could be done in a much more clever way, but it will work for example purposes
4526                steps = self.xmax-self.xmin
4527                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
4528                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
4529                    break
4530            return xnew
4531   
4532    global tsum,nsum
4533    tsum = 0.
4534    nsum = 0
4535   
4536    def getMDparms(item,pfx,parmDict,varyList):
4537        parmDict[pfx+'MDaxis'] = item['axis']
4538        parmDict[pfx+'MDval'] = item['Coef'][0]
4539        if item['Coef'][1]:
4540            varyList += [pfx+'MDval',]
4541            limits = item['Coef'][2]
4542            lower.append(limits[0])
4543            upper.append(limits[1])
4544                       
4545    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
4546        parmDict[pfx+'Atype'] = item['atType']
4547        aTypes |= set([item['atType'],]) 
4548        pstr = ['Ax','Ay','Az']
4549        XYZ = [0,0,0]
4550        for i in range(3):
4551            name = pfx+pstr[i]
4552            parmDict[name] = item['Pos'][0][i]
4553            XYZ[i] = parmDict[name]
4554            if item['Pos'][1][i]:
4555                varyList += [name,]
4556                limits = item['Pos'][2][i]
4557                lower.append(limits[0])
4558                upper.append(limits[1])
4559        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
4560           
4561    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
4562        parmDict[mfx+'MolCent'] = item['MolCent']
4563        parmDict[mfx+'RBId'] = item['RBId']
4564        pstr = ['Px','Py','Pz']
4565        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
4566        for i in range(3):
4567            name = mfx+pstr[i]
4568            parmDict[name] = item['Pos'][0][i]
4569            if item['Pos'][1][i]:
4570                varyList += [name,]
4571                limits = item['Pos'][2][i]
4572                lower.append(limits[0])
4573                upper.append(limits[1])
4574        AV = item['Ori'][0]
4575        A = AV[0]
4576        V = AV[1:]
4577        for i in range(