source: trunk/GSASIImath.py @ 3157

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

add Lorentzian term 'Z' to CW profile instrument parameters - constant term
change sig-q function to d*sig-q from sig-q/d2 for TOF profile instrument parameters
add filter to SumDialog? for PWDR & IMG data

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