source: trunk/GSASIImath.py @ 952

Last change on this file since 952 was 952, checked in by vondreele, 9 years ago

fixes to MC/SA stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 92.9 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2013-06-17 15:14:00 +0000 (Mon, 17 Jun 2013) $
5# $Author: vondreele $
6# $Revision: 952 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 952 2013-06-17 15:14:00Z vondreele $
9########### SVN repository information ###################
10'''
11*GSASIImath: computation module*
12================================
13
14Routines for least-squares minimization and other stuff
15
16'''
17import sys
18import os
19import os.path as ospath
20import random as rn
21import numpy as np
22import numpy.linalg as nl
23import numpy.ma as ma
24import cPickle
25import time
26import math
27import copy
28import GSASIIpath
29GSASIIpath.SetVersionNumber("$Revision: 952 $")
30import GSASIIElem as G2el
31import GSASIIlattice as G2lat
32import GSASIIspc as G2spc
33import numpy.fft as fft
34
35sind = lambda x: np.sin(x*np.pi/180.)
36cosd = lambda x: np.cos(x*np.pi/180.)
37tand = lambda x: np.tan(x*np.pi/180.)
38asind = lambda x: 180.*np.arcsin(x)/np.pi
39acosd = lambda x: 180.*np.arccos(x)/np.pi
40atand = lambda x: 180.*np.arctan(x)/np.pi
41atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
42
43def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
44   
45    """
46    Minimize the sum of squares of a set of equations.
47
48    ::
49   
50                    Nobs
51        x = arg min(sum(func(y)**2,axis=0))
52                    y=0
53
54    :param function func: callable method or function
55        should take at least one (possibly length N vector) argument and
56        returns M floating point numbers.
57    :param np.ndarray x0: The starting estimate for the minimization of length N
58    :param function Hess: callable method or function
59        A required function or method to compute the weighted vector and Hessian for func.
60        It must be a symmetric NxN array
61    :param tuple args: Any extra arguments to func are placed in this tuple.
62    :param float ftol: Relative error desired in the sum of squares.
63    :param float xtol: Relative error desired in the approximate solution.
64    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
65        until other limits are met (ftol, xtol)
66
67    :returns: (x,cov_x,infodict) where
68
69      * x : ndarray
70        The solution (or the result of the last iteration for an unsuccessful
71        call).
72      * cov_x : ndarray
73        Uses the fjac and ipvt optional outputs to construct an
74        estimate of the jacobian around the solution.  ``None`` if a
75        singular matrix encountered (indicates very flat curvature in
76        some direction).  This matrix must be multiplied by the
77        residual standard deviation to get the covariance of the
78        parameter estimates -- see curve_fit.
79      * infodict : dict
80        a dictionary of optional outputs with the keys:
81
82         * 'fvec' : the function evaluated at the output
83         * 'num cyc':
84         * 'nfev':
85         * 'lamMax':
86         * 'psing':
87           
88    """
89               
90    x0 = np.array(x0, ndmin=1)      #might be redundant?
91    n = len(x0)
92    if type(args) != type(()):
93        args = (args,)
94       
95    icycle = 0
96    One = np.ones((n,n))
97    lam = 0.001
98    lamMax = lam
99    nfev = 0
100    while icycle < maxcyc:
101        lamMax = max(lamMax,lam)
102        M = func(x0,*args)
103        nfev += 1
104        chisq0 = np.sum(M**2)
105        Yvec,Amat = Hess(x0,*args)
106        Adiag = np.sqrt(np.diag(Amat))
107        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
108        if np.any(psing):                #hard singularity in matrix
109            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
110        Anorm = np.outer(Adiag,Adiag)
111        Yvec /= Adiag
112        Amat /= Anorm       
113        while True:
114            Lam = np.eye(Amat.shape[0])*lam
115            Amatlam = Amat*(One+Lam)
116            try:
117                Xvec = nl.solve(Amatlam,Yvec)
118            except nl.LinAlgError:
119                print 'ouch #1'
120                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
121                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
122            Xvec /= Adiag
123            M2 = func(x0+Xvec,*args)
124            nfev += 1
125            chisq1 = np.sum(M2**2)
126            if chisq1 > chisq0:
127                lam *= 10.
128            else:
129                x0 += Xvec
130                lam /= 10.
131                break
132        if (chisq0-chisq1)/chisq0 < ftol:
133            break
134        icycle += 1
135    M = func(x0,*args)
136    nfev += 1
137    Yvec,Amat = Hess(x0,*args)
138    try:
139        Bmat = nl.inv(Amat)
140        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
141    except nl.LinAlgError:
142        print 'ouch #2 linear algebra error in LS'
143        psing = []
144        if maxcyc:
145            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
146        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
147   
148def getVCov(varyNames,varyList,covMatrix):
149    '''obtain variance-covariance terms for a set of variables. NB: the varyList
150    and covMatrix were saved by the last least squares refineemnt so they must match
151   
152    :param list varyNames: variable names to find v-cov matric for
153    :param list varyList: full list of all variables in v-cov matrix
154    :param nparray covMatrix: full variance-covariance matrix from the last
155     least squares refinement
156   
157    :returns: nparray vcov: variance-covariance matrix for the variables given
158     in varyNames
159   
160    '''
161    vcov = np.zeros((len(varyNames),len(varyNames)))
162    for i1,name1 in enumerate(varyNames):
163        for i2,name2 in enumerate(varyNames):
164            try:
165                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
166            except ValueError:
167                vcov[i1][i2] = 0.0
168    return vcov
169
170def FindAtomIndexByIDs(atomData,IDs,Draw=True):
171    '''finds the set of atom array indices for a list of atom IDs. Will search
172    either the Atom table or the drawAtom table.
173   
174    :param list atomData: Atom or drawAtom table containting coordinates, etc.
175    :param list IDs: atom IDs to be found
176    :param bool Draw: True if drawAtom table to be searched; False if Atom table
177    is searched
178   
179    :returns: list indx: atom (or drawAtom) indices
180   
181    '''
182    indx = []
183    for i,atom in enumerate(atomData):
184        if Draw and atom[-3] in IDs:
185            indx.append(i)
186        elif atom[-1] in IDs:
187            indx.append(i)
188    return indx
189
190def FillAtomLookUp(atomData):
191    '''create a dictionary of atom indexes with atom IDs as keys
192   
193    :param list atomData: Atom table to be used
194   
195    :returns: dict atomLookUp: dictionary of atom indexes with atom IDs as keys
196   
197    '''
198    atomLookUp = {}
199    for iatm,atom in enumerate(atomData):
200        atomLookUp[atom[-1]] = iatm
201    return atomLookUp
202
203def GetAtomsById(atomData,atomLookUp,IdList):
204    '''gets a list of atoms from Atom table that match a set of atom IDs
205   
206    :param list atomData: Atom table to be used
207    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
208    :param list IdList: atom IDs to be found
209   
210    :returns: list atoms: list of atoms found
211   
212    '''
213    atoms = []
214    for id in IdList:
215        atoms.append(atomData[atomLookUp[id]])
216    return atoms
217   
218def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
219    '''gets atom parameters for atoms using atom IDs
220   
221    :param list atomData: Atom table to be used
222    :param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
223    :param list IdList: atom IDs to be found
224    :param int itemLoc: pointer to desired 1st item in an atom table entry
225    :param int numItems: number of items to be retrieved
226   
227    :returns: type name: description
228   
229    '''
230    Items = []
231    if not isinstance(IdList,list):
232        IdList = [IdList,]
233    for id in IdList:
234        if numItems == 1:
235            Items.append(atomData[atomLookUp[id]][itemLoc])
236        else:
237            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
238    return Items
239   
240def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
241    '''default doc string
242   
243    :param type name: description
244   
245    :returns: type name: description
246   
247    '''
248    pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
249    dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
250    XYZ = []
251    for ind in indx:
252        names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
253        dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
254        XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
255    return XYZ
256
257def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
258    '''default doc string
259   
260    :param type name: description
261   
262    :returns: type name: description
263   
264    '''
265    for atom in atomData:
266        XYZ = np.inner(Amat,atom[cx:cx+3])
267        if atom[cia] == 'A':
268            UIJ = atom[cia+2:cia+8]
269               
270def TLS2Uij(xyz,g,Amat,rbObj):    #not used anywhere, but could be?
271    '''default doc string
272   
273    :param type name: description
274   
275    :returns: type name: description
276   
277    '''
278    TLStype,TLS = rbObj['ThermalMotion'][:2]
279    Tmat = np.zeros((3,3))
280    Lmat = np.zeros((3,3))
281    Smat = np.zeros((3,3))
282    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
283        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
284    if 'T' in TLStype:
285        Tmat = G2lat.U6toUij(TLS[:6])
286    if 'L' in TLStype:
287        Lmat = G2lat.U6toUij(TLS[6:12])
288    if 'S' in TLStype:
289        Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
290    XYZ = np.inner(Amat,xyz)
291    Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
292    Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
293    beta = np.inner(np.inner(g,Umat),g)
294    return G2lat.UijtoU6(beta)*gvec   
295       
296def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj):    #not used anywhere, but could be?
297    '''default doc string
298   
299    :param type name: description
300   
301    :returns: type name: description
302   
303    '''
304    cx,ct,cs,cia = atPtrs
305    TLStype,TLS = rbObj['ThermalMotion'][:2]
306    Tmat = np.zeros((3,3))
307    Lmat = np.zeros((3,3))
308    Smat = np.zeros((3,3))
309    G,g = G2lat.A2Gmat(Amat)
310    gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
311    if 'T' in TLStype:
312        Tmat = G2lat.U6toUij(TLS[:6])
313    if 'L' in TLStype:
314        Lmat = G2lat.U6toUij(TLS[6:12])
315    if 'S' in TLStype:
316        Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
317    for atom in atomData:
318        XYZ = np.inner(Amat,atom[cx:cx+3])
319        Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
320        if 'U' in TSLtype:
321            atom[cia+1] = TLS[0]
322            atom[cia] = 'I'
323        else:
324            atom[cia] = 'A'
325            Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
326            beta = np.inner(np.inner(g,Umat),g)
327            atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
328
329def GetXYZDist(xyz,XYZ,Amat):
330    '''gets distance from position xyz to all XYZ, xyz & XYZ are np.array
331        and are in crystal coordinates; Amat is crystal to Cart matrix
332   
333    :param type name: description
334   
335    :returns: type name: description
336   
337    '''
338    return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
339
340def getAtomXYZ(atoms,cx):
341    '''default doc string
342   
343    :param type name: description
344   
345    :returns: type name: description
346   
347    '''
348    XYZ = []
349    for atom in atoms:
350        XYZ.append(atom[cx:cx+3])
351    return np.array(XYZ)
352
353def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
354    '''default doc string
355   
356    :param type name: description
357   
358    :returns: type name: description
359   
360    '''
361    ''' returns crystal coordinates for atoms described by RBObj
362    '''
363    RBRes = RBData[RBType][RBObj['RBId']]
364    if RBType == 'Vector':
365        vecs = RBRes['rbVect']
366        mags = RBRes['VectMag']
367        Cart = np.zeros_like(vecs[0])
368        for vec,mag in zip(vecs,mags):
369            Cart += vec*mag
370    elif RBType == 'Residue':
371        Cart = np.array(RBRes['rbXYZ'])
372        for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
373            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
374            for ride in seq[3]:
375                Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
376    XYZ = np.zeros_like(Cart)
377    for i,xyz in enumerate(Cart):
378        X = prodQVQ(RBObj['Orient'][0],xyz)
379        XYZ[i] = np.inner(Bmat,X)+RBObj['Orig'][0]
380    return XYZ,Cart
381
382def UpdateMCSAxyz(Bmat,MCSA):
383    '''default doc string
384   
385    :param type name: description
386   
387    :returns: type name: description
388   
389    '''
390    xyz = []
391    atTypes = []
392    iatm = 0
393    for model in MCSA['Models'][1:]:        #skip the MD model
394        if model['Type'] == 'Atom':
395            xyz.append(model['Pos'][0])
396            atTypes.append(model['atType'])
397            iatm += 1
398        else:
399            RBRes = MCSA['rbData'][model['Type']][model['RBId']]
400            Pos = np.array(model['Pos'][0])
401            Qori = np.array(model['Ori'][0])
402            if model['Type'] == 'Vector':
403                vecs = RBRes['rbVect']
404                mags = RBRes['VectMag']
405                Cart = np.zeros_like(vecs[0])
406                for vec,mag in zip(vecs,mags):
407                    Cart += vec*mag
408            elif model['Type'] == 'Residue':
409                Cart = np.array(RBRes['rbXYZ'])
410                for itor,seq in enumerate(RBRes['rbSeq']):
411                    QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
412                    for ride in seq[3]:
413                        Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
414            if model['MolCent'][1]:
415                Cart -= model['MolCent'][0]
416            for i,x in enumerate(Cart):
417                xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
418                atType = RBRes['rbTypes'][i]
419                atTypes.append(atType)
420                iatm += 1
421    return np.array(xyz),atTypes
422   
423def SetMolCent(model,RBData):
424    '''default doc string
425   
426    :param type name: description
427   
428    :returns: type name: description
429   
430    '''
431    rideList = []
432    RBRes = RBData[model['Type']][model['RBId']]
433    if model['Type'] == 'Vector':
434        vecs = RBRes['rbVect']
435        mags = RBRes['VectMag']
436        Cart = np.zeros_like(vecs[0])
437        for vec,mag in zip(vecs,mags):
438            Cart += vec*mag
439    elif model['Type'] == 'Residue':
440        Cart = np.array(RBRes['rbXYZ'])
441        for seq in RBRes['rbSeq']:
442            rideList += seq[3]
443    centList = set(range(len(Cart)))-set(rideList)
444    cent = np.zeros(3)
445    for i in centList:
446        cent += Cart[i]
447    model['MolCent'][0] = cent/len(centList)   
448   
449def UpdateRBUIJ(Bmat,Cart,RBObj):
450    '''default doc string
451   
452    :param type name: description
453   
454    :returns: type name: description
455   
456    '''
457    ''' returns atom I/A, Uiso or UIJ for atoms at XYZ as described by RBObj
458    '''
459    TLStype,TLS = RBObj['ThermalMotion'][:2]
460    T = np.zeros(6)
461    L = np.zeros(6)
462    S = np.zeros(8)
463    if 'T' in TLStype:
464        T = TLS[:6]
465    if 'L' in TLStype:
466        L = np.array(TLS[6:12])*(np.pi/180.)**2
467    if 'S' in TLStype:
468        S = np.array(TLS[12:])*(np.pi/180.)
469    g = nl.inv(np.inner(Bmat,Bmat))
470    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
471        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
472    Uout = []
473    Q = RBObj['Orient'][0]
474    for X in Cart:
475        X = prodQVQ(Q,X)
476        if 'U' in TLStype:
477            Uout.append(['I',TLS[0],0,0,0,0,0,0])
478        elif not 'N' in TLStype:
479            U = [0,0,0,0,0,0]
480            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])
481            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])
482            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])
483            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]+  \
484                S[4]*X[0]-S[5]*X[1]-(S[6]+S[7])*X[2]
485            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]+  \
486                S[3]*X[2]-S[2]*X[0]+S[6]*X[1]
487            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]+  \
488                S[0]*X[1]-S[1]*X[2]+S[7]*X[0]
489            Umat = G2lat.U6toUij(U)
490            beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
491            Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
492        else:
493            Uout.append(['N',])
494    return Uout
495   
496def GetSHCoeff(pId,parmDict,SHkeys):
497    '''default doc string
498   
499    :param type name: description
500   
501    :returns: type name: description
502   
503    '''
504    SHCoeff = {}
505    for shkey in SHkeys:
506        shname = str(pId)+'::'+shkey
507        SHCoeff[shkey] = parmDict[shname]
508    return SHCoeff
509       
510def getMass(generalData):
511    '''default doc string
512   
513    :param type name: description
514   
515    :returns: type name: description
516   
517    '''
518    'Computes mass of unit cell contents'
519    mass = 0.
520    for i,elem in enumerate(generalData['AtomTypes']):
521        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
522    return mass   
523
524def getDensity(generalData):
525    '''default doc string
526   
527    :param type name: description
528   
529    :returns: type name: description
530   
531    '''
532    mass = getMass(generalData)
533    Volume = generalData['Cell'][7]
534    density = mass/(0.6022137*Volume)
535    return density,Volume/mass
536   
537def getWave(Parms):
538    '''default doc string
539   
540    :param type name: description
541   
542    :returns: type name: description
543   
544    '''
545    try:
546        return Parms['Lam'][1]
547    except KeyError:
548        return Parms['Lam1'][1]
549   
550################################################################################
551##### distance, angle, planes, torsion stuff stuff
552################################################################################
553
554def getSyXYZ(XYZ,ops,SGData):
555    '''default doc string
556   
557    :param type name: description
558   
559    :returns: type name: description
560   
561    '''
562    XYZout = np.zeros_like(XYZ)
563    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
564        if op == '1':
565            XYZout[i] = xyz
566        else:
567            oprs = op.split('+')
568            unit = [0,0,0]
569            if oprs[1]:
570                unit = np.array(list(eval(oprs[1])))
571            syop =int(oprs[0])
572            inv = syop/abs(syop)
573            syop *= inv
574            cent = syop/100
575            syop %= 100
576            syop -= 1
577            M,T = SGData['SGOps'][syop]
578            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
579    return XYZout
580   
581def getRestDist(XYZ,Amat):
582    '''default doc string
583   
584    :param type name: description
585   
586    :returns: type name: description
587   
588    '''
589    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
590   
591def getRestDeriv(Func,XYZ,Amat,ops,SGData):
592    '''default doc string
593   
594    :param type name: description
595   
596    :returns: type name: description
597   
598    '''
599    deriv = np.zeros((len(XYZ),3))
600    dx = 0.00001
601    for j,xyz in enumerate(XYZ):
602        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
603            XYZ[j] -= x
604            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
605            XYZ[j] += 2*x
606            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
607            XYZ[j] -= x
608            deriv[j][i] = (d1-d2)/(2*dx)
609    return deriv.flatten()
610
611def getRestAngle(XYZ,Amat):
612    '''default doc string
613   
614    :param type name: description
615   
616    :returns: type name: description
617   
618    '''
619   
620    def calcVec(Ox,Tx,Amat):
621        return np.inner(Amat,(Tx-Ox))
622
623    VecA = calcVec(XYZ[1],XYZ[0],Amat)
624    VecA /= np.sqrt(np.sum(VecA**2))
625    VecB = calcVec(XYZ[1],XYZ[2],Amat)
626    VecB /= np.sqrt(np.sum(VecB**2))
627    edge = VecB-VecA
628    edge = np.sum(edge**2)
629    angle = (2.-edge)/2.
630    angle = max(angle,-1.)
631    return acosd(angle)
632   
633def getRestPlane(XYZ,Amat):
634    '''default doc string
635   
636    :param type name: description
637   
638    :returns: type name: description
639   
640    '''
641    sumXYZ = np.zeros(3)
642    for xyz in XYZ:
643        sumXYZ += xyz
644    sumXYZ /= len(XYZ)
645    XYZ = np.array(XYZ)-sumXYZ
646    XYZ = np.inner(Amat,XYZ).T
647    Zmat = np.zeros((3,3))
648    for i,xyz in enumerate(XYZ):
649        Zmat += np.outer(xyz.T,xyz)
650    Evec,Emat = nl.eig(Zmat)
651    Evec = np.sqrt(Evec)/(len(XYZ)-3)
652    Order = np.argsort(Evec)
653    return Evec[Order[0]]
654   
655def getRestChiral(XYZ,Amat):   
656    '''default doc string
657   
658    :param type name: description
659   
660    :returns: type name: description
661   
662    '''
663    VecA = np.empty((3,3))   
664    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
665    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
666    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
667    return nl.det(VecA)
668   
669def getRestTorsion(XYZ,Amat):
670    '''default doc string
671   
672    :param type name: description
673   
674    :returns: type name: description
675   
676    '''
677    VecA = np.empty((3,3))
678    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
679    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
680    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
681    D = nl.det(VecA)
682    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
683    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
684    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
685    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
686    Ang = 1.0
687    if abs(P12) < 1.0 and abs(P23) < 1.0:
688        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
689    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
690    return TOR
691   
692def calcTorsionEnergy(TOR,Coeff=[]):
693    '''default doc string
694   
695    :param type name: description
696   
697    :returns: type name: description
698   
699    '''
700    sum = 0.
701    if len(Coeff):
702        cof = np.reshape(Coeff,(3,3)).T
703        delt = TOR-cof[1]
704        delt = np.where(delt<-180.,delt+360.,delt)
705        delt = np.where(delt>180.,delt-360.,delt)
706        term = -cof[2]*delt**2
707        val = cof[0]*np.exp(term/1000.0)
708        pMax = cof[0][np.argmin(val)]
709        Eval = np.sum(val)
710        sum = Eval-pMax
711    return sum,Eval
712
713def getTorsionDeriv(XYZ,Amat,Coeff):
714    '''default doc string
715   
716    :param type name: description
717   
718    :returns: type name: description
719   
720    '''
721    deriv = np.zeros((len(XYZ),3))
722    dx = 0.00001
723    for j,xyz in enumerate(XYZ):
724        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
725            XYZ[j] -= x
726            tor = getRestTorsion(XYZ,Amat)
727            p,d1 = calcTorsionEnergy(tor,Coeff)
728            XYZ[j] += 2*x
729            tor = getRestTorsion(XYZ,Amat)
730            p,d2 = calcTorsionEnergy(tor,Coeff)           
731            XYZ[j] -= x
732            deriv[j][i] = (d2-d1)/(2*dx)
733    return deriv.flatten()
734
735def getRestRama(XYZ,Amat):
736    '''Computes a pair of torsion angles in a 5 atom string
737   
738    :param nparray XYZ: crystallographic coordinates of 5 atoms
739    :param nparray Amat: crystal to cartesian transformation matrix
740   
741    :returns: list (phi,psi) two torsion angles in degrees
742   
743    '''
744    phi = getRestTorsion(XYZ[:5],Amat)
745    psi = getRestTorsion(XYZ[1:],Amat)
746    return phi,psi
747   
748def calcRamaEnergy(phi,psi,Coeff=[]):
749    '''Computes pseudo potential energy from a pair of torsion angles and a
750     numerical description of the potential energy surface. Used to create
751     penalty function in LS refinement.
752     
753     Eval(phi,psi) = C[0]*exp(-V/1000) where
754     
755     V = -C[3]*(phi-C[1])**2-C[4]*(psi-C[2])**2-2*(phi-C[1])*(psi-C[2])
756   
757    :param float phi: first torsion angle
758    :param float psi: second torsion angle
759    :param list Coeff: pseudo potential coefficients
760   
761    :returns: list (sum,Eval): pseudo-potential difference from minimum
762     & value; sum is used for penalty function.
763   
764    '''
765    sum = 0.
766    Eval = 0.
767    if len(Coeff):
768        cof = Coeff.T
769        dPhi = phi-cof[1]
770        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
771        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
772        dPsi = psi-cof[2]
773        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
774        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
775        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
776        val = cof[0]*np.exp(val/1000.)
777        pMax = cof[0][np.argmin(val)]
778        Eval = np.sum(val)
779        sum = Eval-pMax
780    return sum,Eval
781
782def getRamaDeriv(XYZ,Amat,Coeff):
783    '''Computes numerical derivatives of torsion angle pair pseudo potential
784    with respect of crystallographic atom coordinates of the 5 atom sequence
785   
786    :param nparray XYZ: crystallographic coordinates of 5 atoms
787    :param nparray Amat: crystal to cartesian transformation matrix
788    :param list Coeff: pseudo potential coefficients
789   
790    :returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
791     crystallographic xyz coordinates.
792   
793    '''
794    deriv = np.zeros((len(XYZ),3))
795    dx = 0.00001
796    for j,xyz in enumerate(XYZ):
797        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
798            XYZ[j] -= x
799            phi,psi = getRestRama(XYZ,Amat)
800            p,d1 = calcRamaEnergy(phi,psi,Coeff)
801            XYZ[j] += 2*x
802            phi,psi = getRestRama(XYZ,Amat)
803            p,d2 = calcRamaEnergy(phi,psi,Coeff)
804            XYZ[j] -= x
805            deriv[j][i] = (d2-d1)/(2*dx)
806    return deriv.flatten()
807
808def getRestPolefig(ODFln,SamSym,Grid):
809    '''default doc string
810   
811    :param type name: description
812   
813    :returns: type name: description
814   
815    '''
816    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
817    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
818    R = np.where(R <= 1.,2.*atand(R),0.0)
819    Z = np.zeros_like(R)
820    Z = G2lat.polfcal(ODFln,SamSym,R,P)
821    Z = np.reshape(Z,(Grid,Grid))
822    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
823
824def getRestPolefigDerv(HKL,Grid,SHCoeff):
825    '''default doc string
826   
827    :param type name: description
828   
829    :returns: type name: description
830   
831    '''
832    pass
833       
834def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
835    '''default doc string
836   
837    :param type name: description
838   
839    :returns: type name: description
840   
841    '''
842    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
843        TxT = inv*(np.inner(M,Tx)+T)+C+U
844        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
845       
846    inv = Top/abs(Top)
847    cent = abs(Top)/100
848    op = abs(Top)%100-1
849    M,T = SGData['SGOps'][op]
850    C = SGData['SGCen'][cent]
851    dx = .00001
852    deriv = np.zeros(6)
853    for i in [0,1,2]:
854        Oxyz[i] -= dx
855        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
856        Oxyz[i] += 2*dx
857        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
858        Oxyz[i] -= dx
859        Txyz[i] -= dx
860        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
861        Txyz[i] += 2*dx
862        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
863        Txyz[i] -= dx
864    return deriv
865   
866def getAngSig(VA,VB,Amat,SGData,covData={}):
867    '''default doc string
868   
869    :param type name: description
870   
871    :returns: type name: description
872   
873    '''
874    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
875        TxT = inv*(np.inner(M,Tx)+T)+C
876        TxT = G2spc.MoveToUnitCell(TxT)+U
877        return np.inner(Amat,(TxT-Ox))
878       
879    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
880        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
881        VecA /= np.sqrt(np.sum(VecA**2))
882        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
883        VecB /= np.sqrt(np.sum(VecB**2))
884        edge = VecB-VecA
885        edge = np.sum(edge**2)
886        angle = (2.-edge)/2.
887        angle = max(angle,-1.)
888        return acosd(angle)
889       
890    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
891    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
892    invA = invB = 1
893    invA = TopA/abs(TopA)
894    invB = TopB/abs(TopB)
895    centA = abs(TopA)/100
896    centB = abs(TopB)/100
897    opA = abs(TopA)%100-1
898    opB = abs(TopB)%100-1
899    MA,TA = SGData['SGOps'][opA]
900    MB,TB = SGData['SGOps'][opB]
901    CA = SGData['SGCen'][centA]
902    CB = SGData['SGCen'][centB]
903    if 'covMatrix' in covData:
904        covMatrix = covData['covMatrix']
905        varyList = covData['varyList']
906        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
907        dx = .00001
908        dadx = np.zeros(9)
909        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
910        for i in [0,1,2]:
911            OxA[i] -= dx
912            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
913            OxA[i] += 2*dx
914            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
915            OxA[i] -= dx
916           
917            TxA[i] -= dx
918            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
919            TxA[i] += 2*dx
920            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
921            TxA[i] -= dx
922           
923            TxB[i] -= dx
924            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
925            TxB[i] += 2*dx
926            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
927            TxB[i] -= dx
928           
929        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
930        if sigAng < 0.01:
931            sigAng = 0.0
932        return Ang,sigAng
933    else:
934        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
935
936def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
937    '''default doc string
938   
939    :param type name: description
940   
941    :returns: type name: description
942   
943    '''
944    def calcDist(Atoms,SyOps,Amat):
945        XYZ = []
946        for i,atom in enumerate(Atoms):
947            Inv,M,T,C,U = SyOps[i]
948            XYZ.append(np.array(atom[1:4]))
949            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
950            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
951        V1 = XYZ[1]-XYZ[0]
952        return np.sqrt(np.sum(V1**2))
953       
954    Inv = []
955    SyOps = []
956    names = []
957    for i,atom in enumerate(Oatoms):
958        names += atom[-1]
959        Op,unit = Atoms[i][-1]
960        inv = Op/abs(Op)
961        m,t = SGData['SGOps'][abs(Op)%100-1]
962        c = SGData['SGCen'][abs(Op)/100]
963        SyOps.append([inv,m,t,c,unit])
964    Dist = calcDist(Oatoms,SyOps,Amat)
965   
966    sig = -0.001
967    if 'covMatrix' in covData:
968        parmNames = []
969        dx = .00001
970        dadx = np.zeros(6)
971        for i in range(6):
972            ia = i/3
973            ix = i%3
974            Oatoms[ia][ix+1] += dx
975            a0 = calcDist(Oatoms,SyOps,Amat)
976            Oatoms[ia][ix+1] -= 2*dx
977            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
978        covMatrix = covData['covMatrix']
979        varyList = covData['varyList']
980        DistVcov = getVCov(names,varyList,covMatrix)
981        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
982        if sig < 0.001:
983            sig = -0.001
984   
985    return Dist,sig
986
987def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
988    '''default doc string
989   
990    :param type name: description
991   
992    :returns: type name: description
993   
994    '''
995
996    def calcAngle(Atoms,SyOps,Amat):
997        XYZ = []
998        for i,atom in enumerate(Atoms):
999            Inv,M,T,C,U = SyOps[i]
1000            XYZ.append(np.array(atom[1:4]))
1001            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1002            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1003        V1 = XYZ[1]-XYZ[0]
1004        V1 /= np.sqrt(np.sum(V1**2))
1005        V2 = XYZ[1]-XYZ[2]
1006        V2 /= np.sqrt(np.sum(V2**2))
1007        V3 = V2-V1
1008        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1009        return acosd(cang)
1010
1011    Inv = []
1012    SyOps = []
1013    names = []
1014    for i,atom in enumerate(Oatoms):
1015        names += atom[-1]
1016        Op,unit = Atoms[i][-1]
1017        inv = Op/abs(Op)
1018        m,t = SGData['SGOps'][abs(Op)%100-1]
1019        c = SGData['SGCen'][abs(Op)/100]
1020        SyOps.append([inv,m,t,c,unit])
1021    Angle = calcAngle(Oatoms,SyOps,Amat)
1022   
1023    sig = -0.01
1024    if 'covMatrix' in covData:
1025        parmNames = []
1026        dx = .00001
1027        dadx = np.zeros(9)
1028        for i in range(9):
1029            ia = i/3
1030            ix = i%3
1031            Oatoms[ia][ix+1] += dx
1032            a0 = calcAngle(Oatoms,SyOps,Amat)
1033            Oatoms[ia][ix+1] -= 2*dx
1034            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1035        covMatrix = covData['covMatrix']
1036        varyList = covData['varyList']
1037        AngVcov = getVCov(names,varyList,covMatrix)
1038        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1039        if sig < 0.01:
1040            sig = -0.01
1041   
1042    return Angle,sig
1043
1044def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
1045    '''default doc string
1046   
1047    :param type name: description
1048   
1049    :returns: type name: description
1050   
1051    '''
1052
1053    def calcTorsion(Atoms,SyOps,Amat):
1054       
1055        XYZ = []
1056        for i,atom in enumerate(Atoms):
1057            Inv,M,T,C,U = SyOps[i]
1058            XYZ.append(np.array(atom[1:4]))
1059            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1060            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1061        V1 = XYZ[1]-XYZ[0]
1062        V2 = XYZ[2]-XYZ[1]
1063        V3 = XYZ[3]-XYZ[2]
1064        V1 /= np.sqrt(np.sum(V1**2))
1065        V2 /= np.sqrt(np.sum(V2**2))
1066        V3 /= np.sqrt(np.sum(V3**2))
1067        M = np.array([V1,V2,V3])
1068        D = nl.det(M)
1069        Ang = 1.0
1070        P12 = np.dot(V1,V2)
1071        P13 = np.dot(V1,V3)
1072        P23 = np.dot(V2,V3)
1073        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1074        return Tors
1075           
1076    Inv = []
1077    SyOps = []
1078    names = []
1079    for i,atom in enumerate(Oatoms):
1080        names += atom[-1]
1081        Op,unit = Atoms[i][-1]
1082        inv = Op/abs(Op)
1083        m,t = SGData['SGOps'][abs(Op)%100-1]
1084        c = SGData['SGCen'][abs(Op)/100]
1085        SyOps.append([inv,m,t,c,unit])
1086    Tors = calcTorsion(Oatoms,SyOps,Amat)
1087   
1088    sig = -0.01
1089    if 'covMatrix' in covData:
1090        parmNames = []
1091        dx = .00001
1092        dadx = np.zeros(12)
1093        for i in range(12):
1094            ia = i/3
1095            ix = i%3
1096            Oatoms[ia][ix+1] -= dx
1097            a0 = calcTorsion(Oatoms,SyOps,Amat)
1098            Oatoms[ia][ix+1] += 2*dx
1099            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1100            Oatoms[ia][ix+1] -= dx           
1101        covMatrix = covData['covMatrix']
1102        varyList = covData['varyList']
1103        TorVcov = getVCov(names,varyList,covMatrix)
1104        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
1105        if sig < 0.01:
1106            sig = -0.01
1107   
1108    return Tors,sig
1109       
1110def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
1111    '''default doc string
1112   
1113    :param type name: description
1114   
1115    :returns: type name: description
1116   
1117    '''
1118
1119    def calcDist(Atoms,SyOps,Amat):
1120        XYZ = []
1121        for i,atom in enumerate(Atoms):
1122            Inv,M,T,C,U = SyOps[i]
1123            XYZ.append(np.array(atom[1:4]))
1124            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1125            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1126        V1 = XYZ[1]-XYZ[0]
1127        return np.sqrt(np.sum(V1**2))
1128       
1129    def calcAngle(Atoms,SyOps,Amat):
1130        XYZ = []
1131        for i,atom in enumerate(Atoms):
1132            Inv,M,T,C,U = SyOps[i]
1133            XYZ.append(np.array(atom[1:4]))
1134            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1135            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1136        V1 = XYZ[1]-XYZ[0]
1137        V1 /= np.sqrt(np.sum(V1**2))
1138        V2 = XYZ[1]-XYZ[2]
1139        V2 /= np.sqrt(np.sum(V2**2))
1140        V3 = V2-V1
1141        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1142        return acosd(cang)
1143
1144    def calcTorsion(Atoms,SyOps,Amat):
1145       
1146        XYZ = []
1147        for i,atom in enumerate(Atoms):
1148            Inv,M,T,C,U = SyOps[i]
1149            XYZ.append(np.array(atom[1:4]))
1150            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1151            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1152        V1 = XYZ[1]-XYZ[0]
1153        V2 = XYZ[2]-XYZ[1]
1154        V3 = XYZ[3]-XYZ[2]
1155        V1 /= np.sqrt(np.sum(V1**2))
1156        V2 /= np.sqrt(np.sum(V2**2))
1157        V3 /= np.sqrt(np.sum(V3**2))
1158        M = np.array([V1,V2,V3])
1159        D = nl.det(M)
1160        Ang = 1.0
1161        P12 = np.dot(V1,V2)
1162        P13 = np.dot(V1,V3)
1163        P23 = np.dot(V2,V3)
1164        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
1165        return Tors
1166           
1167    Inv = []
1168    SyOps = []
1169    names = []
1170    for i,atom in enumerate(Oatoms):
1171        names += atom[-1]
1172        Op,unit = Atoms[i][-1]
1173        inv = Op/abs(Op)
1174        m,t = SGData['SGOps'][abs(Op)%100-1]
1175        c = SGData['SGCen'][abs(Op)/100]
1176        SyOps.append([inv,m,t,c,unit])
1177    M = len(Oatoms)
1178    if M == 2:
1179        Val = calcDist(Oatoms,SyOps,Amat)
1180    elif M == 3:
1181        Val = calcAngle(Oatoms,SyOps,Amat)
1182    else:
1183        Val = calcTorsion(Oatoms,SyOps,Amat)
1184   
1185    sigVals = [-0.001,-0.01,-0.01]
1186    sig = sigVals[M-3]
1187    if 'covMatrix' in covData:
1188        parmNames = []
1189        dx = .00001
1190        N = M*3
1191        dadx = np.zeros(N)
1192        for i in range(N):
1193            ia = i/3
1194            ix = i%3
1195            Oatoms[ia][ix+1] += dx
1196            if M == 2:
1197                a0 = calcDist(Oatoms,SyOps,Amat)
1198            elif M == 3:
1199                a0 = calcAngle(Oatoms,SyOps,Amat)
1200            else:
1201                a0 = calcTorsion(Oatoms,SyOps,Amat)
1202            Oatoms[ia][ix+1] -= 2*dx
1203            if M == 2:
1204                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1205            elif M == 3:
1206                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
1207            else:
1208                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1209        covMatrix = covData['covMatrix']
1210        varyList = covData['varyList']
1211        Vcov = getVCov(names,varyList,covMatrix)
1212        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
1213        if sig < sigVals[M-3]:
1214            sig = sigVals[M-3]
1215   
1216    return Val,sig
1217       
1218def ValEsd(value,esd=0,nTZ=False):
1219    '''Format a floating point number with a given level of precision or
1220    with in crystallographic format with a "esd", as value(esd). If esd is
1221    negative the number is formatted with the level of significant figures
1222    appropriate if abs(esd) were the esd, but the esd is not included.
1223    if the esd is zero, approximately 6 significant figures are printed.
1224    nTZ=True causes "extra" zeros to be removed after the decimal place.
1225    for example:
1226
1227      * "1.235(3)" for value=1.2346 & esd=0.003
1228      * "1.235(3)e4" for value=12346. & esd=30
1229      * "1.235(3)e6" for value=0.12346e7 & esd=3000
1230      * "1.235" for value=1.2346 & esd=-0.003
1231      * "1.240" for value=1.2395 & esd=-0.003
1232      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
1233      * "1.23460" for value=1.2346 & esd=0.0
1234
1235    :param float value: number to be formatted
1236    :param float esd: uncertainty or if esd < 0, specifies level of
1237      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
1238    :param bool nTZ: True to remove trailing zeros (default is False)
1239    :returns: value(esd) or value as a string
1240
1241    '''
1242    # Note: this routine is Python 3 compatible -- I think
1243    if esd != 0:
1244        # transform the esd to a one or two digit integer
1245        l = math.log10(abs(esd)) % 1
1246        # cut off of 19 1.9==>(19) but 1.95==>(2) N.B. log10(1.95) = 0.2900...
1247        if l < 0.290034611362518: l+= 1.       
1248        intesd = int(round(10**l)) # esd as integer
1249        # determine the number of digits offset for the esd
1250        esdoff = int(round(math.log10(intesd*1./abs(esd))))
1251    else:
1252        esdoff = 5
1253    valoff = 0
1254    if esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
1255        # where the digit offset is to the left of the decimal place or where too many
1256        # digits are needed
1257        if abs(value) > 1:
1258            valoff = int(math.log10(abs(value)))
1259        else:
1260            valoff = int(math.log10(abs(value))-0.9999999)
1261    if esd != 0:
1262        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
1263    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
1264        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
1265    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
1266        extra = -math.log10(abs(value))
1267        if extra > 0: extra += 1
1268        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
1269    if esd > 0:
1270        out += ("({:d})").format(intesd)  # add the esd
1271    elif nTZ and '.' in out:
1272        out = out.rstrip('0')  # strip digits to right of decimal
1273    if valoff != 0:
1274        out += ("e{:d}").format(valoff) # add an exponent, when needed
1275    return out
1276
1277################################################################################
1278##### Fourier & charge flip stuff
1279################################################################################
1280
1281def adjHKLmax(SGData,Hmax):
1282    '''default doc string
1283   
1284    :param type name: description
1285   
1286    :returns: type name: description
1287   
1288    '''
1289    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1290        Hmax[0] = ((Hmax[0]+3)/6)*6
1291        Hmax[1] = ((Hmax[1]+3)/6)*6
1292        Hmax[2] = ((Hmax[2]+1)/4)*4
1293    else:
1294        Hmax[0] = ((Hmax[0]+2)/4)*4
1295        Hmax[1] = ((Hmax[1]+2)/4)*4
1296        Hmax[2] = ((Hmax[2]+1)/4)*4
1297
1298def OmitMap(data,reflData):
1299    '''default doc string
1300   
1301    :param type name: description
1302   
1303    :returns: type name: description
1304   
1305    '''
1306    generalData = data['General']
1307    if not generalData['Map']['MapType']:
1308        print '**** ERROR - Fourier map not defined'
1309        return
1310    mapData = generalData['Map']
1311    dmin = mapData['Resolution']
1312    SGData = generalData['SGData']
1313    cell = generalData['Cell'][1:8]       
1314    A = G2lat.cell2A(cell[:6])
1315    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1316    adjHKLmax(SGData,Hmax)
1317    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1318    time0 = time.time()
1319    for ref in reflData:
1320        if ref[4] >= dmin:
1321            Fosq,Fcsq,ph = ref[8:11]
1322            for i,hkl in enumerate(ref[11]):
1323                hkl = np.asarray(hkl,dtype='i')
1324                dp = 360.*ref[12][i]
1325                a = cosd(ph+dp)
1326                b = sind(ph+dp)
1327                phasep = complex(a,b)
1328                phasem = complex(a,-b)
1329                F = np.sqrt(Fosq)
1330                h,k,l = hkl+Hmax
1331                Fhkl[h,k,l] = F*phasep
1332                h,k,l = -hkl+Hmax
1333                Fhkl[h,k,l] = F*phasem
1334    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1335    print 'NB: this is just an Fobs map for now - under development'
1336    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1337    print rho.shape
1338    mapData['rho'] = np.real(rho)
1339    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1340    return mapData
1341
1342def FourierMap(data,reflData):
1343    '''default doc string
1344   
1345    :param type name: description
1346   
1347    :returns: type name: description
1348   
1349    '''
1350    generalData = data['General']
1351    if not generalData['Map']['MapType']:
1352        print '**** ERROR - Fourier map not defined'
1353        return
1354    mapData = generalData['Map']
1355    dmin = mapData['Resolution']
1356    SGData = generalData['SGData']
1357    cell = generalData['Cell'][1:8]       
1358    A = G2lat.cell2A(cell[:6])
1359    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1360    adjHKLmax(SGData,Hmax)
1361    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1362#    Fhkl[0,0,0] = generalData['F000X']
1363    time0 = time.time()
1364    for ref in reflData:
1365        if ref[4] >= dmin:
1366            Fosq,Fcsq,ph = ref[8:11]
1367            for i,hkl in enumerate(ref[11]):
1368                hkl = np.asarray(hkl,dtype='i')
1369                dp = 360.*ref[12][i]
1370                a = cosd(ph+dp)
1371                b = sind(ph+dp)
1372                phasep = complex(a,b)
1373                phasem = complex(a,-b)
1374                if 'Fobs' in mapData['MapType']:
1375                    F = np.sqrt(Fosq)
1376                    h,k,l = hkl+Hmax
1377                    Fhkl[h,k,l] = F*phasep
1378                    h,k,l = -hkl+Hmax
1379                    Fhkl[h,k,l] = F*phasem
1380                elif 'Fcalc' in mapData['MapType']:
1381                    F = np.sqrt(Fcsq)
1382                    h,k,l = hkl+Hmax
1383                    Fhkl[h,k,l] = F*phasep
1384                    h,k,l = -hkl+Hmax
1385                    Fhkl[h,k,l] = F*phasem
1386                elif 'delt-F' in mapData['MapType']:
1387                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
1388                    h,k,l = hkl+Hmax
1389                    Fhkl[h,k,l] = dF*phasep
1390                    h,k,l = -hkl+Hmax
1391                    Fhkl[h,k,l] = dF*phasem
1392                elif '2*Fo-Fc' in mapData['MapType']:
1393                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
1394                    h,k,l = hkl+Hmax
1395                    Fhkl[h,k,l] = F*phasep
1396                    h,k,l = -hkl+Hmax
1397                    Fhkl[h,k,l] = F*phasem
1398                elif 'Patterson' in mapData['MapType']:
1399                    h,k,l = hkl+Hmax
1400                    Fhkl[h,k,l] = complex(Fosq,0.)
1401                    h,k,l = -hkl+Hmax
1402                    Fhkl[h,k,l] = complex(Fosq,0.)
1403    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1404    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1405    mapData['rho'] = np.real(rho)
1406    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1407    return mapData
1408   
1409# map printing for testing purposes
1410def printRho(SGLaue,rho,rhoMax):                         
1411    '''default doc string
1412   
1413    :param type name: description
1414   
1415    :returns: type name: description
1416   
1417    '''
1418    dim = len(rho.shape)
1419    if dim == 2:
1420        ix,jy = rho.shape
1421        for j in range(jy):
1422            line = ''
1423            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1424                line += (jy-j)*'  '
1425            for i in range(ix):
1426                r = int(100*rho[i,j]/rhoMax)
1427                line += '%4d'%(r)
1428            print line+'\n'
1429    else:
1430        ix,jy,kz = rho.shape
1431        for k in range(kz):
1432            print 'k = ',k
1433            for j in range(jy):
1434                line = ''
1435                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1436                    line += (jy-j)*'  '
1437                for i in range(ix):
1438                    r = int(100*rho[i,j,k]/rhoMax)
1439                    line += '%4d'%(r)
1440                print line+'\n'
1441## keep this
1442               
1443def findOffset(SGData,A,Fhkl):   
1444    '''default doc string
1445   
1446    :param type name: description
1447   
1448    :returns: type name: description
1449   
1450    '''
1451    if SGData['SpGrp'] == 'P 1':
1452        return [0,0,0]   
1453    hklShape = Fhkl.shape
1454    hklHalf = np.array(hklShape)/2
1455    sortHKL = np.argsort(Fhkl.flatten())
1456    Fdict = {}
1457    for hkl in sortHKL:
1458        HKL = np.unravel_index(hkl,hklShape)
1459        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
1460        if F == 0.:
1461            break
1462        Fdict['%.6f'%(np.absolute(F))] = hkl
1463    Flist = np.flipud(np.sort(Fdict.keys()))
1464    F = str(1.e6)
1465    i = 0
1466    DH = []
1467    Dphi = []
1468    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1469    for F in Flist:
1470        hkl = np.unravel_index(Fdict[F],hklShape)
1471        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
1472        Uniq = np.array(Uniq,dtype='i')
1473        Phi = np.array(Phi)
1474        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
1475        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
1476        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
1477        ang0 = np.angle(Fh0,deg=True)/360.
1478        for H,phi in zip(Uniq,Phi)[1:]:
1479            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
1480            dH = H-hkl
1481            dang = ang-ang0
1482            if np.any(np.abs(dH)-Hmax > 0):    #keep low order DHs
1483                continue
1484            DH.append(dH)
1485            Dphi.append((dang+.5) % 1.0)
1486        if i > 20 or len(DH) > 30:
1487            break
1488        i += 1
1489    DH = np.array(DH)
1490    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1491    Dphi = np.array(Dphi)
1492    steps = np.array(hklShape)
1493    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
1494    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
1495    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
1496    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
1497    hist,bins = np.histogram(Mmap,bins=1000)
1498#    for i,item in enumerate(hist[:10]):
1499#        print item,bins[i]
1500    chisq = np.min(Mmap)
1501    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1502    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
1503#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
1504    return DX
1505   
1506def ChargeFlip(data,reflData,pgbar):
1507    '''default doc string
1508   
1509    :param type name: description
1510   
1511    :returns: type name: description
1512   
1513    '''
1514    generalData = data['General']
1515    mapData = generalData['Map']
1516    flipData = generalData['Flip']
1517    FFtable = {}
1518    if 'None' not in flipData['Norm element']:
1519        normElem = flipData['Norm element'].upper()
1520        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1521        for ff in FFs:
1522            if ff['Symbol'] == normElem:
1523                FFtable.update(ff)
1524    dmin = flipData['Resolution']
1525    SGData = generalData['SGData']
1526    cell = generalData['Cell'][1:8]       
1527    A = G2lat.cell2A(cell[:6])
1528    Vol = cell[6]
1529    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1530    adjHKLmax(SGData,Hmax)
1531    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
1532    time0 = time.time()
1533    for ref in reflData:
1534        dsp = ref[4]
1535        if dsp >= dmin:
1536            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
1537            if FFtable:
1538                SQ = 0.25/dsp**2
1539                ff *= G2el.ScatFac(FFtable,SQ)[0]
1540            if ref[8] > 0.:         #use only +ve Fobs**2
1541                E = np.sqrt(ref[8])/ff
1542            else:
1543                E = 0.
1544            ph = ref[10]
1545            ph = rn.uniform(0.,360.)
1546            for i,hkl in enumerate(ref[11]):
1547                hkl = np.asarray(hkl,dtype='i')
1548                dp = 360.*ref[12][i]
1549                a = cosd(ph+dp)
1550                b = sind(ph+dp)
1551                phasep = complex(a,b)
1552                phasem = complex(a,-b)
1553                h,k,l = hkl+Hmax
1554                Ehkl[h,k,l] = E*phasep
1555                h,k,l = -hkl+Hmax       #Friedel pair refl.
1556                Ehkl[h,k,l] = E*phasem
1557#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
1558    CEhkl = copy.copy(Ehkl)
1559    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1560    Emask = ma.getmask(MEhkl)
1561    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1562    Ncyc = 0
1563    old = np.seterr(all='raise')
1564    while True:       
1565        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1566        CEsig = np.std(CErho)
1567        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1568        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem! make 20. adjustible
1569        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1570        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
1571        phase = CFhkl/np.absolute(CFhkl)
1572        CEhkl = np.absolute(Ehkl)*phase
1573        Ncyc += 1
1574        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1575        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1576        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1577        if Rcf < 5.:
1578            break
1579        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
1580        if not GoOn or Ncyc > 10000:
1581            break
1582    np.seterr(**old)
1583    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
1584    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
1585    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
1586    roll = findOffset(SGData,A,CEhkl)
1587       
1588    mapData['Rcf'] = Rcf
1589    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1590    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1591    mapData['rollMap'] = [0,0,0]
1592    return mapData
1593   
1594def SearchMap(data):
1595    '''Does a search of a density map for peaks meeting the criterion of peak
1596    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
1597    mapData is data['General']['mapData']; the map is also in mapData.
1598
1599    :param data: the phase data structure
1600    :returns: (peaks,mags,dzeros) where
1601        * peaks : ndarray
1602            x,y,z positions of the peaks found in the map
1603        * mags : ndarray
1604            the magnitudes of the peaks
1605        * dzeros : ndarray
1606            the distance of the peaks from  the unit cell origin
1607
1608    '''       
1609    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1610   
1611    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
1612   
1613    def noDuplicate(xyz,peaks,Amat):
1614        XYZ = np.inner(Amat,xyz)
1615        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1616            print ' Peak',xyz,' <0.5A from another peak'
1617            return False
1618        return True
1619                           
1620    def fixSpecialPos(xyz,SGData,Amat):
1621        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
1622        X = []
1623        xyzs = [equiv[0] for equiv in equivs]
1624        for x in xyzs:
1625            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
1626                X.append(x)
1627        if len(X) > 1:
1628            return np.average(X,axis=0)
1629        else:
1630            return xyz
1631       
1632    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
1633        Mag,x0,y0,z0,sig = parms
1634        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
1635#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
1636        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
1637       
1638    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
1639        Mag,x0,y0,z0,sig = parms
1640        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1641        return M
1642       
1643    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
1644        Mag,x0,y0,z0,sig = parms
1645        dMdv = np.zeros(([5,]+list(rX.shape)))
1646        delt = .01
1647        for i in range(5):
1648            parms[i] -= delt
1649            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1650            parms[i] += 2.*delt
1651            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1652            parms[i] -= delt
1653            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
1654        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1655        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
1656        dMdv = np.reshape(dMdv,(5,rX.size))
1657        Hess = np.inner(dMdv,dMdv)
1658       
1659        return Vec,Hess
1660       
1661    generalData = data['General']
1662    phaseName = generalData['Name']
1663    SGData = generalData['SGData']
1664    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1665    drawingData = data['Drawing']
1666    peaks = []
1667    mags = []
1668    dzeros = []
1669    try:
1670        mapData = generalData['Map']
1671        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
1672        rho = copy.copy(mapData['rho'])     #don't mess up original
1673        mapHalf = np.array(rho.shape)/2
1674        res = mapData['Resolution']
1675        incre = np.array(rho.shape,dtype=np.float)
1676        step = max(1.0,1./res)+1
1677        steps = np.array(3*[step,])
1678    except KeyError:
1679        print '**** ERROR - Fourier map not defined'
1680        return peaks,mags
1681    rhoMask = ma.array(rho,mask=(rho<contLevel))
1682    indices = (-1,0,1)
1683    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1684    for roll in rolls:
1685        if np.any(roll):
1686            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
1687    indx = np.transpose(rhoMask.nonzero())
1688    peaks = indx/incre
1689    mags = rhoMask[rhoMask.nonzero()]
1690    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
1691        rho = rollMap(rho,ind)
1692        rMM = mapHalf-steps
1693        rMP = mapHalf+steps+1
1694        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1695        peakInt = np.sum(rhoPeak)*res**3
1696        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1697        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
1698        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
1699            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
1700        x1 = result[0]
1701        if not np.any(x1 < 0):
1702            mag = x1[0]
1703            peak = (np.array(x1[1:4])-ind)/incre
1704        peak = fixSpecialPos(peak,SGData,Amat)
1705        rho = rollMap(rho,-ind)       
1706    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
1707    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
1708   
1709def sortArray(data,pos,reverse=False):
1710    '''data is a list of items
1711    sort by pos in list; reverse if True
1712    '''
1713    T = []
1714    for i,M in enumerate(data):
1715        T.append((M[pos],i))
1716    D = dict(zip(T,data))
1717    T.sort()
1718    if reverse:
1719        T.reverse()
1720    X = []
1721    for key in T:
1722        X.append(D[key])
1723    return X
1724
1725def PeaksEquiv(data,Ind):
1726    '''Find the equivalent map peaks for those selected. Works on the
1727    contents of data['Map Peaks'].
1728
1729    :param data: the phase data structure
1730    :param list Ind: list of selected peak indices
1731    :returns: augmented list of peaks including those related by symmetry to the
1732      ones in Ind
1733
1734    '''       
1735    def Duplicate(xyz,peaks,Amat):
1736        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1737            return True
1738        return False
1739                           
1740    generalData = data['General']
1741    cell = generalData['Cell'][1:7]
1742    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1743    A = G2lat.cell2A(cell)
1744    SGData = generalData['SGData']
1745    mapPeaks = data['Map Peaks']
1746    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
1747    Indx = {}
1748    for ind in Ind:
1749        xyz = np.array(mapPeaks[ind][1:4])
1750        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
1751#        for x in xyzs: print x
1752        for jnd,xyz in enumerate(XYZ):       
1753            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
1754    Ind = []
1755    for ind in Indx:
1756        if Indx[ind]:
1757            Ind.append(ind)
1758    return Ind
1759               
1760def PeaksUnique(data,Ind):
1761    '''Finds the symmetry unique set of peaks from those selected. Works on the
1762    contents of data['Map Peaks'].
1763
1764    :param data: the phase data structure
1765    :param list Ind: list of selected peak indices
1766    :returns: the list of symmetry unique peaks from among those given in Ind
1767
1768    '''       
1769#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
1770
1771    def noDuplicate(xyz,peaks,Amat):
1772        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1773            return False
1774        return True
1775                           
1776    generalData = data['General']
1777    cell = generalData['Cell'][1:7]
1778    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1779    A = G2lat.cell2A(cell)
1780    SGData = generalData['SGData']
1781    mapPeaks = data['Map Peaks']
1782    Indx = {}
1783    XYZ = {}
1784    for ind in Ind:
1785        XYZ[ind] = np.array(mapPeaks[ind][1:4])
1786        Indx[ind] = True
1787    for ind in Ind:
1788        if Indx[ind]:
1789            xyz = XYZ[ind]
1790            for jnd in Ind:
1791                if ind != jnd and Indx[jnd]:                       
1792                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
1793                    xyzs = np.array([equiv[0] for equiv in Equiv])
1794                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
1795    Ind = []
1796    for ind in Indx:
1797        if Indx[ind]:
1798            Ind.append(ind)
1799    return Ind
1800   
1801################################################################################
1802##### single peak fitting profile fxn stuff
1803################################################################################
1804
1805def getCWsig(ins,pos):
1806    '''default doc string
1807   
1808    :param type name: description
1809   
1810    :returns: type name: description
1811   
1812    '''
1813    tp = tand(pos/2.0)
1814    return ins['U']*tp**2+ins['V']*tp+ins['W']
1815   
1816def getCWsigDeriv(pos):
1817    '''default doc string
1818   
1819    :param type name: description
1820   
1821    :returns: type name: description
1822   
1823    '''
1824    tp = tand(pos/2.0)
1825    return tp**2,tp,1.0
1826   
1827def getCWgam(ins,pos):
1828    '''default doc string
1829   
1830    :param type name: description
1831   
1832    :returns: type name: description
1833   
1834    '''
1835    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
1836   
1837def getCWgamDeriv(pos):
1838    '''default doc string
1839   
1840    :param type name: description
1841   
1842    :returns: type name: description
1843   
1844    '''
1845    return 1./cosd(pos/2.0),tand(pos/2.0)
1846   
1847def getTOFsig(ins,dsp):
1848    '''default doc string
1849   
1850    :param type name: description
1851   
1852    :returns: type name: description
1853   
1854    '''
1855    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-q']*dsp
1856   
1857def getTOFsigDeriv(dsp):
1858    '''default doc string
1859   
1860    :param type name: description
1861   
1862    :returns: type name: description
1863   
1864    '''
1865    return 1.0,dsp**2,dsp
1866   
1867def getTOFgamma(ins,dsp):
1868    '''default doc string
1869   
1870    :param type name: description
1871   
1872    :returns: type name: description
1873   
1874    '''
1875    return ins['X']*dsp+ins['Y']*dsp**2
1876   
1877def getTOFgammaDeriv(dsp):
1878    '''default doc string
1879   
1880    :param type name: description
1881   
1882    :returns: type name: description
1883   
1884    '''
1885    return dsp,dsp**2
1886   
1887def getTOFbeta(ins,dsp):
1888    '''default doc string
1889   
1890    :param type name: description
1891   
1892    :returns: type name: description
1893   
1894    '''
1895    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp
1896   
1897def getTOFbetaDeriv(dsp):
1898    '''default doc string
1899   
1900    :param type name: description
1901   
1902    :returns: type name: description
1903   
1904    '''
1905    return 1.0,1./dsp**4,1./dsp
1906   
1907def getTOFalpha(ins,dsp):
1908    '''default doc string
1909   
1910    :param type name: description
1911   
1912    :returns: type name: description
1913   
1914    '''
1915    return ins['alpha']/dsp
1916   
1917def getTOFalphaDeriv(dsp):
1918    '''default doc string
1919   
1920    :param type name: description
1921   
1922    :returns: type name: description
1923   
1924    '''
1925    return 1./dsp
1926   
1927def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
1928    '''default doc string
1929   
1930    :param type name: description
1931   
1932    :returns: type name: description
1933   
1934    '''
1935    ind = 0
1936    if useFit:
1937        ind = 1
1938    ins = {}
1939    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
1940        for x in ['U','V','W','X','Y']:
1941            ins[x] = Parms[x][ind]
1942        if ifQ:                              #qplot - convert back to 2-theta
1943            pos = 2.0*asind(pos*wave/(4*math.pi))
1944        sig = getCWsig(ins,pos)
1945        gam = getCWgam(ins,pos)           
1946        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
1947    else:
1948        if ifQ:
1949            dsp = 2.*np.pi/pos
1950            pos = Parms['difC']*dsp
1951        else:
1952            dsp = pos/Parms['difC'][1]
1953        if 'Pdabc' in Parms2:
1954            for x in ['sig-0','sig-1','sig-q','X','Y']:
1955                ins[x] = Parms[x][ind]
1956            Pdabc = Parms2['Pdabc'].T
1957            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1958            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1959        else:
1960            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-q','X','Y']:
1961                ins[x] = Parms[x][ind]
1962            alp = getTOFalpha(ins,dsp)
1963            bet = getTOFbeta(ins,dsp)
1964        sig = getTOFsig(ins,dsp)
1965        gam = getTOFgamma(ins,dsp)
1966        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
1967    return XY
1968   
1969################################################################################
1970##### MC/SA stuff
1971################################################################################
1972
1973#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
1974# Original Author: Travis Oliphant 2002
1975# Bug-fixes in 2006 by Tim Leslie
1976
1977
1978import numpy
1979from numpy import asarray, tan, exp, ones, squeeze, sign, \
1980     all, log, sqrt, pi, shape, array, minimum, where
1981from numpy import random
1982
1983__all__ = ['anneal']
1984
1985_double_min = numpy.finfo(float).min
1986_double_max = numpy.finfo(float).max
1987class base_schedule(object):
1988    def __init__(self):
1989        self.dwell = 20
1990        self.learn_rate = 0.5
1991        self.lower = -10
1992        self.upper = 10
1993        self.Ninit = 50
1994        self.accepted = 0
1995        self.tests = 0
1996        self.feval = 0
1997        self.k = 0
1998        self.T = None
1999
2000    def init(self, **options):
2001        self.__dict__.update(options)
2002        self.lower = asarray(self.lower)
2003        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2004        self.upper = asarray(self.upper)
2005        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2006        self.k = 0
2007        self.accepted = 0
2008        self.feval = 0
2009        self.tests = 0
2010
2011    def getstart_temp(self, best_state):
2012        """ Find a matching starting temperature and starting parameters vector
2013        i.e. find x0 such that func(x0) = T0.
2014
2015        Parameters
2016        ----------
2017        best_state : _state
2018            A _state object to store the function value and x0 found.
2019
2020        returns
2021        -------
2022        x0 : array
2023            The starting parameters vector.
2024        """
2025
2026        assert(not self.dims is None)
2027        lrange = self.lower
2028        urange = self.upper
2029        fmax = _double_min
2030        fmin = _double_max
2031        for _ in range(self.Ninit):
2032            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2033            fval = self.func(x0, *self.args)
2034            self.feval += 1
2035            if fval > fmax:
2036                fmax = fval
2037            if fval < fmin:
2038                fmin = fval
2039                best_state.cost = fval
2040                best_state.x = array(x0)
2041
2042        self.T0 = (fmax-fmin)*1.5
2043        return best_state.x
2044
2045    def accept_test(self, dE):
2046        T = self.T
2047        self.tests += 1
2048        if dE < 0:
2049            self.accepted += 1
2050            return 1
2051        p = exp(-dE*1.0/self.boltzmann/T)
2052        if (p > random.uniform(0.0, 1.0)):
2053            self.accepted += 1
2054            return 1
2055        return 0
2056
2057    def update_guess(self, x0):
2058        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2059
2060    def update_temp(self, x0):
2061        pass
2062
2063
2064#  A schedule due to Lester Ingber
2065class fast_sa(base_schedule):
2066    def init(self, **options):
2067        self.__dict__.update(options)
2068        if self.m is None:
2069            self.m = 1.0
2070        if self.n is None:
2071            self.n = 1.0
2072        self.c = self.m * exp(-self.n * self.quench)
2073
2074    def update_guess(self, x0):
2075        x0 = asarray(x0)
2076        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2077        T = self.T
2078        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0
2079        xc = y*(self.upper - self.lower)/2.0+self.lower
2080        return xc
2081#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2082#        xc = y*(self.upper - self.lower)
2083#        xnew = x0 + xc
2084#        return xnew
2085
2086    def update_temp(self):
2087        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2088        self.k += 1
2089        return
2090
2091class cauchy_sa(base_schedule):
2092#    def update_guess(self, x0):
2093#        x0 = asarray(x0)
2094#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2095#        xc = self.learn_rate * self.T * tan(numbers)
2096#        xnew = x0 + xc
2097#        return xnew
2098
2099    def update_temp(self):
2100        self.T = self.T0/(1+self.k)
2101        self.k += 1
2102        return
2103
2104class boltzmann_sa(base_schedule):
2105#    def update_guess(self, x0):
2106#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2107#        x0 = asarray(x0)
2108#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2109#
2110#        xnew = x0 + xc*std*self.learn_rate
2111#        return xnew
2112
2113    def update_temp(self):
2114        self.k += 1
2115        self.T = self.T0 / log(self.k+1.0)
2116        return
2117
2118class log_sa(base_schedule):
2119
2120    def init(self,**options):
2121        self.__dict__.update(options)
2122       
2123#    def update_guess(self,x0):
2124#        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2125       
2126    def update_temp(self):
2127        self.k += 1
2128        self.T = self.T0*self.slope**self.k
2129       
2130class Tremayne_sa(base_schedule):   #needs fixing for two steps
2131
2132    def init(self,**options):
2133        self.__dict__.update(options)
2134
2135#    def update_guess(self,x0):
2136#        x0 = np.asarray(x0)
2137#        u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
2138#        xnew = x0+u
2139#        return xnew       
2140   
2141    def update_temp(self):
2142        self.k += 1
2143        self.T = self.T0*self.slope**k
2144   
2145class _state(object):
2146    def __init__(self):
2147        self.x = None
2148        self.cost = None
2149
2150# TODO:
2151#     allow for general annealing temperature profile
2152#     in that case use update given by alpha and omega and
2153#     variation of all previous updates and temperature?
2154
2155# Simulated annealing
2156
2157def anneal(func, x0, args=(), schedule='fast', full_output=0,
2158           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2159           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2160           lower=-100, upper=100, dwell=50, slope=0.9,dlg=None):
2161    """Minimize a function using simulated annealing.
2162
2163    Schedule is a schedule class implementing the annealing schedule.
2164    Available ones are 'fast', 'cauchy', 'boltzmann'
2165
2166    :param callable func: f(x, \*args)
2167        Function to be optimized.
2168    :param ndarray x0:
2169        Initial guess.
2170    :param tuple args:
2171        Extra parameters to `func`.
2172    :param base_schedule schedule:
2173        Annealing schedule to use (a class).
2174    :param bool full_output:
2175        Whether to return optional outputs.
2176    :param float T0:
2177        Initial Temperature (estimated as 1.2 times the largest
2178        cost-function deviation over random points in the range).
2179    :param float Tf:
2180        Final goal temperature.
2181    :param int maxeval:
2182        Maximum function evaluations.
2183    :param int maxaccept:
2184        Maximum changes to accept.
2185    :param int maxiter:
2186        Maximum cooling iterations.
2187    :param float learn_rate:
2188        Scale constant for adjusting guesses.
2189    :param float boltzmann:
2190        Boltzmann constant in acceptance test
2191        (increase for less stringent test at each temperature).
2192    :param float feps:
2193        Stopping relative error tolerance for the function value in
2194        last four coolings.
2195    :param float quench,m,n:
2196        Parameters to alter fast_sa schedule.
2197    :param float/ndarray lower,upper:
2198        Lower and upper bounds on `x`.
2199    :param int dwell:
2200        The number of times to search the space at each temperature.
2201    :param float slope:
2202        Parameter for log schedule
2203
2204    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2205
2206     * xmin (ndarray): Point giving smallest value found.
2207     * Jmin (float): Minimum value of function found.
2208     * T (float): Final temperature.
2209     * feval (int): Number of function evaluations.
2210     * iters (int): Number of cooling iterations.
2211     * accept (int): Number of tests accepted.
2212     * retval (int): Flag indicating stopping condition:
2213
2214              *  0: Points no longer changing
2215              *  1: Cooled to final temperature
2216              *  2: Maximum function evaluations
2217              *  3: Maximum cooling iterations reached
2218              *  4: Maximum accepted query locations reached
2219              *  5: Final point not the minimum amongst encountered points
2220
2221    *Notes*:
2222    Simulated annealing is a random algorithm which uses no derivative
2223    information from the function being optimized. In practice it has
2224    been more useful in discrete optimization than continuous
2225    optimization, as there are usually better algorithms for continuous
2226    optimization problems.
2227
2228    Some experimentation by trying the difference temperature
2229    schedules and altering their parameters is likely required to
2230    obtain good performance.
2231
2232    The randomness in the algorithm comes from random sampling in numpy.
2233    To obtain the same results you can call numpy.random.seed with the
2234    same seed immediately before calling scipy.optimize.anneal.
2235
2236    We give a brief description of how the three temperature schedules
2237    generate new points and vary their temperature. Temperatures are
2238    only updated with iterations in the outer loop. The inner loop is
2239    over xrange(dwell), and new points are generated for
2240    every iteration in the inner loop. (Though whether the proposed
2241    new points are accepted is probabilistic.)
2242
2243    For readability, let d denote the dimension of the inputs to func.
2244    Also, let x_old denote the previous state, and k denote the
2245    iteration number of the outer loop. All other variables not
2246    defined below are input variables to scipy.optimize.anneal itself.
2247
2248    In the 'fast' schedule the updates are ::
2249
2250        u ~ Uniform(0, 1, size=d)
2251        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2252        xc = y * (upper - lower)
2253        x_new = x_old + xc
2254
2255        c = n * exp(-n * quench)
2256        T_new = T0 * exp(-c * k**quench)
2257
2258
2259    In the 'cauchy' schedule the updates are ::
2260
2261        u ~ Uniform(-pi/2, pi/2, size=d)
2262        xc = learn_rate * T * tan(u)
2263        x_new = x_old + xc
2264
2265        T_new = T0 / (1+k)
2266
2267    In the 'boltzmann' schedule the updates are ::
2268
2269        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2270        y ~ Normal(0, std, size=d)
2271        x_new = x_old + learn_rate * y
2272
2273        T_new = T0 / log(1+k)
2274
2275    """
2276    x0 = asarray(x0)
2277    lower = asarray(lower)
2278    upper = asarray(upper)
2279
2280    schedule = eval(schedule+'_sa()')
2281    #   initialize the schedule
2282    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2283                  learn_rate=learn_rate, lower=lower, upper=upper,
2284                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2285
2286    current_state, last_state, best_state = _state(), _state(), _state()
2287    if T0 is None:
2288        x0 = schedule.getstart_temp(best_state)
2289    else:
2290        best_state.x = None
2291        best_state.cost = numpy.Inf
2292
2293    last_state.x = asarray(x0).copy()
2294    fval = func(x0,*args)
2295    schedule.feval += 1
2296    last_state.cost = fval
2297    if last_state.cost < best_state.cost:
2298        best_state.cost = fval
2299        best_state.x = asarray(x0).copy()
2300    schedule.T = schedule.T0
2301    fqueue = [100, 300, 500, 700]
2302    iters = 0
2303    while 1:
2304        for n in xrange(dwell):
2305            current_state.x = schedule.update_guess(last_state.x)
2306            current_state.cost = func(current_state.x,*args)
2307            schedule.feval += 1
2308
2309            dE = current_state.cost - last_state.cost
2310            if schedule.accept_test(dE):
2311                last_state.x = current_state.x.copy()
2312                last_state.cost = current_state.cost
2313                if last_state.cost < best_state.cost:
2314                    best_state.x = last_state.x.copy()
2315                    best_state.cost = last_state.cost
2316        if dlg:
2317            GoOn = dlg.Update(best_state.cost*100,
2318                newmsg='%s%8.3f\n%s%8.3f%s'%('Temperature =',schedule.T,'MC/SA Residual =',best_state.cost*100,'%'))[0]
2319            if not GoOn:
2320                break
2321        schedule.update_temp()
2322        iters += 1
2323        # Stopping conditions
2324        # 0) last saved values of f from each cooling step
2325        #     are all very similar (effectively cooled)
2326        # 1) Tf is set and we are below it
2327        # 2) maxeval is set and we are past it
2328        # 3) maxiter is set and we are past it
2329        # 4) maxaccept is set and we are past it
2330
2331        fqueue.append(squeeze(last_state.cost))
2332        fqueue.pop(0)
2333        af = asarray(fqueue)*1.0
2334        if all(abs((af-af[0])/af[0]) < feps):
2335            retval = 0
2336            if abs(af[-1]-best_state.cost) > feps*10:
2337                retval = 5
2338                print "Warning: Cooled to %f at %s but this is not" \
2339                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2340                      + " the smallest point found."
2341            break
2342        if (Tf is not None) and (schedule.T < Tf):
2343            retval = 1
2344            break
2345        if (maxeval is not None) and (schedule.feval > maxeval):
2346            retval = 2
2347            break
2348        if (iters > maxiter):
2349            print "Warning: Maximum number of iterations exceeded."
2350            retval = 3
2351            break
2352        if (maxaccept is not None) and (schedule.accepted > maxaccept):
2353            retval = 4
2354            break
2355
2356    if full_output:
2357        return best_state.x, best_state.cost, schedule.T, \
2358               schedule.feval, iters, schedule.accepted, retval
2359    else:
2360        return best_state.x, retval
2361
2362
2363def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
2364    '''default doc string
2365   
2366    :param type name: description
2367   
2368    :returns: type name: description
2369   
2370    '''
2371    gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
2372   
2373    def getMDparms(item,pfx,parmDict,varyList):
2374        parmDict[pfx+'MDaxis'] = item['axis']
2375        parmDict[pfx+'MDval'] = item['Coef'][0]
2376        if item['Coef'][1]:
2377            varyList += [pfx+'MDval',]
2378            limits = item['Coef'][2]
2379            lower.append(limits[0])
2380            upper.append(limits[1])
2381                       
2382    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
2383        parmDict[pfx+'Atype'] = item['atType']
2384        aTypes |= set([item['atType'],]) 
2385        pstr = ['Ax','Ay','Az']
2386        XYZ = [0,0,0]
2387        for i in range(3):
2388            name = pfx+pstr[i]
2389            parmDict[name] = item['Pos'][0][i]
2390            XYZ[i] = parmDict[name]
2391            if item['Pos'][1][i]:
2392                varyList += [name,]
2393                limits = item['Pos'][2][i]
2394                lower.append(limits[0])
2395                upper.append(limits[1])
2396        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
2397           
2398    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
2399        parmDict[mfx+'MolCent'] = item['MolCent']
2400        parmDict[mfx+'RBId'] = item['RBId']
2401        pstr = ['Px','Py','Pz']
2402        ostr = ['Qa','Qi','Qj','Qk']
2403        for i in range(3):
2404            name = mfx+pstr[i]
2405            parmDict[name] = item['Pos'][0][i]
2406            if item['Pos'][1][i]:
2407                varyList += [name,]
2408                limits = item['Pos'][2][i]
2409                lower.append(limits[0])
2410                upper.append(limits[1])
2411        for i in range(4):
2412            name = mfx+ostr[i]
2413            parmDict[name] = item['Ori'][0][i]
2414            if item['Ovar'] == 'AV' and i:
2415                varyList += [name,]
2416                limits = item['Ori'][2][i]
2417                lower.append(limits[0])
2418                upper.append(limits[1])
2419            elif item['Ovar'] == 'A' and not i:
2420                varyList += [name,]
2421                limits = item['Ori'][2][i]
2422                lower.append(limits[0])
2423                upper.append(limits[1])
2424        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
2425            for i in range(len(item['Tor'][0])):
2426                name = mfx+'Tor'+str(i)
2427                parmDict[name] = item['Tor'][0][i]
2428                if item['Tor'][1][i]:
2429                    varyList += [name,]
2430                    limits = item['Tor'][2][i]
2431                    lower.append(limits[0])
2432                    upper.append(limits[1])
2433        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2434        aTypes |= set(atypes)
2435        atNo += len(atypes)
2436        return atNo
2437       
2438    def GetAtomM(Xdata,SGData):
2439        Mdata = []
2440        for xyz in Xdata.T:
2441            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
2442        return np.array(Mdata)
2443
2444    def GetAtomTX(RBdata,parmDict):
2445        'Needs a doc string'
2446        Bmat = parmDict['Bmat']
2447        atNo = parmDict['atNo']
2448        nfixAt = parmDict['nfixAt']
2449        Tdata = atNo*[' ',]
2450        Xdata = np.zeros((3,atNo))
2451        keys = {':Atype':Tdata,':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
2452        for iatm in range(nfixAt):
2453            for key in keys:
2454                parm = ':'+str(iatm)+key
2455                if parm in parmDict:
2456                    keys[key][iatm] = parmDict[parm]
2457        iatm = nfixAt
2458        for iObj in range(parmDict['nObj']):
2459            pfx = str(iObj)+':'
2460            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2461                if parmDict[pfx+'Type'] == 'Vector':
2462                    RBId = parmDict[pfx+'RBId']
2463                    RBRes = RBdata['Vector'][RBId]
2464                    aTypes = RBRes['rbTypes']
2465                    vecs = RBRes['rbVect']
2466                    mags = RBRes['VectMag']
2467                    Cart = np.zeros_like(vecs[0])
2468                    for vec,mag in zip(vecs,mags):
2469                        Cart += vec*mag
2470                elif parmDict[pfx+'Type'] == 'Residue':
2471                    RBId = parmDict[pfx+'RBId']
2472                    RBRes = RBdata['Residue'][RBId]
2473                    aTypes = RBRes['rbTypes']
2474                    Cart = np.array(RBRes['rbXYZ'])
2475                    for itor,seq in enumerate(RBRes['rbSeq']):
2476                        tName = pfx+'Tor'+str(itor)
2477                        QuatA = AVdeg2Q(parmDict[tName],Cart[seq[0]]-Cart[seq[1]])
2478                        for ride in seq[3]:
2479                            Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
2480                if parmDict[pfx+'MolCent'][1]:
2481                    Cart -= parmDict[pfx+'MolCent'][0]
2482                Qori = np.array([parmDict[pfx+'Qa'],parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2483                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
2484                for i,x in enumerate(Cart):
2485                    X = np.inner(Bmat,prodQVQ(Qori,x))+Pos
2486                    for j in range(3):
2487                        Xdata[j][iatm] = X[j]
2488                    Tdata[iatm] = aTypes[i]
2489                    iatm += 1
2490            elif parmDict[pfx+'Type'] == 'Atom':
2491                atNo = parmDict[pfx+'atNo']
2492                for key in keys:
2493                    parm = pfx+key[1:]              #remove extra ':'
2494                    if parm in parmDict:
2495                        keys[key][atNo] = parmDict[parm]
2496            else:
2497                continue        #skips March Dollase
2498        return Tdata,Xdata
2499   
2500    def calcMDcorr(MDval,MDaxis,Uniq,G):
2501        sumMD = 0
2502        for H in Uniq:           
2503            cosP,sinP = G2lat.CosSinAngle(H,MDaxis,G)
2504            A = 1.0/np.sqrt((MDval*cosP)**2+sinP**2/MDval)
2505            sumMD += A**3
2506        return sumMD
2507       
2508    def mcsaCalc(values,refList,rcov,ifInv,RBdata,varyList,parmDict):
2509#        ''' Compute structure factors for all h,k,l for phase
2510#        input:
2511#            refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
2512#            ParmDict:
2513#        puts result F^2 in each ref[8] in refList
2514#        '''       
2515        twopi = 2.0*np.pi
2516        parmDict.update(dict(zip(varyList,values)))
2517        Tdata,Xdata = GetAtomTX(RBdata,parmDict)
2518        Mdata = parmDict['Mdata']
2519        FF = np.zeros(len(Tdata))
2520        MDval = parmDict['0:MDval']
2521        MDaxis = parmDict['0:MDaxis']
2522        Gmat = parmDict['Gmat']
2523        Srefs = np.zeros(len(refList))
2524        sumFcsq = 0.
2525        for refl in refList:
2526            fbs = 0.
2527            H = refl[:3]
2528            for i,El in enumerate(Tdata):
2529                FF[i] = refl[7][El]
2530            Uniq = refl[8]
2531            phi = refl[9]
2532            phase = twopi*(np.inner(Uniq,(Xdata.T))+phi[:,np.newaxis])
2533            sinp = np.sin(phase)
2534            cosp = np.cos(phase)
2535            occ = Mdata/len(Uniq)
2536            fa = np.array(FF*occ*cosp)
2537            fas = np.sum(fa)
2538            if not ifInv:
2539                fb = np.array(FF*occ*sinp)
2540                fbs = np.sum(fb)
2541            fcsq = (fas**2+fbs**2)*refl[3]      #*calcMDcorr(MDval,MDaxis,Uniq,Gmat)
2542            sumFcsq += fcsq
2543            refl[5] = fcsq
2544        scale = (parmDict['sumFosq']/sumFcsq)
2545        for iref,refl in enumerate(refList):
2546            refl[5] *= scale
2547            Srefs[iref] = refl[4]-refl[5]
2548        M = np.inner(Srefs,np.inner(rcov,Srefs))
2549        return M/parmDict['sumFosq']**2
2550   
2551    sq8ln2 = np.sqrt(8*np.log(2))
2552    sq2pi = np.sqrt(2*np.pi)
2553    sq4pi = np.sqrt(4*np.pi)
2554    generalData = data['General']
2555    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2556    Gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])[0]
2557    SGData = generalData['SGData']
2558    fixAtoms = data['Atoms']                       #if any
2559    cx,ct,cs = generalData['AtomPtrs'][:3]
2560    aTypes = set([])
2561    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
2562    varyList = []
2563    atNo = 0
2564    for atm in fixAtoms:
2565        pfx = ':'+str(atNo)+':'
2566        parmDict[pfx+'Atype'] = atm[ct]
2567        aTypes |= set([atm[ct],])
2568        pstr = ['Ax','Ay','Az']
2569        parmDict[pfx+'Amul'] = atm[cs+1]
2570        for i in range(3):
2571            name = pfx+pstr[i]
2572            parmDict[name] = atm[cx+i]
2573        atNo += 1
2574    parmDict['nfixAt'] = len(fixAtoms)       
2575    MCSA = generalData['MCSA controls']
2576    reflName = MCSA['Data source']
2577    phaseName = generalData['Name']
2578    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
2579    upper = []
2580    lower = []
2581    for i,item in enumerate(MCSAObjs):
2582        mfx = str(i)+':'
2583        parmDict[mfx+'Type'] = item['Type']
2584        if item['Type'] == 'MD':
2585            getMDparms(item,mfx,parmDict,varyList)
2586        elif item['Type'] == 'Atom':
2587            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
2588            parmDict[mfx+'atNo'] = atNo
2589            atNo += 1
2590        elif item['Type'] in ['Residue','Vector']:
2591            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
2592    parmDict['atNo'] = atNo                 #total no. of atoms
2593    parmDict['nObj'] = len(MCSAObjs)
2594    Tdata,Xdata = GetAtomTX(RBdata,parmDict)
2595    parmDict['Mdata'] = GetAtomM(Xdata,SGData)
2596    FFtables = G2el.GetFFtable(aTypes)
2597    refs = []
2598    sumFosq = 0
2599    if 'PWDR' in reflName:
2600        for ref in reflData:
2601            h,k,l,m,d,pos,sig,gam,f = ref[:9]
2602            if d >= MCSA['dmin']:
2603                sig = gamFW(sig,gam)/sq8ln2        #--> sig from FWHM
2604                SQ = 0.25/d**2
2605                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2606                FFs = G2el.getFFvalues(FFtables,SQ)
2607                refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi])
2608                sumFosq += f*m
2609        nRef = len(refs)
2610        rcov = np.zeros((nRef,nRef))
2611        for iref,refI in enumerate(refs):
2612            rcov[iref][iref] = 1./(sq4pi*refI[6])
2613            for jref,refJ in enumerate(refs[:iref]):
2614                t1 = refI[6]**2+refJ[6]**2
2615                t2 = (refJ[5]-refI[5])**2/(2.*t1)
2616                if t2 > 10.:
2617                    rcov[iref][jref] = 0.
2618                else:
2619                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
2620        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2621        Rdiag = np.sqrt(np.diag(rcov))
2622        Rnorm = np.outer(Rdiag,Rdiag)
2623        rcov /= Rnorm
2624    elif 'Pawley' in reflName:
2625        covMatrix = covData['covMatrix']
2626        vList = covData['varyList']
2627        for iref,refI in enumerate(reflData):
2628            h,k,l,m,d,v,f,s = refI
2629            if d >= MCSA['dmin'] and v:       #skip unrefined ones
2630                SQ = 0.25/d**2
2631                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2632                FFs = G2el.getFFvalues(FFtables,SQ)
2633                refs.append([h,k,l,m,f*m,iref,0.,FFs,Uniq,phi])
2634                sumFosq += f*m
2635        nRef = len(refs)
2636        pfx = str(data['pId'])+'::PWLref:'
2637        rcov = np.zeros((nRef,nRef))       
2638        for iref,refI in enumerate(refs):
2639            I = refI[5]
2640            nameI = pfx+str(I)
2641            if nameI in vList:
2642                Iindx = vList.index(nameI)
2643                rcov[iref][iref] = covMatrix[Iindx][Iindx]
2644                for jref,refJ in enumerate(refs[:iref]):
2645                    J = refJ[5]
2646                    nameJ = pfx+str(J)
2647                    try:
2648                        rcov[iref][jref] = covMatrix[vList.index(nameI)][vList.index(nameJ)]
2649                    except ValueError:
2650                        rcov[iref][jref] = rcov[iref][jref-1]
2651            else:
2652                rcov[iref] = rcov[iref-1]
2653                rcov[iref][iref] = rcov[iref-1][iref-1]
2654        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2655        Rdiag = np.sqrt(np.diag(rcov))
2656        Rnorm = np.outer(Rdiag,Rdiag)
2657        rcov /= Rnorm
2658    elif 'HKLF' in reflName:
2659        for ref in reflData:
2660            [h,k,l,m,d],f = ref[:5],ref[6]
2661            if d >= MCSA['dmin']:
2662                SQ = 0.25/d**2
2663                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2664                FFs = G2el.getFFvalues(FFtables,SQ)
2665                refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi])
2666                sumFosq += f*m
2667        rcov = np.identity(len(refs))
2668    parmDict['sumFosq'] = sumFosq
2669    x0 = [parmDict[val] for val in varyList]
2670    ifInv = SGData['SGInv']
2671    results = anneal(mcsaCalc,x0,args=(refs,rcov,ifInv,RBdata,varyList,parmDict), 
2672        schedule=MCSA['Algorithm'], full_output=True,
2673        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
2674        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
2675        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
2676        lower=lower, upper=upper, slope=MCSA['log slope'],dlg=pgbar)
2677    return [False,results[1],results[2],results[0],varyList]
2678
2679       
2680################################################################################
2681##### Quaternion stuff
2682################################################################################
2683
2684def prodQQ(QA,QB):
2685    ''' Grassman quaternion product
2686        QA,QB quaternions; q=r+ai+bj+ck
2687    '''
2688    D = np.zeros(4)
2689    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
2690    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
2691    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
2692    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
2693    return D
2694   
2695def normQ(QA):
2696    ''' get length of quaternion & normalize it
2697        q=r+ai+bj+ck
2698    '''
2699    n = np.sqrt(np.sum(np.array(QA)**2))
2700    return QA/n
2701   
2702def invQ(Q):
2703    '''
2704        get inverse of quaternion
2705        q=r+ai+bj+ck; q* = r-ai-bj-ck
2706    '''
2707    return Q*np.array([1,-1,-1,-1])
2708   
2709def prodQVQ(Q,V):
2710    """
2711    compute the quaternion vector rotation qvq-1 = v'
2712    q=r+ai+bj+ck
2713    """
2714    VP = np.zeros(3)
2715    T2 = Q[0]*Q[1]
2716    T3 = Q[0]*Q[2]
2717    T4 = Q[0]*Q[3]
2718    T5 = -Q[1]*Q[1]
2719    T6 = Q[1]*Q[2]
2720    T7 = Q[1]*Q[3]
2721    T8 = -Q[2]*Q[2]
2722    T9 = Q[2]*Q[3]
2723    T10 = -Q[3]*Q[3]
2724    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
2725    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
2726    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2] 
2727    return VP   
2728   
2729def Q2Mat(Q):
2730    ''' make rotation matrix from quaternion
2731        q=r+ai+bj+ck
2732    '''
2733    QN = normQ(Q)
2734    aa = QN[0]**2
2735    ab = QN[0]*QN[1]
2736    ac = QN[0]*QN[2]
2737    ad = QN[0]*QN[3]
2738    bb = QN[1]**2
2739    bc = QN[1]*QN[2]
2740    bd = QN[1]*QN[3]
2741    cc = QN[2]**2
2742    cd = QN[2]*QN[3]
2743    dd = QN[3]**2
2744    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
2745        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
2746        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
2747    return np.array(M)
2748   
2749def AV2Q(A,V):
2750    ''' convert angle (radians) & vector to quaternion
2751        q=r+ai+bj+ck
2752    '''
2753    Q = np.zeros(4)
2754    d = np.sqrt(np.sum(np.array(V)**2))
2755    if d:
2756        V /= d
2757    else:
2758        V = np.array([0.,0.,1.])
2759    p = A/2.
2760    Q[0] = np.cos(p)
2761    Q[1:4] = V*np.sin(p)
2762    return Q
2763   
2764def AVdeg2Q(A,V):
2765    ''' convert angle (degrees) & vector to quaternion
2766        q=r+ai+bj+ck
2767    '''
2768    Q = np.zeros(4)
2769    d = np.sqrt(np.sum(np.array(V)**2))
2770    if d:
2771        V /= d
2772    else:
2773        V = np.array([0.,0.,1.])
2774    p = A/2.
2775    Q[0] = cosd(p)
2776    Q[1:4] = V*sind(p)
2777    return Q
2778   
2779def Q2AVdeg(Q):
2780    ''' convert quaternion to angle (degrees 0-360) & normalized vector
2781        q=r+ai+bj+ck
2782    '''
2783    A = 2.*acosd(Q[0])
2784    V = np.array(Q[1:])
2785    d = np.sqrt(np.sum(V**2))
2786    if d:
2787        V /= d
2788    else:
2789        A = 0.
2790        V = np.array([0.,0.,0.])
2791    return A,V
2792   
2793def Q2AV(Q):
2794    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
2795        q=r+ai+bj+ck
2796    '''
2797    A = 2.*np.arccos(Q[0])
2798    V = np.array(Q[1:])
2799    d = np.sqrt(np.sum(V**2))
2800    if d:
2801        V /= d
2802    else:
2803        A = 0.
2804        V = np.array([0.,0.,0.])
2805    return A,V
2806   
2807def makeQuat(A,B,C):
2808    ''' Make quaternion from rotation of A vector to B vector about C axis
2809
2810        :param np.array A,B,C: Cartesian 3-vectors
2811        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
2812    '''
2813
2814    V1 = np.cross(A,C)
2815    V2 = np.cross(B,C)
2816    if nl.norm(V1)*nl.norm(V2):
2817        V1 /= nl.norm(V1)
2818        V2 /= nl.norm(V2)
2819        V3 = np.cross(V1,V2)
2820    else:
2821        V3 = np.zeros(3)
2822    Q = np.array([0.,0.,0.,1.])
2823    D = 0.
2824    if nl.norm(V3):
2825        V3 /= nl.norm(V3)
2826        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
2827        D = np.arccos(D1)/2.0
2828        V1 = C-V3
2829        V2 = C+V3
2830        DM = nl.norm(V1)
2831        DP = nl.norm(V2)
2832        S = np.sin(D)
2833        Q[0] = np.cos(D)
2834        Q[1:] = V3*S
2835        D *= 2.
2836        if DM > DP:
2837            D *= -1.
2838    return Q,D
2839   
2840if __name__ == "__main__":
2841    from numpy import cos
2842    # minimum expected at ~-0.195
2843    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
2844    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
2845    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
2846    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
2847
2848    # minimum expected at ~[-0.195, -0.1]
2849    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
2850    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='cauchy')
2851    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='fast')
2852    print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='boltzmann')
Note: See TracBrowser for help on using the repository browser.