source: trunk/GSASIImath.py @ 3164

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

fix neutron resonant scattering factor calcs. - neglected constant term
put tyr/except around import pypowder in G2pwd & G2math - to allow testing

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