Changeset 1930 for trunk/GSASIImath.py


Ignore:
Timestamp:
Jul 16, 2015 1:37:20 PM (7 years ago)
Author:
vondreele
Message:

allow exclude atoms from H-atom position calcs.
hydrogen add complete(?) - needs testing
new hydrogen update implemented

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r1929 r1930  
    401401        if j != Orig:
    402402            if AtNames[j] != notName:
    403                 Neigh.append([AtNames[j],dist[j]])
     403                Neigh.append([AtNames[j],dist[j],True])
    404404                Ids.append(Atoms[j][cia+8])
    405405    return Neigh,[OId,Ids]
    406406   
    407407def AddHydrogens(AtLookUp,General,Atoms,AddHydId):
     408   
     409    def getTransMat(RXYZ,OXYZ,TXYZ,Amat):
     410        Vec = np.inner(Amat,np.array([OXYZ-TXYZ[0],RXYZ-TXYZ[0]])).T           
     411        Vec /= np.sqrt(np.sum(Vec**2,axis=1))[:,np.newaxis]
     412        Mat2 = np.cross(Vec[0],Vec[1])      #UxV
     413        Mat2 /= np.sqrt(np.sum(Mat2**2))
     414        Mat3 = np.cross(Mat2,Vec[0])        #(UxV)xU
     415        return nl.inv(np.array([Vec[0],Mat2,Mat3]))       
    408416   
    409417    cx,ct,cs,cia = General['AtomPtrs']
     
    439447            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
    440448            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]))
     449            iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
    447450            a = 0.96*cosd(70.5)
    448451            b = 0.96*sind(70.5)
     
    460463            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
    461464            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]))
    468             print 'add 2 H'
    469             return []
     465            iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
     466            a = 0.93*cosd(60.)
     467            b = 0.93*sind(60.)
     468            Hpos = [[a,b,0],[a,-b,0]]
     469            Hpos = np.inner(Bmat,np.inner(iMat,Hpos).T).T+OXYZ
     470            return Hpos
    470471    else:   #2 bonds
    471472        if 'C' in Oatom[ct]:
     
    479480            Ratom = GetAtomsById(Atoms,AtLookUp,[AddHydId[2],])[0]
    480481            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]))
     482            iMat = getTransMat(RXYZ,OXYZ,TXYZ,Amat)
    487483            a = 0.82*cosd(70.5)
    488484            b = 0.82*sind(70.5)
Note: See TracChangeset for help on using the changeset viewer.