source: trunk/GSASIImath.py @ 3888

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

provide contour plot (not finished yet) on structure drawing

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