source: trunk/GSASIImath.py @ 3861

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

another try at incomm. mag. str. fctr.

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