Changeset 1927


Ignore:
Timestamp:
Jul 14, 2015 3:55:03 PM (7 years ago)
Author:
vondreele
Message:

more hydrogen add work

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1926 r1927  
    376376    return XYZ
    377377
    378 def FindNeighbors(phase,FrstName,AtNames):
     378def FindNeighbors(phase,FrstName,AtNames,notName=''):
    379379    General = phase['General']
    380380    cx,ct,cs,cia = General['AtomPtrs']
     
    400400    for j in IndB[0]:
    401401        if j != Orig:
    402             Neigh.append([AtNames[j],dist[j]])
    403             Ids.append(Atoms[j][cia+8])
     402            if AtNames[j] != notName:
     403                Neigh.append([AtNames[j],dist[j]])
     404                Ids.append(Atoms[j][cia+8])
    404405    return Neigh,[OId,Ids]
    405406   
    406407def AddHydrogens(AtLookUp,General,Atoms,AddHydId):
     408   
    407409    cx,ct,cs,cia = General['AtomPtrs']
    408410    Cell = General['Cell'][1:7]
    409411    Amat,Bmat = G2lat.cell2AB(Cell)
    410     nBonds = AddHydId[2]+len(AddHydId[1])
     412    nBonds = AddHydId[-1]+len(AddHydId[1])
    411413    Oatom = GetAtomsById(Atoms,AtLookUp,[AddHydId[0],])[0]
    412414    OXYZ = np.array(Oatom[cx:cx+3])
     
    415417    DX = np.inner(Amat,TXYZ-OXYZ).T
    416418    if nBonds == 4:
    417         if AddHydId[2] == 1:
     419        if AddHydId[-1] == 1:
    418420            Vec = TXYZ-OXYZ
    419421            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2,axis=0))
     
    422424            Hpos = OXYZ-0.98*np.inner(Bmat,Vec).T/Len
    423425            return [Hpos,]
    424         elif AddHydId[2] == 2:
    425             print 'add 2 H'
    426             return []
     426        elif AddHydId[-1] == 2:
     427            Vec = np.inner(Amat,TXYZ-OXYZ).T
     428            Vec[0] += Vec[1]            #U - along bisector
     429            Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,np.newaxis]
     430            Mat2 = np.cross(Vec[0],Vec[1])      #UxV
     431            Mat2 /= np.sqrt(np.sum(Mat2**2))
     432            Mat3 = np.cross(Mat2,Vec[0])        #(UxV)xU
     433            iMat = nl.inv(np.array([Vec[0],Mat2,Mat3]))
     434            Hpos = np.array([[-0.97*cosd(54.75),0.97*sind(54.75),0.],
     435                [-0.97*cosd(54.75),-0.97*sind(54.75),0.]])
     436            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
     437            return Hpos
    427438        else:
    428             print 'add 3 H'
    429             return []           
     439            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
     440            RXYZ = np.array(Ratom[cx:cx+3])
     441            Vec = np.inner(Amat,np.array([OXYZ-TXYZ[0],RXYZ-TXYZ[0]])).T           
     442            Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,np.newaxis]
     443            Mat2 = np.cross(Vec[0],Vec[1])      #UxV
     444            Mat2 /= np.sqrt(np.sum(Mat2**2))
     445            Mat3 = np.cross(Mat2,Vec[0])        #(UxV)xU
     446            iMat = nl.inv(np.array([Vec[0],Mat2,Mat3]))
     447            a = 0.96*cosd(70.5)
     448            b = 0.96*sind(70.5)
     449            Hpos = np.array([[a,0.,-b],[a,-b*cosd(30.),0.5*b],[a,b*cosd(30.),0.5*b]])
     450            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
     451            return Hpos           
    430452    elif nBonds == 3:
    431         if AddHydId[2] == 1:
     453        if AddHydId[-1] == 1:
    432454            Vec = np.sum(TXYZ-OXYZ,axis=0)               
    433455            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2))
     
    435457            Hpos = OXYZ+Vec
    436458            return [Hpos,]
    437         elif AddHydId[2] == 2:
     459        elif AddHydId[-1] == 2:
     460            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
     461            RXYZ = np.array(Ratom[cx:cx+3])
     462            Vec = np.inner(Amat,np.array([OXYZ-TXYZ[0],RXYZ-TXYZ[0]])).T           
     463            Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,np.newaxis]
     464            Mat2 = np.cross(Vec[0],Vec[1])      #UxV
     465            Mat2 /= np.sqrt(np.sum(Mat2**2))
     466            Mat3 = np.cross(Mat2,Vec[0])        #(UxV)xU
     467            iMat = nl.inv(np.array([Vec[0],Mat2,Mat3]))
    438468            print 'add 2 H'
    439469            return []
    440470    else:   #2 bonds
    441471        if 'C' in Oatom[ct]:
    442             print 'sp atom',Oatom[ct-1]
    443             return []
     472            Vec = TXYZ[0]-OXYZ
     473            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2))
     474            Vec = -0.93*Vec/Len
     475            Hpos = OXYZ+Vec
     476            return [Hpos,]
    444477        elif 'O' in Oatom[ct]:
    445             #idea - construct ring at 0.82 from O atom & find high spot?
    446             print 'sp3 atom ',Oatom[ct-1]
    447             print 'add 1 H'
    448             return []
     478            mapData = General['Map']
     479            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
     480            RXYZ = np.array(Ratom[cx:cx+3])
     481            Vec = np.inner(Amat,np.array([OXYZ-TXYZ[0],RXYZ-TXYZ[0]])).T           
     482            Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,np.newaxis]
     483            Mat2 = np.cross(Vec[0],Vec[1])      #UxV
     484            Mat2 /= np.sqrt(np.sum(Mat2**2))
     485            Mat3 = np.cross(Mat2,Vec[0])        #(UxV)xU
     486            iMat = nl.inv(np.array([Vec[0],Mat2,Mat3]))
     487            a = 0.82*cosd(70.5)
     488            b = 0.82*sind(70.5)
     489            azm = np.arange(0.,360.,5.)
     490            Hpos = np.array([[a,b*cosd(x),b*sind(x)] for x in azm])
     491            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
     492            return Hpos
    449493    return []
    450    
    451    
    452        
    453    
    454494       
    455495def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
  • trunk/GSASIIphsGUI.py

    r1926 r1927  
    15721572                    nH = 4-len(neigh[1][0])
    15731573                bonds = dict(neigh[1][0])
     1574                nextName = ''
     1575                if len(bonds) == 1:
     1576                    nextName = bonds.keys()[0]
    15741577                for bond in bonds:
    15751578                    if 'C' in neigh[0]:
     
    15831586                        nH -= 1
    15841587                        break
     1588                nextneigh = []
     1589                if nextName:
     1590                    nextneigh = G2mth.FindNeighbors(data,nextName,AtNames,notName=neigh[0])
     1591                    neigh[1][1].append(nextneigh[1][1][0])
    15851592                neigh[2] = max(0,nH)  #set expected no. H's needed
    15861593                AddHydIds.append(neigh[1][1])
     
    15911598                    Nat = len(atomData)
    15921599                    Neigh = dlg.GetData()
    1593                     for ineigh,neigh in enumerate(Neigh):
     1600                    mapData = generalData['Map']
     1601                    mapError = False
     1602                    for ineigh,neigh in enumerate(Neigh):                       
    15941603                        AddHydIds[ineigh].append(neigh[2])
     1604                        if 'O' in neigh[0] and not len(mapData['rho']) or not 'delt-F' in mapData['MapType']:
     1605                            mapError = True
     1606                            continue                           
    15951607                        Hxyz = G2mth.AddHydrogens(AtLookUp,generalData,atomData,AddHydIds[ineigh])
    15961608                        for X in Hxyz:
    15971609                            Nat += 1
    15981610                            AtomAdd(X[0],X[1],X[2],'H','H(%d)'%(Nat))
    1599 
     1611                if mapError:
     1612                    G2frame.ErrorDialog('Add H atom error','Adding O-H atoms requires delt-F map')
    16001613                SetupGeneral()
    16011614                FillAtomsGrid(Atoms)
Note: See TracChangeset for help on using the changeset viewer.