Changeset 4952 for trunk/GSASIIpwd.py

Timestamp:

Jun 11, 2021 3:43:57 PM (2 years ago)

Author:

toby

Message:

switch to native fullrmc input modes rather than convert files

File:

• trunk/GSASIIpwd.py (modified) (5 diffs)

trunk/GSASIIpwd.py

@@ -3017 +3017 @@
 # otherwise build it 
 ENGINE = Engine(path=None)
+Ebc = ENGINE.boundaryConditions
 if not ENGINE.is_engine(engineFileName) or FRESH_START:
     ## create structure
@@ -3024 +3025 @@
                                  unitcellBC = cell,
                                  supercell  = supercell)
-    rmax = min( [ENGINE.boundaryConditions.get_a(), ENGINE.boundaryConditions.get_b(), ENGINE.boundaryConditions.get_c()] ) /2.
+    rmax = min( [Ebc.get_a(), Ebc.get_b(), Ebc.get_c()] ) /2.
 '''    
     import atmdata
-    rundata += '# conversion factors (may be needed)\n'
-    rundata += '    sumCiBi2 = 0.\n'
-    for elem in Phase['General']['AtomTypes']:
-        rundata += '    Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
-        rundata += '    sumCiBi2 += (Ci*{})**2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
+    # rundata += '    # conversion factors (may be needed)\n'
+    # rundata += '    sumCiBi2 = 0.\n'
+    # for elem in Phase['General']['AtomTypes']:
+    #     rundata += '    Ci = ENGINE.numberOfAtomsPerElement["{}"]/len(ENGINE.allElements)\n'.format(elem)
+    #     rundata += '    sumCiBi2 += (Ci*{})**2\n'.format(atmdata.AtmBlens[elem+'_']['SL'][0])
     rundata += '    rho0 = len(ENGINE.allNames)/ENGINE.volume\n'
     # settings that require a new Engine
@@ -3043 +3044 @@
             rundata += '    GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
             if filDat[3] == 0:
-                rundata += '''    # read and xform G(r) as defined in RMCProfile
+                #rundata += '''    # read and xform G(r) as defined in RMCProfile
     # see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n'''
-                rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
-                rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
+                #rundata += '    GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n'
+                #rundata += '    GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
+                rundata += '    # G(r) as defined in RMCProfile\n'
+                rundata += '    GofR = fullrmc.Constraints.RadialDistributionConstraints.RadialDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt
             elif filDat[3] == 1:
                 rundata += '    # This is G(r) as defined in PDFFIT\n'
@@ -3060 +3063 @@
             rundata += '    SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0]
             if filDat[3] == 0:
-                rundata += '    # Read & xform F(Q) as defined in RMCProfile to S(Q)-1\n'
-                rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
+                rundata += '    # F(Q) as defined in RMCProfile\n'
+                #rundata += '    SOQ[1] *= 1 / sumCiBi2\n'
+                if filDat[4]:
+                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=rmax)\n'
+                rundata += '    SofQ = fullrmc.Constraints.StructureFactorConstraints.NormalizedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
             elif filDat[3] == 1:
-                rundata += '    # This is S(Q) as defined in PDFFIT\n'
+                rundata += '    # S(Q) as defined in PDFFIT\n'
                 rundata += '    SOQ[1] -= 1\n'
-            if filDat[4]:
-                rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=rmax)\n'
-            rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
+                if filDat[4]:
+                    rundata += '    SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=rmax)\n'
+                rundata += '    SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt
+            else:
+                raise ValueError('Invalid S(Q) type: '+str(filDat[3]))
             rundata += '    ENGINE.add_constraints([SofQ])\n'
         else:
@@ -3096 +3104 @@
 else:
     ENGINE = ENGINE.load(path=engineFileName)
-    rmax = min( [ENGINE.boundaryConditions.get_a(), ENGINE.boundaryConditions.get_b(), ENGINE.boundaryConditions.get_c()] ) /2.
+    rmax = min( [Ebc.get_a(), Ebc.get_b(), Ebc.get_c()] ) /2.
 
 ENGINE.set_log_file(os.path.join(dirName,prefix))