Changeset 4455 for trunk


Ignore:
Timestamp:
Jun 1, 2020 7:16:10 PM (4 years ago)
Author:
toby
Message:

fix and speed up the "Complete Molecule" command; w/o for mult calc bug

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIconstrGUI.py

    r4431 r4455  
    3838import GSASIIdataGUI as G2gd
    3939import GSASIIctrlGUI as G2G
     40import GSASIIfiles as G2fl
    4041import GSASIIplot as G2plt
    4142import GSASIIobj as G2obj
    … …  
    25112512            rd.errors = ""
    25122513            if not rd.ContentsValidator(filename):   # Report error
    2513                 G2fil.G2Print("Warning: File {} has a validation error".format(filename))
     2514                G2fl.G2Print("Warning: File {} has a validation error".format(filename))
    25142515                return
    25152516            if len(rd.selections) > 1:
    … …  
    25212522            try:
    25222523                flag = rd.Reader(filename)
    2523             except:
    2524                 G2fil.G2Print("Warning: read of file {} failed".format(filename))
     2524            except Exception as msg:
     2525                G2fl.G2Print("Warning: read of file {} failed\n{}".format(
     2526                    filename,rd.errors))
     2527                if GSASIIpath.GetConfigValue('debug'):
     2528                    print(msg)
     2529                    import traceback
     2530                    print (traceback.format_exc())
     2531                    GSASIIpath.IPyBreak()
    25252532                return
    25262533

  • trunk/GSASIIphsGUI.py

    r4454 r4455  
    12241224                DisAglCtls['BondRadii'].append(dis)
    12251225    return DisAglCtls['AtomTypes'],DisAglCtls['BondRadii']
     1226
     1227def FindCoordinationByLabel(data):
     1228    '''Map out molecular connectivity by determining the atoms bonded
     1229    to each atom, by label. The atoms boced to each atom in the asymmetric
     1230    unit is determined and returned in a dict. Works best
     1231    '''
     1232    generalData = data['General']
     1233    cx,ct,cs,cia = generalData['AtomPtrs']
     1234    atomTypes,radii = getAtomRadii(data)
     1235    SGData = generalData['SGData']
     1236    cellArray = G2lat.CellBlock(1)
     1237    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     1238
     1239    neighborArray = {}
     1240    error = ''
     1241    for atomA in data['Atoms']:
     1242        lblA = atomA[0]
     1243        if lblA in neighborArray:
     1244            if error: error += ', '
     1245            error += lblA
     1246        else:
     1247            neighborArray[lblA] = []
     1248        xyzA = np.array(atomA[cx:cx+3])
     1249        indA = atomTypes.index(atomA[ct])
     1250        for atomB in data['Atoms']:
     1251            indB = atomTypes.index(atomB[ct])
     1252            sumR = data['Drawing']['radiusFactor']*(radii[indA]+radii[indB])
     1253            symAtms = [atm[0] for atm in  G2spc.GenAtom(np.array(atomB[cx:cx+3]),SGData,False,6*[0],True)]
     1254            symCellAtms = np.concatenate([cellArray+i for i in symAtms])
     1255            dists = np.sqrt(np.sum(np.inner(Amat,symCellAtms-xyzA)**2,axis=0))
     1256            if np.any(np.logical_and(dists < sumR, dists != 0)):
     1257                if atomB[0] not in neighborArray[lblA]:
     1258                    neighborArray[lblA].append(atomB[0])
     1259    if error:
     1260        print('Warning, duplicated atom labels:',error)
     1261    return neighborArray
    12261262   
    1227 def FindCoordination(ind,data,cmx=0,targets=None):
    1228     'Find atoms coordinating atom ind, spe'
     1263# def FindCoordination(ind,data,cmx=0,targets=None):
     1264#     'Find atoms coordinating atom ind, somewhat faster version'
     1265#     time1 = time.time()
     1266#     generalData = data['General']
     1267#     Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     1268#     atomTypes,radii = getAtomRadii(data)
     1269#     atomData = data['Drawing']['Atoms']
     1270#     numAtoms = len(atomData)
     1271#     cx,ct,cs,ci = data['Drawing']['atomPtrs']
     1272#     cij = ci+2
     1273#     SGData = generalData['SGData']
     1274#     cellArray = G2lat.CellBlock(1)
     1275   
     1276#     newAtomList = []
     1277#     atomA = atomData[ind]
     1278#     xyzA = np.array(atomA[cx:cx+3])
     1279#     indA = atomTypes.index(atomA[ct])
     1280#     for atomB in atomData:
     1281#         if targets and atomB[ct] not in targets:
     1282#             continue
     1283#         indB = atomTypes.index(atomB[ct])
     1284#         sumR = radii[indA]+radii[indB]
     1285#         xyzB = np.array(atomB[cx:cx+3])
     1286#         Uij = atomB[cs+5:cs+5+6]
     1287#         for item in G2spc.GenAtom(xyzB,SGData,False,Uij,True):
     1288#             atom = copy.copy(atomB)
     1289#             atom[cx:cx+3] = item[0]
     1290#             Opr = abs(item[2])%100
     1291#             M = SGData['SGOps'][Opr-1][0]
     1292#             if cmx:
     1293#                 opNum = G2spc.GetOpNum(item[2],SGData)
     1294#                 mom = np.array(atom[cmx:cmx+3])
     1295#                 if SGData['SGGray']:
     1296#                     atom[cmx:cmx+3] = np.inner(mom,M)*nl.det(M)
     1297#                 else:   
     1298#                     atom[cmx:cmx+3] = np.inner(mom,M)*nl.det(M)*SpnFlp[opNum-1]
     1299#             atom[cs-1] = str(item[2])+'+'
     1300#             atom[cs+5:cs+5+6] = item[1]
     1301#             posInAllCells = cellArray+np.array(atom[cx:cx+3])
     1302#             dists = np.sqrt(np.sum(np.inner(Amat,posInAllCells-xyzA)**2,axis=0))
     1303#             bonded = np.logical_and(dists < data['Drawing']['radiusFactor']*sumR, dists !=0)
     1304#             for xyz in posInAllCells[bonded]:
     1305#                 if True in [np.allclose(np.array(xyz),np.array(atom[cx:cx+3]),atol=0.0002) for atom in atomData]: continue
     1306#                 C = xyz-atom[cx:cx+3]+item[3]
     1307#                 newAtom = atom[:]
     1308#                 newAtom[cx:cx+3] = xyz
     1309#                 newAtom[cs-1] += str(int(round(C[0])))+','+str(int(round(C[1])))+','+str(int(round(C[2])))
     1310#                 newAtomList.append(newAtom)
     1311#     print ('Search time: %.2fs'%(time.time()-time1))
     1312#     return newAtomList
     1313
     1314def FindCoordination(ind,data,neighborArray,coordsArray,cmx=0,targets=None):
     1315    'Find atoms coordinating atom ind, speed-up version'
    12291316    generalData = data['General']
    12301317    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
    … …  
    12341321    SGData = generalData['SGData']
    12351322    cellArray = G2lat.CellBlock(1)
    1236    
    12371323    newAtomList = []
    12381324    atomA = atomData[ind]
    12391325    xyzA = np.array(atomA[cx:cx+3])
     1326    lblA = atomA[0]
    12401327    indA = atomTypes.index(atomA[ct])
    12411328    for atomB in atomData:
    1242         if targets and atomB[ct] not in targets:
    1243             continue
     1329        if targets and atomB[ct] not in targets: continue
     1330        if atomB[0] not in neighborArray[lblA]: continue
    12441331        indB = atomTypes.index(atomB[ct])
    1245         sumR = radii[indA]+radii[indB]
     1332        sumR = data['Drawing']['radiusFactor']*(radii[indA]+radii[indB])
    12461333        xyzB = np.array(atomB[cx:cx+3])
    12471334        Uij = atomB[cs+5:cs+5+6]
     1335        coords = []
     1336        symMisc = []
    12481337        for item in G2spc.GenAtom(xyzB,SGData,False,Uij,True):
    1249             atom = copy.copy(atomB)
    1250             atom[cx:cx+3] = item[0]
     1338            coords.append(item[0])
     1339            symMisc.append(item[1:4])
     1340        symCoords = np.array(coords)
     1341        dists = np.sqrt(np.sum(np.inner(Amat,
     1342                np.array([symCoords+i-xyzA for i in cellArray])
     1343                                           )**2,axis=0))
     1344        for icell,isym in np.argwhere(np.logical_and(dists < sumR, dists != 0.0)):
     1345            xyz = symCoords[isym] + cellArray[icell]
     1346            item = [None]+list(symMisc[isym])
     1347            atom = copy.deepcopy(atomB)
     1348            atom[cx:cx+3] = xyz
    12511349            Opr = abs(item[2])%100
    12521350            M = SGData['SGOps'][Opr-1][0]
    … …  
    12601358            atom[cs-1] = str(item[2])+'+'
    12611359            atom[cs+5:cs+5+6] = item[1]
    1262             posInAllCells = cellArray+np.array(atom[cx:cx+3])
    1263             dists = np.sqrt(np.sum(np.inner(Amat,posInAllCells-xyzA)**2,axis=0))
    1264             bonded = np.logical_and(dists < data['Drawing']['radiusFactor']*sumR, dists !=0)
    1265             for xyz in posInAllCells[bonded]:
    1266                 if True in [np.allclose(np.array(xyz),np.array(atom[cx:cx+3]),atol=0.0002) for atom in atomData]: continue
    1267                 C = xyz-atom[cx:cx+3]+item[3]
    1268                 newAtom = atom[:]
    1269                 newAtom[cx:cx+3] = xyz
    1270                 newAtom[cs-1] += str(int(round(C[0])))+','+str(int(round(C[1])))+','+str(int(round(C[2])))
    1271                 newAtomList.append(newAtom)
     1360            # have we already found an atom at this site?
     1361            if np.any(np.all(np.isclose(xyz,coordsArray,atol=0.0002),axis=1)): continue
     1362            # are we going to add it already?
     1363            if True in [np.allclose(np.array(xyz),np.array(a[cx:cx+3]),atol=0.0002) for a in newAtomList]: continue
     1364            C = xyz-symCoords[isym]+item[3]
     1365            atom[cs-1] += str(int(round(C[0])))+','+str(int(round(C[1])))+','+str(int(round(C[2])))
     1366            newAtomList.append(atom)
    12721367    return newAtomList
    12731368
    … …  
    78547949           
    78557950    def FillMolecule(event):
     7951        '''This is called by the Complete Molecule command. It adds a layer of bonded atoms
     7952        of the selected types for all selected atoms in the Draw Atoms table. If the number
     7953        of repetitions is greater than one, the added atoms (other than H atoms, which are assumed
     7954        to only have one bond) are then searched for the next surrounding layer of bonded atoms.
     7955        '''
    78567956        indx = getAtomSelections(drawAtoms)
    78577957        if not indx: return
    … …  
    78627962        dlg = wx.Dialog(G2frame,wx.ID_ANY,'Addition criteria',
    78637963            pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE)
     7964        dlg.CenterOnParent()
    78647965        mainSizer = wx.BoxSizer(wx.VERTICAL)
    78657966        mainSizer.Add(wx.StaticText(dlg,wx.ID_ANY,'Molecular completion parameters'),0,WACV)
    … …  
    78677968        topSizer.Add(wx.StaticText(dlg,wx.ID_ANY,'Max # of repetitions: '),0,WACV)
    78687969        choices = [1,2,5,10,50]
    7869         params = {'maxrep':5, 'maxatm':1000}
     7970        params = {'maxrep':10, 'maxatm':1000}
    78707971        topSizer.Add(G2G.EnumSelector(dlg,params,'maxrep',
    78717972                        [str(i) for i in choices],choices))
    … …  
    79138014        if 'Mx' in colLabels:
    79148015            cmx = colLabels.index('Mx')
     8016        SGData = generalData['SGData']
     8017        cellArray = G2lat.CellBlock(1)
     8018        neighborArray = FindCoordinationByLabel(data)
     8019
     8020        time1 = time.time()
    79158021        pgbar = wx.ProgressDialog('Fill molecular coordination',
    7916                                     'Passes done=0 %0',params['maxrep'],
     8022            'Passes done=0 %0',params['maxrep'],
     8023            parent=G2frame,
    79178024            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_CAN_ABORT)
    79188025        screenSize = wx.ClientDisplayRect()
    … …  
    79278034        for rep in range(params['maxrep']):
    79288035            startlen = len(data['Drawing']['Atoms'])
     8036            coordsArray = np.array([a[cx:cx+3] for a in atomData])
    79298037            addedAtoms = []
    79308038            for Ind,ind in enumerate(indx):
    7931                 addedAtoms += FindCoordination(ind,data,cmx,targets)
     8039                addedAtoms += FindCoordination(ind,data,neighborArray,coordsArray,cmx,targets)
    79328040                GoOn = pgbar.Update(rep+1,
    79338041                    newmsg='Passes done={} atom #{} of {}'
    79348042                                        .format(rep+1,Ind+1,len(indx)))
    7935                 wx.Yield()
    79368043                if not GoOn[0]: break
    79378044            if not GoOn[0]: break
    7938             print('pass {} processed {} atoms adding {}'.format(
    7939                 rep+1,len(indx),len(addedAtoms)))
     8045            print('pass {} processed {} atoms adding {}; Search time: {:.2f}s'.format(
     8046                rep+1,len(indx),len(addedAtoms),time.time()-time1))
     8047            time1 = time.time()           
     8048
    79408049            if len(addedAtoms) == 0: break
    79418050            added += len(addedAtoms)
    79428051            if added > params['maxatm']: break
    79438052            data['Drawing']['Atoms'] += addedAtoms
    7944             indx = list(range(startlen,startlen+len(addedAtoms)))
     8053            # atoms to search over (omit H)
     8054            indx = [i+startlen for i in range(len(addedAtoms)) if addedAtoms[i][ct] != 'H'] 
    79458055            UpdateDrawAtoms()
    79468056            G2plt.PlotStructure(G2frame,data)
    7947             #GSASIIpath.IPyBreak()
    79488057        pgbar.Destroy()   
    79498058        UpdateDrawAtoms()
    79508059        drawAtoms.ClearSelection()
    79518060        G2plt.PlotStructure(G2frame,data)
    7952         #GSASIIpath.IPyBreak()
    79538061       
    79548062    def FillCoordSphere(event):

  • trunk/GSASIIspc.py

    r4441 r4455  
    33243324                        Isym += 2**(jx-1)
    33253325    if Isym == 1073741824: Isym = 0
    3326     Mult = len(SGData['SGOps'])*Ncen*inv//Jdup
    3327          
     3326    try:
     3327        Mult = len(SGData['SGOps'])*Ncen*inv//Jdup
     3328    except: # patch because Jdup is not getting incremented for most atoms!
     3329        Mult = 0
     3330       
    33283331    return GetKNsym(str(Isym)),Mult,Ndup,dupDir
    33293332   
Note: See TracChangeset for help on using the changeset viewer.