source: trunk/GSASIImath.py @ 3243

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

G2ctrlGUI: comment about PyGridTableBase? - python 3?
move VALIDPROTEIN to General Compute menu from Atoms/Compute?
trap the too many atoms in box for proein validation
deal with gap in protein chain for protein validation
G2restGUI" replace all wx.ALIGN_CENTER_VERTICAL with WACV
revamp grid displays to put scroll bars on grids & make them fit frame better

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