source: trunk/GSASIImath.py @ 3853

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

incomm.mag. mods

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