source: trunk/GSASIImath.py @ 3864

Last change on this file since 3864 was 3864, checked in by vondreele, 3 years ago

incommensurate mgnetic str. fctr. calcs.latest attempt

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