Changeset 808


Ignore:
Timestamp:
Dec 5, 2012 10:06:52 AM (9 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.