# source:trunk/GSASIImath.py@3187

Last change on this file since 3187 was 3187, checked in by toby, 5 years ago

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