Changeset 4425

Timestamp:

May 17, 2020 4:31:21 PM (3 years ago)

Author:

vondreele

Message:

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

Location:

trunk

Files:

• GSASIImath.py (modified) (2 diffs)
• GSASIIphsGUI.py (modified) (1 diff)
• GSASIIpwd.py (modified) (5 diffs)
• imports/G2pdf_gr.py (modified) (3 diffs)

trunk/GSASIImath.py

@@ -732 +732 @@
                 IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR[:,iA],0.))
                 for iU in IndB[0]:
-                    if AtNames[iA] != notName:
+                    if AtNames[iA%len(AtNames)] != notName:
                         unit = Units[iU]
                         if np.any(unit):
@@ -739 +739 @@
                             Topstr = ' +(%4d)'%(Top)
                         if Short:
-                            Neigh.append([AtNames[iA],dist[iU],True])
+                            Neigh.append([AtNames[iA%len(AtNames)],dist[iU],True])
                         else:
                             Neigh.append([AtNames[iA]+Topstr,atTypes[iA],dist[iU],dx[iU]])

trunk/GSASIIphsGUI.py

@@ -5551 +5551 @@
                                     NamesT = []
                                     for item in rdfDict:
+                                        if 'va' in item:
+                                            continue
                                         if 'rdf' not in item and 'g(r)' in title:
                                             Ylab = r'$\mathsf{G(r),\AA^{-1}}$'

trunk/GSASIIpwd.py

@@ -2207 +2207 @@
     return fname
 
-def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):
-    
-    def findDup(Atoms):
-        Dup = []
-        Fracs = []
-        for iat1,at1 in enumerate(Atoms):
-            if any([at1[0] in dup for dup in Dup]):
-                continue
-            else:
-                Dup.append([at1[0],])
-                Fracs.append([at1[6],])
-            for iat2,at2 in enumerate(Atoms[(iat1+1):]):
-                if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
-                    Dup[-1] += [at2[0],]
-                    Fracs[-1]+= [at2[6],]
-        return Dup,Fracs
+def findDup(Atoms):
+    Dup = []
+    Fracs = []
+    for iat1,at1 in enumerate(Atoms):
+        if any([at1[0] in dup for dup in Dup]):
+            continue
+        else:
+            Dup.append([at1[0],])
+            Fracs.append([at1[6],])
+        for iat2,at2 in enumerate(Atoms[(iat1+1):]):
+            if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001:
+                Dup[-1] += [at2[0],]
+                Fracs[-1]+= [at2[6],]
+    return Dup,Fracs
+
+def MakeRMC6f(PWDdata,Name,Phase,RMCPdict):    
     
     Meta = RMCPdict['metadata']
@@ -2753 +2753 @@
     generalData = Phase['General']
     Atseq = RMCPdict['atSeq']
+    Dups,Fracs = findDup(Phase['Atoms'])
+    Sfracs = [np.cumsum(fracs) for fracs in Fracs]
     Supercell = RMCPdict['SuperCell']
     Cell = generalData['Cell'][1:7]
@@ -2761 +2763 @@
     newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=False)
     Atoms = newPhase['Atoms']
-    XYZ = np.array([atom[3:6] for atom in Atoms]).T
+
+    Natm = np.core.defchararray.count(np.array(Atcodes),'+')    #no. atoms in original unit cell
+    Natm = np.count_nonzero(Natm-1)
+    Atoms = newPhase['Atoms']
+    Satoms = G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms,5),4),3)
+    Datoms = [[atom for atom in Satoms if atom[0] in dup] for dup in Dups]
+    Natoms = []
+    for idup,dup in enumerate(Dups):
+        ldup = len(dup)
+        datoms = Datoms[idup]
+        natm = len(datoms)
+        i = 0
+        while i < natm:
+            atoms = datoms[i:i+ldup]
+            try:
+                atom = atoms[np.searchsorted(Sfracs[idup],rand.random())]
+                Natoms.append(atom)
+            except IndexError:      #what about vacancies?
+                if 'Va' not in Atseq:
+                    Atseq.append('Va')
+                    RMCPdict['aTypes']['Va'] = 0.0
+                atom = atoms[0]
+                atom[1] = 'Va'
+                Natoms.append(atom)
+            i += ldup
+
+
+
+
+
+    XYZ = np.array([atom[3:6] for atom in Natoms]).T
     XYZptp = np.array([ma.ptp(XYZ[0]),ma.ptp(XYZ[1]),ma.ptp(XYZ[2])])/2.
     Cell = newPhase['General']['Cell'][1:7]
