source: trunk/GSASIImath.py @ 4185

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

mag math again

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