source: trunk/GSASIImath.py @ 4054

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

changes to incommensurate magnetism math

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