source: trunk/GSASIImath.py @ 3837

Last change on this file since 3837 was 3837, checked in by vondreele, 4 years ago

more on incomm. mag. str. fctr.

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