Changeset 4448 for trunk/GSASIIphsGUI.py


Ignore:
Timestamp:
May 29, 2020 9:36:54 AM (3 years ago)
Author:
toby
Message:

add complete molecule (slow) & fix use of radii

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIphsGUI.py

    r4447 r4448  
    205205    '''
    206206    def __init__(self,parent,general,drawing,indx):
    207         wx.Dialog.__init__(self,parent,wx.ID_ANY,'Setup phase transformation',
     207        wx.Dialog.__init__(self,parent,wx.ID_ANY,'Supply sphere info',
    208208            pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE)
    209209        self.panel = wx.Panel(self)         #just a dummy - gets destroyed in Draw!
    … …  
    214214        self.centers = []
    215215        self.atomTypes = [[item,False] for item in self.General['AtomTypes']]
    216        
     216        self.CenterOnParent()
    217217        self.Draw()
    218218       
    … …  
    11861186    cell = generalData['Cell'][1:7]
    11871187    FindBondsDrawCell(data,cell)
     1188
     1189def getAtomRadii(data):
     1190    '''Get radii for atoms, using generalData['DisAglCtls']['BondRadii']
     1191    to override generalData['BondRadii'] when present. Fix to make sure
     1192    that all elements in generalData are present in DisAglCtls.
     1193    '''
     1194    generalData = data['General']
     1195    if 'DisAglCtls' not in generalData:
     1196        return generalData['AtomTypes'],generalData['BondRadii']
     1197    if 'BondRadii' not in generalData['DisAglCtls']:
     1198        return generalData['AtomTypes'],generalData['BondRadii']
     1199    DisAglCtls = generalData['DisAglCtls']
     1200    if len(generalData['BondRadii']) != len(DisAglCtls['BondRadii']):
     1201        for typ,dis in zip(generalData['AtomTypes'],generalData['BondRadii']):
     1202            if typ not in DisAglCtls['AtomTypes']:
     1203                DisAglCtls['AtomTypes'].append(typ)
     1204                DisAglCtls['AngleRadii'].append(dis)
     1205                DisAglCtls['BondRadii'].append(dis)
     1206    return DisAglCtls['AtomTypes'],DisAglCtls['BondRadii']
    11881207   
     1208def FindCoordination(ind,data,cmx=0,targets=None):
     1209    'Find atoms coordinating atom ind, spe'
     1210    generalData = data['General']
     1211    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     1212    atomTypes,radii = getAtomRadii(data)
     1213    atomData = data['Drawing']['Atoms']
     1214    numAtoms = len(atomData)
     1215    cx,ct,cs,ci = data['Drawing']['atomPtrs']
     1216    cij = ci+2
     1217    SGData = generalData['SGData']
     1218    cellArray = G2lat.CellBlock(1)
     1219   
     1220    newAtomList = []
     1221    atomA = atomData[ind]
     1222    xyzA = np.array(atomA[cx:cx+3])
     1223    indA = atomTypes.index(atomA[ct])
     1224    for atomB in atomData:
     1225        if targets and atomB[ct] not in targets:
     1226            continue
     1227        indB = atomTypes.index(atomB[ct])
     1228        sumR = radii[indA]+radii[indB]
     1229        xyzB = np.array(atomB[cx:cx+3])
     1230        Uij = atomB[cs+5:cs+5+6]
     1231        for item in G2spc.GenAtom(xyzB,SGData,False,Uij,True):
     1232            atom = copy.copy(atomB)
     1233            atom[cx:cx+3] = item[0]
     1234            Opr = abs(item[2])%100
     1235            M = SGData['SGOps'][Opr-1][0]
     1236            if cmx:
     1237                opNum = G2spc.GetOpNum(item[2],SGData)
     1238                mom = np.array(atom[cmx:cmx+3])
     1239                if SGData['SGGray']:
     1240                    atom[cmx:cmx+3] = np.inner(mom,M)*nl.det(M)
     1241                else:   
     1242                    atom[cmx:cmx+3] = np.inner(mom,M)*nl.det(M)*SpnFlp[opNum-1]
     1243            atom[cs-1] = str(item[2])+'+'
     1244            atom[cs+5:cs+5+6] = item[1]
     1245            posInAllCells = cellArray+np.array(atom[cx:cx+3])
     1246            dists = np.sqrt(np.sum(np.inner(Amat,posInAllCells-xyzA)**2,axis=0))
     1247            bonded = np.logical_and(dists < data['Drawing']['radiusFactor']*sumR, dists !=0)
     1248            for xyz in posInAllCells[bonded]:
     1249                if True in [np.allclose(np.array(xyz),np.array(atom[cx:cx+3]),atol=0.0002) for atom in atomData]: continue
     1250                C = xyz-atom[cx:cx+3]+item[3]
     1251                newAtom = atom[:]
     1252                newAtom[cx:cx+3] = xyz
     1253                newAtom[cs-1] += str(int(round(C[0])))+','+str(int(round(C[1])))+','+str(int(round(C[2])))
     1254                newAtomList.append(newAtom)
     1255    return newAtomList
     1256
    11891257def FindBondsDrawCell(data,cell):   
    11901258    '''uses numpy & masks - very fast even for proteins!
    … …  
    11971265    generalData = data['General']
    11981266    Amat,Bmat = G2lat.cell2AB(cell)
    1199     radii = generalData['BondRadii']
    1200 #    if generalData.get('DisAglCtls',{}):
    1201 #        radii = generalData['DisAglCtls']['BondRadii']
    1202     atomTypes = generalData['AtomTypes']
     1267    atomTypes,radii = getAtomRadii(data)
    12031268    try:
    12041269        indH = atomTypes.index('H')
    … …  
    76157680        if 'Mx' in colLabels:
    76167681            cmx = colLabels.index('Mx')
    7617         generalData = data['General']
    76187682        SGData = generalData['SGData']
    76197683        SpnFlp = SGData.get('SpnFlp',[])
    … …  
    77287792        G2plt.PlotStructure(G2frame,data)
    77297793           
     7794    def FillMolecule(event):
     7795        indx = getAtomSelections(drawAtoms)
     7796        if not indx: return
     7797        generalData = data['General']
     7798        Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     7799        atomTypes,radii = getAtomRadii(data)
     7800       
     7801        dlg = wx.Dialog(G2frame,wx.ID_ANY,'Addition criteria',
     7802            pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE)
     7803        mainSizer = wx.BoxSizer(wx.VERTICAL)
     7804        mainSizer.Add(wx.StaticText(dlg,wx.ID_ANY,'Molecular completion parameters'),0,WACV)
     7805        topSizer = wx.BoxSizer(wx.HORIZONTAL)
     7806        topSizer.Add(wx.StaticText(dlg,wx.ID_ANY,'Max # of repetitions: '),0,WACV)
     7807        choices = [1,2,5,10,50]
     7808        params = {'maxrep':5, 'maxatm':1000}
     7809        topSizer.Add(G2G.EnumSelector(dlg,params,'maxrep',
     7810                        [str(i) for i in choices],choices))
     7811        mainSizer.Add(topSizer,0,WACV)
     7812        topSizer = wx.BoxSizer(wx.HORIZONTAL)
     7813        topSizer.Add(wx.StaticText(dlg,wx.ID_ANY,'Max # of added atoms: '),0,WACV)
     7814        choices = [100,500,1000,5000,10000]
     7815        topSizer.Add(G2G.EnumSelector(dlg,params,'maxatm',
     7816                        [str(i) for i in choices],choices))
     7817        mainSizer.Add(topSizer,0,WACV)
     7818        mainSizer.Add(wx.StaticText(dlg,wx.ID_ANY,'Atom types to add:'),0,WACV)
     7819        atSizer = wx.BoxSizer(wx.HORIZONTAL)
     7820        for i,item in enumerate(atomTypes):
     7821            params[item] = True
     7822            atm = G2G.G2CheckBox(dlg,item,params,item)
     7823            atSizer.Add(atm,0,WACV)
     7824        mainSizer.Add(atSizer,0,WACV)
     7825       
     7826        OkBtn = wx.Button(dlg,-1,"Ok")
     7827        OkBtn.Bind(wx.EVT_BUTTON, lambda x: dlg.EndModal(wx.ID_OK))
     7828        cancelBtn = wx.Button(dlg,-1,"Cancel")
     7829        cancelBtn.Bind(wx.EVT_BUTTON, lambda x: dlg.EndModal(wx.ID_CANCEL))
     7830        btnSizer = wx.BoxSizer(wx.HORIZONTAL)
     7831        btnSizer.Add((20,20),1)
     7832        btnSizer.Add(OkBtn)
     7833        btnSizer.Add((20,20),1)
     7834        btnSizer.Add(cancelBtn)
     7835        btnSizer.Add((20,20),1)
     7836       
     7837        mainSizer.Add(btnSizer,0,wx.EXPAND|wx.BOTTOM|wx.TOP, 10)
     7838        dlg.SetSizer(mainSizer)
     7839        dlg.Fit()
     7840        if dlg.ShowModal() != wx.ID_OK:
     7841            dlg.Destroy()
     7842            return
     7843        dlg.Destroy()
     7844       
     7845        try:
     7846            indH = atomTypes.index('H')
     7847            radii[indH] = 0.5
     7848        except:
     7849            pass
     7850        atomData = data['Drawing']['Atoms']
     7851        numAtoms = len(atomData)
     7852        cx,ct,cs,ci = data['Drawing']['atomPtrs']
     7853        cij = ci+2
     7854        colLabels = [drawAtoms.GetColLabelValue(c) for c in range(drawAtoms.GetNumberCols())]
     7855        cmx = 0
     7856        if 'Mx' in colLabels:
     7857            cmx = colLabels.index('Mx')
     7858        SGData = generalData['SGData']
     7859        cellArray = G2lat.CellBlock(1)
     7860        nind = len(indx)
     7861        pgbar = wx.ProgressDialog('Fill molecular coordination',
     7862                                    'Passes done=0 %0',params['maxrep'],
     7863            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_CAN_ABORT)
     7864        screenSize = wx.ClientDisplayRect()
     7865        Size = pgbar.GetSize()
     7866        if 50 < Size[0] < 500: # sanity check on size, since this fails w/Win & wx3.0
     7867            pgbar.SetSize((int(Size[0]*1.2),Size[1])) # increase size a bit along x
     7868            pgbar.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
     7869        pgbar.Raise()
     7870        wx.Yield()
     7871        added = 0
     7872        targets = [item for item in atomTypes if params[item]]
     7873        for rep in range(params['maxrep']):
     7874            startlen = len(data['Drawing']['Atoms'])
     7875            addedAtoms = []
     7876            for Ind,ind in enumerate(indx):
     7877                addedAtoms += FindCoordination(ind,data,cmx,targets)
     7878                GoOn = pgbar.Update(rep+1,
     7879                    newmsg='Passes done={} atom #{} of {}'
     7880                                        .format(rep+1,Ind+1,len(indx)))
     7881                wx.Yield()
     7882                if not GoOn[0]: break
     7883            if not GoOn[0]: break
     7884            print('pass {} processed {} atoms adding {}'.format(
     7885                rep+1,len(indx),len(addedAtoms)))
     7886            if len(addedAtoms) == 0: break
     7887            added += len(addedAtoms)
     7888            if added > params['maxatm']: break
     7889            data['Drawing']['Atoms'] += addedAtoms
     7890            indx = list(range(startlen,startlen+len(addedAtoms)))
     7891            UpdateDrawAtoms()
     7892            G2plt.PlotStructure(G2frame,data)
     7893            #GSASIIpath.IPyBreak()
     7894        pgbar.Destroy()   
     7895        UpdateDrawAtoms()
     7896        drawAtoms.ClearSelection()
     7897        G2plt.PlotStructure(G2frame,data)
     7898        #GSASIIpath.IPyBreak()
     7899       
    77307900    def FillCoordSphere(event):
    77317901        indx = getAtomSelections(drawAtoms)
    … …  
    77337903        generalData = data['General']
    77347904        Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
    7735         radii = generalData['BondRadii']
    7736         atomTypes = generalData['AtomTypes']
     7905        atomTypes,radii = getAtomRadii(data)
    77377906        try:
    77387907            indH = atomTypes.index('H')
    77397908            radii[indH] = 0.5
    77407909        except:
    7741             pass           
     7910            pass
    77427911        indx.sort()
    77437912        atomData = data['Drawing']['Atoms']
    … …  
    77457914        cx,ct,cs,ci = data['Drawing']['atomPtrs']
    77467915        cij = ci+2
     7916        colLabels = [drawAtoms.GetColLabelValue(c) for c in range(drawAtoms.GetNumberCols())]
     7917        cmx = 0
     7918        if 'Mx' in colLabels:
     7919            cmx = colLabels.index('Mx')
    77477920        SGData = generalData['SGData']
    77487921        cellArray = G2lat.CellBlock(1)
    … …  
    77567929            pgbar.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
    77577930        for Ind,ind in enumerate(indx):
    7758             atomA = atomData[ind]
    7759             xyzA = np.array(atomA[cx:cx+3])
    7760             indA = atomTypes.index(atomA[ct])
    7761             for atomB in atomData[:numAtoms]:
    7762                 indB = atomTypes.index(atomB[ct])
    7763                 sumR = radii[indA]+radii[indB]
    7764                 xyzB = np.array(atomB[cx:cx+3])
    7765                 for xyz in cellArray+xyzB:
    7766                     dist = np.sqrt(np.sum(np.inner(Amat,xyz-xyzA)**2))
    7767                     if 0 < dist <= data['Drawing']['radiusFactor']*sumR:
    7768                         if noDuplicate(xyz,atomData):
    7769                             oprB = atomB[cs-1]
    7770                             C = xyz-xyzB
    7771                             newOp = '1+'+str(int(round(C[0])))+','+str(int(round(C[1])))+','+str(int(round(C[2])))
    7772                             newAtom = atomB[:]
    7773                             newAtom[cx:cx+3] = xyz
    7774                             newAtom[cs-1] = G2spc.StringOpsProd(oprB,newOp,SGData)
    7775                             atomData.append(newAtom[:cij+9])  #not SS stuff
     7931            atomData += FindCoordination(ind,data,cmx)
    77767932            GoOn = pgbar.Update(Ind,newmsg='Atoms done=%d'%(Ind))
    7777             if not GoOn[0]:
    7778                 break
     7933            if not GoOn[0]: break
    77797934        pgbar.Destroy()   
    77807935        data['Drawing']['Atoms'] = atomData
    … …  
    1167411829        G2frame.Bind(wx.EVT_MENU, DrawAtomColor, id=G2G.wxID_DRAWATOMCOLOR)
    1167511830        G2frame.Bind(wx.EVT_MENU, ResetAtomColors, id=G2G.wxID_DRAWATOMRESETCOLOR)
    11676         #G2frame.Bind(wx.EVT_MENU, OnEditAtomRadii, id=G2G.wxID_DRWAEDITRADII)           # TODO: removed until fixed
     11831        G2frame.Bind(wx.EVT_MENU, OnEditAtomRadii, id=G2G.wxID_DRWAEDITRADII)
    1167711832        G2frame.Bind(wx.EVT_MENU, SetViewPoint, id=G2G.wxID_DRAWVIEWPOINT)
    1167811833        G2frame.Bind(wx.EVT_MENU, AddSymEquiv, id=G2G.wxID_DRAWADDEQUIV)
    … …  
    1169011845        G2frame.Bind(wx.EVT_MENU, OnRestraint, id=G2G.wxID_DRAWRESTRCHIRAL)
    1169111846        G2frame.Bind(wx.EVT_MENU, OnDefineRB, id=G2G.wxID_DRAWDEFINERB)
     11847        G2frame.Bind(wx.EVT_MENU, FillMolecule, id=G2G.wxID_DRAWADDMOLECULE)
     11848       
    1169211849        # RB Models
    1169311850        FillSelectPageMenu(TabSelectionIdDict, G2frame.dataWindow.RigidBodiesMenu)
Note: See TracChangeset for help on using the changeset viewer.