source: trunk/GSASIImath.py @ 3045

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

plot both protein vlidators (#2 still in erroe)
remove config value for choice - not needed

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