source: trunk/GSASIImath.py @ 3777

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

fix wx.SAVE --> wx.FD_SAVE bug
work on incommensurate magnetic structure factor calcs. close but no cigar yet

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