Changeset 4366 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Mar 13, 2020 5:44:04 AM (3 years ago)
Author:
vondreele
Message:

Add plotting of bond distance & bond angle distributions found in potential energy restraint in RMCProfile

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r4347 r4366  
    23832383            fl.write('  > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search'])
    23842384            for angle in RMCPdict['Potentials']['Angles']:
    2385                 fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f Ang %.2f Ang\n'%
     2385                fl.write('  > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'%
    23862386                    (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5]))
    23872387    if RMCPdict['useBVS']:
     
    25032503    fl.close()
    25042504    return fname
     2505
     2506def GetRMCBonds(general,RMCPdict,Atoms,bondList):
     2507    bondDist = []
     2508    Cell = general['Cell'][1:7]
     2509    Supercell =  RMCPdict['SuperCell']
     2510    Trans = np.eye(3)*np.array(Supercell)
     2511    Cell = G2lat.TransformCell(Cell,Trans)[:6]
     2512    Amat,Bmat = G2lat.cell2AB(Cell)
     2513    indices = (-1,0,1)
     2514    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
     2515    for bonds in bondList:
     2516        Oxyz = np.array(Atoms[bonds[0]][1:])
     2517        Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]])       
     2518        Dx = np.array([Txyz-Oxyz+unit for unit in Units])
     2519        Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2))
     2520        for dx in Dx.T:
     2521            bondDist.append(np.min(dx))
     2522    return np.array(bondDist)
     2523   
     2524def GetRMCAngles(general,RMCPdict,Atoms,angleList):
     2525    bondAngles = []
     2526    Cell = general['Cell'][1:7]
     2527    Supercell =  RMCPdict['SuperCell']
     2528    Trans = np.eye(3)*np.array(Supercell)
     2529    Cell = G2lat.TransformCell(Cell,Trans)[:6]
     2530    Amat,Bmat = G2lat.cell2AB(Cell)
     2531    indices = (-1,0,1)
     2532    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
     2533    for angle in angleList:
     2534        Oxyz = np.array(Atoms[angle[0]][1:])
     2535        TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]])       
     2536        TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]])       
     2537        DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat)
     2538        DAx = np.sqrt(np.sum(DAxV**2,axis=2))
     2539        DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat)
     2540        DBx = np.sqrt(np.sum(DBxV**2,axis=2))
     2541        iDAx = np.argmin(DAx,axis=0)
     2542        iDBx = np.argmin(DBx,axis=0)
     2543        for i,[iA,iB] in enumerate(zip(iDAx,iDBx)):
     2544            DAv = DAxV[iA,i]/DAx[iA,i]
     2545            DBv = DBxV[iB,i]/DBx[iB,i]
     2546            bondAngles.append(npacosd(np.sum(DAv*DBv)))
     2547    return np.array(bondAngles)
    25052548   
    25062549################################################################################
Note: See TracChangeset for help on using the changeset viewer.