Changeset 808


Ignore:
Timestamp:
Dec 5, 2012 10:06:52 AM (10 years ago)
Author:
vondreele
Message:

Add Babinet modification to form factors - works OK
Use atom Ids for restraints
add 'all' in atom constraint editing routines
New add routines for bond restraints (angle add not done yet)
FindAtomIndexByIDs, FillAtomLookUp?, GetAtomsByID, & GetAtomItemsById? all now in GSASIImath.py
removed from GSASIIphsGUI.py
change phase tick mark colors
get GSASIIstruct.py to recognize protein atoms & do a protein calc.

Location:
trunk
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIgrid.py

    r807 r808  
    1414import sys
    1515import numpy as np
     16import numpy.ma as ma
    1617import os.path
    1718import wx.html        # could postpone this for quicker startup
    … …  
    2122import GSASIImath as G2mth
    2223import GSASIIIO as G2IO
     24import GSASIIlattice as G2lat
    2325import GSASIIplot as G2plt
    2426import GSASIIpwdGUI as G2pdG
    … …  
    2628import GSASIIphsGUI as G2phG
    2729import GSASIIstruct as G2str
     30import GSASIIspc as G2spc
    2831import GSASIImapvars as G2mv
    2932
    … …  
    13851388    phaseList.sort()
    13861389    phaseAtNames = {}
     1390    phaseAtTypes = {}
     1391    TypeList = []
    13871392    for item in phaseList:
    13881393        Split = item.split(':')
    13891394        if Split[2][:2] in ['AU','Af','dA']:
    1390             phaseAtNames[item] = AtomDict[int(Split[0])][int(Split[3])][0]
     1395            Id = int(Split[0])
     1396            phaseAtNames[item] = AtomDict[Id][int(Split[3])][0]
     1397            phaseAtTypes[item] = AtomDict[Id][int(Split[3])][1]
     1398            if phaseAtTypes[item] not in TypeList:
     1399                TypeList.append(phaseAtTypes[item])
    13911400        else:
    13921401            phaseAtNames[item] = ''
     1402            phaseAtTypes[item] = ''
    13931403           
    13941404    hapVary,hapDict,controlDict = G2str.GetHistogramPhaseData(Phases,Histograms,Print=False)
    … …  
    14521462        if page[1] == 'phs':
    14531463            atchoice = [item+' for '+phaseAtNames[item] for item in varList]
     1464            atchoice += [FrstVarb+' for all']
     1465            atchoice += [FrstVarb+' for all '+atype for atype in TypeList]
    14541466            dlg = wx.MultiChoiceDialog(G2frame,'Select more variables:'+legend,
    14551467                'Constrain '+FrstVarb+' and...',atchoice)
    … …  
    14601472        if dlg.ShowModal() == wx.ID_OK:
    14611473            sel = dlg.GetSelections()
    1462             for x in sel:
    1463                 varbs.append(varList[x])
     1474            try:
     1475                for x in sel:
     1476                    varbs.append(varList[x])
     1477            except IndexError:      # one of the 'all' chosen - supercedes any individual selection
     1478                varbs = [FrstVarb,]
     1479                Atypes = []
     1480                for x in sel:
     1481                    item = atchoice[x]
     1482                    if 'all' in item:
     1483                        Atypes.append(item.split('all')[1].strip())
     1484                if '' in Atypes:
     1485                    varbs += varList
     1486                else:
     1487                    for item in varList:
     1488                        if phaseAtTypes[item] in Atypes:
     1489                            varbs.append(item) 
    14641490        dlg.Destroy()
    14651491        if len(varbs) > 1:
    … …  
    18761902    if 'Chiral' not in restrData:
    18771903        restrData['Chiral'] = {'wtFactor':1.0,'Volumes':[],'Use':True}
     1904    General = phasedata['General']
     1905    Cell = General['Cell'][1:7]          #skip flag & volume   
     1906    Amat,Bmat = G2lat.cell2AB(Cell)
     1907    SGData = General['SGData']
     1908    cx,ct = General['AtomPtrs'][:2]
     1909    Atoms = phasedata['Atoms']
     1910    AtLookUp = G2mth.FillAtomLookUp(Atoms)
     1911    Names = ['all '+ name for name in General['AtomTypes']]
     1912    iBeg = len(Names)
     1913    Types = [name for name in General['AtomTypes']]
     1914    Coords = [ [] for type in Types]
     1915    Ids = [ 0 for type in Types]
     1916    Names += [atom[ct-1] for atom in Atoms]
     1917    Types += [atom[ct] for atom in Atoms]
     1918    Coords += [atom[cx:cx+3] for atom in Atoms]
     1919    Ids += [atom[-1] for atom in Atoms]
    18781920   
    18791921    def OnSelectPhase(event):
    … …  
    18891931        page = G2frame.dataDisplay.GetSelection()
    18901932        if 'Bond' in G2frame.dataDisplay.GetPageText(page):
    1891             AddBondRestraint()
     1933            bondRestData = restrData['Bond']
     1934            AddBondRestraint(bondRestData)
    18921935        elif 'Angle' in G2frame.dataDisplay.GetPageText(page):
    1893             AddAngleRestraint()
     1936            angleRestData = restrData['Angle']
     1937            AddAngleRestraint(angleRestData)
    18941938        elif 'Plane' in G2frame.dataDisplay.GetPageText(page):
    18951939            AddPlaneRestraint()
    … …  
    18971941            AddChiralRestraint()
    18981942           
    1899     def AddBondRestraint():
    1900         General = phasedata['General']
    1901         cx,ct = General['AtomPtrs'][:2]
    1902         Atoms = phasedata['Atoms']
    1903         radiiDict = dict(zip(General['AtomTypes'],General['BondRadii']))
    1904         Names = ['all '+ name for name in General['AtomTypes']]
    1905         Types = [name for name in General['AtomTypes']]
    1906         Names += [atom[ct-1] for atom in Atoms]
    1907         Types += [atom[ct] for atom in Atoms]
     1943    def AddBondRestraint(bondRestData):
     1944        Radii = dict(zip(General['AtomTypes'],General['BondRadii']))
    19081945        Lists = {'origin':[],'target':[]}
    19091946        for listName in ['origin','target']:
    1910             dlg = wx.MultiChoiceDialog(G2frame,'Bond restraint'+listName+' for '+General['Name'],
     1947            dlg = wx.MultiChoiceDialog(G2frame,'Bond restraint '+listName+' for '+General['Name'],
    19111948                    'Select bond restraint '+listName+' atoms',Names)
    19121949            if dlg.ShowModal() == wx.ID_OK:
    19131950                sel = dlg.GetSelections()
    19141951                for x in sel:
    1915                     Lists[listName].append([Names[x],Types[x]])
    1916         SGData = General['SGData']
    1917         Bonds = restrData['Bond']['Bonds']
    1918 
    1919     def AddAngleRestraint():
    1920         print 'Angle restraint'
     1952                    if 'all' in Names[x]:
     1953                        allType = Types[x]
     1954                        for name,Type,coords,id in zip(Names,Types,Coords,Ids):
     1955                            if Type == allType and 'all' not in name:
     1956                                Lists[listName].append([id,Type,coords])
     1957                    else:
     1958                        Lists[listName].append([Ids[x],Types[x],Coords[x],])
     1959        Factor = .85
     1960        indices = (-1,0,1)
     1961        Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
     1962        origAtoms = Lists['origin']
     1963        targAtoms = Lists['target']
     1964        for Oid,Otype,Ocoord in origAtoms:
     1965            for Tid,Ttype,Tcoord in targAtoms:
     1966                result = G2spc.GenAtom(Tcoord,SGData,False,Move=False)
     1967                BsumR = (Radii[Otype]+Radii[Ttype])*Factor
     1968                for Txyz,Top,Tunit in result:
     1969                    Dx = (Txyz-np.array(Ocoord))+Units
     1970                    dx = np.inner(Amat,Dx)
     1971                    dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
     1972                    IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
     1973                    if np.any(IndB):
     1974                        for indb in IndB:
     1975                            for i in range(len(indb)):
     1976                                unit = Units[indb][i]+Tunit
     1977                                if np.any(unit):
     1978                                    Topstr = '%d+%d,%d,%d'%(Top,unit[0],unit[1],unit[2])
     1979                                else:
     1980                                    Topstr = str(Top)
     1981                                bondRestData['Bonds'].append([[Oid,Tid],['1',Topstr], \
     1982                                    ma.getdata(dist[indb])[i],1.54,0.01])
     1983        UpdateBondRestr(bondRestData)               
     1984
     1985    def AddAngleRestraint(angleRestData):
     1986        Radii = dict(zip(General['AtomTypes'],General['AngleRadii']))
     1987        origAtoms = []
     1988        targAtoms = [[Ids[x],Types[x],Coords[x]] for x in enumerate(Names[iBeg:])]
     1989        dlg = wx.MultiChoiceDialog(G2frame,'Select atom B for angle A-B-C for '+General['Name'],
     1990                'Select angle restraint origin atoms',Names)
     1991        if dlg.ShowModal() == wx.ID_OK:
     1992            sel = dlg.GetSelections()
     1993            for x in sel:
     1994                if 'all' in Names[x]:
     1995                    allType = Types[x]
     1996                    for name,Type,coords,id in zip(Names,Types,Coords,Ids):
     1997                        if Type == allType and 'all' not in name:
     1998                            origAtoms.append([id,Type,coords])
     1999                else:
     2000                    origAtoms.append([Ids[x],Types[x],Coords[x]])
     2001
     2002        Factor = .85
     2003        indices = (-1,0,1)
     2004        Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
     2005        for Oid,Otype,Ocoord in origAtoms:
     2006            IndBlist = []
     2007            Dist = []
     2008            Vect = []
     2009            VectA = []
     2010            angles = []
     2011            for Tid,Ttype,Tcoord in targAtoms:
     2012                result = G2spc.GenAtom(Tcoord,SGData,False,Move=False)
     2013                BsumR = (Radii[Otype]+Radii[Ttype])*Factor
     2014                AsumR = (Radii[Otype]+Radii[Ttype])*Factor
     2015                for Txyz,Top,Tunit in result:
     2016                    Dx = (Txyz-np.array(Ocoord))+Units
     2017                    dx = np.inner(Amat,Dx)
     2018                    dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
     2019                    IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
     2020                    if np.any(IndB):
     2021                        for indb in IndB:
     2022                            for i in range(len(indb)):
     2023                                if str(dx.T[indb][i]) not in IndBlist:
     2024                                    IndBlist.append(str(dx.T[indb][i]))
     2025                                    unit = Units[indb][i]+Tunit
     2026                                if np.any(unit):
     2027                                    Topstr = '%d+%d,%d,%d'%(Top,unit[0],unit[1],unit[2])
     2028                                else:
     2029                                    Topstr = str(Top)
     2030                                    tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
     2031                                    Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig])
     2032                                    if (Dist[-1][-1]-AsumR) <= 0.:
     2033                                        Vect.append(dx.T[indb][i]/Dist[-1][-2])
     2034                                        VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
     2035                                    else:
     2036                                        Vect.append([0.,0.,0.])
     2037                                        VectA.append([])
     2038            Vect = np.array(Vect)
     2039            angles = np.zeros((len(Vect),len(Vect)))
     2040            angsig = np.zeros((len(Vect),len(Vect)))
     2041            for i,veca in enumerate(Vect):
     2042                if np.any(veca):
     2043                    for j,vecb in enumerate(Vect):
     2044                        if np.any(vecb):
     2045                            angles[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData)[0]
     2046
     2047
    19212048
    19222049    def AddPlaneRestraint():
    … …  
    20102137
    20112138        bondList = bondRestData['Bonds']
     2139        if len(bondList[0]) == 6:   #patch
     2140            bondList = bondRestData['Bonds'] = []
    20122141        if len(bondList):
    20132142            table = []
    … …  
    20152144            colLabels = ['A+SymOp  B+SymOp','d-calc','d-obs','esd']
    20162145            Types = [wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,3',]
    2017             for i,[atoms,ops,indx,dcalc,dobs,esd] in enumerate(bondList):
     2146            for i,[indx,ops,dcalc,dobs,esd] in enumerate(bondList):
     2147                atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
    20182148                table.append([atoms[0]+'+ ('+ops[0]+')  '+atoms[1]+'+ ('+ops[1]+')',dcalc,dobs,esd])
    20192149                rowLabels.append(str(i))
    … …  
    21042234            colLabels = ['A+SymOp  B+SymOp  C+SymOp','calc','obs','esd']
    21052235            Types = [wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,2',]
    2106             for i,[atoms,ops,indx,dcalc,dobs,esd] in enumerate(angleList):
     2236            for i,[indx,ops,dcalc,dobs,esd] in enumerate(angleList):
     2237                atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
    21072238                table.append([atoms[0]+'+ ('+ops[0]+')  '+atoms[1]+'+ ('+ops[1]+')  '+atoms[2]+ \
    21082239                '+ ('+ops[2]+')',dcalc,dobs,esd])
    … …  
    21602291            colLabels = ['atom+SymOp','calc','obs','esd']
    21612292            Types = [wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,2',]
    2162             for i,[atoms,ops,indx,dcalc,dobs,esd] in enumerate(planeList):
     2293            for i,[indx,ops,dcalc,dobs,esd] in enumerate(planeList):
     2294                atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
    21632295                atString = ''
    21642296                for a,atom in enumerate(atoms):
    … …  
    22132345            colLabels = ['O+SymOp  A+SymOp  B+SymOp  C+SymOp','calc','obs','esd']
    22142346            Types = [wg.GRID_VALUE_STRING,]+3*[wg.GRID_VALUE_FLOAT+':10,2',]
    2215             for i,[atoms,ops,indx,dcalc,dobs,esd] in enumerate(volumeList):
     2347            for i,[indx,ops,dcalc,dobs,esd] in enumerate(volumeList):
     2348                atoms = G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,ct-1)
    22162349                table.append([atoms[0]+'+ ('+ops[0]+') '+atoms[1]+'+ ('+ops[1]+') '+atoms[2]+ \
    22172350                '+ ('+ops[2]+') '+atoms[3]+'+ ('+ops[3]+')',dcalc,dobs,esd])

  • trunk/GSASIImath.py

    r805 r808  
    155155    return vcov
    156156
     157def FindAtomIndexByIDs(atomData,IDs,Draw=True):
     158    indx = []
     159    for i,atom in enumerate(atomData):
     160        if Draw and atom[-3] in IDs:
     161            indx.append(i)
     162        elif atom[-1] in IDs:
     163            indx.append(i)
     164    return indx
     165
     166def FillAtomLookUp(atomData):
     167    atomLookUp = {}
     168    for iatm,atom in enumerate(atomData):
     169        atomLookUp[atom[-1]] = iatm
     170    return atomLookUp
     171
     172def GetAtomsById(atomData,atomLookUp,IdList):
     173    atoms = []
     174    for id in IdList:
     175        atoms.append(atomData[atomLookUp[id]])
     176    return atoms
     177   
     178def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
     179    Items = []
     180    for id in IdList:
     181        if numItems == 1:
     182            Items.append(atomData[atomLookUp[id]][itemLoc])
     183        else:
     184            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
     185    return Items
     186       
    157187def getMass(generalData):
    158188    mass = 0.

  • trunk/GSASIIphsGUI.py

    r807 r808  
    368368            evt.Skip()
    369369
    370 def FindAtomIndexByIDs(atomData,IDs,Draw=True):
    371     indx = []
    372     for i,atom in enumerate(atomData):
    373         if Draw and atom[-3] in IDs:
    374             indx.append(i)
    375         elif atom[-1] in IDs:
    376             indx.append(i)
    377     return indx
    378        
    379370def UpdatePhaseData(G2frame,Item,data,oldPage):
    380371
    … …  
    15461537            SGData = generalData['SGData']
    15471538            dlg = SymOpDialog(G2frame,SGData,True,True)
     1539            New = False
    15481540            try:
    15491541                if dlg.ShowModal() == wx.ID_OK:
    … …  
    16141606                DisAglData['pId'] = data['pId']
    16151607                DisAglData['covData'] = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.root, 'Covariance'))
    1616             G2str.DistAngle(DisAglCtls,DisAglData)
     1608            try:
     1609                G2str.DistAngle(DisAglCtls,DisAglData)
     1610            except KeyError:        # inside DistAngle for missing atom types in DisAglCtls
     1611                print '**** ERROR - try again but do "Reset" to fill in missing atom types ****'
    16171612                       
    16181613################################################################################
    … …  
    17161711        IDs = [ID,]
    17171712        atomData = drawingData['Atoms']
    1718         indx = FindAtomIndexByIDs(atomData,IDs)
     1713        indx = G2mth.FindAtomIndexByIDs(atomData,IDs)
    17191714        for ind in indx:
    17201715            atomData[ind] = MakeDrawAtom(atom,atomData[ind])
    … …  
    17341729        atIndx = []
    17351730        for item in indx:
    1736             atNames.append(atomData[item][ct-1])
    17371731            atXYZ.append(np.array(atomData[item][cx:cx+3]))
    17381732            atSymOp.append(atomData[item][cs-1])
    … …  
    17461740                restData[PhaseName]['Bond'] = bondData
    17471741            dist = G2mth.getRestDist(atXYZ,Amat)
    1748             bondData['Bonds'].append([atNames,atSymOp,atIndx,dist,1.54,0.01])
     1742            bondData['Bonds'].append([atIndx,atSymOp,dist,1.54,0.01])
    17491743        elif event.GetId() == G2gd.wxID_DRAWRESTRANGLE and len(indx) == 3:
    17501744            try:
    … …  
    17551749                restData[PhaseName]['Angle'] = angleData
    17561750            angle = G2mth.getRestAngle(atXYZ,Amat)
    1757             angleData['Angles'].append([atNames,atSymOp,atIndx,angle,109.5,1.0])           
     1751            angleData['Angles'].append([atIndx,atSymOp,angle,109.5,1.0])           
    17581752        elif event.GetId() == G2gd.wxID_DRAWRESTRPLANE and len(indx) > 3:
    17591753            try:
    … …  
    17641758                restData[PhaseName]['Plane'] = planeData
    17651759            plane = G2mth.getRestPlane(atXYZ,Amat)
    1766             planeData['Planes'].append([atNames,atSymOp,atIndx,plane,0.0,0.01])           
     1760            planeData['Planes'].append([atIndx,atSymOp,plane,0.0,0.01])           
    17671761        elif event.GetId() == G2gd.wxID_DRAWRESTRCHIRAL and len(indx) == 4:
    17681762            try:
    … …  
    17731767                restData[PhaseName]['Chiral'] = chiralData
    17741768            volume = G2mth.getRestChiral(atXYZ,Amat)
    1775             chiralData['Volumes'].append([atNames,atSymOp,atIndx,volume,2.5,0.1])           
     1769            chiralData['Volumes'].append([atIndx,atSymOp,volume,2.5,0.1])           
    17761770        else:
    17771771            print '**** ERROR wrong number of atoms selected for this restraint'
    … …  
    24022396        event.StopPropagation()
    24032397       
    2404     def FindAtomIndexByIDs(atomData,IDs,Draw=True):
    2405         indx = []
    2406         for i,atom in enumerate(atomData):
    2407             if Draw and atom[-3] in IDs:
    2408                 indx.append(i)
    2409             elif atom[-1] in IDs:
    2410                 indx.append(i)
    2411         return indx
    2412        
    24132398    def DrawAtomsDeleteByIDs(IDs):
    24142399        atomData = data['Drawing']['Atoms']
    2415         indx = FindAtomIndexByIDs(atomData,IDs)
     2400        indx = G2mth.FindAtomIndexByIDs(atomData,IDs)
    24162401        indx.reverse()
    24172402        for ind in indx:
    … …  
    24272412        elif colName == 'I/A':
    24282413            col = cs
    2429         indx = FindAtomIndexByIDs(atomData,IDs)
     2414        indx = G2mth.FindAtomIndexByIDs(atomData,IDs)
    24302415        for ind in indx:
    24312416            atomData[ind][col] = value
    … …  
    24942479            atom = atomDData[i]
    24952480            xyz.append([i,]+atom[cn:cn+2]+atom[cx:cx+4]) #also gets Sym Op
    2496             id = FindAtomIndexByIDs(atomData,[atom[cid],],False)[0]
     2481            id = G2mth.FindAtomIndexByIDs(atomData,[atom[cid],],False)[0]
    24972482            Oxyz.append([id,]+atomData[id][cx+1:cx+4])
    24982483        DATData['Datoms'] = xyz
    … …  
    31643149            hist = Indx[Obj.GetId()]
    31653150            sourceDict = UseList[hist]
    3166             copyNames = ['Scale','Pref.Ori.','Size','Mustrain','HStrain','Extinction']
     3151            copyNames = ['Scale','Pref.Ori.','Size','Mustrain','HStrain','Extinction','Babinet']
    31673152            copyDict = {}
    31683153            for name in copyNames:
    … …  
    31923177            sourceDict = UseList[hist]
    31933178            copyDict = {}
    3194             copyNames = ['Scale','Pref.Ori.','Size','Mustrain','HStrain','Extinction']
     3179            copyNames = ['Scale','Pref.Ori.','Size','Mustrain','HStrain','Extinction','Babinet']
     3180            babNames = ['BabA','BabU']
    31953181            for name in copyNames:
    31963182                if name in ['Scale','Extinction','HStrain']:
    … …  
    32053191                        for item in SHterms:
    32063192                            SHflags[item] = SHterms[item][1]
    3207                         copyDict[name].append(SHflags)                       
     3193                        copyDict[name].append(SHflags)
     3194                elif name == 'Babinet':
     3195                    copyDict[name] = {}
     3196                    for bab in babNames:
     3197                        copyDict[name][bab] = sourceDict[name][bab][1]                       
    32083198            keyList = ['All',]+UseList.keys()
    32093199            if UseList:
    … …  
    32353225                                        SHterms = copy.copy(sourceDict[name][5])
    32363226                                        for item in SHflags:
    3237                                             SHterms[item][1] = copy.copy(SHflags[item])                                             
     3227                                            SHterms[item][1] = copy.copy(SHflags[item])
     3228                                elif name == 'Babinet':
     3229                                    for bab in babNames:
     3230                                        UseList[item][name][bab][1] = copy.copy(copyDict[name][bab])                                             
    32383231                        wx.CallAfter(UpdateDData)
    32393232                finally:
    … …  
    34483441                pass
    34493442            Obj.SetValue("%.2f"%(UseList[Indx[Obj.GetId()]]['Extinction'][0]))
    3450            
     3443
     3444        def OnBabRef(event):
     3445            Obj = event.GetEventObject()
     3446            item,bab = Indx[Obj.GetId()]
     3447            UseList[item]['Babinet']['Bab'+bab][1] = Obj.GetValue()
     3448           
     3449        def OnBabVal(event):
     3450            Obj = event.GetEventObject()
     3451            item,bab = Indx[Obj.GetId()]
     3452            try:
     3453                val = float(Obj.GetValue())
     3454                if val >= 0:
     3455                    UseList[item]['Babinet']['Bab'+bab][0] = val
     3456            except ValueError:
     3457                pass
     3458            Obj.SetValue("%.3f"%(UseList[item]['Babinet']['Bab'+bab][0]))
     3459
    34513460        def OnTbarVal(event):
    34523461            Obj = event.GetEventObject()
    … …  
    37463755            return extSizer
    37473756           
     3757        def BabSizer():
     3758            babSizer = wx.BoxSizer(wx.HORIZONTAL)
     3759            for bab in ['A','U']:
     3760                babRef = wx.CheckBox(DData,-1,label=' Babinet '+bab+': ')
     3761                babRef.SetValue(UseList[item]['Babinet']['Bab'+bab][1])
     3762                Indx[babRef.GetId()] = [item,bab]
     3763                babRef.Bind(wx.EVT_CHECKBOX, OnBabRef)
     3764                babSizer.Add(babRef,0,wx.ALIGN_CENTER_VERTICAL)
     3765                babVal = wx.TextCtrl(DData,wx.ID_ANY,
     3766                    '%.3f'%(UseList[item]['Babinet']['Bab'+bab][0]),style=wx.TE_PROCESS_ENTER)
     3767                Indx[babVal.GetId()] = [item,bab]
     3768                babVal.Bind(wx.EVT_TEXT_ENTER,OnBabVal)
     3769                babVal.Bind(wx.EVT_KILL_FOCUS,OnBabVal)
     3770                babSizer.Add(babVal,0,wx.ALIGN_CENTER_VERTICAL)
     3771                babSizer.Add((5,5),0)
     3772            return babSizer
     3773           
    37483774        if DData.GetSizer():
    37493775            DData.GetSizer().Clear(True)
    … …  
    37573783            if 'Use' not in UseList[item]:      #patch
    37583784                UseList[item]['Use'] = True
     3785            if 'Babinet' not in UseList[item]:
     3786                UseList[item]['Babinet'] = {'BabA':[0.0,False],'BabU':[0.0,False]}
    37593787            showSizer = wx.BoxSizer(wx.HORIZONTAL)
    37603788            showData = wx.CheckBox(DData,-1,label=' Show '+item)
    … …  
    38543882                mainSizer.Add(ExtSizer())
    38553883                mainSizer.Add((0,5),0)
     3884                mainSizer.Add(BabSizer())
     3885                mainSizer.Add((0,5),0)
    38563886            elif item[:4] == 'HKLF' and UseList[item]['Show']:
    38573887                mainSizer.Add((0,5),0)               
    38583888                mainSizer.Add(SCExtSizer())
     3889                mainSizer.Add((0,5),0)
     3890                mainSizer.Add(BabSizer())
    38593891                mainSizer.Add((0,5),0)
    38603892                pass
    … …  
    38893921                        histoName = TextList[i]
    38903922                        UseList[histoName] = {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
     3923                            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
    38913924                            'Extinction':['Lorentzian','Secondary Type I',
    38923925                            {'Tbar':0.0,'Eg':[0.0,False],'Es':[0.0,False],'Ep':[0.0,False]},]}                       
    … …  
    39443977                                NShkl*[0.01,],NShkl*[False,]],
    39453978                            'HStrain':[NDij*[0.0,],NDij*[False,]],                         
    3946                             'Extinction':[0.0,False]}
     3979                            'Extinction':[0.0,False],'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}}
    39473980                        refList = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,pId,'Reflection Lists'))
    39483981                        refList[generalData['Name']] = []                       

  • trunk/GSASIIplot.py

    r807 r808  
    798798        elif G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists'] or \
    799799            'PWDR' in G2frame.PatternTree.GetItemText(PickId):
     800            refColors=['b','r','c','g','m','k']
    800801            Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists'))
    801802            for pId,phase in enumerate(Phases):
    … …  
    806807                pos = G2frame.refOffset-pId*Ymax*G2frame.refDelt*np.ones_like(peak)
    807808                if G2frame.qPlot:
    808                     Plot.plot(2*np.pi/peak.T[0],pos,colors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
     809                    Plot.plot(2*np.pi/peak.T[0],pos,refColors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
    809810                else:
    810                     Plot.plot(peak.T[1],pos,colors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
     811                    Plot.plot(peak.T[1],pos,refColors[pId%6]+'|',mew=1,ms=8,picker=3.,label=phase)
    811812            if len(Phases):
    812813                handles,legends = Plot.get_legend_handles_labels()  #got double entries in the legends for some reason

  • trunk/GSASIIpwd.py

    r803 r808  
    522522            iD += 1       
    523523        except KeyError:
     524            break
     525        except ValueError:
     526            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
    524527            break       
    525528    return yb
    … …  
    616619        except KeyError:
    617620            break
     621        except ValueError:
     622            print '**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****'
     623            break       
    618624    return dydb,dyddb,dydpk
    619625

  • trunk/GSASIIpwdGUI.py

    r802 r808  
    675675    def inst2data(inst,ref,data):
    676676        for item in data:
    677             data[item] = [data[item][0],inst[item],ref[item]]
     677            try:
     678                data[item] = [data[item][0],inst[item],ref[item]]
     679            except KeyError:
     680                pass        #skip 'Polariz.' for N-data
    678681        return data
    679682       
    … …  
    722725                RefreshInstrumentGrid(event,doAnyway=True)          #to get peaks updated
    723726                UpdateInstrumentGrid(G2frame,data)
     727                G2plt.PlotPeakWidths(G2frame)
    724728        finally:
    725729            dlg.Destroy()
    … …  
    749753        RefreshInstrumentGrid(event,doAnyway=True)          #to get peaks updated
    750754        UpdateInstrumentGrid(G2frame,data)
     755        G2plt.PlotPeakWidths(G2frame)
    751756       
    752757    def OnInstFlagCopy(event):
    … …  
    769774            for item in copyList:
    770775                Id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,item)
    771                 instData = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'))
     776                instData = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'))[0]
    772777                if len(data) == len(instData) and instType == instData['Type'][0]:   #don't mix data types or lam & lam1/lam2 parms!
    773                    
    774                     instData[2] = copy.copy(flags)
     778                    for item in instData:
     779                        instData[item][2] = copy.copy(flags[item])
    775780                else:
    776781                    print item+' not copied - instrument parameters not commensurate'
    … …  
    796801            for item in copyList:
    797802                Id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,item)
    798                 instData = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'))
     803                instData = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'))[0]
    799804                if len(data) == len(instData) and instType == instData['Type'][0]:  #don't mix data types or lam & lam1/lam2 parms!
    800                     for i,item in enumerate(data[1:]):                  #skip default values in tuple
    801                         instData[i+1][:-1] = copy.copy(item[:-1])       #skip azimuth at end
     805                    instData.update(data)
    802806                else:
    803807                    print item+' not copied - instrument parameters not commensurate'
    … …  
    863867        Obj.SetValue(fmt%(value))
    864868        data = updateData(insVal,insRef)
     869        G2plt.PlotPeakWidths(G2frame)
    865870       
    866871    def OnItemRef(event):

  • trunk/GSASIIspc.py

    r762 r808  
    988988        'P m 3 n','P n 3 m',),
    989989    'Im3m':('I 2 3','I 21 3','I m -3','I a -3', 'I 4 3 2','I 41 3 2',
    990         'I -4 3 m', 'I -4 3 d','I m -3 m','I a -3 d',),
     990        'I -4 3 m', 'I -4 3 d','I m -3 m','I m 3 m','I a -3 d',),
    991991    'Fm3m':('F 2 3','F m -3','F d -3','F 4 3 2','F 41 3 2','F -4 3 m',
    992         'F -4 3 c','F m -3 m','F m -3 c','F d -3 m','F d -3 c',),
     992        'F -4 3 c','F m -3 m','F m 3 m','F m -3 c','F d -3 m','F d -3 c',),
    993993}
    994994'A few non-standard space groups for test use'

  • trunk/GSASIIstruct.py

    r803 r808  
    543543               
    544544    def PrintAtoms(General,Atoms):
     545        cx,ct,cs,cia = General['AtomPtrs']
    545546        print >>pFile,'\n Atoms:'
    546547        line = '   name    type  refine?   x         y         z    '+ \
    … …  
    549550            line += '   Mx     My     Mz'
    550551        elif General['Type'] == 'macromolecular':
    551             line = ' res no  residue  chain '+line
     552            line = ' res no residue chain'+line
    552553        print >>pFile,line
    553554        if General['Type'] == 'nuclear':
    554555            print >>pFile,135*'-'
    555556            for i,at in enumerate(Atoms):
    556                 line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
    557                     '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
    558                 if at[9] == 'I':
    559                     line += '%8.4f'%(at[10])+48*' '
     557                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
     558                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
     559                if at[cia] == 'I':
     560                    line += '%8.4f'%(at[cia+1])+48*' '
    560561                else:
    561562                    line += 8*' '
    562563                    for j in range(6):
    563                         line += '%8.4f'%(at[11+j])
     564                        line += '%8.4f'%(at[cia+1+j])
     565                print >>pFile,line
     566        elif General['Type'] == 'macromolecular':
     567            print >>pFile,135*'-'           
     568            for i,at in enumerate(Atoms):
     569                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
     570                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
     571                if at[cia] == 'I':
     572                    line += '%8.4f'%(at[cia+1])+48*' '
     573                else:
     574                    line += 8*' '
     575                    for j in range(6):
     576                        line += '%8.4f'%(at[cia+1+j])
    564577                print >>pFile,line
    565578       
    … …  
    639652        Natoms[pfx] = 0
    640653        if Atoms and not General.get('doPawley'):
    641             if General['Type'] == 'nuclear':
     654            cx,ct,cs,cia = General['AtomPtrs']
     655            if General['Type'] in ['nuclear','macromolecular']:
    642656                Natoms[pfx] = len(Atoms)
    643657                for i,at in enumerate(Atoms):
    644658                    atomIndx[at[-1]] = [pfx,i]      #lookup table for restraints
    645                     phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
    646                         pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
     659                    phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
     660                        pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
    647661                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
    648                         pfx+'AI/A:'+str(i):at[9],})
    649                     if at[9] == 'I':
    650                         phaseDict[pfx+'AUiso:'+str(i)] = at[10]
     662                        pfx+'AI/A:'+str(i):at[cia],})
     663                    if at[cia] == 'I':
     664                        phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
    651665                    else:
    652                         phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
    653                             pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
    654                     if 'F' in at[2]:
     666                        phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
     667                            pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
     668                    if 'F' in at[ct+1]:
    655669                        phaseVary.append(pfx+'Afrac:'+str(i))
    656                     if 'X' in at[2]:
    657                         xId,xCoef = G2spc.GetCSxinel(at[7])
     670                    if 'X' in at[ct+1]:
     671                        xId,xCoef = G2spc.GetCSxinel(at[cs])
    658672                        names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
    659673                        equivs = [[],[],[]]
    … …  
    667681                                for eqv in equiv[1:]:
    668682                                    G2mv.StoreEquivalence(name,(eqv,))
    669                     if 'U' in at[2]:
     683                    if 'U' in at[ct+1]:
    670684                        if at[9] == 'I':
    671685                            phaseVary.append(pfx+'AUiso:'+str(i))
    672686                        else:
    673                             uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
     687                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
    674688                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
    675689                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
    … …  
    10851099        print >>pFile,ptstr
    10861100   
     1101    def PrintBabinet(hapData):
     1102        print >>pFile,'\n Babinet form factor modification: '
     1103        ptlbls = ' names :'
     1104        ptstr =  ' values:'
     1105        varstr = ' refine:'
     1106        for name in ['BabA','BabU']:
     1107            ptlbls += '%12s' % (name)
     1108            ptstr += '%12.6f' % (hapData[name][0])
     1109            varstr += '%12s' % (str(hapData[name][1]))
     1110        print >>pFile,ptlbls
     1111        print >>pFile,ptstr
     1112        print >>pFile,varstr
     1113       
     1114   
    10871115    hapDict = {}
    10881116    hapVary = []
    … …  
    11441172                for item in ['Mustrain','Size']:
    11451173                    controlDict[pfx+item+'Type'] = hapData[item][0]
    1146                     hapDict[pfx+item+':mx'] = hapData[item][1][2]
     1174                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
    11471175                    if hapData[item][2][2]:
    1148                         hapVary.append(pfx+item+':mx')
     1176                        hapVary.append(pfx+item+';mx')
    11491177                    if hapData[item][0] in ['isotropic','uniaxial']:
    1150                         hapDict[pfx+item+':i'] = hapData[item][1][0]
     1178                        hapDict[pfx+item+';i'] = hapData[item][1][0]
    11511179                        if hapData[item][2][0]:
    1152                             hapVary.append(pfx+item+':i')
     1180                            hapVary.append(pfx+item+';i')
    11531181                        if hapData[item][0] == 'uniaxial':
    11541182                            controlDict[pfx+item+'Axis'] = hapData[item][3]
    1155                             hapDict[pfx+item+':a'] = hapData[item][1][1]
     1183                            hapDict[pfx+item+';a'] = hapData[item][1][1]
    11561184                            if hapData[item][2][1]:
    1157                                 hapVary.append(pfx+item+':a')
     1185                                hapVary.append(pfx+item+';a')
    11581186                    else:       #generalized for mustrain or ellipsoidal for size
    11591187                        Nterms = len(hapData[item][4])
    … …  
    11701198                            if hapData[item][5][i]:
    11711199                                hapVary.append(pfx+item+sfx)
     1200                for bab in ['BabA','BabU']:
     1201                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
     1202                    if hapData['Babinet'][bab][1]:
     1203                        hapVary.append(pfx+bab)
    11721204                               
    11731205                if Print:
    … …  
    11851217                    PrintMuStrain(hapData['Mustrain'],SGData)
    11861218                    PrintHStrain(hapData['HStrain'],SGData)
     1219                    PrintBabinet(hapData['Babinet'])
    11871220                HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
    11881221                refList = []
    … …  
    12101243                if hapData['Scale'][1]:
    12111244                    hapVary.append(pfx+'Scale')
     1245                               
    12121246                extApprox,extType,extParms = hapData['Extinction']
    12131247                controlDict[pfx+'EType'] = extType
    … …  
    12251259                    if extParms[eKey][1]:
    12261260                        hapVary.append(pfx+eKey)
     1261                for bab in ['BabA','BabU']:
     1262                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
     1263                    if hapData['Babinet'][bab][1]:
     1264                        hapVary.append(pfx+bab)
    12271265                if Print:
    12281266                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
    … …  
    12341272                        text += ' %4s : %10.3g Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
    12351273                    print >>pFile,text
     1274                    PrintBabinet(hapData['Babinet'])
    12361275                Histogram['Reflection Lists'] = phase       
    12371276               
    … …  
    13431382            print >>pFile,sigstr
    13441383       
    1345     def PrintSHPOAndSig(pdx,hapData,POsig):
     1384    def PrintSHPOAndSig(pfx,hapData,POsig):
    13461385        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
    13471386        ptlbls = ' names :'
    … …  
    13581397        print >>pFile,ptstr
    13591398        print >>pFile,sigstr
     1399       
     1400    def PrintBabinetAndSig(pfx,hapData,BabSig):
     1401        print >>pFile,'\n Babinet form factor modification: '
     1402        ptlbls = ' names :'
     1403        ptstr =  ' values:'
     1404        sigstr = ' sig   :'
     1405        for item in hapData:
     1406            ptlbls += '%12s'%(item)
     1407            ptstr += '%12.3f'%(hapData[item][0])
     1408            if pfx+item in BabSig:
     1409                sigstr += '%12.3f'%(BabSig[pfx+item])
     1410            else:
     1411                sigstr += 12*' '
     1412        print >>pFile,ptlbls
     1413        print >>pFile,ptstr
     1414        print >>pFile,sigstr
    13601415   
    13611416    PhFrExtPOSig = {}
    13621417    SizeMuStrSig = {}
    13631418    ScalExtSig = {}
     1419    BabSig = {}
    13641420    wtFrSum = {}
    13651421    for phase in Phases:
    … …  
    14001456                    pfx+'HStrain':{}})                 
    14011457                for item in ['Mustrain','Size']:
    1402                     hapData[item][1][2] = parmDict[pfx+item+':mx']
     1458                    hapData[item][1][2] = parmDict[pfx+item+';mx']
    14031459                    hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
    1404                     if pfx+item+':mx' in sigDict:
    1405                         SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+':mx']
     1460                    if pfx+item+';mx' in sigDict:
     1461                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
    14061462                    if hapData[item][0] in ['isotropic','uniaxial']:                   
    1407                         hapData[item][1][0] = parmDict[pfx+item+':i']
     1463                        hapData[item][1][0] = parmDict[pfx+item+';i']
    14081464                        if item == 'Size':
    14091465                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
    1410                         if pfx+item+':i' in sigDict:
    1411                             SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+':i']
     1466                        if pfx+item+';i' in sigDict:
     1467                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
    14121468                        if hapData[item][0] == 'uniaxial':
    1413                             hapData[item][1][1] = parmDict[pfx+item+':a']
     1469                            hapData[item][1][1] = parmDict[pfx+item+';a']
    14141470                            if item == 'Size':
    14151471                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
    1416                             if pfx+item+':a' in sigDict:
    1417                                 SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+':a']
     1472                            if pfx+item+';a' in sigDict:
     1473                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
    14181474                    else:       #generalized for mustrain or ellipsoidal for size
    14191475                        Nterms = len(hapData[item][4])
    … …  
    14281484                    if pfx+name in sigDict:
    14291485                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
     1486                for name in ['BabA','BabU']:
     1487                    hapData['Babinet'][name][0] = parmDict[pfx+name]
     1488                    if pfx+name in sigDict:
     1489                        BabSig[pfx+name] = sigDict[pfx+name]               
    14301490               
    14311491            elif 'HKLF' in histogram:
    … …  
    14351495                        if pfx+item in sigDict:
    14361496                            ScalExtSig[pfx+item] = sigDict[pfx+item]
     1497                for name in ['BabA','BabU']:
     1498                    hapData['Babinet'][name][0] = parmDict[pfx+name]
     1499                    if pfx+name in sigDict:
     1500                        BabSig[pfx+name] = sigDict[pfx+name]               
    14371501
    14381502    if Print:
    … …  
    14741538                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
    14751539                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
     1540                    PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
    14761541                   
    14771542                elif 'HKLF' in histogram:
    … …  
    14811546                    if 'Scale' in ScalExtSig:
    14821547                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
     1548                    PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
    14831549
    14841550# fix after it runs!               
    … …  
    19712037    twopi = 2.0*np.pi
    19722038    twopisq = 2.0*np.pi**2
     2039    phfx = pfx.split(':')[0]+hfx
    19732040    ast = np.sqrt(np.diag(G))
    19742041    Mast = twopisq*np.multiply.outer(ast,ast)
    … …  
    19892056        H = refl[:3]
    19902057        SQ = 1./(2.*refl[4])**2
     2058        SQfactor = 4.0*SQ*twopisq
     2059        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
    19912060        if not len(refl[-1]):                #no form factors
    19922061            if 'N' in parmDict[hfx+'Type']:
    … …  
    19962065        for i,El in enumerate(Tdata):
    19972066            FF[i] = refl[-1][El]           
    1998         SQfactor = 4.0*SQ*twopisq
    19992067        Uniq = refl[11]
    20002068        phi = refl[12]
    … …  
    20082076        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    20092077        Tcorr = Tiso*Tuij
    2010         fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
     2078        fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
    20112079        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
    20122080        if not SGData['SGInv']:
    2013             fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
     2081            fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
    20142082            fbs = np.sum(np.sum(fb,axis=1),axis=1)
    20152083        fasq = fas**2
    … …  
    20222090    twopi = 2.0*np.pi
    20232091    twopisq = 2.0*np.pi**2
     2092    phfx = pfx.split(':')[0]+hfx
    20242093    ast = np.sqrt(np.diag(G))
    20252094    Mast = twopisq*np.multiply.outer(ast,ast)
    … …  
    20422111    dFdui = np.zeros((len(refList),len(Mdata)))
    20432112    dFdua = np.zeros((len(refList),len(Mdata),6))
     2113    dFdbab = np.zeros((len(refList),2))
    20442114    for iref,refl in enumerate(refList):
    20452115        H = np.array(refl[:3])
    20462116        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
     2117        SQfactor = 8.0*SQ*np.pi**2
     2118        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
     2119        Bab = parmDict[phfx+'BabA']*dBabdA
    20472120        for i,El in enumerate(Tdata):           
    20482121            FF[i] = refl[-1][El]           
    2049         SQfactor = 8.0*SQ*np.pi**2
    20502122        Uniq = refl[11]
    20512123        phi = refl[12]
    … …  
    20612133        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    20622134        Tcorr = Tiso*Tuij
    2063         fot = (FF+FP)*occ*Tcorr
     2135        fot = (FF+FP-Bab)*occ*Tcorr
    20642136        fotp = FPP*occ*Tcorr
    20652137        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
    … …  
    20752147        dfadui = np.sum(-SQfactor*fa,axis=2)
    20762148        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
     2149        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
    20772150        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
    20782151        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
    … …  
    20802153        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
    20812154        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
     2155        dFdbab[iref] = np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
    20822156        if not SGData['SGInv']:
    20832157            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
    … …  
    20852159            dfbdui = np.sum(-SQfactor*fb,axis=2)
    20862160            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
     2161            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
    20872162            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
    20882163            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
    20892164            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
    20902165            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
     2166            dFdbab[iref] += np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
    20912167        #loop over atoms - each dict entry is list of derivatives for all the reflections
    20922168    for i in range(len(Mdata)):     
    … …  
    21022178        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
    21032179        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     2180        dFdvDict[pfx+'BabA'] = dFdbab.T[0]
     2181        dFdvDict[pfx+'BabU'] = dFdbab.T[1]
    21042182    return dFdvDict
    21052183       
    … …  
    23162394    #crystallite size
    23172395    if calcControls[phfx+'SizeType'] == 'isotropic':
    2318         Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
     2396        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
    23192397    elif calcControls[phfx+'SizeType'] == 'uniaxial':
    23202398        H = np.array(refl[:3])
    23212399        P = np.array(calcControls[phfx+'SizeAxis'])
    23222400        cosP,sinP = G2lat.CosSinAngle(H,P,G)
    2323         Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh)
    2324         Sgam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2)
     2401        Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
     2402        Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
    23252403    else:           #ellipsoidal crystallites
    23262404        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
    … …  
    23302408    #microstrain               
    23312409    if calcControls[phfx+'MustrainType'] == 'isotropic':
    2332         Mgam = 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
     2410        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
    23332411    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
    23342412        H = np.array(refl[:3])
    23352413        P = np.array(calcControls[phfx+'MustrainAxis'])
    23362414        cosP,sinP = G2lat.CosSinAngle(H,P,G)
    2337         Si = parmDict[phfx+'Mustrain:i']
    2338         Sa = parmDict[phfx+'Mustrain:a']
     2415        Si = parmDict[phfx+'Mustrain;i']
     2416        Sa = parmDict[phfx+'Mustrain;a']
    23392417        Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
    23402418    else:       #generalized - P.W. Stephens model
    … …  
    23442422            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
    23452423        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
    2346     gam = Sgam*parmDict[phfx+'Size:mx']+Mgam*parmDict[phfx+'Mustrain:mx']
    2347     sig = (Sgam*(1.-parmDict[phfx+'Size:mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain:mx']))**2
     2424    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
     2425    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
    23482426    sig /= ateln2
    23492427    return sig,gam
    … …  
    23562434    #crystallite size derivatives
    23572435    if calcControls[phfx+'SizeType'] == 'isotropic':
    2358         Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
    2359         gamDict[phfx+'Size:i'] = -900.*wave*parmDict[phfx+'Size:mx']/(np.pi*costh)
    2360         sigDict[phfx+'Size:i'] = -1800.*Sgam*wave*(1.-parmDict[phfx+'Size:mx'])**2/(np.pi*costh*ateln2)
     2436        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
     2437        gamDict[phfx+'Size;i'] = -900.*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
     2438        sigDict[phfx+'Size;i'] = -1800.*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
    23612439    elif calcControls[phfx+'SizeType'] == 'uniaxial':
    23622440        H = np.array(refl[:3])
    23632441        P = np.array(calcControls[phfx+'SizeAxis'])
    23642442        cosP,sinP = G2lat.CosSinAngle(H,P,G)
    2365         Si = parmDict[phfx+'Size:i']
    2366         Sa = parmDict[phfx+'Size:a']
     2443        Si = parmDict[phfx+'Size;i']
     2444        Sa = parmDict[phfx+'Size;a']
    23672445        gami = (1.8*wave/np.pi)/(Si*Sa)
    23682446        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
    … …  
    23712449        dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si)
    23722450        dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa)
    2373         gamDict[phfx+'Size:i'] = dsi*parmDict[phfx+'Size:mx']
    2374         gamDict[phfx+'Size:a'] = dsa*parmDict[phfx+'Size:mx']
    2375         sigDict[phfx+'Size:i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size:mx'])**2/ateln2
    2376         sigDict[phfx+'Size:a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size:mx'])**2/ateln2
     2451        gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
     2452        gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
     2453        sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
     2454        sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
    23772455    else:           #ellipsoidal crystallites
    23782456        const = 1.8*wave/(np.pi*costh)
    … …  
    23822460        Sgam = 1.8*wave/(np.pi*costh*lenR)
    23832461        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
    2384             gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size:mx']
    2385             sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size:mx'])**2/ateln2
    2386     gamDict[phfx+'Size:mx'] = Sgam
    2387     sigDict[phfx+'Size:mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size:mx'])/ateln2
     2462            gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
     2463            sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
     2464    gamDict[phfx+'Size;mx'] = Sgam
     2465    sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
    23882466           
    23892467    #microstrain derivatives               
    23902468    if calcControls[phfx+'MustrainType'] == 'isotropic':
    2391         Mgam = 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
    2392         gamDict[phfx+'Mustrain:i'] =  0.018*tanth*parmDict[phfx+'Mustrain:mx']/np.pi
    2393         sigDict[phfx+'Mustrain:i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain:mx'])**2/(np.pi*ateln2)       
     2469        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
     2470        gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
     2471        sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
    23942472    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
    23952473        H = np.array(refl[:3])
    23962474        P = np.array(calcControls[phfx+'MustrainAxis'])
    23972475        cosP,sinP = G2lat.CosSinAngle(H,P,G)
    2398         Si = parmDict[phfx+'Mustrain:i']
    2399         Sa = parmDict[phfx+'Mustrain:a']
     2476        Si = parmDict[phfx+'Mustrain;i']
     2477        Sa = parmDict[phfx+'Mustrain;a']
    24002478        gami = 0.018*Si*Sa*tanth/np.pi
    24012479        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
    … …  
    24032481        dsi = -gami*Si*cosP**2/sqtrm**3
    24042482        dsa = -gami*Sa*sinP**2/sqtrm**3
    2405         gamDict[phfx+'Mustrain:i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain:mx']
    2406         gamDict[phfx+'Mustrain:a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain:mx']
    2407         sigDict[phfx+'Mustrain:i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain:mx'])**2/ateln2
    2408         sigDict[phfx+'Mustrain:a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain:mx'])**2/ateln2       
     2483        gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
     2484        gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
     2485        sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
     2486        sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
    24092487    else:       #generalized - P.W. Stephens model
    24102488        pwrs = calcControls[phfx+'MuPwrs']
    … …  
    24142492            term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
    24152493            sum += parmDict[phfx+'Mustrain:'+str(i)]*term
    2416             gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain:mx']
     2494            gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx']
    24172495            sigDict[phfx+'Mustrain:'+str(i)] = \
    2418                 2.*const*term*(1.-parmDict[phfx+'Mustrain:mx'])**2/ateln2
     2496                2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
    24192497        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
    24202498        for i in range(len(pwrs)):
    24212499            sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
    2422     gamDict[phfx+'Mustrain:mx'] = Mgam
    2423     sigDict[phfx+'Mustrain:mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain:mx'])/ateln2
     2500    gamDict[phfx+'Mustrain;mx'] = Mgam
     2501    sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
    24242502    return sigDict,gamDict
    24252503       
    … …  
    27062784            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
    27072785        elif SGData['SGLaue'] in ['4/m','4/mmm']:
    2708 #            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
    27092786            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
    27102787        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
    2711 #            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
    27122788            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
    27132789        elif SGData['SGLaue'] in ['3R', '3mR']:
    27142790            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
    27152791        elif SGData['SGLaue'] in ['m3m','m3']:
    2716 #            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
    27172792            return [[pfx+'A0',dpdA[0]]]
    27182793    # create a list of dependent variables and set up a dictionary to hold their derivatives
    … …  
    27892864                dMdpk = np.zeros(shape=(6,lenBF))
    27902865                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
    2791                 for i in range(1,5):
     2866                for i in range(5):
    27922867                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[13]*refl[9]*dMdipk[i]
    2793                 dMdpk[0] += 100.*cw[iBeg:iFin]*refl[13]*refl[9]*dMdipk[0]
    27942868                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
    27952869                if Ka2:
    … …  
    28012875                        dMdpk2 = np.zeros(shape=(6,lenBF2))
    28022876                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
    2803                         for i in range(1,5):
     2877                        for i in range(5):
    28042878                            dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[13]*refl[9]*kRatio*dMdipk2[i]
    2805                         dMdpk2[0] = 100.*cw[iBeg2:iFin2]*refl[13]*refl[9]*kRatio*dMdipk2[0]
    28062879                        dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[13]*dMdipk2[0]
    28072880                        dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
    … …  
    29012974                        if Ka2:
    29022975                            depDerivDict[name][iBeg2:iFin2] += sigDict[name]*dervDict2['sig']
    2903                                                
     2976                for name in ['BabA','BabU']:
     2977                    if phfx+name in varylist:
     2978                        dMdv[varylist.index(phfx+name)][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
     2979                        if Ka2:
     2980                            dMdv[varylist.index(phfx+name)][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]
     2981                    elif phfx+name in dependentVars:                   
     2982                        depDerivDict[phfx+name][iBeg:iFin] += dFdvDict[pfx+name][iref]*dervDict['int']*cw[iBeg:iFin]
     2983                        if Ka2:
     2984                            depDerivDict[phfx+name][iBeg2:iFin2] += dFdvDict[pfx+name][iref]*dervDict2['int']*cw[iBeg2:iFin2]                 
    29042985            elif 'T' in calcControls[hfx+'histType']:
    29052986                print 'TOF Undefined at present'

  • trunk/imports/G2pwd_fxye.py

    r803 r808  
    3232            if i==0: # first line is always a comment
    3333                continue
    34             if i==1 and line[:4].lower() == 'inst' and ':' in line:
     34#            if i==1 and line[:4].lower() == 'inst' and ':' in line:
     35            if i==1 and line[:4].lower() == 'inst':
    3536                # 2nd line is optional instrument parameter file
    3637                continue
Note: See TracChangeset for help on using the changeset viewer.