source: trunk/GSASIImath.py @ 3029

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

implement correctly errat for protein validation; erratv2 not right

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