Changeset 4455 for trunk/GSASIIphsGUI.py
 Timestamp:
 Jun 1, 2020 7:16:10 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIphsGUI.py
r4454 r4455 1224 1224 DisAglCtls['BondRadii'].append(dis) 1225 1225 return DisAglCtls['AtomTypes'],DisAglCtls['BondRadii'] 1226 1227 def 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,symCellAtmsxyzA)**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 1226 1262 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'][Opr1][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[opNum1] 1299 # atom[cs1] = 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,posInAllCellsxyzA)**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 = xyzatom[cx:cx+3]+item[3] 1307 # newAtom = atom[:] 1308 # newAtom[cx:cx+3] = xyz 1309 # newAtom[cs1] += 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 1314 def FindCoordination(ind,data,neighborArray,coordsArray,cmx=0,targets=None): 1315 'Find atoms coordinating atom ind, speedup version' 1229 1316 generalData = data['General'] 1230 1317 Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7]) … … 1234 1321 SGData = generalData['SGData'] 1235 1322 cellArray = G2lat.CellBlock(1) 1236 1237 1323 newAtomList = [] 1238 1324 atomA = atomData[ind] 1239 1325 xyzA = np.array(atomA[cx:cx+3]) 1326 lblA = atomA[0] 1240 1327 indA = atomTypes.index(atomA[ct]) 1241 1328 for atomB in atomData: 1242 if targets and atomB[ct] not in targets: 1243 1329 if targets and atomB[ct] not in targets: continue 1330 if atomB[0] not in neighborArray[lblA]: continue 1244 1331 indB = atomTypes.index(atomB[ct]) 1245 sumR = radii[indA]+radii[indB]1332 sumR = data['Drawing']['radiusFactor']*(radii[indA]+radii[indB]) 1246 1333 xyzB = np.array(atomB[cx:cx+3]) 1247 1334 Uij = atomB[cs+5:cs+5+6] 1335 coords = [] 1336 symMisc = [] 1248 1337 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+ixyzA 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 1251 1349 Opr = abs(item[2])%100 1252 1350 M = SGData['SGOps'][Opr1][0] … … 1260 1358 atom[cs1] = str(item[2])+'+' 1261 1359 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,posInAllCellsxyzA)**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 = xyzatom[cx:cx+3]+item[3] 1268 newAtom = atom[:] 1269 newAtom[cx:cx+3] = xyz 1270 newAtom[cs1] += 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 = xyzsymCoords[isym]+item[3] 1365 atom[cs1] += str(int(round(C[0])))+','+str(int(round(C[1])))+','+str(int(round(C[2]))) 1366 newAtomList.append(atom) 1272 1367 return newAtomList 1273 1368 … … 7854 7949 7855 7950 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 ''' 7856 7956 indx = getAtomSelections(drawAtoms) 7857 7957 if not indx: return … … 7862 7962 dlg = wx.Dialog(G2frame,wx.ID_ANY,'Addition criteria', 7863 7963 pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE) 7964 dlg.CenterOnParent() 7864 7965 mainSizer = wx.BoxSizer(wx.VERTICAL) 7865 7966 mainSizer.Add(wx.StaticText(dlg,wx.ID_ANY,'Molecular completion parameters'),0,WACV) … … 7867 7968 topSizer.Add(wx.StaticText(dlg,wx.ID_ANY,'Max # of repetitions: '),0,WACV) 7868 7969 choices = [1,2,5,10,50] 7869 params = {'maxrep': 5, 'maxatm':1000}7970 params = {'maxrep':10, 'maxatm':1000} 7870 7971 topSizer.Add(G2G.EnumSelector(dlg,params,'maxrep', 7871 7972 [str(i) for i in choices],choices)) … … 7913 8014 if 'Mx' in colLabels: 7914 8015 cmx = colLabels.index('Mx') 8016 SGData = generalData['SGData'] 8017 cellArray = G2lat.CellBlock(1) 8018 neighborArray = FindCoordinationByLabel(data) 8019 8020 time1 = time.time() 7915 8021 pgbar = wx.ProgressDialog('Fill molecular coordination', 7916 'Passes done=0 %0',params['maxrep'], 8022 'Passes done=0 %0',params['maxrep'], 8023 parent=G2frame, 7917 8024 style = wx.PD_ELAPSED_TIMEwx.PD_AUTO_HIDEwx.PD_CAN_ABORT) 7918 8025 screenSize = wx.ClientDisplayRect() … … 7927 8034 for rep in range(params['maxrep']): 7928 8035 startlen = len(data['Drawing']['Atoms']) 8036 coordsArray = np.array([a[cx:cx+3] for a in atomData]) 7929 8037 addedAtoms = [] 7930 8038 for Ind,ind in enumerate(indx): 7931 addedAtoms += FindCoordination(ind,data, cmx,targets)8039 addedAtoms += FindCoordination(ind,data,neighborArray,coordsArray,cmx,targets) 7932 8040 GoOn = pgbar.Update(rep+1, 7933 8041 newmsg='Passes done={} atom #{} of {}' 7934 8042 .format(rep+1,Ind+1,len(indx))) 7935 wx.Yield()7936 8043 if not GoOn[0]: break 7937 8044 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 7940 8049 if len(addedAtoms) == 0: break 7941 8050 added += len(addedAtoms) 7942 8051 if added > params['maxatm']: break 7943 8052 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'] 7945 8055 UpdateDrawAtoms() 7946 8056 G2plt.PlotStructure(G2frame,data) 7947 #GSASIIpath.IPyBreak()7948 8057 pgbar.Destroy() 7949 8058 UpdateDrawAtoms() 7950 8059 drawAtoms.ClearSelection() 7951 8060 G2plt.PlotStructure(G2frame,data) 7952 #GSASIIpath.IPyBreak()7953 8061 7954 8062 def FillCoordSphere(event):
Note: See TracChangeset
for help on using the changeset viewer.