Skip to content

Commit f243fb3

Browse files
vondreelevondreele
committed
fix printing of texture & spin RB penalty function results
1 parent 40a18c2 commit f243fb3

2 files changed

Lines changed: 148 additions & 125 deletions

File tree

‎GSASII/GSASIIstrIO.py‎

Lines changed: 147 additions & 124 deletions
Original file line numberDiff line numberDiff line change
@@ -1936,9 +1936,6 @@ def MakeRBSphHarm(rbKey,phaseVary,phaseDict):
19361936
(Vec[0],Vec[1],Vec[2],maxH,vRef))
19371937
if seqHistName is not None:
19381938
PrintTexture(textureData)
1939-
if name in RestraintDict:
1940-
PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1941-
textureData,RestraintDict[name],pFile)
19421939

19431940
elif PawleyRef:
19441941
if Print:
@@ -2045,124 +2042,6 @@ def cellFill(pfx,SGData,parmDict,sigDict):
20452042
sigA = [0,0,0,0,0,0]
20462043
return A,sigA
20472044

2048-
def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
2049-
'''Documents Restraint settings in .lst file
2050-
2051-
TODO: pass in parmDict to evaluate general restraints
2052-
'''
2053-
if phaseRest:
2054-
Amat = G2lat.cell2AB(cell)[0]
2055-
cx,ct,cs = AtPtrs[:3]
2056-
names = G2obj.restraintNames
2057-
for name,rest in names:
2058-
if name not in phaseRest:
2059-
continue
2060-
itemRest = phaseRest[name]
2061-
if rest in itemRest and itemRest[rest] and itemRest['Use']:
2062-
pFile.write('\n %s restraint weight factor %10.3f Use: %s\n'%(name,itemRest['wtFactor'],str(itemRest['Use'])))
2063-
if name in ['Bond','Angle','Plane','Chiral']:
2064-
pFile.write(' calc obs sig delt/sig atoms(symOp): \n')
2065-
for indx,ops,obs,esd in itemRest[rest]:
2066-
try:
2067-
AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
2068-
AtName = ''
2069-
for i,Aname in enumerate(AtNames):
2070-
AtName += Aname
2071-
if ops[i] == '1':
2072-
AtName += '-'
2073-
else:
2074-
AtName += '+('+ops[i]+')-'
2075-
XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
2076-
XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2077-
if name == 'Bond':
2078-
calc = G2mth.getRestDist(XYZ,Amat)
2079-
elif name == 'Angle':
2080-
calc = G2mth.getRestAngle(XYZ,Amat)
2081-
elif name == 'Plane':
2082-
calc = G2mth.getRestPlane(XYZ,Amat)
2083-
elif name == 'Chiral':
2084-
calc = G2mth.getRestChiral(XYZ,Amat)
2085-
pFile.write(' %9.3f %9.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1]))
2086-
except KeyError:
2087-
del itemRest[rest]
2088-
elif name in ['Torsion','Rama']:
2089-
pFile.write(' atoms(symOp) calc obs sig delt/sig torsions: \n')
2090-
coeffDict = itemRest['Coeff']
2091-
for indx,ops,cofName,esd in itemRest[rest]:
2092-
AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
2093-
AtName = ''
2094-
for i,Aname in enumerate(AtNames):
2095-
AtName += Aname+'+('+ops[i]+')-'
2096-
XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
2097-
XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2098-
if name == 'Torsion':
2099-
tor = G2mth.getRestTorsion(XYZ,Amat)
2100-
restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
2101-
pFile.write(' %8.3f %8.3f %.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1]))
2102-
else:
2103-
phi,psi = G2mth.getRestRama(XYZ,Amat)
2104-
restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
2105-
pFile.write(' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,phi,psi,AtName[:-1]))
2106-
elif name == 'ChemComp':
2107-
pFile.write(' atoms mul*frac factor prod\n')
2108-
for indx,factors,obs,esd in itemRest[rest]:
2109-
try:
2110-
atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
2111-
mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
2112-
frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
2113-
mulfrac = mul*frac
2114-
calcs = mul*frac*factors
2115-
for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
2116-
pFile.write(' %10s %8.3f %8.3f %8.3f\n'%(atom,mf,fr,clc))
2117-
pFile.write(' Sum: calc: %8.3f obs: %8.3f esd: %8.3f\n'%(np.sum(calcs),obs,esd))
2118-
except KeyError:
2119-
del itemRest[rest]
2120-
elif name == 'Moments':
2121-
for item in itemRest['Moments']:
2122-
nMom = len(item[0])
2123-
AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,item[0],ct-1)
2124-
moms = G2mth.GetAtomItemsById(Atoms,AtLookup,item[0],cx+4,3)
2125-
obs = 0.
2126-
pFile.write(nMom*' Atom calc'+' obs esd\n')
2127-
line = ''
2128-
for i,mom in enumerate(moms):
2129-
Mom = G2mth.GetMag(mom,cell)
2130-
line += ' %s %8.3f'%(AtNames[i],Mom)
2131-
obs += Mom
2132-
obs /= nMom
2133-
line += '%8.3f %6.3f\n'%(obs,item[-1])
2134-
pFile.write(line)
2135-
elif name == 'Texture' and textureData['Order']:
2136-
SHCoef = textureData['SH Coeff'][1]
2137-
shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2138-
SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2139-
pFile.write (' HKL grid neg esd sum neg num neg use unit? unit esd \n')
2140-
for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
2141-
phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
2142-
ODFln = G2lat.Flnh(SHCoef,phi,beta,SGData)
2143-
R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
2144-
Z = ma.masked_greater(Z,0.0)
2145-
num = ma.count(Z)
2146-
sum = 0
2147-
if num:
2148-
sum = np.sum(Z)
2149-
pFile.write (' %d %d %d %d %8.3f %8.3f %8d %s %8.3f\n'%(hkl[0],hkl[1],hkl[2],grid,esd1,sum,num,str(ifesd2),esd2))
2150-
elif name =='SpinRB':
2151-
continue
2152-
elif name == 'General':
2153-
pFile.write(' target sig obs expression variables \n')
2154-
for expObj,target,esd in itemRest[rest]:
2155-
val = '?'
2156-
#calcobj = G2obj.ExpressionCalcObj(expObj)
2157-
# need to get parmDict to evaluate the value
2158-
#calcobj.SetupCalc(parmDict)
2159-
#val = ' {:8.3g} '.format(calcobj.EvalExpression())
2160-
txt = ''
2161-
for i in expObj.assgnVars:
2162-
if txt: txt += '; '
2163-
txt += str(i) + '=' + str(expObj.assgnVars[i])
2164-
if len(txt) > 80: txt = txt[:77]+'...'
2165-
pFile.write(f'{target:8.5g}{esd:6.3g}{val:>8s} {expObj.expression:17s} {txt}\n')
21662045

21672046
def SummRestraints(restraintDict):
21682047
'Summarize number of Restraints in a single line'
@@ -2453,8 +2332,150 @@ def PrintSHtextureAndSig(textureData,SHtextureSig):
24532332
iFin = min(iBeg+10,nCoeff)
24542333
pFile.write(' Texture index J = %.3f(%d)'%(Tindx,int(1000*np.sqrt(Tvar))))
24552334

2335+
def PrintRestraints(cell,AtPtrs,phaseRest):
2336+
'''Documents Restraint settings in .lst file
2337+
2338+
TODO: pass in parmDict to evaluate general restraints
2339+
'''
2340+
if phaseRest:
2341+
Amat = G2lat.cell2AB(cell)[0]
2342+
cx,ct,cs = AtPtrs[:3]
2343+
names = G2obj.restraintNames
2344+
for name,rest in names:
2345+
if name not in phaseRest:
2346+
continue
2347+
itemRest = phaseRest[name]
2348+
if rest in itemRest and itemRest[rest] and itemRest['Use']:
2349+
pFile.write('\n %s restraint weight factor %10.3f Use: %s\n'%(name,itemRest['wtFactor'],str(itemRest['Use'])))
2350+
if name in ['Bond','Angle','Plane','Chiral']:
2351+
pFile.write(' calc obs sig delt/sig atoms(symOp): \n')
2352+
for indx,ops,obs,esd in itemRest[rest]:
2353+
try:
2354+
AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
2355+
AtName = ''
2356+
for i,Aname in enumerate(AtNames):
2357+
AtName += Aname
2358+
if ops[i] == '1':
2359+
AtName += '-'
2360+
else:
2361+
AtName += '+('+ops[i]+')-'
2362+
XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
2363+
XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2364+
if name == 'Bond':
2365+
calc = G2mth.getRestDist(XYZ,Amat)
2366+
elif name == 'Angle':
2367+
calc = G2mth.getRestAngle(XYZ,Amat)
2368+
elif name == 'Plane':
2369+
calc = G2mth.getRestPlane(XYZ,Amat)
2370+
elif name == 'Chiral':
2371+
calc = G2mth.getRestChiral(XYZ,Amat)
2372+
pFile.write(' %9.3f %9.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1]))
2373+
except KeyError:
2374+
del itemRest[rest]
2375+
elif name in ['Torsion','Rama']:
2376+
pFile.write(' atoms(symOp) calc obs sig delt/sig torsions: \n')
2377+
coeffDict = itemRest['Coeff']
2378+
for indx,ops,cofName,esd in itemRest[rest]:
2379+
AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
2380+
AtName = ''
2381+
for i,Aname in enumerate(AtNames):
2382+
AtName += Aname+'+('+ops[i]+')-'
2383+
XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
2384+
XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2385+
if name == 'Torsion':
2386+
tor = G2mth.getRestTorsion(XYZ,Amat)
2387+
restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
2388+
pFile.write(' %8.3f %8.3f %.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1]))
2389+
else:
2390+
phi,psi = G2mth.getRestRama(XYZ,Amat)
2391+
restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])
2392+
pFile.write(' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,phi,psi,AtName[:-1]))
2393+
elif name == 'ChemComp':
2394+
pFile.write(' atoms mul*frac factor prod\n')
2395+
for indx,factors,obs,esd in itemRest[rest]:
2396+
try:
2397+
atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
2398+
mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
2399+
frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
2400+
mulfrac = mul*frac
2401+
calcs = mul*frac*factors
2402+
for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
2403+
pFile.write(' %10s %8.3f %8.3f %8.3f\n'%(atom,mf,fr,clc))
2404+
pFile.write(' Sum: calc: %8.3f obs: %8.3f esd: %8.3f\n'%(np.sum(calcs),obs,esd))
2405+
except KeyError:
2406+
del itemRest[rest]
2407+
elif name == 'Moments':
2408+
for item in itemRest['Moments']:
2409+
nMom = len(item[0])
2410+
AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,item[0],ct-1)
2411+
moms = G2mth.GetAtomItemsById(Atoms,AtLookup,item[0],cx+4,3)
2412+
obs = 0.
2413+
pFile.write(nMom*' Atom calc'+' obs esd\n')
2414+
line = ''
2415+
for i,mom in enumerate(moms):
2416+
Mom = G2mth.GetMag(mom,cell)
2417+
line += ' %s %8.3f'%(AtNames[i],Mom)
2418+
obs += Mom
2419+
obs /= nMom
2420+
line += '%8.3f %6.3f\n'%(obs,item[-1])
2421+
pFile.write(line)
2422+
elif name == 'Texture' and textureData['Order']:
2423+
SHCoef = textureData['SH Coeff'][1]
2424+
shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2425+
SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2426+
pFile.write (' HKL grid neg esd sum neg num neg use unit? unit esd \n')
2427+
for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
2428+
phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
2429+
ODFln = G2lat.Flnh(SHCoef,phi,beta,SGData)
2430+
R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
2431+
Z = ma.masked_greater(Z,0.0)
2432+
num = ma.count(Z)
2433+
sum = 0
2434+
if num:
2435+
sum = np.sum(Z)
2436+
pFile.write (' %d %d %d %d %8.3f %8.3f %8d %s %8.3f\n'%(hkl[0],hkl[1],hkl[2],grid,esd1,sum,num,str(ifesd2),esd2))
2437+
elif name =='SpinRB':
2438+
RBData = Phases[phase]['RBModels']
2439+
SpinRB = RBData.get('Spin',[])
2440+
Npts = 20
2441+
PSI,GAM = np.mgrid[0:2*Npts,0:Npts] #[azm,pol]
2442+
PSI = PSI.flatten()*360./(2*Npts) #azimuth 0-360 ncl
2443+
GAM = GAM.flatten()*180./Npts #polar 0-180 incl
2444+
pFile.write (' atom shell neg esd sum neg num neg \n')
2445+
for ires,[iAt,esd,ifUnit,esdU] in enumerate(itemRest[rest]):
2446+
for irb,RBObj in enumerate(SpinRB):
2447+
if iAt == RBObj['Ids'][0]:
2448+
atom = Atoms[AtLookup[iAt]]
2449+
for ish in range(len(RBObj['RBId'])): #do this by shell
2450+
SH = np.zeros(2*Npts*Npts)
2451+
for item in RBObj['SHC'][ish]:
2452+
iRb = SRBIds.index(RBObj['RBId'][ish])
2453+
pName = '%d::RBSSh;%d;%s:%d:%d'%(pId,ish,item,AtLookup[iAt],iRb)
2454+
SH += np.array(G2lat.KslCalc(item,PSI,GAM))*parmDict[pName]
2455+
SH = ma.masked_greater(SH,0.0)
2456+
num = ma.count(SH)
2457+
sum = 0
2458+
if num:
2459+
sum = np.sum(SH)
2460+
pFile.write (' %s %d %8.3f %8.3f %6d\n'%(atom[ct-1],ish,esd,sum,num))
2461+
2462+
elif name == 'General':
2463+
pFile.write(' target sig obs expression variables \n')
2464+
for expObj,target,esd in itemRest[rest]:
2465+
val = '?'
2466+
#calcobj = G2obj.ExpressionCalcObj(expObj)
2467+
# need to get parmDict to evaluate the value
2468+
#calcobj.SetupCalc(parmDict)
2469+
#val = ' {:8.3g} '.format(calcobj.EvalExpression())
2470+
txt = ''
2471+
for i in expObj.assgnVars:
2472+
if txt: txt += '; '
2473+
txt += str(i) + '=' + str(expObj.assgnVars[i])
2474+
if len(txt) > 80: txt = txt[:77]+'...'
2475+
pFile.write(f'{target:8.5g}{esd:6.3g}{val:>8s} {expObj.expression:17s} {txt}\n')
2476+
24562477
##########################################################################
2457-
# SetPhaseData starts here
2478+
#### SetPhaseData starts here
24582479
if pFile: pFile.write('\n Phases:\n')
24592480
for phase in Phases:
24602481
if pFile: pFile.write(' Result for phase: %s\n'%phase)
@@ -2640,6 +2661,9 @@ def PrintSHtextureAndSig(textureData,SHtextureSig):
26402661
if pfx+name in sigDict:
26412662
wavesSig[name] = sigDict[pfx+name]
26422663

2664+
2665+
2666+
26432667
Deformations = Phase.get('Deformations',{})
26442668
for iAt in Deformations:
26452669
if iAt < 0:
@@ -2685,8 +2709,7 @@ def PrintSHtextureAndSig(textureData,SHtextureSig):
26852709
SHtextureSig[name] = sigDict[aname]
26862710
PrintSHtextureAndSig(textureData,SHtextureSig)
26872711
if phase in RestraintDict and not Phase['General'].get('doPawley'):
2688-
PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2689-
textureData,RestraintDict[phase],pFile)
2712+
PrintRestraints(cell[1:7],General['AtomPtrs'],RestraintDict[phase])
26902713

26912714
def SetISOmodes(parmDict,sigDict,Phases,pFile=None):
26922715
'''After a refinement, sets the values for the ISODISTORT modes into

‎GSASII/GSASIIstrMath.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -788,10 +788,10 @@ def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
788788
pVals.append(sh)
789789
pWt.append(wt/esd**2)
790790
pWsum[name] += wt*(-sh/esd)**2
791-
pWnum[name] += 1
792791
else:
793792
pVals.append(0.0)
794793
pWt.append(0.0)
794+
pWnum[name] += ma.count(ma.masked_greater(SH,0.0))
795795

796796
elif name == 'Texture':
797797
SHkeys = list(textureData['SH Coeff'][1].keys())

0 commit comments

Comments
 (0)