source: trunk/GSASIImath.py @ 3103

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

show seq. result in structure drawing - steps thru histograms changing structure at each step

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