source: branch/2frame/GSASIImath.py @ 2937

Last change on this file since 2937 was 2937, checked in by vondreele, 4 years ago

put a self.SetAutoLayout?(True) in SetDataSize?
cleanup anneal - remove boltzmann & cauchy option & parameters - fundamentlly flawed
add new function (makeTsched) to generate anneal T schedule - show plot if coeff changed
add a final 'L-BFGS-B' refinement of each anneal result - gives much improved results
show plot of reflection covariances for MC/SA runs
cleanup MC/SA GUI stuff
cleanup Index list stuff - fix bug as well

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