@@ -2780 +2812 @@
     if RMCPdict['byMolec']:
         NPM = RMCPdict['Natoms']
-        for iat,atom in enumerate(Atoms):
+        for iat,atom in enumerate(Natoms):
             XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian;residue = 'RMC'
             fl.write('ATOM  %5d %-4s RMC%6d%12.3f%8.3f%8.3f  1.00  0.00          %2s\n'%(       
@@ -2787 +2819 @@
     else:
         for atm in Atseq:
-            for iat,atom in enumerate(Atoms):
+            for iat,atom in enumerate(Natoms):
                 if atom[1] == atm:
                     XYZ = np.inner(A,np.array(atom[3:6])-XYZptp)    #shift origin to middle & make Cartesian

trunk/imports/G2pdf_gr.py

@@ -23 +23 @@
 GSASIIpath.SetVersionNumber("$Revision$")
 
+class txt_FSQReaderClass(G2obj.ImportPDFData):
+    'Routines to import S(Q) data from a .fq file'
+    def __init__(self):
+        super(self.__class__,self).__init__( # fancy way to self-reference
+            extensionlist=('.fq','.sq'),
+            strictExtension=False,
+            formatName = 'q (A-1) step S(Q) data',
+            longFormatName = 'q (A-1) stepped S(Q) PDF data from pdfGet or GSAS-II'
+            )
+
+    # Validate the contents -- make sure we only have valid lines
+    def ContentsValidator(self, filename):
+        'Look through the file for expected types of lines in a valid r-step file'
+        filepointer = open(filename,'r')
+        Ndata = 0
+        for i,S in enumerate(filepointer):
+            if '#' in S[:2]:
+                break
+            if len(S.split()) != 2:
+                break
+        for i,S in enumerate(filepointer):            
+            vals = S.split()
+            if len(vals) >= 2:
+                try:
+                    data = [float(val) for val in vals]
+                    Ndata += 1
+                except ValueError:
+                    pass
+        if not Ndata:     
+            self.errors = 'No 2 or more column numeric data found'
+            filepointer.close()
+            return False
+        filepointer.close()
+        return True # no errors encountered
+
+    def Reader(self,filename,ParentFrame=None, **unused):
+        print ('Read a q-step text file')
+        x = []
+        y = []
+        ifData = False
+        filepointer = open(filename,'r')
+        for i,S in enumerate(filepointer):
+            if not ifData:
+                if len(S) == 1:     #skip blank line
+                    continue
+                if len(S.split()) != 2:
+                    continue
+                self.comments.append(S[:-1])
+                if '#' in S[:2]:
+                    ifData = True
+            else:
+                vals = S.split()
+                if len(vals) >= 2:
+                    try:
+                        data = [float(val) for val in vals]
+                        x.append(float(data[0]))
+                        y.append(float(data[1]))
+                    except ValueError:
+                        msg = 'Error in line '+str(i+1)
+                        print (msg)
+                        continue
+        self.pdfdata = np.array([
+            np.array(x), # x-axis values q
+            np.array(y), # pdf f(q))
+            ])
+        self.pdfentry[0] = filename
+        self.pdfentry[2] = 1 # xy file only has one bank
+        self.idstring = ospath.basename(filename)
+
+        return True
+
 class txt_PDFReaderClass(G2obj.ImportPDFData):
     'Routines to import PDF G(R) data from a .gr file'
@@ -39 +110 @@
         Ndata = 0
         for i,S in enumerate(filepointer):
-            if '#L' in S[:2]:
+            if '#' in S[:2]:
                 break
         for i,S in enumerate(filepointer):            
@@ -67 +138 @@
                     continue
                 self.comments.append(S[:-1])
-                if '#L' in S[:2]:
+                if '#' in S[:2]:
                     ifData = True
             else: