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

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

tweaks of basinhopping
replace more TextCtrl? with ValidatedTxtCtrl? in phaseGUI

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 182.8 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2017-07-12 20:54:39 +0000 (Wed, 12 Jul 2017) $
5# $Author: vondreele $
6# $Revision: 2930 $
7# $URL: branch/2frame/GSASIImath.py $
8# $Id: GSASIImath.py 2930 2017-07-12 20:54:39Z 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: 2930 $")
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, tan, exp, squeeze, sign, \
3886     all, log, pi, 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.learn_rate = 0.5
3897        self.lower = -10
3898        self.upper = 10
3899        self.Ninit = 50
3900        self.accepted = 0
3901        self.tests = 0
3902        self.feval = 0
3903        self.k = 0
3904        self.T = None
3905
3906    def init(self, **options):
3907        self.__dict__.update(options)
3908        self.lower = asarray(self.lower)
3909        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
3910        self.upper = asarray(self.upper)
3911        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
3912        self.k = 0
3913        self.accepted = 0
3914        self.feval = 0
3915        self.tests = 0
3916
3917    def getstart_temp(self, best_state):
3918        """ Find a matching starting temperature and starting parameters vector
3919        i.e. find x0 such that func(x0) = T0.
3920
3921        Parameters
3922        ----------
3923        best_state : _state
3924            A _state object to store the function value and x0 found.
3925
3926        returns
3927        -------
3928        x0 : array
3929            The starting parameters vector.
3930        """
3931
3932        assert(not self.dims is None)
3933        lrange = self.lower
3934        urange = self.upper
3935        fmax = _double_min
3936        fmin = _double_max
3937        for _ in range(self.Ninit):
3938            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
3939            fval = self.func(x0, *self.args)
3940            self.feval += 1
3941            if fval > fmax:
3942                fmax = fval
3943            if fval < fmin:
3944                fmin = fval
3945                best_state.cost = fval
3946                best_state.x = array(x0)
3947
3948        self.T0 = (fmax-fmin)*1.5
3949        return best_state.x
3950       
3951    def set_range(self,x0,frac):
3952        delrange = frac*(self.upper-self.lower)
3953        self.upper = x0+delrange
3954        self.lower = x0-delrange
3955
3956    def accept_test(self, dE):
3957        T = self.T
3958        self.tests += 1
3959        if dE < 0:
3960            self.accepted += 1
3961            return 1
3962        p = exp(-dE*1.0/self.boltzmann/T)
3963        if (p > random.uniform(0.0, 1.0)):
3964            self.accepted += 1
3965            return 1
3966        return 0
3967
3968    def update_guess(self, x0):
3969        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
3970
3971    def update_temp(self, x0):
3972        pass
3973
3974
3975#  A schedule due to Lester Ingber modified to use bounds - OK
3976class fast_sa(base_schedule):
3977    def init(self, **options):
3978        self.__dict__.update(options)
3979        if self.m is None:
3980            self.m = 1.0
3981        if self.n is None:
3982            self.n = 1.0
3983        self.c = self.m * exp(-self.n * self.quench)
3984
3985    def update_guess(self, x0):
3986        x0 = asarray(x0)
3987        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
3988        T = self.T
3989        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
3990        xnew = xc*(self.upper - self.lower)+self.lower
3991        return xnew
3992#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
3993#        xc = y*(self.upper - self.lower)
3994#        xnew = x0 + xc
3995#        return xnew
3996
3997    def update_temp(self):
3998        self.T = self.T0*exp(-self.c * self.k**(self.quench))
3999        self.k += 1
4000        return
4001
4002class cauchy_sa(base_schedule):     #modified to use bounds - not good
4003    def update_guess(self, x0):
4004        x0 = asarray(x0)
4005        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
4006        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
4007        xnew = xc*(self.upper - self.lower)+self.lower
4008        return xnew
4009#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
4010#        xc = self.learn_rate * self.T * tan(numbers)
4011#        xnew = x0 + xc
4012#        return xnew
4013
4014    def update_temp(self):
4015        self.T = self.T0/(1+self.k)
4016        self.k += 1
4017        return
4018
4019class boltzmann_sa(base_schedule):
4020#    def update_guess(self, x0):
4021#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
4022#        x0 = asarray(x0)
4023#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
4024#
4025#        xnew = x0 + xc*std*self.learn_rate
4026#        return xnew
4027
4028    def update_temp(self):
4029        self.k += 1
4030        self.T = self.T0 / log(self.k+1.0)
4031        return
4032
4033class log_sa(base_schedule):        #OK
4034
4035    def init(self,**options):
4036        self.__dict__.update(options)
4037       
4038    def update_guess(self,x0):     #same as default #TODO - is this a reasonable update procedure?
4039        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
4040       
4041    def update_temp(self):
4042        self.k += 1
4043        self.T = self.T0*self.slope**self.k
4044       
4045class _state(object):
4046    def __init__(self):
4047        self.x = None
4048        self.cost = None
4049
4050# TODO:
4051#     allow for general annealing temperature profile
4052#     in that case use update given by alpha and omega and
4053#     variation of all previous updates and temperature?
4054
4055# Simulated annealing   #TODO - should we switch to scipy basinhopping?
4056
4057def anneal(func, x0, args=(), schedule='fast', full_output=0,
4058           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
4059           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
4060           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
4061           ranRange=0.10,autoRan=False,dlg=None):
4062    """Minimize a function using simulated annealing.
4063
4064    Schedule is a schedule class implementing the annealing schedule.
4065    Available ones are 'fast', 'cauchy', 'boltzmann'
4066
4067    :param callable func: f(x, \*args)
4068        Function to be optimized.
4069    :param ndarray x0:
4070        Initial guess.
4071    :param tuple args:
4072        Extra parameters to `func`.
4073    :param base_schedule schedule:
4074        Annealing schedule to use (a class).
4075    :param bool full_output:
4076        Whether to return optional outputs.
4077    :param float T0:
4078        Initial Temperature (estimated as 1.2 times the largest
4079        cost-function deviation over random points in the range).
4080    :param float Tf:
4081        Final goal temperature.
4082    :param int maxeval:
4083        Maximum function evaluations.
4084    :param int maxaccept:
4085        Maximum changes to accept.
4086    :param int maxiter:
4087        Maximum cooling iterations.
4088    :param float learn_rate:
4089        Scale constant for adjusting guesses.
4090    :param float boltzmann:
4091        Boltzmann constant in acceptance test
4092        (increase for less stringent test at each temperature).
4093    :param float feps:
4094        Stopping relative error tolerance for the function value in
4095        last four coolings.
4096    :param float quench,m,n:
4097        Parameters to alter fast_sa schedule.
4098    :param float/ndarray lower,upper:
4099        Lower and upper bounds on `x`.
4100    :param int dwell:
4101        The number of times to search the space at each temperature.
4102    :param float slope:
4103        Parameter for log schedule
4104    :param bool ranStart=False:
4105        True for set 10% of ranges about x
4106
4107    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
4108
4109     * xmin (ndarray): Point giving smallest value found.
4110     * Jmin (float): Minimum value of function found.
4111     * T (float): Final temperature.
4112     * feval (int): Number of function evaluations.
4113     * iters (int): Number of cooling iterations.
4114     * accept (int): Number of tests accepted.
4115     * retval (int): Flag indicating stopping condition:
4116
4117              *  0: Points no longer changing
4118              *  1: Cooled to final temperature
4119              *  2: Maximum function evaluations
4120              *  3: Maximum cooling iterations reached
4121              *  4: Maximum accepted query locations reached
4122              *  5: Final point not the minimum amongst encountered points
4123
4124    *Notes*:
4125    Simulated annealing is a random algorithm which uses no derivative
4126    information from the function being optimized. In practice it has
4127    been more useful in discrete optimization than continuous
4128    optimization, as there are usually better algorithms for continuous
4129    optimization problems.
4130
4131    Some experimentation by trying the difference temperature
4132    schedules and altering their parameters is likely required to
4133    obtain good performance.
4134
4135    The randomness in the algorithm comes from random sampling in numpy.
4136    To obtain the same results you can call numpy.random.seed with the
4137    same seed immediately before calling scipy.optimize.anneal.
4138
4139    We give a brief description of how the three temperature schedules
4140    generate new points and vary their temperature. Temperatures are
4141    only updated with iterations in the outer loop. The inner loop is
4142    over xrange(dwell), and new points are generated for
4143    every iteration in the inner loop. (Though whether the proposed
4144    new points are accepted is probabilistic.)
4145
4146    For readability, let d denote the dimension of the inputs to func.
4147    Also, let x_old denote the previous state, and k denote the
4148    iteration number of the outer loop. All other variables not
4149    defined below are input variables to scipy.optimize.anneal itself.
4150
4151    In the 'fast' schedule the updates are ::
4152
4153        u ~ Uniform(0, 1, size=d)
4154        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
4155        xc = y * (upper - lower)
4156        x_new = x_old + xc
4157
4158        c = n * exp(-n * quench)
4159        T_new = T0 * exp(-c * k**quench)
4160
4161
4162    In the 'cauchy' schedule the updates are ::
4163
4164        u ~ Uniform(-pi/2, pi/2, size=d)
4165        xc = learn_rate * T * tan(u)
4166        x_new = x_old + xc
4167
4168        T_new = T0 / (1+k)
4169
4170    In the 'boltzmann' schedule the updates are ::
4171
4172        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
4173        y ~ Normal(0, std, size=d)
4174        x_new = x_old + learn_rate * y
4175
4176        T_new = T0 / log(1+k)
4177
4178    """
4179    x0 = asarray(x0)
4180    lower = asarray(lower)
4181    upper = asarray(upper)
4182
4183    schedule = eval(schedule+'_sa()')
4184    #   initialize the schedule
4185    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
4186                  learn_rate=learn_rate, lower=lower, upper=upper,
4187                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
4188
4189    current_state, last_state, best_state = _state(), _state(), _state()
4190    if ranStart:
4191        schedule.set_range(x0,ranRange)
4192    if T0 is None:
4193        x0 = schedule.getstart_temp(best_state)
4194    else:
4195        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
4196        best_state.x = None
4197        best_state.cost = numpy.Inf
4198
4199    last_state.x = asarray(x0).copy()
4200    fval = func(x0,*args)
4201    schedule.feval += 1
4202    last_state.cost = fval
4203    if last_state.cost < best_state.cost:
4204        best_state.cost = fval
4205        best_state.x = asarray(x0).copy()
4206    schedule.T = schedule.T0
4207    fqueue = [100, 300, 500, 700]
4208    iters = 1
4209    keepGoing = True
4210    bestn = 0
4211    while keepGoing:
4212        retval = 0
4213        for n in xrange(dwell):
4214            current_state.x = schedule.update_guess(last_state.x)
4215            current_state.cost = func(current_state.x,*args)
4216            schedule.feval += 1
4217
4218            dE = current_state.cost - last_state.cost
4219            if schedule.accept_test(dE):
4220                last_state.x = current_state.x.copy()
4221                last_state.cost = current_state.cost
4222                if last_state.cost < best_state.cost:
4223                    best_state.x = last_state.x.copy()
4224                    best_state.cost = last_state.cost
4225                    bestn = n
4226                    if best_state.cost < 1.0 and autoRan:
4227                        schedule.set_range(x0,best_state.cost/2.)                       
4228        if dlg:
4229            GoOn = dlg.Update(min(100.,best_state.cost*100),
4230                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
4231                    'Best trial:',bestn,  \
4232                    'MC/SA Residual:',best_state.cost*100,'%', \
4233                    ))[0]
4234            if not GoOn:
4235                best_state.x = last_state.x.copy()
4236                best_state.cost = last_state.cost
4237                retval = 5
4238        schedule.update_temp()
4239        iters += 1
4240        # Stopping conditions
4241        # 0) last saved values of f from each cooling step
4242        #     are all very similar (effectively cooled)
4243        # 1) Tf is set and we are below it
4244        # 2) maxeval is set and we are past it
4245        # 3) maxiter is set and we are past it
4246        # 4) maxaccept is set and we are past it
4247        # 5) user canceled run via progress bar
4248
4249        fqueue.append(squeeze(last_state.cost))
4250        fqueue.pop(0)
4251        af = asarray(fqueue)*1.0
4252        if retval == 5:
4253            print ' User terminated run; incomplete MC/SA'
4254            keepGoing = False
4255            break
4256        if all(abs((af-af[0])/af[0]) < feps):
4257            retval = 0
4258            if abs(af[-1]-best_state.cost) > feps*10:
4259                retval = 5
4260#                print "Warning: Cooled to %f at %s but this is not" \
4261#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
4262#                      + " the smallest point found."
4263            break
4264        if (Tf is not None) and (schedule.T < Tf):
4265            retval = 1
4266            break
4267        if (maxeval is not None) and (schedule.feval > maxeval):
4268            retval = 2
4269            break
4270        if (iters > maxiter):
4271            print "Warning: Maximum number of iterations exceeded."
4272            retval = 3
4273            break
4274        if (maxaccept is not None) and (schedule.accepted > maxaccept):
4275            retval = 4
4276            break
4277
4278    if full_output:
4279        return best_state.x, best_state.cost, schedule.T, \
4280               schedule.feval, iters, schedule.accepted, retval
4281    else:
4282        return best_state.x, retval
4283
4284def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
4285    outlist = []
4286    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
4287    for n in range(iCyc):
4288        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
4289        outlist.append(result[0])
4290        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
4291    out_q.put(outlist)
4292
4293def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
4294    import multiprocessing as mp
4295   
4296    nprocs = mp.cpu_count()
4297    out_q = mp.Queue()
4298    procs = []
4299    iCyc = np.zeros(nprocs)
4300    for i in range(nCyc):
4301        iCyc[i%nprocs] += 1
4302    for i in range(nprocs):
4303        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
4304        procs.append(p)
4305        p.start()
4306    resultlist = []
4307    for i in range(nprocs):
4308        resultlist += out_q.get()
4309    for p in procs:
4310        p.join()
4311    return resultlist
4312
4313def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
4314    '''default doc string
4315   
4316    :param type name: description
4317   
4318    :returns: type name: description
4319    '''
4320   
4321    class RandomDisplacementBounds(object):
4322        """random displacement with bounds"""
4323        def __init__(self, xmin, xmax, stepsize=0.5):
4324            self.xmin = xmin
4325            self.xmax = xmax
4326            self.stepsize = stepsize
4327   
4328        def __call__(self, x):
4329            """take a random step but ensure the new position is within the bounds"""
4330            while True:
4331                # this could be done in a much more clever way, but it will work for example purposes
4332                steps = self.xmax-self.xmin
4333                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
4334                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
4335                    break
4336            return xnew
4337   
4338    global tsum
4339    tsum = 0.
4340   
4341    def getMDparms(item,pfx,parmDict,varyList):
4342        parmDict[pfx+'MDaxis'] = item['axis']
4343        parmDict[pfx+'MDval'] = item['Coef'][0]
4344        if item['Coef'][1]:
4345            varyList += [pfx+'MDval',]
4346            limits = item['Coef'][2]
4347            lower.append(limits[0])
4348            upper.append(limits[1])
4349                       
4350    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
4351        parmDict[pfx+'Atype'] = item['atType']
4352        aTypes |= set([item['atType'],]) 
4353        pstr = ['Ax','Ay','Az']
4354        XYZ = [0,0,0]
4355        for i in range(3):
4356            name = pfx+pstr[i]
4357            parmDict[name] = item['Pos'][0][i]
4358            XYZ[i] = parmDict[name]
4359            if item['Pos'][1][i]:
4360                varyList += [name,]
4361                limits = item['Pos'][2][i]
4362                lower.append(limits[0])
4363                upper.append(limits[1])
4364        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
4365           
4366    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
4367        parmDict[mfx+'MolCent'] = item['MolCent']
4368        parmDict[mfx+'RBId'] = item['RBId']
4369        pstr = ['Px','Py','Pz']
4370        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
4371        for i in range(3):
4372            name = mfx+pstr[i]
4373            parmDict[name] = item['Pos'][0][i]
4374            if item['Pos'][1][i]:
4375                varyList += [name,]
4376                limits = item['Pos'][2][i]
4377                lower.append(limits[0])
4378                upper.append(limits[1])
4379        AV = item['Ori'][0]
4380        A = AV[0]
4381        V = AV[1:]
4382        for i in range(4):
4383            name = mfx+ostr[i]
4384            if i:
4385                parmDict[name] = V[i-1]
4386            else:
4387                parmDict[name] = A
4388            if item['Ovar'] == 'AV':
4389                varyList += [name,]
4390                limits = item['Ori'][2][i]
4391                lower.append(limits[0])
4392                upper.append(limits[1])
4393            elif item['Ovar'] == 'A' and not i:
4394                varyList += [name,]
4395                limits = item['Ori'][2][i]
4396                lower.append(limits[0])
4397                upper.append(limits[1])
4398        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
4399            for i in range(len(item['Tor'][0])):
4400                name = mfx+'Tor'+str(i)
4401                parmDict[name] = item['Tor'][0][i]
4402                if item['Tor'][1][i]:
4403                    varyList += [name,]
4404                    limits = item['Tor'][2][i]
4405                    lower.append(limits[0])
4406                    upper.append(limits[1])
4407        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
4408        aTypes |= set(atypes)
4409        atNo += len(atypes)
4410        return atNo
4411               
4412    def GetAtomM(Xdata,SGData):
4413        Mdata = []
4414        for xyz in Xdata:
4415            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
4416        return np.array(Mdata)
4417       
4418    def GetAtomT(RBdata,parmDict):
4419        'Needs a doc string'
4420        atNo = parmDict['atNo']
4421        nfixAt = parmDict['nfixAt']
4422        Tdata = atNo*[' ',]
4423        for iatm in range(nfixAt):
4424            parm = ':'+str(iatm)+':Atype'
4425            if parm in parmDict:
4426                Tdata[iatm] = aTypes.index(parmDict[parm])
4427        iatm = nfixAt
4428        for iObj in range(parmDict['nObj']):
4429            pfx = str(iObj)+':'
4430            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4431                if parmDict[pfx+'Type'] == 'Vector':
4432                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
4433                    nAtm = len(RBRes['rbVect'][0])
4434                else:       #Residue
4435                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
4436                    nAtm = len(RBRes['rbXYZ'])
4437                for i in range(nAtm):
4438                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
4439                    iatm += 1
4440            elif parmDict[pfx+'Type'] == 'Atom':
4441                atNo = parmDict[pfx+'atNo']
4442                parm = pfx+'Atype'              #remove extra ':'
4443                if parm in parmDict:
4444                    Tdata[atNo] = aTypes.index(parmDict[parm])
4445                iatm += 1
4446            else:
4447                continue        #skips March Dollase
4448        return Tdata
4449       
4450    def GetAtomX(RBdata,parmDict):
4451        'Needs a doc string'
4452        Bmat = parmDict['Bmat']
4453        atNo = parmDict['atNo']
4454        nfixAt = parmDict['nfixAt']
4455        Xdata = np.zeros((3,atNo))
4456        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
4457        for iatm in range(nfixAt):
4458            for key in keys:
4459                parm = ':'+str(iatm)+key
4460                if parm in parmDict:
4461                    keys[key][iatm] = parmDict[parm]
4462        iatm = nfixAt
4463        for iObj in range(parmDict['nObj']):
4464            pfx = str(iObj)+':'
4465            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4466                if parmDict[pfx+'Type'] == 'Vector':
4467                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
4468                    vecs = RBRes['rbVect']
4469                    mags = RBRes['VectMag']
4470                    Cart = np.zeros_like(vecs[0])
4471                    for vec,mag in zip(vecs,mags):
4472                        Cart += vec*mag
4473                elif parmDict[pfx+'Type'] == 'Residue':
4474                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
4475                    Cart = np.array(RBRes['rbXYZ'])
4476                    for itor,seq in enumerate(RBRes['rbSeq']):
4477                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
4478                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
4479                if parmDict[pfx+'MolCent'][1]:
4480                    Cart -= parmDict[pfx+'MolCent'][0]
4481                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
4482                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
4483                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
4484                iatm += len(Cart)
4485            elif parmDict[pfx+'Type'] == 'Atom':
4486                atNo = parmDict[pfx+'atNo']
4487                for key in keys:
4488                    parm = pfx+key[1:]              #remove extra ':'
4489                    if parm in parmDict:
4490                        keys[key][atNo] = parmDict[parm]
4491                iatm += 1
4492            else:
4493                continue        #skips March Dollase
4494        return Xdata.T
4495       
4496    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
4497        allX = np.inner(Xdata,SGM)+SGT
4498        allT = np.repeat(Tdata,allX.shape[1])
4499        allM = np.repeat(Mdata,allX.shape[1])
4500        allX = np.reshape(allX,(-1,3))
4501        return allT,allM,allX
4502
4503    def getAllX(Xdata,SGM,SGT):
4504        allX = np.inner(Xdata,SGM)+SGT
4505        allX = np.reshape(allX,(-1,3))
4506        return allX
4507       
4508    def normQuaternions(RBdata,parmDict,varyList,values):
4509        for iObj in range(parmDict['nObj']):
4510            pfx = str(iObj)+':'
4511            if parmDict[pfx+'Type'] in ['Vector','Residue']:
4512                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
4513                A,V = Q2AVdeg(Qori)
4514                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
4515                    if i:
4516                        parmDict[pfx+name] = V[i-1]
4517                    else:
4518                        parmDict[pfx+name] = A
4519       
4520    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
4521        ''' Compute structure factors for all h,k,l for phase
4522        input:
4523            refList: [ref] where each ref = h,k,l,m,d,...
4524            rcov:   array[nref,nref] covariance terms between Fo^2 values
4525            ifInv:  bool True if centrosymmetric
4526            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
4527            RBdata: [dict] rigid body dictionary
4528            varyList: [list] names of varied parameters in MC/SA (not used here)           
4529            ParmDict: [dict] problem parameters
4530        puts result F^2 in each ref[5] in refList
4531        returns:
4532            delt-F*rcov*delt-F/sum(Fo^2)
4533        '''   
4534           
4535        global tsum
4536        t0 = time.time()
4537        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
4538        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
4539        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
4540        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
4541        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
4542        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
4543        refList[5] = Aterm
4544        if not ifInv:
4545            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
4546        if len(cosTable):        #apply MD correction
4547            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
4548        sumFcsq = np.sum(refList[5])
4549        scale = parmDict['sumFosq']/sumFcsq
4550        refList[5] *= scale
4551        refList[6] = refList[4]-refList[5]
4552        M = np.inner(refList[6],np.inner(rcov,refList[6]))
4553        tsum += (time.time()-t0)
4554        return M/np.sum(ref