source: trunk/GSASIImath.py @ 3003

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

more development of protein validation - still not correct
substitute old "paired" & "Paired-r" continuous colors for the installed one (16 discrete colors in mpl 2.0.X
trap use of tick_label in barh (not in mpl <2.0)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 182.5 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2017-08-12 21:23:31 +0000 (Sat, 12 Aug 2017) $
5# $Author: vondreele $
6# $Revision: 3003 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 3003 2017-08-12 21:23:31Z 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: 3003 $")
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, skipping main chain C & N; make cartesian
2602    for atom in Atoms:
2603        if atom[1] in resNames:
2604            if atom[4].strip() in ['S','Se']:
2605                atom[4] = 'O'
2606#            if atom[3].strip() in ['C','N']:    #skip main chain C & N atoms
2607#                continue
2608            cartAtoms.append(atom[:cx+3])
2609            cartAtoms[-1][cx:cx+3] = np.inner(Amat,cartAtoms[-1][cx:cx+3])
2610    XYZ = np.array([atom[cx:cx+3] for atom in cartAtoms])
2611    xyzmin = np.array([np.min(XYZ.T[i]) for i in [0,1,2]])
2612    xyzmax = np.array([np.max(XYZ.T[i]) for i in [0,1,2]])
2613    nbox = list(np.array(np.ceil((xyzmax-xyzmin)/4.),dtype=int))+[15,]
2614    Boxes = np.zeros(nbox,dtype=int)
2615    iBox = np.array([np.trunc((XYZ.T[i]-xyzmin[i])/4.) for i in [0,1,2]],dtype=int).T
2616    for ib,box in enumerate(iBox):  #put in a try for too many atoms in box (IndexError)?
2617        Boxes[box[0],box[1],box[2],0] += 1
2618        Boxes[box[0],box[1],box[2],Boxes[box[0],box[1],box[2],0]] = ib
2619    #Box content checks with errat.f ibox1 array
2620    indices = (-1,0,1)
2621    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices]) 
2622    dsmax = 3.75**2
2623    if old:
2624        dsmax = 3.5**2
2625    chains = []
2626    resIntAct = []
2627    chainIntAct = []
2628    res = []
2629    resNames = []
2630    resname = []
2631    newChain = True
2632    intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2633    for ia,atom in enumerate(cartAtoms):
2634        if atom[2] not in chains:   #get chain id & save residue sequence from last chain
2635            chains.append(atom[2])
2636            if len(resIntAct):
2637                resIntAct.append(sumintact(intact))
2638                chainIntAct.append(resIntAct)
2639                resNames += resname
2640                res = []
2641                resname = []
2642                resIntAct = []
2643                intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2644                newChain = True
2645        if atom[0] not in res:  #new residue, get residue no.
2646            if len(res): print res[-1],resname[-1],intact
2647            res.append(atom[0])
2648            resname.append('%s-%s%s'%(atom[2],atom[0],atom[1]))
2649            if not newChain:
2650                resIntAct.append(sumintact(intact))
2651            intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
2652            newChain = False
2653        ibox = iBox[ia]         #box location of atom
2654        tgts = []
2655        for unit in Units:      #assemble list of all possible target atoms
2656            jbox = ibox+unit
2657            if np.all(jbox>=0) and np.all((jbox-nbox[:3])<0):               
2658                tgts += list(Boxes[jbox[0],jbox[1],jbox[2]])
2659        tgts = list(set(tgts))
2660        tgts = [tgt for tgt in tgts if atom[1:3] != cartAtoms[tgt][1:3]]    #exclude same residue
2661        if atom[3].strip() in ['C',]:
2662            tgts = [tgt for tgt in tgts if cartAtoms[tgt][3].strip() not in ['N',]]  #exclude main chain C-N
2663        if atom[3].strip() in ['N',]:
2664            tgts = [tgt for tgt in tgts if cartAtoms[tgt][3].strip() not in ['C',]]  #exclude main chain N-C
2665        tgts = [tgt for tgt in tgts if np.sum((XYZ[ia]-XYZ[tgt])**2) <= dsmax]
2666        for tgt in tgts:
2667            dsqt = np.sqrt(np.sum((XYZ[ia]-XYZ[tgt])**2))
2668            mult = 1.0
2669            if dsqt > 3.25 and not old:
2670                mult = 2.*(3.75-dsqt)
2671            intype = atom[4].strip()+cartAtoms[tgt][4].strip()
2672            intact[intype] += mult
2673    resNames += resname
2674    resIntAct.append(sumintact(intact))
2675    chainIntAct.append(resIntAct)
2676    chainProb = []
2677    for ich,chn in enumerate(chains):
2678        IntAct = chainIntAct[ich]
2679        nRes = len(IntAct)
2680        Probs = [0.,0.,0.,0.]
2681        for i in range(4,nRes-4):
2682            print i,IntAct[i]
2683            mtrx = np.zeros(5)
2684            summ = 0.
2685            for j in range(i-4,i+5):
2686                summ += np.sum(np.array(IntAct[j].values()))
2687                if old:
2688                    mtrx[0] += IntAct[j]['CC']
2689                    mtrx[1] += IntAct[j]['CO']
2690                    mtrx[2] += IntAct[j]['NN']
2691                    mtrx[3] += IntAct[j]['NO']
2692                    mtrx[4] += IntAct[j]['OO']
2693                else:
2694                    mtrx[0] += IntAct[j]['CC']
2695                    mtrx[1] += IntAct[j]['CN']
2696                    mtrx[2] += IntAct[j]['CO']
2697                    mtrx[3] += IntAct[j]['NN']
2698                    mtrx[4] += IntAct[j]['NO']
2699            print mtrx,np.sum(mtrx)
2700            print summ
2701            mtrx /= summ
2702            if old:
2703                mtrx -= avg_old
2704                prob = np.inner(np.inner(mtrx,b1_old),mtrx)
2705            else:
2706                mtrx -= avg
2707                prob = np.inner(np.inner(mtrx,b1),mtrx)
2708            print i, mtrx
2709            Probs.append(prob)
2710        Probs += 4*[0.,]
2711        chainProb += Probs
2712    return resNames,chainProb
2713   
2714################################################################################
2715##### Texture fitting stuff
2716################################################################################
2717
2718def FitTexture(General,Gangls,refData,keyList,pgbar):
2719    import pytexture as ptx
2720    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2721   
2722    def printSpHarm(textureData,SHtextureSig):
2723        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
2724        names = ['omega','chi','phi']
2725        namstr = '  names :'
2726        ptstr =  '  values:'
2727        sigstr = '  esds  :'
2728        for name in names:
2729            namstr += '%12s'%('Sample '+name)
2730            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2731            if 'Sample '+name in SHtextureSig:
2732                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2733            else:
2734                sigstr += 12*' '
2735        print namstr
2736        print ptstr
2737        print sigstr
2738        print '\n Texture coefficients:'
2739        SHcoeff = textureData['SH Coeff'][1]
2740        SHkeys = SHcoeff.keys()
2741        nCoeff = len(SHcoeff)
2742        nBlock = nCoeff/10+1
2743        iBeg = 0
2744        iFin = min(iBeg+10,nCoeff)
2745        for block in range(nBlock):
2746            namstr = '  names :'
2747            ptstr =  '  values:'
2748            sigstr = '  esds  :'
2749            for name in SHkeys[iBeg:iFin]:
2750                if 'C' in name:
2751                    namstr += '%12s'%(name)
2752                    ptstr += '%12.3f'%(SHcoeff[name])
2753                    if name in SHtextureSig:
2754                        sigstr += '%12.3f'%(SHtextureSig[name])
2755                    else:
2756                        sigstr += 12*' '
2757            print namstr
2758            print ptstr
2759            print sigstr
2760            iBeg += 10
2761            iFin = min(iBeg+10,nCoeff)
2762           
2763    def Dict2Values(parmdict, varylist):
2764        '''Use before call to leastsq to setup list of values for the parameters
2765        in parmdict, as selected by key in varylist'''
2766        return [parmdict[key] for key in varylist] 
2767       
2768    def Values2Dict(parmdict, varylist, values):
2769        ''' Use after call to leastsq to update the parameter dictionary with
2770        values corresponding to keys in varylist'''
2771        parmdict.update(zip(varylist,values))
2772       
2773    def errSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2774        parmDict.update(zip(varyList,values))
2775        Mat = np.empty(0)
2776        sumObs = 0
2777        Sangls = [parmDict['Sample '+'omega'],parmDict['Sample '+'chi'],parmDict['Sample '+'phi']]
2778        for hist in Gangls.keys():
2779            Refs = refData[hist]
2780            Refs[:,5] = np.where(Refs[:,5]>0.,Refs[:,5],0.)
2781            wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2782#            wt = 1./np.max(Refs[:,4],.25)
2783            sumObs += np.sum(wt*Refs[:,5])
2784            Refs[:,6] = 1.
2785            H = Refs[:,:3]
2786            phi,beta = G2lat.CrsAng(H,cell,SGData)
2787            psi,gam,x,x = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2788            for item in parmDict:
2789                if 'C' in item:
2790                    L,M,N = eval(item.strip('C'))
2791                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2792                    Ksl,x,x = G2lat.GetKsl(L,M,shModel,psi,gam)
2793                    Lnorm = G2lat.Lnorm(L)
2794                    Refs[:,6] += parmDict[item]*Lnorm*Kcl*Ksl
2795            mat = wt*(Refs[:,5]-Refs[:,6])
2796            Mat = np.concatenate((Mat,mat))
2797        sumD = np.sum(np.abs(Mat))
2798        R = min(100.,100.*sumD/sumObs)
2799        pgbar.Update(R,newmsg='Residual = %5.2f'%(R))
2800        print ' Residual: %.3f%%'%(R)
2801        return Mat
2802       
2803    def dervSpHarm(values,SGData,cell,Gangls,shModel,refData,parmDict,varyList,pgbar):
2804        Mat = np.empty(0)
2805        Sangls = [parmDict['Sample omega'],parmDict['Sample chi'],parmDict['Sample phi']]
2806        for hist in Gangls.keys():
2807            mat = np.zeros((len(varyList),len(refData[hist])))
2808            Refs = refData[hist]
2809            H = Refs[:,:3]
2810            wt = 1./np.sqrt(np.fmax(Refs[:,4],.25))
2811#            wt = 1./np.max(Refs[:,4],.25)
2812            phi,beta = G2lat.CrsAng(H,cell,SGData)
2813            psi,gam,dPdA,dGdA = G2lat.SamAng(Refs[:,3]/2.,Gangls[hist],Sangls,False) #assume not Bragg-Brentano!
2814            for j,item in enumerate(varyList):
2815                if 'C' in item:
2816                    L,M,N = eval(item.strip('C'))
2817                    Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
2818                    Ksl,dKdp,dKdg = G2lat.GetKsl(L,M,shModel,psi,gam)
2819                    Lnorm = G2lat.Lnorm(L)
2820                    mat[j] = -wt*Lnorm*Kcl*Ksl
2821                    for k,itema in enumerate(['Sample omega','Sample chi','Sample phi']):
2822                        try:
2823                            l = varyList.index(itema)
2824                            mat[l] -= parmDict[item]*wt*Lnorm*Kcl*(dKdp*dPdA[k]+dKdg*dGdA[k])
2825                        except ValueError:
2826                            pass
2827            if len(Mat):
2828                Mat = np.concatenate((Mat,mat.T))
2829            else:
2830                Mat = mat.T
2831        print 'deriv'
2832        return Mat
2833
2834    print ' Fit texture for '+General['Name']
2835    SGData = General['SGData']
2836    cell = General['Cell'][1:7]
2837    Texture = General['SH Texture']
2838    if not Texture['Order']:
2839        return 'No spherical harmonics coefficients'
2840    varyList = []
2841    parmDict = copy.copy(Texture['SH Coeff'][1])
2842    for item in ['Sample omega','Sample chi','Sample phi']:
2843        parmDict[item] = Texture[item][1]
2844        if Texture[item][0]:
2845            varyList.append(item)
2846    if Texture['SH Coeff'][0]:
2847        varyList += Texture['SH Coeff'][1].keys()
2848    while True:
2849        begin = time.time()
2850        values =  np.array(Dict2Values(parmDict, varyList))
2851        result = so.leastsq(errSpHarm,values,Dfun=dervSpHarm,full_output=True,ftol=1.e-6,
2852            args=(SGData,cell,Gangls,Texture['Model'],refData,parmDict,varyList,pgbar))
2853        ncyc = int(result[2]['nfev']/2)
2854        if ncyc:
2855            runtime = time.time()-begin   
2856            chisq = np.sum(result[2]['fvec']**2)
2857            Values2Dict(parmDict, varyList, result[0])
2858            GOF = chisq/(len(result[2]['fvec'])-len(varyList))       #reduced chi^2
2859            print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',len(result[2]['fvec']),' Number of parameters: ',len(varyList)
2860            print 'refinement time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)
2861            try:
2862                sig = np.sqrt(np.diag(result[1])*GOF)
2863                if np.any(np.isnan(sig)):
2864                    print '*** Least squares aborted - some invalid esds possible ***'
2865                break                   #refinement succeeded - finish up!
2866            except ValueError:          #result[1] is None on singular matrix
2867                print '**** Refinement failed - singular matrix ****'
2868                return None
2869        else:
2870            break
2871   
2872    if ncyc:
2873        for parm in parmDict:
2874            if 'C' in parm:
2875                Texture['SH Coeff'][1][parm] = parmDict[parm]
2876            else:
2877                Texture[parm][1] = parmDict[parm] 
2878        sigDict = dict(zip(varyList,sig))
2879        printSpHarm(Texture,sigDict)
2880       
2881    return None
2882   
2883################################################################################
2884##### Fourier & charge flip stuff
2885################################################################################
2886
2887def adjHKLmax(SGData,Hmax):
2888    '''default doc string
2889   
2890    :param type name: description
2891   
2892    :returns: type name: description
2893   
2894    '''
2895    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
2896        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
2897        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
2898        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
2899    else:
2900        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
2901        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
2902        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
2903
2904def OmitMap(data,reflDict,pgbar=None):
2905    '''default doc string
2906   
2907    :param type name: description
2908   
2909    :returns: type name: description
2910   
2911    '''
2912    generalData = data['General']
2913    if not generalData['Map']['MapType']:
2914        print '**** ERROR - Fourier map not defined'
2915        return
2916    mapData = generalData['Map']
2917    dmin = mapData['Resolution']
2918    SGData = generalData['SGData']
2919    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2920    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2921    cell = generalData['Cell'][1:8]       
2922    A = G2lat.cell2A(cell[:6])
2923    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2924    adjHKLmax(SGData,Hmax)
2925    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2926    time0 = time.time()
2927    for iref,ref in enumerate(reflDict['RefList']):
2928        if ref[4] >= dmin:
2929            Fosq,Fcsq,ph = ref[8:11]
2930            Uniq = np.inner(ref[:3],SGMT)
2931            Phi = np.inner(ref[:3],SGT)
2932            for i,hkl in enumerate(Uniq):        #uses uniq
2933                hkl = np.asarray(hkl,dtype='i')
2934                dp = 360.*Phi[i]                #and phi
2935                a = cosd(ph+dp)
2936                b = sind(ph+dp)
2937                phasep = complex(a,b)
2938                phasem = complex(a,-b)
2939                if '2Fo-Fc' in mapData['MapType']:
2940                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
2941                else:
2942                    F = np.sqrt(Fosq)
2943                h,k,l = hkl+Hmax
2944                Fhkl[h,k,l] = F*phasep
2945                h,k,l = -hkl+Hmax
2946                Fhkl[h,k,l] = F*phasem
2947    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
2948    M = np.mgrid[0:4,0:4,0:4]
2949    blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
2950    iBeg = blkIds*rho0.shape/4
2951    iFin = (blkIds+1)*rho0.shape/4
2952    rho_omit = np.zeros_like(rho0)
2953    nBlk = 0
2954    for iB,iF in zip(iBeg,iFin):
2955        rho1 = np.copy(rho0)
2956        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
2957        Fnew = fft.ifftshift(fft.ifftn(rho1))
2958        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
2959        phase = Fnew/np.absolute(Fnew)
2960        OFhkl = np.absolute(Fhkl)*phase
2961        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
2962        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]])
2963        nBlk += 1
2964        pgbar.Update(nBlk)
2965    mapData['rho'] = np.real(rho_omit)/cell[6]
2966    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2967    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2968    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
2969    return mapData
2970   
2971def FourierMap(data,reflDict):
2972    '''default doc string
2973   
2974    :param type name: description
2975   
2976    :returns: type name: description
2977   
2978    '''
2979    generalData = data['General']
2980    mapData = generalData['Map']
2981    dmin = mapData['Resolution']
2982    SGData = generalData['SGData']
2983    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2984    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2985    cell = generalData['Cell'][1:8]       
2986    A = G2lat.cell2A(cell[:6])
2987    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
2988    adjHKLmax(SGData,Hmax)
2989    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
2990#    Fhkl[0,0,0] = generalData['F000X']
2991    time0 = time.time()
2992    for iref,ref in enumerate(reflDict['RefList']):
2993        if ref[4] > dmin:
2994            Fosq,Fcsq,ph = ref[8:11]
2995            Uniq = np.inner(ref[:3],SGMT)
2996            Phi = np.inner(ref[:3],SGT)
2997            for i,hkl in enumerate(Uniq):        #uses uniq
2998                hkl = np.asarray(hkl,dtype='i')
2999                dp = 360.*Phi[i]                #and phi
3000                a = cosd(ph+dp)
3001                b = sind(ph+dp)
3002                phasep = complex(a,b)
3003                phasem = complex(a,-b)
3004                if 'Fobs' in mapData['MapType']:
3005                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
3006                    h,k,l = hkl+Hmax
3007                    Fhkl[h,k,l] = F*phasep
3008                    h,k,l = -hkl+Hmax
3009                    Fhkl[h,k,l] = F*phasem
3010                elif 'Fcalc' in mapData['MapType']:
3011                    F = np.sqrt(Fcsq)
3012                    h,k,l = hkl+Hmax
3013                    Fhkl[h,k,l] = F*phasep
3014                    h,k,l = -hkl+Hmax
3015                    Fhkl[h,k,l] = F*phasem
3016                elif 'delt-F' in mapData['MapType']:
3017                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
3018                    h,k,l = hkl+Hmax
3019                    Fhkl[h,k,l] = dF*phasep
3020                    h,k,l = -hkl+Hmax
3021                    Fhkl[h,k,l] = dF*phasem
3022                elif '2*Fo-Fc' in mapData['MapType']:
3023                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
3024                    h,k,l = hkl+Hmax
3025                    Fhkl[h,k,l] = F*phasep
3026                    h,k,l = -hkl+Hmax
3027                    Fhkl[h,k,l] = F*phasem
3028                elif 'Patterson' in mapData['MapType']:
3029                    h,k,l = hkl+Hmax
3030                    Fhkl[h,k,l] = complex(Fosq,0.)
3031                    h,k,l = -hkl+Hmax
3032                    Fhkl[h,k,l] = complex(Fosq,0.)
3033    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
3034    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
3035    mapData['Type'] = reflDict['Type']
3036    mapData['rho'] = np.real(rho)
3037    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3038    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3039   
3040def Fourier4DMap(data,reflDict):
3041    '''default doc string
3042   
3043    :param type name: description
3044   
3045    :returns: type name: description
3046   
3047    '''
3048    generalData = data['General']
3049    map4DData = generalData['4DmapData']
3050    mapData = generalData['Map']
3051    dmin = mapData['Resolution']
3052    SGData = generalData['SGData']
3053    SSGData = generalData['SSGData']
3054    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3055    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3056    cell = generalData['Cell'][1:8]       
3057    A = G2lat.cell2A(cell[:6])
3058    maxM = 4
3059    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
3060    adjHKLmax(SGData,Hmax)
3061    Hmax = np.asarray(Hmax,dtype='i')+1
3062    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
3063    time0 = time.time()
3064    for iref,ref in enumerate(reflDict['RefList']):
3065        if ref[5] > dmin:
3066            Fosq,Fcsq,ph = ref[9:12]
3067            Fosq = np.where(Fosq>0.,Fosq,0.)    #can't use Fo^2 < 0
3068            Uniq = np.inner(ref[:4],SSGMT)
3069            Phi = np.inner(ref[:4],SSGT)
3070            for i,hkl in enumerate(Uniq):        #uses uniq
3071                hkl = np.asarray(hkl,dtype='i')
3072                dp = 360.*Phi[i]                #and phi
3073                a = cosd(ph+dp)
3074                b = sind(ph+dp)
3075                phasep = complex(a,b)
3076                phasem = complex(a,-b)
3077                if 'Fobs' in mapData['MapType']:
3078                    F = np.sqrt(Fosq)
3079                    h,k,l,m = hkl+Hmax
3080                    Fhkl[h,k,l,m] = F*phasep
3081                    h,k,l,m = -hkl+Hmax
3082                    Fhkl[h,k,l,m] = F*phasem
3083                elif 'Fcalc' in mapData['MapType']:
3084                    F = np.sqrt(Fcsq)
3085                    h,k,l,m = hkl+Hmax
3086                    Fhkl[h,k,l,m] = F*phasep
3087                    h,k,l,m = -hkl+Hmax
3088                    Fhkl[h,k,l,m] = F*phasem                   
3089                elif 'delt-F' in mapData['MapType']:
3090                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
3091                    h,k,l,m = hkl+Hmax
3092                    Fhkl[h,k,l,m] = dF*phasep
3093                    h,k,l,m = -hkl+Hmax
3094                    Fhkl[h,k,l,m] = dF*phasem
3095    SSrho = fft.fftn(fft.fftshift(Fhkl))/cell[6]          #4D map
3096    rho = fft.fftn(fft.fftshift(Fhkl[:,:,:,maxM+1]))/cell[6]    #3D map
3097    map4DData['rho'] = np.real(SSrho)
3098    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3099    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3100    map4DData['Type'] = reflDict['Type']
3101    mapData['Type'] = reflDict['Type']
3102    mapData['rho'] = np.real(rho)
3103    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3104    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3105    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
3106
3107# map printing for testing purposes
3108def printRho(SGLaue,rho,rhoMax):                         
3109    '''default doc string
3110   
3111    :param type name: description
3112   
3113    :returns: type name: description
3114   
3115    '''
3116    dim = len(rho.shape)
3117    if dim == 2:
3118        ix,jy = rho.shape
3119        for j in range(jy):
3120            line = ''
3121            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3122                line += (jy-j)*'  '
3123            for i in range(ix):
3124                r = int(100*rho[i,j]/rhoMax)
3125                line += '%4d'%(r)
3126            print line+'\n'
3127    else:
3128        ix,jy,kz = rho.shape
3129        for k in range(kz):
3130            print 'k = ',k
3131            for j in range(jy):
3132                line = ''
3133                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
3134                    line += (jy-j)*'  '
3135                for i in range(ix):
3136                    r = int(100*rho[i,j,k]/rhoMax)
3137                    line += '%4d'%(r)
3138                print line+'\n'
3139## keep this
3140               
3141def findOffset(SGData,A,Fhkl):   
3142    '''default doc string
3143   
3144    :param type name: description
3145   
3146    :returns: type name: description
3147   
3148    '''
3149    if SGData['SpGrp'] == 'P 1':
3150        return [0,0,0]   
3151    hklShape = Fhkl.shape
3152    hklHalf = np.array(hklShape)/2
3153    sortHKL = np.argsort(Fhkl.flatten())
3154    Fdict = {}
3155    for hkl in sortHKL:
3156        HKL = np.unravel_index(hkl,hklShape)
3157        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
3158        if F == 0.:
3159            break
3160        Fdict['%.6f'%(np.absolute(F))] = hkl
3161    Flist = np.flipud(np.sort(Fdict.keys()))
3162    F = str(1.e6)
3163    i = 0
3164    DH = []
3165    Dphi = []
3166    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3167    for F in Flist:
3168        hkl = np.unravel_index(Fdict[F],hklShape)
3169        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
3170            continue
3171        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
3172        Uniq = np.array(Uniq,dtype='i')
3173        Phi = np.array(Phi)
3174        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
3175        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
3176        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
3177        ang0 = np.angle(Fh0,deg=True)/360.
3178        for H,phi in zip(Uniq,Phi)[1:]:
3179            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
3180            dH = H-hkl
3181            dang = ang-ang0
3182            DH.append(dH)
3183            Dphi.append((dang+.5) % 1.0)
3184        if i > 20 or len(DH) > 30:
3185            break
3186        i += 1
3187    DH = np.array(DH)
3188    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
3189    Dphi = np.array(Dphi)
3190    steps = np.array(hklShape)
3191    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
3192    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
3193    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
3194    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3195    hist,bins = np.histogram(Mmap,bins=1000)
3196#    for i,item in enumerate(hist[:10]):
3197#        print item,bins[i]
3198    chisq = np.min(Mmap)
3199    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3200    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
3201#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
3202    return DX
3203   
3204def ChargeFlip(data,reflDict,pgbar):
3205    '''default doc string
3206   
3207    :param type name: description
3208   
3209    :returns: type name: description
3210   
3211    '''
3212    generalData = data['General']
3213    mapData = generalData['Map']
3214    flipData = generalData['Flip']
3215    FFtable = {}
3216    if 'None' not in flipData['Norm element']:
3217        normElem = flipData['Norm element'].upper()
3218        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3219        for ff in FFs:
3220            if ff['Symbol'] == normElem:
3221                FFtable.update(ff)
3222    dmin = flipData['Resolution']
3223    SGData = generalData['SGData']
3224    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
3225    SGT = np.array([ops[1] for ops in SGData['SGOps']])
3226    cell = generalData['Cell'][1:8]       
3227    A = G2lat.cell2A(cell[:6])
3228    Vol = cell[6]
3229    im = 0
3230    if generalData['Type'] in ['modulated','magnetic',]:
3231        im = 1
3232    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
3233    adjHKLmax(SGData,Hmax)
3234    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3235    time0 = time.time()
3236    for iref,ref in enumerate(reflDict['RefList']):
3237        dsp = ref[4+im]
3238        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
3239            continue
3240        if dsp > dmin:
3241            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3242            if FFtable:
3243                SQ = 0.25/dsp**2
3244                ff *= G2el.ScatFac(FFtable,SQ)[0]
3245            if ref[8+im] > 0.:         #use only +ve Fobs**2
3246                E = np.sqrt(ref[8+im])/ff
3247            else:
3248                E = 0.
3249            ph = ref[10]
3250            ph = rn.uniform(0.,360.)
3251            Uniq = np.inner(ref[:3],SGMT)
3252            Phi = np.inner(ref[:3],SGT)
3253            for i,hkl in enumerate(Uniq):        #uses uniq
3254                hkl = np.asarray(hkl,dtype='i')
3255                dp = 360.*Phi[i]                #and phi
3256                a = cosd(ph+dp)
3257                b = sind(ph+dp)
3258                phasep = complex(a,b)
3259                phasem = complex(a,-b)
3260                h,k,l = hkl+Hmax
3261                Ehkl[h,k,l] = E*phasep
3262                h,k,l = -hkl+Hmax
3263                Ehkl[h,k,l] = E*phasem
3264#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3265    testHKL = np.array(flipData['testHKL'])+Hmax
3266    CEhkl = copy.copy(Ehkl)
3267    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3268    Emask = ma.getmask(MEhkl)
3269    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3270    Ncyc = 0
3271    old = np.seterr(all='raise')
3272    twophases = []
3273    while True:       
3274        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3275        CEsig = np.std(CErho)
3276        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3277        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3278        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3279        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3280        phase = CFhkl/np.absolute(CFhkl)
3281        twophases.append([np.angle(phase[h,k,l]) for h,k,l in testHKL])
3282        CEhkl = np.absolute(Ehkl)*phase
3283        Ncyc += 1
3284        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3285        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3286        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3287        if Rcf < 5.:
3288            break
3289        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3290        if not GoOn or Ncyc > 10000:
3291            break
3292    np.seterr(**old)
3293    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
3294    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.  #? to get on same scale as e-map
3295    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
3296    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3297       
3298    mapData['Rcf'] = Rcf
3299    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3300    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3301    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3302    mapData['Type'] = reflDict['Type']
3303    return mapData,twophases
3304   
3305def findSSOffset(SGData,SSGData,A,Fhklm):   
3306    '''default doc string
3307   
3308    :param type name: description
3309   
3310    :returns: type name: description
3311   
3312    '''
3313    if SGData['SpGrp'] == 'P 1':
3314        return [0,0,0,0]   
3315    hklmShape = Fhklm.shape
3316    hklmHalf = np.array(hklmShape)/2
3317    sortHKLM = np.argsort(Fhklm.flatten())
3318    Fdict = {}
3319    for hklm in sortHKLM:
3320        HKLM = np.unravel_index(hklm,hklmShape)
3321        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
3322        if F == 0.:
3323            break
3324        Fdict['%.6f'%(np.absolute(F))] = hklm
3325    Flist = np.flipud(np.sort(Fdict.keys()))
3326    F = str(1.e6)
3327    i = 0
3328    DH = []
3329    Dphi = []
3330    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3331    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3332    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
3333    for F in Flist:
3334        hklm = np.unravel_index(Fdict[F],hklmShape)
3335        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
3336            continue
3337        Uniq = np.inner(hklm-hklmHalf,SSGMT)
3338        Phi = np.inner(hklm-hklmHalf,SSGT)
3339        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
3340        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
3341        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
3342        ang0 = np.angle(Fh0,deg=True)/360.
3343        for H,phi in zip(Uniq,Phi)[1:]:
3344            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
3345            dH = H-hklm
3346            dang = ang-ang0
3347            DH.append(dH)
3348            Dphi.append((dang+.5) % 1.0)
3349        if i > 20 or len(DH) > 30:
3350            break
3351        i += 1
3352    DH = np.array(DH)
3353    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
3354    Dphi = np.array(Dphi)
3355    steps = np.array(hklmShape)
3356    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]]
3357    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
3358    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
3359    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
3360    hist,bins = np.histogram(Mmap,bins=1000)
3361#    for i,item in enumerate(hist[:10]):
3362#        print item,bins[i]
3363    chisq = np.min(Mmap)
3364    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
3365    print ' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
3366#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
3367    return DX
3368   
3369def SSChargeFlip(data,reflDict,pgbar):
3370    '''default doc string
3371   
3372    :param type name: description
3373   
3374    :returns: type name: description
3375   
3376    '''
3377    generalData = data['General']
3378    mapData = generalData['Map']
3379    map4DData = {}
3380    flipData = generalData['Flip']
3381    FFtable = {}
3382    if 'None' not in flipData['Norm element']:
3383        normElem = flipData['Norm element'].upper()
3384        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
3385        for ff in FFs:
3386            if ff['Symbol'] == normElem:
3387                FFtable.update(ff)
3388    dmin = flipData['Resolution']
3389    SGData = generalData['SGData']
3390    SSGData = generalData['SSGData']
3391    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
3392    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
3393    cell = generalData['Cell'][1:8]       
3394    A = G2lat.cell2A(cell[:6])
3395    Vol = cell[6]
3396    maxM = 4
3397    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
3398    adjHKLmax(SGData,Hmax)
3399    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
3400    time0 = time.time()
3401    for iref,ref in enumerate(reflDict['RefList']):
3402        dsp = ref[5]
3403        if dsp > dmin:
3404            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
3405            if FFtable:
3406                SQ = 0.25/dsp**2
3407                ff *= G2el.ScatFac(FFtable,SQ)[0]
3408            if ref[9] > 0.:         #use only +ve Fobs**2
3409                E = np.sqrt(ref[9])/ff
3410            else:
3411                E = 0.
3412            ph = ref[11]
3413            ph = rn.uniform(0.,360.)
3414            Uniq = np.inner(ref[:4],SSGMT)
3415            Phi = np.inner(ref[:4],SSGT)
3416            for i,hklm in enumerate(Uniq):        #uses uniq
3417                hklm = np.asarray(hklm,dtype='i')
3418                dp = 360.*Phi[i]                #and phi
3419                a = cosd(ph+dp)
3420                b = sind(ph+dp)
3421                phasep = complex(a,b)
3422                phasem = complex(a,-b)
3423                h,k,l,m = hklm+Hmax
3424                Ehkl[h,k,l,m] = E*phasep
3425                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
3426                Ehkl[h,k,l,m] = E*phasem
3427#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
3428    CEhkl = copy.copy(Ehkl)
3429    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
3430    Emask = ma.getmask(MEhkl)
3431    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
3432    Ncyc = 0
3433    old = np.seterr(all='raise')
3434    while True:       
3435        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
3436        CEsig = np.std(CErho)
3437        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
3438        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
3439        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
3440        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
3441        phase = CFhkl/np.absolute(CFhkl)
3442        CEhkl = np.absolute(Ehkl)*phase
3443        Ncyc += 1
3444        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
3445        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
3446        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
3447        if Rcf < 5.:
3448            break
3449        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
3450        if not GoOn or Ncyc > 10000:
3451            break
3452    np.seterr(**old)
3453    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
3454    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))/10.    #? to get on same scale as e-map
3455    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))/10.                  #? ditto
3456    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
3457    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
3458
3459    mapData['Rcf'] = Rcf
3460    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3461    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
3462    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
3463    mapData['Type'] = reflDict['Type']
3464
3465    map4DData['Rcf'] = Rcf
3466    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))
3467    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
3468    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
3469    map4DData['Type'] = reflDict['Type']
3470    return mapData,map4DData
3471   
3472def getRho(xyz,mapData):
3473    ''' get scattering density at a point by 8-point interpolation
3474    param xyz: coordinate to be probed
3475    param: mapData: dict of map data
3476   
3477    :returns: density at xyz
3478    '''
3479    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3480    if not len(mapData):
3481        return 0.0
3482    rho = copy.copy(mapData['rho'])     #don't mess up original
3483    if not len(rho):
3484        return 0.0
3485    mapShape = np.array(rho.shape)
3486    mapStep = 1./mapShape
3487    X = np.array(xyz)%1.    #get into unit cell
3488    I = np.array(X*mapShape,dtype='int')
3489    D = X-I*mapStep         #position inside map cell
3490    D12 = D[0]*D[1]
3491    D13 = D[0]*D[2]
3492    D23 = D[1]*D[2]
3493    D123 = np.prod(D)
3494    Rho = rollMap(rho,-I)    #shifts map so point is in corner
3495    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]+  \
3496        Rho[1,1,1]*D123+Rho[0,1,1]*(D23-D123)+Rho[1,0,1]*(D13-D123)+Rho[1,1,0]*(D12-D123)+  \
3497        Rho[0,0,0]*(D12+D13+D23-D123)-Rho[0,0,1]*(D13+D23-D123)-    \
3498        Rho[0,1,0]*(D23+D12-D123)-Rho[1,0,0]*(D13+D12-D123)   
3499    return R
3500       
3501def SearchMap(generalData,drawingData,Neg=False):
3502    '''Does a search of a density map for peaks meeting the criterion of peak
3503    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
3504    mapData is data['General']['mapData']; the map is also in mapData.
3505
3506    :param generalData: the phase data structure; includes the map
3507    :param drawingData: the drawing data structure
3508    :param Neg:  if True then search for negative peaks (i.e. H-atoms & neutron data)
3509
3510    :returns: (peaks,mags,dzeros) where
3511
3512        * peaks : ndarray
3513          x,y,z positions of the peaks found in the map
3514        * mags : ndarray
3515          the magnitudes of the peaks
3516        * dzeros : ndarray
3517          the distance of the peaks from  the unit cell origin
3518        * dcent : ndarray
3519          the distance of the peaks from  the unit cell center
3520
3521    '''       
3522    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
3523   
3524    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
3525   
3526#    def noDuplicate(xyz,peaks,Amat):
3527#        XYZ = np.inner(Amat,xyz)
3528#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3529#            print ' Peak',xyz,' <0.5A from another peak'
3530#            return False
3531#        return True
3532#                           
3533    def fixSpecialPos(xyz,SGData,Amat):
3534        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
3535        X = []
3536        xyzs = [equiv[0] for equiv in equivs]
3537        for x in xyzs:
3538            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
3539                X.append(x)
3540        if len(X) > 1:
3541            return np.average(X,axis=0)
3542        else:
3543            return xyz
3544       
3545    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
3546        Mag,x0,y0,z0,sig = parms
3547        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
3548#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
3549        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
3550       
3551    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
3552        Mag,x0,y0,z0,sig = parms
3553        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3554        return M
3555       
3556    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
3557        Mag,x0,y0,z0,sig = parms
3558        dMdv = np.zeros(([5,]+list(rX.shape)))
3559        delt = .01
3560        for i in range(5):
3561            parms[i] -= delt
3562            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3563            parms[i] += 2.*delt
3564            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3565            parms[i] -= delt
3566            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
3567        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
3568        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
3569        dMdv = np.reshape(dMdv,(5,rX.size))
3570        Hess = np.inner(dMdv,dMdv)
3571       
3572        return Vec,Hess
3573       
3574    SGData = generalData['SGData']
3575    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3576    peaks = []
3577    mags = []
3578    dzeros = []
3579    dcent = []
3580    try:
3581        mapData = generalData['Map']
3582        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
3583        if Neg:
3584            rho = -copy.copy(mapData['rho'])     #flip +/-
3585        else:
3586            rho = copy.copy(mapData['rho'])     #don't mess up original
3587        mapHalf = np.array(rho.shape)/2
3588        res = mapData['Resolution']
3589        incre = np.array(rho.shape,dtype=np.float)
3590        step = max(1.0,1./res)+1
3591        steps = np.array((3*[step,]),dtype='int32')
3592    except KeyError:
3593        print '**** ERROR - Fourier map not defined'
3594        return peaks,mags
3595    rhoMask = ma.array(rho,mask=(rho<contLevel))
3596    indices = (-1,0,1)
3597    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
3598    for roll in rolls:
3599        if np.any(roll):
3600            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
3601    indx = np.transpose(rhoMask.nonzero())
3602    peaks = indx/incre
3603    mags = rhoMask[rhoMask.nonzero()]
3604    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
3605        rho = rollMap(rho,ind)
3606        rMM = mapHalf-steps
3607        rMP = mapHalf+steps+1
3608        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
3609        peakInt = np.sum(rhoPeak)*res**3
3610        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
3611        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
3612        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
3613            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
3614        x1 = result[0]
3615        if not np.any(x1 < 0):
3616            peak = (np.array(x1[1:4])-ind)/incre
3617        peak = fixSpecialPos(peak,SGData,Amat)
3618        rho = rollMap(rho,-ind)
3619    cent = np.ones(3)*.5     
3620    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
3621    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
3622    if Neg:     #want negative magnitudes for negative peaks
3623        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3624    else:
3625        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
3626   
3627def sortArray(data,pos,reverse=False):
3628    '''data is a list of items
3629    sort by pos in list; reverse if True
3630    '''
3631    T = []
3632    for i,M in enumerate(data):
3633        try:
3634            T.append((M[pos],i))
3635        except IndexError:
3636            return data
3637    D = dict(zip(T,data))
3638    T.sort()
3639    if reverse:
3640        T.reverse()
3641    X = []
3642    for key in T:
3643        X.append(D[key])
3644    return X
3645
3646def PeaksEquiv(data,Ind):
3647    '''Find the equivalent map peaks for those selected. Works on the
3648    contents of data['Map Peaks'].
3649
3650    :param data: the phase data structure
3651    :param list Ind: list of selected peak indices
3652    :returns: augmented list of peaks including those related by symmetry to the
3653      ones in Ind
3654
3655    '''       
3656    def Duplicate(xyz,peaks,Amat):
3657        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3658            return True
3659        return False
3660                           
3661    generalData = data['General']
3662    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3663    SGData = generalData['SGData']
3664    mapPeaks = data['Map Peaks']
3665    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
3666    Indx = {}
3667    for ind in Ind:
3668        xyz = np.array(mapPeaks[ind][1:4])
3669        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
3670        for jnd,xyz in enumerate(XYZ):       
3671            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
3672    Ind = []
3673    for ind in Indx:
3674        if Indx[ind]:
3675            Ind.append(ind)
3676    return Ind
3677               
3678def PeaksUnique(data,Ind):
3679    '''Finds the symmetry unique set of peaks from those selected. Works on the
3680    contents of data['Map Peaks'].
3681
3682    :param data: the phase data structure
3683    :param list Ind: list of selected peak indices
3684    :returns: the list of symmetry unique peaks from among those given in Ind
3685
3686    '''       
3687#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
3688
3689    def noDuplicate(xyz,peaks,Amat):
3690        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
3691            return False
3692        return True
3693                           
3694    generalData = data['General']
3695    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3696    SGData = generalData['SGData']
3697    mapPeaks = data['Map Peaks']
3698    Indx = {}
3699    XYZ = {}
3700    for ind in Ind:
3701        XYZ[ind] = np.array(mapPeaks[ind][1:4])
3702        Indx[ind] = True
3703    for ind in Ind:
3704        if Indx[ind]:
3705            xyz = XYZ[ind]
3706            for jnd in Ind:
3707                if ind != jnd and Indx[jnd]:                       
3708                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
3709                    xyzs = np.array([equiv[0] for equiv in Equiv])
3710                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
3711    Ind = []
3712    for ind in Indx:
3713        if Indx[ind]:
3714            Ind.append(ind)
3715    return Ind
3716   
3717################################################################################
3718##### single peak fitting profile fxn stuff
3719################################################################################
3720
3721def getCWsig(ins,pos):
3722    '''get CW peak profile sigma^2
3723   
3724    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
3725      as values only
3726    :param float pos: 2-theta of peak
3727    :returns: float getCWsig: peak sigma^2
3728   
3729    '''
3730    tp = tand(pos/2.0)
3731    return ins['U']*tp**2+ins['V']*tp+ins['W']
3732   
3733def getCWsigDeriv(pos):
3734    '''get derivatives of CW peak profile sigma^2 wrt U,V, & W
3735   
3736    :param float pos: 2-theta of peak
3737   
3738    :returns: list getCWsigDeriv: d(sig^2)/dU, d(sig)/dV & d(sig)/dW
3739   
3740    '''
3741    tp = tand(pos/2.0)
3742    return tp**2,tp,1.0
3743   
3744def getCWgam(ins,pos):
3745    '''get CW peak profile gamma
3746   
3747    :param dict ins: instrument parameters with at least 'X' & 'Y'
3748      as values only
3749    :param float pos: 2-theta of peak
3750    :returns: float getCWgam: peak gamma
3751   
3752    '''
3753    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
3754   
3755def getCWgamDeriv(pos):
3756    '''get derivatives of CW peak profile gamma wrt X & Y
3757   
3758    :param float pos: 2-theta of peak
3759   
3760    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
3761   
3762    '''
3763    return 1./cosd(pos/2.0),tand(pos/2.0)
3764   
3765def getTOFsig(ins,dsp):
3766    '''get TOF peak profile sigma^2
3767   
3768    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
3769      as values only
3770    :param float dsp: d-spacing of peak
3771   
3772    :returns: float getTOFsig: peak sigma^2
3773   
3774    '''
3775    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
3776   
3777def getTOFsigDeriv(dsp):
3778    '''get derivatives of TOF peak profile sigma^2 wrt sig-0, sig-1, & sig-q
3779   
3780    :param float dsp: d-spacing of peak
3781   
3782    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
3783   
3784    '''
3785    return 1.0,dsp**2,dsp**4,1./dsp**2
3786   
3787def getTOFgamma(ins,dsp):
3788    '''get TOF peak profile gamma
3789   
3790    :param dict ins: instrument parameters with at least 'X' & 'Y'
3791      as values only
3792    :param float dsp: d-spacing of peak
3793   
3794    :returns: float getTOFgamma: peak gamma
3795   
3796    '''
3797    return ins['X']*dsp+ins['Y']*dsp**2
3798   
3799def getTOFgammaDeriv(dsp):
3800    '''get derivatives of TOF peak profile gamma wrt X & Y
3801   
3802    :param float dsp: d-spacing of peak
3803   
3804    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
3805   
3806    '''
3807    return dsp,dsp**2
3808   
3809def getTOFbeta(ins,dsp):
3810    '''get TOF peak profile beta
3811   
3812    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
3813      as values only
3814    :param float dsp: d-spacing of peak
3815   
3816    :returns: float getTOFbeta: peak beat
3817   
3818    '''
3819    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
3820   
3821def getTOFbetaDeriv(dsp):
3822    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
3823   
3824    :param float dsp: d-spacing of peak
3825   
3826    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
3827   
3828    '''
3829    return 1.0,1./dsp**4,1./dsp**2
3830   
3831def getTOFalpha(ins,dsp):
3832    '''get TOF peak profile alpha
3833   
3834    :param dict ins: instrument parameters with at least 'alpha'
3835      as values only
3836    :param float dsp: d-spacing of peak
3837   
3838    :returns: flaot getTOFalpha: peak alpha
3839   
3840    '''
3841    return ins['alpha']/dsp
3842   
3843def getTOFalphaDeriv(dsp):
3844    '''get derivatives of TOF peak profile beta wrt alpha
3845   
3846    :param float dsp: d-spacing of peak
3847   
3848    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
3849   
3850    '''
3851    return 1./dsp
3852   
3853def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
3854    '''set starting peak parameters for single peak fits from plot selection or auto selection
3855   
3856    :param dict Parms: instrument parameters dictionary
3857    :param dict Parms2: table lookup for TOF profile coefficients
3858    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
3859    :param float mag: peak top magnitude from pick
3860    :param bool ifQ: True if pos in Q
3861    :param bool useFit: True if use fitted CW Parms values (not defaults)
3862   
3863    :returns: list XY: peak list entry:
3864        for CW: [pos,0,mag,1,sig,0,gam,0]
3865        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3866        NB: mag refinement set by default, all others off
3867   
3868    '''
3869    ind = 0
3870    if useFit:
3871        ind = 1
3872    ins = {}
3873    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
3874        for x in ['U','V','W','X','Y']:
3875            ins[x] = Parms[x][ind]
3876        if ifQ:                              #qplot - convert back to 2-theta
3877            pos = 2.0*asind(pos*getWave(Parms)/(4*math.pi))
3878        sig = getCWsig(ins,pos)
3879        gam = getCWgam(ins,pos)           
3880        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
3881    else:
3882        if ifQ:
3883            dsp = 2.*np.pi/pos
3884            pos = Parms['difC']*dsp
3885        else:
3886            dsp = pos/Parms['difC'][1]
3887        if 'Pdabc' in Parms2:
3888            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
3889                ins[x] = Parms[x][ind]
3890            Pdabc = Parms2['Pdabc'].T
3891            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
3892            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
3893        else:
3894            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
3895                ins[x] = Parms[x][ind]
3896            alp = getTOFalpha(ins,dsp)
3897            bet = getTOFbeta(ins,dsp)
3898        sig = getTOFsig(ins,dsp)
3899        gam = getTOFgamma(ins,dsp)
3900        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
3901    return XY
3902   
3903################################################################################
3904##### MC/SA stuff
3905################################################################################
3906
3907#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
3908# Original Author: Travis Oliphant 2002
3909# Bug-fixes in 2006 by Tim Leslie
3910
3911
3912import numpy
3913from numpy import asarray, exp, squeeze, sign, \
3914     all, shape, array, where
3915from numpy import random
3916
3917#__all__ = ['anneal']
3918
3919_double_min = numpy.finfo(float).min
3920_double_max = numpy.finfo(float).max
3921class base_schedule(object):
3922    def __init__(self):
3923        self.dwell = 20
3924        self.lower = 0.
3925        self.upper = 1.
3926        self.Ninit = 50
3927        self.accepted = 0
3928        self.tests = 0
3929        self.feval = 0
3930        self.k = 0
3931        self.T = None
3932
3933    def init(self, **options):
3934        self.__dict__.update(options)
3935        self.lower = asarray(self.lower)
3936        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
3937        self.upper = asarray(self.upper)
3938        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
3939        self.k = 0
3940        self.accepted = 0
3941        self.feval = 0
3942        self.tests = 0
3943
3944    def getstart_temp(self, best_state):
3945        """ Find a matching starting temperature and starting parameters vector
3946        i.e. find x0 such that func(x0) = T0.
3947
3948        :param best_state: _state
3949            A _state object to store the function value and x0 found.
3950
3951        :returns: x0 : array
3952            The starting parameters vector.
3953        """
3954
3955        assert(not self.dims is None)
3956        lrange = self.lower
3957        urange = self.upper
3958        fmax = _double_min
3959        fmin = _double_max
3960        for _ in range(self.Ninit):
3961            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
3962            fval = self.func(x0, *self.args)
3963            self.feval += 1
3964            if fval > fmax:
3965                fmax = fval
3966            if fval < fmin:
3967                fmin = fval
3968                best_state.cost = fval
3969                best_state.x = array(x0)
3970
3971        self.T0 = (fmax-fmin)*1.5
3972        return best_state.x
3973       
3974    def set_range(self,x0,frac):
3975        delrange = frac*(self.upper-self.lower)
3976        self.upper = x0+delrange
3977        self.lower = x0-delrange
3978
3979    def accept_test(self, dE):
3980        T = self.T
3981        self.tests += 1
3982        if dE < 0:
3983            self.accepted += 1
3984            return 1
3985        p = exp(-dE*1.0/T)
3986        if (p > random.uniform(0.0, 1.0)):
3987            self.accepted += 1
3988            return 1
3989        return 0
3990
3991    def update_guess(self, x0):
3992        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
3993
3994    def update_temp(self, x0):
3995        pass
3996
3997class fast_sa(base_schedule):
3998    def init(self, **options):
3999        self.__dict__.update(options)
4000
4001    def update_guess(self, x0):
4002        x0 = asarray(x0)
4003        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4004        T = self.T
4005        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
4006        xnew = xc*(self.upper - self.lower)+self.lower
4007        return xnew
4008
4009    def update_temp(self):
4010        self.T = self.T0*exp(-self.c * self.k**(self.quench))
4011        self.k += 1
4012        return
4013
4014class log_sa(base_schedule):        #OK
4015
4016    def init(self,**options):
4017        self.__dict__.update(options)
4018       
4019    def update_guess(self,x0):     #same as default #TODO - is this a reasonable update procedure?
4020        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
4021        T = self.T
4022        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
4023        xnew = xc*(self.upper - self.lower)+self.lower
4024        return xnew
4025       
4026    def update_temp(self):
4027        self.k += 1
4028        self.T = self.T0*self.slope**self.k
4029       
4030class _state(object):
4031    def __init__(self):
4032        self.x = None
4033        self.cost = None
4034
4035def makeTsched(data):
4036    if data['Algorithm'] == 'fast':
4037        sched = fast_sa()
4038        sched.quench = data['fast parms'][0]
4039        sched.c = data['fast parms'][1]
4040    elif data['Algorithm'] == 'log':
4041        sched = log_sa()
4042        sched.slope = data['log slope']
4043    sched.T0 = data['Annealing'][0]
4044    if not sched.T0:
4045        sched.T0 = 50.
4046    Tf = data['Annealing'][1]
4047    if not Tf:
4048        Tf = 0.001
4049    Tsched = [sched.T0,]
4050    while Tsched[-1] > Tf:
4051        sched.update_temp()
4052        Tsched.append(sched.T)
4053    return Tsched[1:]
4054   
4055def anneal(func, x0, args=(), schedule='fast', 
4056           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
4057           feps=1e-6, quench=1.0, c=1.0,
4058           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
4059           ranRange=0.10,autoRan=False,dlg=None):
4060    """Minimize a function using simulated annealing.
4061
4062    Schedule is a schedule class implementing the annealing schedule.
4063    Available ones are 'fast', 'cauchy', 'boltzmann'
4064
4065    :param callable func: f(x, \*args)
4066        Function to be optimized.
4067    :param ndarray x0:
4068        Initial guess.
4069    :param tuple args:
4070        Extra parameters to `func`.
4071    :param base_schedule schedule:
4072        Annealing schedule to use (a class).
4073    :param float T0:
4074        Initial Temperature (estimated as 1.2 times the largest
4075        cost-function deviation over random points in the range).
4076    :param float Tf:
4077        Final goal temperature.
4078    :param int maxeval:
4079        Maximum function evaluations.
4080    :param int maxaccept:
4081        Maximum changes to accept.
4082    :param int maxiter:
4083        Maximum cooling iterations.
4084    :param float feps:
4085        Stopping relative error tolerance for the function value in
4086        last four coolings.
4087    :param float quench,c:
4088        Parameters to alter fast_sa schedule.
4089    :param float/ndarray lower,upper:
4090        Lower and upper bounds on `x`.
4091    :param int dwell:
4092        The number of times to search the space at each temperature.
4093    :param float slope:
4094        Parameter for log schedule
4095    :param bool ranStart=False:
4096        True for set 10% of ranges about x
4097
4098    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
4099
4100     * xmin (ndarray): Point giving smallest value found.
4101     * Jmin (float): Minimum value of function found.
4102     * T (float): Final temperature.
4103     * feval (int): Number of function evaluations.
4104     * iters (int): Number of cooling iterations.
4105     * accept (int): Number of tests accepted.
4106     * retval (int): Flag indicating stopping condition:
4107
4108              *  0: Points no longer changing
4109              *  1: Cooled to final temperature
4110              *  2: Maximum function evaluations
4111              *  3: Maximum cooling iterations reached
4112              *  4: Maximum accepted query locations reached
4113              *  5: Final point not the minimum amongst encountered points
4114
4115    *Notes*:
4116    Simulated annealing is a random algorithm which uses no derivative
4117    information from the function being optimized. In practice it has
4118    been more useful in discrete optimization than continuous
4119    optimization, as there are usually better algorithms for continuous
4120    optimization problems.
4121
4122    Some experimentation by trying the difference temperature
4123    schedules and altering their parameters is likely required to
4124    obtain good performance.
4125
4126    The randomness in the algorithm comes from random sampling in numpy.
4127    To obtain the same results you can call numpy.random.seed with the
4128    same seed immediately before calling scipy.optimize.anneal.
4129
4130    We give a brief description of how the three temperature schedules
4131    generate new points and vary their temperature. Temperatures are
4132    only updated with iterations in the outer loop. The inner loop is
4133    over xrange(dwell), and new points are generated for
4134    every iteration in the inner loop. (Though whether the proposed
4135    new points are accepted is probabilistic.)
4136
4137    For readability, let d denote the dimension of the inputs to func.
4138    Also, let x_old denote the previous state, and k denote the
4139    iteration number of the outer loop. All other variables not
4140    defined below are input variables to scipy.optimize.anneal itself.
4141
4142    In the 'fast' schedule the updates are ::
4143
4144        u ~ Uniform(0, 1, size=d)
4145        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
4146        xc = y * (upper - lower)
4147        x_new = x_old + xc
4148
4149        T_new = T0 * exp(-c * k**quench)
4150
4151
4152    """
4153    x0 = asarray(x0)
4154    lower = asarray(lower)
4155    upper = asarray(upper)
4156
4157    schedule = eval(schedule+'_sa()')
4158    #   initialize the schedule
4159    schedule.init(dims=shape(x0),func=func,args=args,T0=T0,lower=lower, upper=upper,
4160        c=c, quench=quench, dwell=dwell, slope=slope)
4161
4162    current_state, last_state, best_state = _state(), _state(), _state()
4163    if ranStart:
4164        schedule.set_range(x0,ranRange)
4165    if T0 is None:
4166        x0 = schedule.getstart_temp(best_state)
4167    else:
4168        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
4169        best_state.x = None
4170        best_state.cost = numpy.Inf
4171
4172    last_state.x = asarray(x0).copy()
4173    fval = func(x0,*args)
4174    schedule.feval += 1
4175    last_state.cost = fval
4176    if last_state.cost < best_state.cost:
4177        best_state.cost = fval
4178        best_state.x = asarray(x0).copy()
4179    schedule.T = schedule.T0
4180    fqueue = [100, 300, 500, 700]
4181    iters = 1
4182    keepGoing = True
4183    bestn = 0
4184    while keepGoing:
4185        retval = 0
4186        for n in xrange(dwell):
4187            current_state.x = schedule.update_guess(last_state.x)
4188            current_state.cost = func(current_state.x,*args)
4189            schedule.feval += 1
4190
4191            dE = current_state.cost - last_state.cost
4192            if schedule.accept_test(dE):
4193                last_state.x = current_state.x.copy()
4194                last_state.cost = current_state.cost
4195                if last_state.cost < best_state.cost:
4196                    best_state.x = last_state.x.copy()
4197                    best_state.cost = last_state.cost
4198                    bestn = n
4199                    if best_state.cost < 1.0 and autoRan:
4200                        schedule.set_range(x0,best_state.cost/2.)                       
4201        if dlg:
4202            GoOn = dlg.Update(min(100.,best_state.cost*100),
4203                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
4204                    'Best trial:',bestn,  \
4205                    'MC/SA Residual:',best_state.cost*100,'%', \
4206                    ))[0]
4207            if not GoOn:
4208                best_state.x = last_state.x.copy()
4209                best_state.cost = last_state.cost
4210                retval = 5
4211        schedule.update_temp()
4212        iters += 1
4213        # Stopping conditions
4214        # 0) last saved values of f from each cooling step
4215        #     are all very similar (effectively cooled)
4216        # 1) Tf is set and we are below it
4217        # 2) maxeval is set and we are past it
4218        # 3) maxiter is set and we are past it
4219        # 4) maxaccept is set and we are past it
4220        # 5) user canceled run via progress bar
4221
4222        fqueue.append(squeeze(last_state.cost))
4223        fqueue.pop(0)
4224        af = asarray(fqueue)*1.0
4225        if retval == 5:
4226            print ' User terminated run; incomplete MC/SA'
4227            keepGoing = False
4228            break
4229        if all(abs((af-af[0])/af[0]) < feps):
4230            retval = 0
4231            if abs(af[-1]-best_state.cost) > feps*10:
4232                retval = 5
4233                print " Warning: Cooled to %.4f > selected Tmin %.4f in %d steps"%(squeeze(last_state.cost),Tf,iters-1)
4234            break
4235        if (Tf is not None) and (schedule.T < Tf):
4236#            print ' Minimum T reached in %d steps'%(iters-1)
4237            retval = 1
4238            break
4239        if (maxeval is not None) and (schedule.feval > maxeval):
4240            retval = 2
4241            break
4242        if (iters > maxiter):
4243            print  " Warning: Maximum number of iterations exceeded."
4244            retval = 3
4245            break
4246        if (maxaccept is not None) and (schedule.accepted > maxaccept):
4247            retval = 4
4248            break
4249
4250    return best_state.x, best_state.cost, schedule.T, \
4251           schedule.feval, iters, schedule.accepted, retval
4252
4253def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,nprocess=-1):
4254    outlist = []
4255    timelist = []
4256    nsflist = []
4257    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
4258    for n in range(iCyc):
4259        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None,False)         #mcsa result,time,rcov
4260        outlist.append(result[0])
4261        timelist.append(result[1])
4262        nsflist.append(result[2])
4263        print ' MC/SA final fit: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
4264    out_q.put(outlist)
4265    out_t.put(timelist)
4266    out_n.put(nsflist)
4267
4268def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData,nprocs):
4269    import multiprocessing as mp
4270   
4271    out_q = mp.Queue()
4272    out_t = mp.Queue()
4273    out_n = mp.Queue()
4274    procs = []
4275    totsftime = 0.
4276    totnsf = 0
4277    iCyc = np.zeros(nprocs)
4278    for i in range(nCyc):
4279        iCyc[i%nprocs] += 1
4280    for i in range(nprocs):
4281        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,out_t,out_n,i))
4282        procs.append(p)
4283        p.start()
4284    resultlist = []
4285    for i in range(nprocs):
4286        resultlist += out_q.get()
4287        totsftime += np.sum(out_t.get())
4288        totnsf += np.sum(out_n.get())
4289    for p in procs:
4290        p.join()
4291    return resultlist,totsftime,totnsf
4292
4293def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar,start=True):
4294    '''default doc string
4295   
4296    :param type name: description
4297   
4298    :returns: type name: description
4299    '''
4300   
4301    class RandomDisplacementBounds(object):
4302        """random displacement with bounds"""
4303        def __init__(self, xmin, xmax, stepsize=0.5):
4304            self.xmin = xmin
4305            self.xmax = xmax
4306            self.stepsize = stepsize
4307   
4308        def __call__(self, x):
4309            """take a random step but ensure the new position is within the bounds"""
4310            while True:
4311                # this could be done in a much more clever way, but it will work for example purposes
4312                steps = self.xmax-self.xmin
4313                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
4314                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
4315                    break
4316            return xnew
4317   
4318    global tsum,nsum
4319    tsum = 0.
4320    nsum = 0
4321   
4322    def getMDparms(item,pfx,parmDict,varyList):
4323        parmDict[pfx+'MDaxis'] = item['axis']
4324        parmDict[pfx+'MDval'] = item['Coef'][0]
4325        if item['Coef'][1]:
4326            varyList += [pfx+'MDval',]
4327            limits = item['Coef'][2]
4328            lower.append(limits[0])
4329            upper.append(limits[1])
4330                       
4331    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
4332        parmDict[pfx+'Atype'] = item['atType']
4333        aTypes |= set([item['atType'],]) 
4334        pstr = ['Ax','Ay','Az']
4335        XYZ = [0,0,0]
4336        for i in range(3):
4337            name = pfx+pstr[i]
4338            parmDict[name] = item['Pos'][0][i]
4339            XYZ[i] = parmDict[name]
4340            if item['Pos'][1][i]:
4341                varyList += [name,]
4342                limits = item['Pos'][2][i]
4343                lower.append(limits[0])
4344                upper.append(limits[1])
4345        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
4346           
4347    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
4348        parmDict[mfx+'MolCent'] = item['MolCent']
4349        parmDict[mfx+'RBId'] = item['RBId']
4350        pstr = ['Px','Py','Pz']
4351        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
4352        for i in range(3):
4353            name = mfx+pstr[i]
4354            parmDict[name] = item['Pos'][0][i]
4355            if item['Pos'][1][i]:
4356                varyList += [name,]
4357                limits = item['Pos'][2][i]
4358                lower.append(limits[0])
4359                upper.append(limits[1])
4360        AV = item['Ori'][0]
4361        A = AV[0]
4362        V = AV[1:]
4363        for i in range(4):
4364            name = mfx+ostr[i]
4365            if i:
4366                parmDict[name] = V[i-1]
4367            else:
4368                parmDict[name] = A
4369            if item['Ovar'] == 'AV':
4370                varyList += [name,]
4371                limits = item['Ori'][2][i]
4372                lower.append(limits[0])
4373                upper.append(limits[1])
4374            elif item['Ovar'] == 'A' and not i:
4375                varyList += [name,]
4376                limits = item['Ori'][2][i]
4377                lower.append(limits[0])
4378                upper.append(limits[1])
4379        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
4380            for i in range(len(item['Tor'][0])):
4381                name = mfx+'Tor'+str(i)
4382                parmDict[name] = item['Tor'][0][i]
4383                if item['Tor'][1][i]:
4384                    varyList += [name,]
4385                    limits = item['Tor'][2][i]
4386                    lower.append(limits[0])
4387                    upper.append(limits[1])
4388        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
4389        aTypes |= set(atypes)
4390        atNo += len(atypes)
4391        return atNo
4392               
4393    def GetAtomM(Xdata,SGData):
4394        Mdata = []
4395        for xyz in Xdata:
4396            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
4397        return np.array(Mdata)
4398       
4399    def GetAtomT(RBdata,parmDict):
4400        'Needs a doc string'
4401        atNo = parmDict['atNo']
4402        nfixAt = parmDict['nfixAt']
4403        Tdata = atNo*[' ',]
4404        for iatm in range(nfixAt):
4405            parm = ':'+str(iatm)+':Atype'
4406            if parm in parmDict:
4407                Tdata[iatm] = aTypes.index(parmDict[parm])
4408        iatm = nfixAt
4409        for iObj in range(parmDict['nObj']):
4410            pfx = str(iObj)+':'
4411            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4412                if parmDict[pfx+'Type'] == 'Vector':
4413                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
4414                    nAtm = len(RBRes['rbVect'][0])
4415                else:       #Residue
4416                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
4417                    nAtm = len(RBRes['rbXYZ'])
4418                for i in range(nAtm):
4419                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
4420                    iatm += 1
4421            elif parmDict[pfx+'Type'] == 'Atom':
4422                atNo = parmDict[pfx+'atNo']
4423                parm = pfx+'Atype'              #remove extra ':'
4424                if parm in parmDict:
4425                    Tdata[atNo] = aTypes.index(parmDict[parm])
4426                iatm += 1
4427            else:
4428                continue        #skips March Dollase
4429        return Tdata
4430       
4431    def GetAtomX(RBdata,parmDict):
4432        'Needs a doc string'
4433        Bmat = parmDict['Bmat']
4434        atNo = parmDict['atNo']
4435        nfixAt = parmDict['nfixAt']
4436        Xdata = np.zeros((3,atNo))
4437        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
4438        for iatm in range(nfixAt):
4439            for key in keys:
4440                parm = ':'+str(iatm)+key
4441                if parm in parmDict:
4442                    keys[key][iatm] = parmDict[parm]
4443        iatm = nfixAt
4444        for iObj in range(parmDict['nObj']):
4445            pfx = str(iObj)+':'
4446            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4447                if parmDict[pfx+'Type'] == 'Vector':
4448                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
4449                    vecs = RBRes['rbVect']
4450                    mags = RBRes['VectMag']
4451                    Cart = np.zeros_like(vecs[0])
4452                    for vec,mag in zip(vecs,mags):
4453                        Cart += vec*mag
4454                elif parmDict[pfx+'Type'] == 'Residue':
4455                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
4456                    Cart = np.array(RBRes['rbXYZ'])
4457                    for itor,seq in enumerate(RBRes['rbSeq']):
4458                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
4459                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
4460                if parmDict[pfx+'MolCent'][1]:
4461                    Cart -= parmDict[pfx+'MolCent'][0]
4462                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
4463                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
4464                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
4465                iatm += len(Cart)
4466            elif parmDict[pfx+'Type'] == 'Atom':
4467                atNo = parmDict[pfx+'atNo']
4468                for key in keys:
4469                    parm = pfx+key[1:]              #remove extra ':'
4470                    if parm in parmDict:
4471                        keys[key][atNo] = parmDict[parm]
4472                iatm += 1
4473            else:
4474                continue        #skips March Dollase
4475        return Xdata.T
4476       
4477    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
4478        allX = np.inner(Xdata,SGM)+SGT
4479        allT = np.repeat(Tdata,allX.shape[1])
4480        allM = np.repeat(Mdata,allX.shape[1])
4481        allX = np.reshape(allX,(-1,3))
4482        return allT,allM,allX
4483
4484    def getAllX(Xdata,SGM,SGT):
4485        allX = np.inner(Xdata,SGM)+SGT
4486        allX = np.reshape(allX,(-1,3))
4487        return allX
4488       
4489    def normQuaternions(RBdata,parmDict,varyList,values):
4490        for iObj in range(parmDict['nObj']):
4491            pfx = str(iObj)+':'
4492            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4493                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
4494                A,V = Q2AVdeg(Qori)
4495                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
4496                    if i:
4497                        parmDict[pfx+name] = V[i-1]
4498                    else:
4499                        parmDict[pfx+name] = A
4500       
4501    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
4502        ''' Compute structure factors for all h,k,l for phase
4503        input:
4504            refList: [ref] where each ref = h,k,l,m,d,...
4505            rcov:   array[nref,nref] covariance terms between Fo^2 values
4506            ifInv:  bool True if centrosymmetric
4507            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
4508            RBdata: [dict] rigid body dictionary
4509            varyList: [list] names of varied parameters in MC/SA (not used here)           
4510            ParmDict: [dict] problem parameters
4511        puts result F^2 in each ref[5] in refList
4512        returns:
4513            delt-F*rcov*delt-F/sum(Fo^2)
4514        '''   
4515           
4516        global tsum,nsum
4517        t0 = time.time()
4518        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
4519        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
4520        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
4521        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
4522        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
4523        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
4524        refList[5] = Aterm
4525        if not ifInv:
4526            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
4527        if len(cosTable):        #apply MD correction
4528            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))