source: trunk/GSASIImath.py @ 3867

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

modulated mag moments now ok; not sf calc

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