Changeset 5078


Ignore:
Timestamp:
Nov 13, 2021 11:41:21 AM (7 months ago)
Author:
toby
Message:

start on mmCIF export

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/exports/G2export_CIF.py

    r5052 r5078  
    731731        WriteCIFitem(fp, s)
    732732
     733def WriteAtomsMM(fp, phasedict, phasenam, parmDict, sigDict,
     734                          RBparms={}):
     735    'Write atom positions to CIF using mmCIF items'
     736    # phasedict = self.Phases[phasenam] # pointer to current phase info
     737    General = phasedict['General']
     738    cx,ct,cs,cia = General['AtomPtrs']
     739    GS = G2lat.cell2GS(General['Cell'][1:7])
     740    Amat = G2lat.cell2AB(General['Cell'][1:7])[0]
     741    Atoms = phasedict['Atoms']
     742    cfrac = cx+3
     743    fpfx = str(phasedict['pId'])+'::Afrac:'
     744    for i,at in enumerate(Atoms):
     745        fval = parmDict.get(fpfx+str(i),at[cfrac])
     746        if fval != 0.0:
     747            break
     748    else:
     749        WriteCIFitem(fp, '\n# PHASE HAS NO ATOMS!')
     750        return
     751
     752    WriteCIFitem(fp, '\n# ATOMIC COORDINATES AND DISPLACEMENT PARAMETERS')
     753    WriteCIFitem(fp, 'loop_ '+
     754                 '\n   _atom_site.id'+
     755                 '\n   _atom_site.type_symbol'+
     756                 '\n   _atom_site.label_atom_id'+
     757                 '\n   _atom_site.label_alt_id'+
     758                 '\n   _atom_site.label_comp_id'+
     759                 '\n   _atom_site.label_asym_id'+
     760                 '\n   _atom_site.label_entity_id'+
     761                 '\n   _atom_site.label_seq_id'+
     762                 '\n   _atom_site.pdbx_PDB_ins_code'+
     763                 '\n   _atom_site.fract_x'+
     764                 '\n   _atom_site.fract_y'+
     765                 '\n   _atom_site.fract_z'+
     766                 '\n   _atom_site.occupancy'+
     767#                 '\n   _atom_site_adp_type'+
     768                 '\n   _atom_site.U_iso_or_equiv'
     769#                 '\n   _atom_site_site_symmetry_multiplicity'
     770                 )
     771# _atom_site.group_PDB
     772# _atom_site.Cartn_x
     773# _atom_site.Cartn_y
     774# _atom_site.Cartn_z
     775# _atom_site.Cartn_x_esd
     776# _atom_site.Cartn_y_esd
     777# _atom_site.Cartn_z_esd
     778# _atom_site.occupancy_esd
     779# _atom_site.B_iso_or_equiv_esd
     780# _atom_site.pdbx_formal_charge
     781# _atom_site.auth_seq_id
     782# _atom_site.auth_comp_id
     783# _atom_site.auth_asym_id
     784# _atom_site.auth_atom_id
     785# _atom_site.pdbx_PDB_model_num
     786    varnames = {cx:'Ax',cx+1:'Ay',cx+2:'Az',cx+3:'Afrac',
     787                cia+1:'AUiso',cia+2:'AU11',cia+3:'AU22',cia+4:'AU33',
     788                cia+5:'AU12',cia+6:'AU13',cia+7:'AU23'}
     789
     790    pfx = str(phasedict['pId'])+'::'
     791    # loop over all atoms
     792    naniso = 0
     793    for i,at in enumerate(Atoms):
     794        s = ''
     795        s += PutInCol(str(i),5) # atom number
     796        s += PutInCol(FmtAtomType(at[ct]),4) # type
     797        s += PutInCol(at[ct-1],4) # label_atom_id
     798        s += PutInCol('.',2) # alt_id
     799        s += PutInCol(at[ct-3],4) # comp_id
     800        s += PutInCol(at[ct-2],3) # _asym_id
     801        s += PutInCol(at[ct-4],3) # entity_id
     802        s += PutInCol('?',2) # _seq_id
     803        s += PutInCol('?',2) # pdbx_PDB_ins_code
     804        fval = parmDict.get(fpfx+str(i),at[cfrac])
     805        if fval == 0.0: continue # ignore any atoms that have a occupancy set to 0 (exact)
     806        if at[cia] == 'I':
     807            adp = 'Uiso '
     808        else:
     809            adp = 'Uani '
     810            naniso += 1
     811            t = G2lat.Uij2Ueqv(at[cia+2:cia+8],GS,Amat)[0]
     812            for j in (2,3,4):
     813                var = pfx+varnames[cia+j]+":"+str(i)
     814        for j in (cx,cx+1,cx+2,cx+3):
     815            if j in (cx,cx+1,cx+2):
     816                dig = 11
     817                sigdig = -0.00009
     818            else:
     819                dig = 5
     820                sigdig = -0.09
     821            var = pfx+varnames[j]+":"+str(i)
     822            dvar = pfx+"d"+varnames[j]+":"+str(i)
     823            if dvar not in sigDict:
     824                dvar = var
     825            #print var,(var in parmDict),(var in sigDict)
     826            val = parmDict.get(var,at[j])
     827            sig = sigDict.get(dvar,sigdig)
     828            if dvar in G2mv.GetDependentVars(): # do not include an esd for dependent vars
     829                sig = -abs(sig)
     830            s += PutInCol(G2mth.ValEsd(val,sig),dig)
     831        s += PutInCol(at[cs+1],3)
     832        WriteCIFitem(fp, s)
     833    # save information about rigid bodies
     834#     header = False
     835#     num = 0
     836#     rbAtoms = []
     837#     for irb,RBObj in enumerate(phasedict['RBModels'].get('Residue',[])):
     838#         if not header:
     839#             header = True
     840#             RBheader(fp)
     841#         jrb = RBparms['RBIds']['Residue'].index(RBObj['RBId'])
     842#         rbsx = str(irb)+':'+str(jrb)
     843#         num += 1
     844#         WriteCIFitem(fp,'',str(num))
     845#         RBModel = RBparms['Residue'][RBObj['RBId']]
     846#         SGData = phasedict['General']['SGData']
     847#         Sytsym,Mult = G2spc.SytSym(RBObj['Orig'][0],SGData)[:2]
     848#         s = '''GSAS-II residue rigid body "{}" with {} atoms
     849#   Site symmetry @ origin: {}, multiplicity: {}
     850# '''.format(RBObj['RBname'],len(RBModel['rbTypes']),Sytsym,Mult)
     851#         for i in G2stIO.WriteResRBModel(RBModel):
     852#             s += i
     853#         s += '\n Location:\n'
     854#         for i in G2stIO.WriteRBObjPOAndSig(pfx,'RBR',rbsx,parmDict,sigDict):
     855#             s += i+'\n'
     856#         for i in G2stIO.WriteRBObjTLSAndSig(pfx,'RBR',rbsx,
     857#                         RBObj['ThermalMotion'][0],parmDict,sigDict):
     858#             s += i
     859#         nTors = len(RBObj['Torsions'])
     860#         if nTors:
     861#             for i in G2stIO.WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,
     862#                         nTors):
     863#                 s += i
     864#         WriteCIFitem(fp,'',s.rstrip())
     865       
     866#         pId = phasedict['pId']
     867#         for i in RBObj['Ids']:
     868#             lbl = G2obj.LookupAtomLabel(pId,G2obj.LookupAtomId(pId,i))[0]
     869#             rbAtoms.append('{:7s} 1_555 {:3d} ?'.format(lbl,num))
     870#         #GSASIIpath.IPyBreak()
     871
     872#     for irb,RBObj in enumerate(phasedict['RBModels'].get('Vector',[])):
     873#         if not header:
     874#             header = True
     875#             RBheader(fp)
     876#         jrb = RBparms['RBIds']['Vector'].index(RBObj['RBId'])
     877#         rbsx = str(irb)+':'+str(jrb)
     878#         num += 1
     879#         WriteCIFitem(fp,'',str(num))
     880#         RBModel = RBparms['Vector'][RBObj['RBId']]
     881#         SGData = phasedict['General']['SGData']
     882#         Sytsym,Mult = G2spc.SytSym(RBObj['Orig'][0],SGData)[:2]
     883#         s = '''GSAS-II vector rigid body "{}" with {} atoms
     884#   Site symmetry @ origin: {}, multiplicity: {}
     885# '''.format(RBObj['RBname'],len(RBModel['rbTypes']),Sytsym,Mult)
     886#         for i in G2stIO.WriteVecRBModel(RBModel,sigDict,irb):
     887#             s += i
     888#         s += '\n Location:\n'
     889#         for i in G2stIO.WriteRBObjPOAndSig(pfx,'RBV',rbsx,parmDict,sigDict):
     890#             s += i+'\n'
     891#         for i in G2stIO.WriteRBObjTLSAndSig(pfx,'RBV',rbsx,
     892#                         RBObj['ThermalMotion'][0],parmDict,sigDict):
     893#             s += i
     894#         WriteCIFitem(fp,'',s.rstrip())
     895       
     896#         pId = phasedict['pId']
     897#         for i in RBObj['Ids']:
     898#             lbl = G2obj.LookupAtomLabel(pId,G2obj.LookupAtomId(pId,i))[0]
     899#             rbAtoms.append('{:7s} 1_555 {:3d} ?'.format(lbl,num))
     900
     901#     if rbAtoms:
     902#         WriteCIFitem(fp,'loop_\n    _restr_rigid_body.id'+
     903#             '\n    _restr_rigid_body.atom_site_label\n    _restr_rigid_body.site_symmetry'+
     904#             '\n    _restr_rigid_body.class_id\n    _restr_rigid_body.details')
     905#         for i,l in enumerate(rbAtoms):
     906#             WriteCIFitem(fp,'   {:5d} {}'.format(i+1,l))
     907
    733908# Refactored over here to allow access by GSASIIscriptable.py
    734909def WriteSeqAtomsNuclear(fp, cell, phasedict, phasenam, hist, seqData, RBparms):
     
    11081283                  G2mth.ValEsd(cellmass/Z,-0.09,True))
    11091284
     1285def WriteCompositionMM(fp, phasedict, phasenam, parmDict, quickmode=True, keV=None):
     1286    '''determine the composition for the unit cell, crudely determine Z and
     1287    then compute the composition in formula units.
     1288
     1289    If quickmode is False, then scattering factors are added to the element loop.
     1290
     1291    If keV is specified, then resonant scattering factors are also computed and included.
     1292    '''
     1293    General = phasedict['General']
     1294    Z = General.get('cellZ',0.0)
     1295    cx,ct,cs,cia = General['AtomPtrs']
     1296    Atoms = phasedict['Atoms']
     1297    fpfx = str(phasedict['pId'])+'::Afrac:'
     1298    cfrac = cx+3
     1299    cmult = cs+1
     1300    compDict = {} # combines H,D & T
     1301    sitemultlist = []
     1302    massDict = dict(zip(General['AtomTypes'],General['AtomMass']))
     1303    cellmass = 0
     1304    elmLookup = {}
     1305    for i,at in enumerate(Atoms):
     1306        atype = at[ct].strip()
     1307        if atype.find('-') != -1: atype = atype.split('-')[0]
     1308        if atype.find('+') != -1: atype = atype.split('+')[0]
     1309        atype = atype[0].upper()+atype[1:2].lower() # force case conversion
     1310        if atype == "D" or atype == "D": atype = "H"
     1311        fvar = fpfx+str(i)
     1312        fval = parmDict.get(fvar,at[cfrac])
     1313        mult = at[cmult]
     1314        if not massDict.get(at[ct]):
     1315            print('Error: No mass found for atom type '+at[ct])
     1316            print('Will not compute cell contents for phase '+phasenam)
     1317            return
     1318        cellmass += massDict[at[ct]]*mult*fval
     1319        compDict[atype] = compDict.get(atype,0.0) + mult*fval
     1320        elmLookup[atype] = at[ct].strip()
     1321        if fval == 1: sitemultlist.append(mult)
     1322    if len(compDict.keys()) == 0: return # no elements!
     1323    if Z < 1: # Z has not been computed or set by user
     1324        Z = 1
     1325        if not sitemultlist:
     1326            General['cellZ'] = 1
     1327            return
     1328        for i in range(2,min(sitemultlist)+1):
     1329            for m in sitemultlist:
     1330                if m % i != 0:
     1331                    break
     1332                else:
     1333                    Z = i
     1334        General['cellZ'] = Z # save it
     1335
     1336    if not quickmode:
     1337        FFtable = G2el.GetFFtable(General['AtomTypes'])
     1338        BLtable = G2el.GetBLtable(General)
     1339
     1340    WriteCIFitem(fp, '\nloop_  _atom_site.type_symbol _atom_type.number_in_cell')
     1341    s = '       '
     1342    if not quickmode:
     1343        for j in ('a1','a2','a3','a4','b1','b2','b3','b4','c',2,1):
     1344            if len(s) > 80:
     1345                WriteCIFitem(fp, s)
     1346                s = '       '
     1347            if j==1:
     1348                s += ' _atom_type.scat_source'
     1349            elif j==2:
     1350                s += ' _atom_type.scat_length_neutron'
     1351            else:
     1352                s += ' _atom_type.scat_Cromer_Mann_'
     1353                s += j
     1354        if keV:
     1355            WriteCIFitem(fp, s)
     1356            s = '        _atom_type.scat_dispersion_real _atom_type.scat_dispersion_imag _atom_type_scat_dispersion_source'
     1357        WriteCIFitem(fp, s)
     1358
     1359           
     1360    formula = ''
     1361    for elem in HillSortElements(list(compDict.keys())):
     1362        s = '  '
     1363        elmsym = elmLookup[elem]
     1364        # CIF does not allow underscore in element symbol (https://www.iucr.org/__data/iucr/cifdic_html/1/cif_core.dic/Iatom_type_symbol.html)
     1365        if elmsym.endswith("_"):
     1366             s += PutInCol(elmsym.replace('_',''))
     1367        elif '_' in elmsym:
     1368             s += PutInCol(elmsym.replace('_','~'))
     1369        else:
     1370            s += PutInCol(elmsym,7)
     1371        s += PutInCol(G2mth.ValEsd(compDict[elem],-0.009,True),5)
     1372        if not quickmode:
     1373            for i in 'fa','fb','fc':
     1374                if i != 'fc':
     1375                    for j in range(4):
     1376                        if elmsym in FFtable:
     1377                            val = G2mth.ValEsd(FFtable[elmsym][i][j],-0.0009,True)
     1378                        else:
     1379                            val = '?'
     1380                        s += ' '
     1381                        s += PutInCol(val,9)
     1382                else:
     1383                    if elmsym in FFtable:
     1384                        val = G2mth.ValEsd(FFtable[elmsym][i],-0.0009,True)
     1385                    else:
     1386                        val = '?'
     1387                    s += ' '
     1388                    s += PutInCol(val,9)
     1389            if elmsym in BLtable:
     1390                bldata = BLtable[elmsym]
     1391                #isotope = bldata[0]
     1392                #mass = bldata[1]['Mass']
     1393                if 'BW-LS' in bldata[1]:
     1394                    val = 0
     1395                else:
     1396                    val = G2mth.ValEsd(bldata[1]['SL'][0],-0.0009,True)
     1397            else:
     1398                val = '?'
     1399            s += ' '
     1400            s += PutInCol(val,9)
     1401            WriteCIFitem(fp,s.rstrip())
     1402            WriteCIFitem(fp,'      https://subversion.xray.aps.anl.gov/pyGSAS/trunk/atmdata.py')
     1403            if keV:
     1404                Orbs = G2el.GetXsectionCoeff(elem.split('+')[0].split('-')[0])
     1405                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
     1406                WriteCIFitem(fp,'  {:8.3f}{:8.3f}   https://subversion.xray.aps.anl.gov/pyGSAS/trunk/atmdata.py'.format(FP,FPP))
     1407        else:
     1408            WriteCIFitem(fp,s.rstrip())
     1409        if formula: formula += " "
     1410        formula += elem
     1411        if compDict[elem] == Z: continue
     1412        formula += G2mth.ValEsd(compDict[elem]/Z,-0.009,True)
     1413    WriteCIFitem(fp,  '\n# Note that Z affects _cell_formula.sum and .weight')
     1414    WriteCIFitem(fp,  '_cell.formula_units_Z',str(Z))
     1415    WriteCIFitem(fp,  '_chemical_formula.sum',formula)
     1416    WriteCIFitem(fp,  '_chemical_formula.weight',
     1417                  G2mth.ValEsd(cellmass/Z,-0.09,True))
     1418
    11101419class ExportCIF(G2IO.ExportBaseclass):
    11111420    '''Base class for CIF exports
     
    22142523            # report cell contents                       
    22152524            WriteComposition(self.fp, self.Phases[phasenam], phasenam, self.parmDict, self.quickmode, keV)
     2525            if not self.quickmode and phasedict['General']['Type'] == 'nuclear':      # report distances and angles
     2526                WriteDistances(phasenam)
     2527            if 'Map' in phasedict['General'] and 'minmax' in phasedict['General']['Map']:
     2528                WriteCIFitem(self.fp, '\n# Difference density results')
     2529                MinMax = phasedict['General']['Map']['minmax']
     2530                WriteCIFitem(self.fp, '_refine_diff_density_max',G2mth.ValEsd(MinMax[0],-0.009))
     2531                WriteCIFitem(self.fp, '_refine_diff_density_min',G2mth.ValEsd(MinMax[1],-0.009))
     2532
     2533        def WritePhaseInfoMM(phasenam,quick=True,oneblock=True):
     2534            'Write out the phase information for the selected phase for a macromolecular phase'
     2535            WriteCIFitem(self.fp, '\n# phase info for '+str(phasenam) + ' follows')
     2536            phasedict = self.Phases[phasenam] # pointer to current phase info
     2537            WriteCIFitem(self.fp, '_pd_phase_name', phasenam)
     2538            cellList,cellSig = self.GetCell(phasenam)
     2539            if oneblock:
     2540                pass # temperature should be written when the histogram saved later
     2541            else: # get T set in _SelectPhaseT_CellSelectHist and possibly get new cell params
     2542                T,hRanId = self.CellHistSelection.get(phasedict['ranId'],
     2543                                                          ('?',None))
     2544                try:
     2545                    T = G2mth.ValEsd(T,-1.0)
     2546                except:
     2547                    pass
     2548                WriteCIFitem(self.fp,"_cell_measurement.temp",T)
     2549                for h in self.Histograms:
     2550                    if self.Histograms[h]['ranId'] == hRanId:
     2551                        pId = phasedict['pId']
     2552                        hId = self.Histograms[h]['hId']
     2553                        cellList,cellSig = G2stIO.getCellSU(pId,hId,
     2554                                        phasedict['General']['SGData'],
     2555                                        self.parmDict,
     2556                                        self.OverallParms['Covariance'])
     2557                        break
     2558                else:
     2559                    T = '?'
     2560
     2561            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
     2562            prevsig = 0
     2563            for lbl,defsig,val,sig in zip(cellNames,defsigL,cellList,cellSig):
     2564                if sig:
     2565                    txt = G2mth.ValEsd(val,sig)
     2566                    prevsig = -sig # use this as the significance for next value
     2567                else:
     2568                    txt = G2mth.ValEsd(val,min(defsig,prevsig),True)
     2569                WriteCIFitem(self.fp, '_cell.'+lbl,txt)
     2570               
     2571            density = G2mth.getDensity(phasedict['General'])[0]
     2572            WriteCIFitem(self.fp, '_exptl_crystal.density_diffrn',
     2573                    G2mth.ValEsd(density,-0.001))                   
     2574
     2575            WriteCIFitem(self.fp, '_symmetry.cell_setting',
     2576                         phasedict['General']['SGData']['SGSys'])
     2577
     2578            spacegroup = phasedict['General']['SGData']['SpGrp'].strip()
     2579            # regularize capitalization and remove trailing H/R
     2580            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
     2581            WriteCIFitem(self.fp, '_symmetry.space_group_name_H-M',spacegroup)
     2582
     2583            # generate symmetry operations including centering and center of symmetry
     2584            SymOpList,offsetList,symOpList,G2oprList,G2opcodes = G2spc.AllOps(
     2585                phasedict['General']['SGData'])
     2586            WriteCIFitem(self.fp, 'loop_\n    _space_group.symop_id\n    _space_group.symop_operation_xyz')
     2587            for i,op in enumerate(SymOpList,start=1):
     2588                WriteCIFitem(self.fp, '   {:3d}  {:}'.format(i,op.lower()))
     2589
     2590            # loop over histogram(s) used in this phase
     2591            if not oneblock and not self.quickmode:
     2592                # report pointers to the histograms used in this phase
     2593                histlist = []
     2594                for hist in self.Phases[phasenam]['Histograms']:
     2595                    if self.Phases[phasenam]['Histograms'][hist]['Use']:
     2596                        if phasebyhistDict.get(hist):
     2597                            phasebyhistDict[hist].append(phasenam)
     2598                        else:
     2599                            phasebyhistDict[hist] = [phasenam,]
     2600                        blockid = datablockidDict.get(hist)
     2601                        if not blockid:
     2602                            print("Internal error: no block for data. Phase "+str(
     2603                                phasenam)+" histogram "+str(hist))
     2604                            histlist = []
     2605                            break
     2606                        histlist.append(blockid)
     2607
     2608                if len(histlist) == 0:
     2609                    WriteCIFitem(self.fp, '# Note: phase has no associated data')
     2610
     2611            # report atom params
     2612            try:
     2613                self.labellist
     2614            except AttributeError:
     2615                self.labellist = []
     2616            WriteAtomsMM(self.fp, self.Phases[phasenam], phasenam,
     2617                              self.parmDict, self.sigDict,
     2618                                  self.OverallParms['Rigid bodies'])
     2619
     2620            keV = None             
     2621            if oneblock: # get xray wavelength
     2622                lamlist = []
     2623                for hist in self.Histograms:
     2624                    if 'X' not in self.Histograms[hist]['Instrument Parameters'][0]['Type'][0]:
     2625                        continue
     2626                    for k in ('Lam','Lam1'):
     2627                        if k in self.Histograms[hist]['Instrument Parameters'][0]:
     2628                            lamlist.append(self.Histograms[hist]['Instrument Parameters'][0][k][0])
     2629                            break
     2630                if len(lamlist) == 1:
     2631                    keV = 12.397639/lamlist[0]
     2632
     2633            # report cell contents                       
     2634            WriteCompositionMM(self.fp, self.Phases[phasenam], phasenam, self.parmDict, self.quickmode, keV)
    22162635            if not self.quickmode and phasedict['General']['Type'] == 'nuclear':      # report distances and angles
    22172636                WriteDistances(phasenam)
     
    30113430            #phaseblk = self.Phases[phaseOnly] # pointer to current phase info
    30123431            # report the phase info
    3013             WritePhaseInfo(phaseOnly)
     3432            if self.Phases[phaseOnly]['General']['Type'] == 'macromolecular':
     3433                WritePhaseInfoMM(phaseOnly)
     3434            else:
     3435                WritePhaseInfo(phaseOnly)
    30143436            return
    30153437        elif histOnly: #====Histogram only CIF ================================
     
    32183640            writeCIFtemplate(self.Phases[phasenam]['General'],'phase',phasenam) # write phase template
    32193641            # report the phase info
    3220             WritePhaseInfo(phasenam,False)
     3642            if self.Phases[phasenam]['General']['Type'] == 'macromolecular':
     3643                WritePhaseInfoMM(phasenam,False)
     3644            else:
     3645                WritePhaseInfo(phasenam,False)
    32213646            if hist.startswith("PWDR"):  # this is invoked for single-block CIFs
    32223647                # preferred orientation
     
    35333958                    # report the phase
    35343959                    writeCIFtemplate(self.Phases[phasenam]['General'],'phase',phasenam) # write phase template
    3535                     WritePhaseInfo(phasenam,False,False)
     3960                    if self.Phases[phasenam]['General']['Type'] == 'macromolecular':
     3961                        WritePhaseInfoMM(phasenam,False,False)
     3962                    else:
     3963                        WritePhaseInfo(phasenam,False,False)
    35363964                    # preferred orientation
    35373965                    if self.ifPWDR:
Note: See TracChangeset for help on using the changeset viewer.