source: trunk/GSASIImath.py @ 3964

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

fix RB bug involving zipped lists

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