Changeset 4425 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
May 17, 2020 4:31:21 PM (17 months ago)
Author:
vondreele
Message:

fix generation of fullrmc pdb file for structures with partial substitutions & vacancies.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r4420 r4425  
    22072207    return fname
    22082208
    2209 def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):
    2210    
    2211     def findDup(Atoms):
    2212         Dup = []
    2213         Fracs = []
    2214         for iat1,at1 in enumerate(Atoms):
    2215             if any([at1[0] in dup for dup in Dup]):
    2216                 continue
    2217             else:
    2218                 Dup.append([at1[0],])
    2219                 Fracs.append([at1[6],])
    2220             for iat2,at2 in enumerate(Atoms[(iat1+1):]):
    2221                 if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
    2222                     Dup[-1] += [at2[0],]
    2223                     Fracs[-1]+= [at2[6],]
    2224         return Dup,Fracs
     2209def findDup(Atoms):
     2210    Dup = []
     2211    Fracs = []
     2212    for iat1,at1 in enumerate(Atoms):
     2213        if any([at1[0] in dup for dup in Dup]):
     2214            continue
     2215        else:
     2216            Dup.append([at1[0],])
     2217            Fracs.append([at1[6],])
     2218        for iat2,at2 in enumerate(Atoms[(iat1+1):]):
     2219            if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
     2220                Dup[-1] += [at2[0],]
     2221                Fracs[-1]+= [at2[6],]
     2222    return Dup,Fracs
     2223
     2224def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):   
    22252225   
    22262226    Meta = RMCPdict['metadata']
     
    27532753    generalData = Phase['General']
    27542754    Atseq = RMCPdict['atSeq']
     2755    Dups,Fracs = findDup(Phase['Atoms'])
     2756    Sfracs = [np.cumsum(fracs) for fracs in Fracs]
    27552757    Supercell = RMCPdict['SuperCell']
    27562758    Cell = generalData['Cell'][1:7]
     
    27612763    newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=False)
    27622764    Atoms = newPhase['Atoms']
    2763     XYZ = np.array([atom[3:6] for atom in Atoms]).T
     2765
     2766    Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
     2767    Natm = np.count_nonzero(Natm-1)
     2768    Atoms = newPhase['Atoms']
     2769    Satoms = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms,5),4),3)
     2770    Datoms = [[atom for atom in Satoms if atom[0] in dup] for dup in Dups]
     2771    Natoms = []
     2772    for idup,dup in enumerate(Dups):
     2773        ldup = len(dup)
     2774        datoms = Datoms[idup]
     2775        natm = len(datoms)
     2776        i = 0
     2777        while i < natm:
     2778            atoms = datoms[i:i+ldup]
     2779            try:
     2780                atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
     2781                Natoms.append(atom)
     2782            except IndexError:      #what about vacancies?
     2783                if 'Va' not in Atseq:
     2784                    Atseq.append('Va')
     2785                    RMCPdict['aTypes']['Va'] = 0.0
     2786                atom = atoms[0]
     2787                atom[1] = 'Va'
     2788                Natoms.append(atom)
     2789            i += ldup
     2790
     2791
     2792
     2793
     2794
     2795    XYZ = np.array([atom[3:6] for atom in Natoms]).T
    27642796    XYZptp = np.array([ma.ptp(XYZ[0]),ma.ptp(XYZ[1]),ma.ptp(XYZ[2])])/2.
    27652797    Cell = newPhase['General']['Cell'][1:7]
     
    27802812    if RMCPdict['byMolec']:
    27812813        NPM = RMCPdict['Natoms']
    2782         for iat,atom in enumerate(Atoms):
     2814        for iat,atom in enumerate(Natoms):
    27832815            XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian;residue = 'RMC'
    27842816            fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
     
    27872819    else:
    27882820        for atm in Atseq:
    2789             for iat,atom in enumerate(Atoms):
     2821            for iat,atom in enumerate(Natoms):
    27902822                if atom[1] == atm:
    27912823                    XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian
Note: See TracChangeset for help on using the changeset viewer.