source: branch/2frame/GSASIImath.py @ 2964

Last change on this file since 2964 was 2964, checked in by vondreele, 5 years ago

work on protein validation - not correct yet; gives plot tho.
fix refinement of Dij values in a LeBail? refinement (esp. sequential LeBail?).

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