source: trunk/GSASIImath.py @ 3079

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

implement rigid body reference vectors

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