diff --git a/src/diffpy/pdfgui/control/calculation.py b/src/diffpy/pdfgui/control/calculation.py index 4417be2b..e81fc8ce 100644 --- a/src/diffpy/pdfgui/control/calculation.py +++ b/src/diffpy/pdfgui/control/calculation.py @@ -25,6 +25,9 @@ from diffpy.pdfgui.control.pdfcomponent import PDFComponent from diffpy.pdfgui.utils import safeCPickleDumps, pickle_loads +from diffpy.srreal.pdfcalculator import PDFCalculator, DebyePDFCalculator +from diffpy.srreal.structureadapter import nometa + class Calculation(PDFComponent): """Perform a theoretical computation of PDF from model structure. @@ -39,8 +42,11 @@ class Calculation(PDFComponent): rcalc -- list of r values, this is set after calculation is finished Gcalc -- list of calculated G values stype -- scattering type, 'X' or 'N' + pctype -- PDF calculator type, 'PC' or 'DPC', to use RealSpacePC + or DebyePC. qmax -- maximum value of Q in inverse Angstroms. Termination ripples are ignored for qmax=0. + qmin -- minimum value of Q in inverse Angstroms. Default qmin=0. qdamp -- specifies width of Gaussian damping factor in pdf_obs due to imperfect Q resolution qbroad -- quadratic peak broadening factor related to dataset @@ -63,8 +69,10 @@ def __init__(self, name): self.rcalc = [] self.Gcalc = [] self.stype = 'X' + self.pctype = 'PC' # default use PDF real space calculator # user must specify qmax to get termination ripples self.qmax = 0.0 + self.qmin = 0.0 self.qdamp = 0.001 self.qbroad = 0.0 self.spdiameter = None @@ -135,7 +143,81 @@ def start(self): return def calculate(self): - """do the real calculation + """ + Perform calculation in diffpy.srreal + """ + # clean up old results + self.rcalc = [] + self.Gcalc = [] + + # do the job + if len(self.owner.strucs) == 0: + raise ControlConfigError("No structure is given for calculation") + # make sure parameters are initialized + self.owner.updateParameters() + + # initialize the calculator + if self.pctype == 'PC': # use PDFCalculator + pc = PDFCalculator() + elif self.pctype == 'DPC': # use DebyePDFCalculator + pc = DebyePDFCalculator() + # x-ray or neutron, PC default x-ray + if self.stype == 'N': + pc.scatteringfactortable = "neutron" + + pc.qmax = self.qmax + pc.qmin = self.qmin + pc.qdamp = self.qdamp + pc.qbroad = self.qbroad + pc.rmin = self.rmin + pc.rmax = self.rmax + pc.rstep = self.rstep + pc.scale = self.dscale + + # load structure and disable metadata using the nometa function + # and set any calculator attributes as needed as above + r_list = [] + g_list = [] + self.owner.applyParameters() + for struc in self.owner.strucs: + # pc is for one calculation. the common setting + # pc_temp is for each phase, specific setting + pc_temp = pc.copy() + pc_temp.delta1 = struc.getvar('delta1') + pc_temp.delta2 = struc.getvar('delta2') + if struc.getvar('spdiameter'): + pc_temp.addEnvelope('sphericalshape') + pc_temp.spdiameter = struc.getvar('spdiameter') + if struc.getvar('stepcut'): + pc_temp.addEnvelope('stepcut') + pc_temp.stepcut = struc.getvar('stepcut') + + + ##TODO: pair mask. the baseline is not correct with PDFFIT. + struc.applyCMIPairSelection(pc_temp) + # pc_temp.setTypeMask("Ni","O",True) + + r, g = pc_temp(nometa(struc)) + if struc.getvar('pscale'): + g = g * struc.getvar('pscale') + r_list.append(r) + g_list.append(g) + print("len(r_list)", len(r_list)) + print("len(g_list)", len(g_list)) + + # get results + + self.rcalc = r_list[0].tolist() # r0, r1, r2 are the same, so just use r0 + # sum up multi-phase PDFs + gsum = 0 + for i in range(len(self.owner.strucs)): + gsum += g_list[i] + self.Gcalc = gsum.tolist() + + return + + def calculate_pdffit(self): + """do the real calculation in PDFFIT """ # clean up old results self.rcalc = [] @@ -217,6 +299,11 @@ def writeStr(self): lines.append('stype=X x-ray scattering') elif self.stype == 'N': lines.append('stype=N neutron scattering') + # pctype + if self.pctype == 'PC': + lines.append('pctype=PC real space PDF calculator') + elif self.pctype == 'DPC': + lines.append('pctype=DPC Debye PDF calculator') # dscale if self.dscale: lines.append('dscale=%g' % self.dscale) @@ -226,6 +313,12 @@ def writeStr(self): else: qmax_line = 'qmax=%.2f' % self.qmax lines.append(qmax_line) + # qmin + if self.qmin == 0: + qmin_line = 'qmin=0' + else: + qmin_line = 'qmin=%.2f' % self.qmin + lines.append(qmin_line) # qdamp if isinstance(self.qdamp, float): lines.append('qdamp=%g' % self.qdamp) @@ -257,7 +350,9 @@ def load(self, z, subpath): self.rcalc = config['rcalc'] self.Gcalc = config['Gcalc'] self.stype = config['stype'] + self.pctype = config['pctype'] self.qmax = config['qmax'] + self.qmin = config['qmin'] self.qdamp = config.get('qdamp', config.get('qsig')) self.qbroad = config.get('qbroad', config.get('qalp', 0.0)) self.spdiameter = config.get('spdiameter') @@ -278,7 +373,9 @@ def save(self, z, subpath): 'rcalc' : self.rcalc, 'Gcalc' : self.Gcalc, 'stype' : self.stype, + 'pctype' : self.pctype, 'qmax' : self.qmax, + 'qmin' : self.qmin, 'qdamp' : self.qdamp, 'qbroad' : self.qbroad, 'dscale' : self.dscale, @@ -299,7 +396,7 @@ def copy(self, other=None): # rcalc and Gcalc may be assigned, they get replaced by new lists # after every calculation assign_attributes = ( 'rmin', 'rstep', 'rmax', 'rlen', - 'rcalc', 'Gcalc', 'stype', 'qmax', 'qdamp', + 'rcalc', 'Gcalc', 'stype', 'pctype', 'qmax', 'qmin', 'qdamp', 'qbroad', 'dscale', ) copy_attributes = ( ) for a in assign_attributes: diff --git a/src/diffpy/pdfgui/control/fitstructure.py b/src/diffpy/pdfgui/control/fitstructure.py index 0dc52c8f..9311a648 100644 --- a/src/diffpy/pdfgui/control/fitstructure.py +++ b/src/diffpy/pdfgui/control/fitstructure.py @@ -621,6 +621,25 @@ def applyPairSelection(self, server, phaseidx): server.selectAtomIndex(phaseidx, 'j', idx, jflag) return + def applyCMIPairSelection(self, pc): + """Apply pair selection for CMI calculations of partial PDF. + + pc -- instance of CMI PDFCalculator engine + """ + psf = self.getPairSelectionFlags() + # to begin masking, first remove ALL + pc.setTypeMask("all", "all", False) + #option1: mask by site index. has problem + # for i, iflag in enumerate(psf['firstflags']): + # for j, jflag in enumerate(psf['secondflags']): + # print("i, iflag", i, iflag) + # print("j, jflag", j, jflag) + # pc.setPairMask(i, j, (iflag and jflag)) + #option2: mask by element type + firstelement = str(psf['fixed_pair_string'].split("-")[0]) + secondelement = str(psf['fixed_pair_string'].split("-")[1]) + pc.setTypeMask(firstelement, secondelement, True) + return def getSelectedIndices(self, s): '''Indices of the atoms that match the specified selection string. diff --git a/src/diffpy/pdfgui/control/fitting.py b/src/diffpy/pdfgui/control/fitting.py index 7e42b1dc..7b400855 100644 --- a/src/diffpy/pdfgui/control/fitting.py +++ b/src/diffpy/pdfgui/control/fitting.py @@ -24,6 +24,15 @@ from diffpy.pdfgui.control.controlerrors import ControlValueError from diffpy.pdfgui.utils import safeCPickleDumps, pickle_loads +from diffpy.srfit.pdf import PDFContribution +from diffpy.srfit.fitbase import Profile +from diffpy.srreal.structureadapter import nometa +from diffpy.srfit.fitbase import FitContribution, FitRecipe, FitResults +from diffpy.srfit.pdf import PDFParser +from diffpy.srfit.pdf import PDFGenerator, DebyePDFGenerator +from diffpy.srfit.fitbase import FitRecipe, FitResults +from scipy.optimize.minpack import leastsq + # helper routines to deal with PDFfit2 exceptions def getEngineExceptions(): @@ -124,6 +133,13 @@ def __init__(self, name): # the PDFfit2 server instance. self.server = None + # the CMI server instance. + self.cmipdfgen = None + self.cmiprofile = None + self.cmicontribution = None + self.cmirecipe = None + self.cmiresults = None + # public data members self.step = 0 self.parameters = {} @@ -320,6 +336,9 @@ def getServer(self): # create a new instance of calculation server from diffpy.pdffit2 import PdfFit self.server = PdfFit() + + self.cmiprofile = Profile() + self.__changeStatus(fitStatus=Fitting.CONNECTED) @@ -330,10 +349,90 @@ def configure(self): # make sure parameters are initialized self.updateParameters() + + #long CMI part + self.applyParameters() + + if self.datasets[0].pctype == 'PC': + self.cmipdfgen = PDFGenerator("cmipdfgen") + elif self.datasets[0].pctype == 'DPC': + self.cmipdfgen = DebyePDFGenerator("cmipdfgen") + + self.cmipdfgen.setStructure(self.strucs[0]) + parser = PDFParser() + parser.parseString(self.datasets[0].writeResampledObsStr()) + self.cmiprofile.loadParsedData(parser) + self.cmiprofile.setCalculationRange(xmin = self.datasets[0].fitrmin, + xmax = self.datasets[0].fitrmax, + dx = self.datasets[0].fitrstep) + + self.cmicontribution = FitContribution("cmicontribution") + self.cmicontribution.addProfileGenerator(self.cmipdfgen) + self.cmicontribution.setProfile(self.cmiprofile, xname ="r") + self.cmicontribution.setEquation("scale * cmipdfgen") + + # add qmax, qdamp, qbroad into cmipdfgen + self.cmipdfgen.setQmax(self.datasets[0].qmax) + self.cmipdfgen.qdamp.value = self.datasets[0].qdamp + self.cmipdfgen.qbroad.value = self.datasets[0].qbroad + + self.cmirecipe = FitRecipe() + self.cmirecipe.addContribution(self.cmicontribution) + # self.cmirecipe.addVar(self.cmicontribution.scale, 1.0) + + for index, par in self.parameters.items(): + # clean any refined value + par.refined = None + # self.server.setpar(index, par.initialValue()) # info[0] = init value + var_name = "var" + str(index) + print("var_name") + print(var_name) + print("par.initialValue()") + print(par.initialValue()) + self.cmirecipe.newVar(var_name, par.initialValue()) + # fix if fixed. Note: all parameters are free after self.server.reset(). + # if par.fixed: + # self.server.fixpar(index) + + # phase constrains + for struc in self.strucs: + struc.clearRefined() + for key, var in struc.constraints.items(): + self.cmiConstrain(key, var) + + # data constrains + for dataset in self.datasets: + dataset.clearRefined() + for key, var in dataset.constraints.items(): + self.cmiConstrain(key, var) + # Removed call to pdfrange call, because data were already + # resampled to at fit range. + # + # Pair selection applies only to the current dataset, + # therefore it has to be done here. + # TODO the following comment lines + # nstrucs = len(self.strucs) + # for phaseidx, struc in zip(range(1, nstrucs + 1), self.strucs): + # print(phaseidx) + # print(struc) + # struc.applyCMIPairSelection() + # struc.applyPairSelection(self.server, phaseidx) + + + # turn on printout fithook in each refinement step + self.cmirecipe.fithooks[0].verbose = 3 + + leastsq(self.cmirecipe.residual, self.cmirecipe.values) + self.cmiresults = "\n=============================== CMI RESULTS ==================================\n" + self.cmiresults += str(FitResults(self.cmirecipe)) + self.cmiresults += "============================ END OF CMI RESULTS ==============================\n\n" + + #Long end CMI part + self.server.reset() for struc in self.strucs: struc.clearRefined() - self.server.read_struct_string(struc.initial.writeStr("pdffit") ) + self.server.read_struct_string(struc.initial.writeStr("pdffit")) for key, var in struc.constraints.items(): self.server.constrain(key, var.formula) @@ -376,6 +475,14 @@ def resetStatus(self): """reset status back to initialized""" self.snapshots = [] self.step = 0 + # long + # initialize cmi + self.cmipdfgen = None + self.cmiprofile = None + self.cmicontribution = None + self.cmirecipe = None + self.cmiresults = None + # end long if self.fitStatus == Fitting.INITIALIZED: return # already reset @@ -809,4 +916,118 @@ def _getData(self, id, name, step=-1): else: return self.snapshots[step][index] + # Long new helper function + def cmiConstrain(self, key, var): + """Constrain structure parameters into cmi receipe. + key -- names of parameters, like pscale, lat(n). + var -- var.formula represents names of constrains, like @1, @2 + 1. + """ + key_ref, key_arg = self.__getRef(key) + var_name = self.transVar(var.formula) + + lat = self.cmipdfgen.phase.getLattice() + atoms = self.cmipdfgen.phase.getScatterers() + + if key_ref == 'pscale': + self.cmirecipe.constrain(self.cmicontribution.scale, var_name) + if key_ref == 'lat': + if key_arg == '1': + self.cmirecipe.constrain(lat.a, var_name) + if key_arg == '2': + self.cmirecipe.constrain(lat.b, var_name) + if key_arg == '3': + self.cmirecipe.constrain(lat.c, var_name) + if key_arg == '4': + self.cmirecipe.constrain(lat.alpha, var_name) + if key_arg == '5': + self.cmirecipe.constrain(lat.beta, var_name) + if key_arg == '6': + self.cmirecipe.constrain(lat.gamma, var_name) + + # delta term + if key_ref == 'delta1': + self.cmirecipe.constrain(self.cmipdfgen.delta1, var_name) + if key_ref == 'delta2': + self.cmirecipe.constrain(self.cmipdfgen.delta2, var_name) + + # ADP + ## TODO key_ascii == 'u11(i)', constrain the ith atom's ADP U11. + if key_ref == 'u11': + self.cmirecipe.constrain(atoms[key_arg - 1].U11, var_name) + if key_ref == 'u22': + self.cmirecipe.constrain(atoms[key_arg - 1].U22, var_name) + if key_ref == 'u33': + self.cmirecipe.constrain(atoms[key_arg - 1].U33, var_name) + if key_ref == 'u12': + self.cmirecipe.constrain(atoms[key_arg - 1].U12, var_name) + if key_ref == 'u13': + self.cmirecipe.constrain(atoms[key_arg - 1].U13, var_name) + if key_ref == 'u23': + self.cmirecipe.constrain(atoms[key_arg - 1].U23, var_name) + + # atom positions + if key_ref == 'x': + self.cmirecipe.constrain(atoms[key_arg - 1].x, var_name) + if key_ref == 'y': + self.cmirecipe.constrain(atoms[key_arg - 1].y, var_name) + if key_ref == 'z': + self.cmirecipe.constrain(atoms[key_arg - 1].z, var_name) + + # occupancy + if key_ref == 'occ': + self.cmirecipe.constrain(atoms[key_arg - 1].occupancy, var_name) + + # data parameters + if key_ref == 'qdamp': + self.cmirecipe.constrain(self.cmipdfgen.qdamp, var_name) + if key_ref == 'qbroad': + self.cmirecipe.constrain(self.cmipdfgen.qbroad, var_name) + # TODO how to deal with `dscale`. cmipdfgen don't have `dscale` parameter. + + return + + def transVar(self, str): + # input "@11" + # output "var11" + return str.replace("@", "var") + + def __getRef(self, var_string): + # copy from __getRef in pdffit.py from PDFfit2 package. + # input pscale, output method_string = "pscale", arg_int = None + # input lat(1), output method_string = "lat", arg_int = 1 + # input u11(2), output method_string = "u11", arg_int = 2 + """Return the actual reference to the variable in the var_string. + + This function must be called before trying to actually reference an + internal variable. See the constrain method for an example. + + Raises: + pdffit2.unassignedError if variable is not yet assigned + ValueError if variable index does not exist (e.g. lat(7)) + """ + var_string = _convertCallable(var_string) + arg_int = None + try: + method_string, arg_string = var_string.split("(") + method_string = method_string.strip() + arg_int = int(arg_string.strip(")").strip()) + except ValueError: #There is no arg_string + method_string = var_string.strip() + return method_string, arg_int + + +# Long helper routines +def _convertCallable(var): + """Convert an object to the result of its call when callable. + + var -- string or callable object that returns string + + Return var or var(). + """ + if callable(var): + rv = var() + else: + rv = var + return rv + # End of file diff --git a/src/diffpy/pdfgui/control/pdfdataset.py b/src/diffpy/pdfgui/control/pdfdataset.py index 48fb53e0..e71ddfed 100644 --- a/src/diffpy/pdfgui/control/pdfdataset.py +++ b/src/diffpy/pdfgui/control/pdfdataset.py @@ -35,8 +35,11 @@ class PDFDataSet(PDFComponent): drobs -- list of standard deviations of robs dGobs -- list of standard deviations of Gobs stype -- scattering type, 'X' or 'N' + pctype -- PDF calculator type, 'PC' or 'DPC', to use RealSpacePC + or DebyePC. qmax -- maximum value of Q in inverse Angstroms. Termination ripples are neglected for qmax=0. + qmin -- minimum value of Q in inverse Angstroms. Default qmin=0. qdamp -- specifies width of Gaussian damping factor in pdf_obs due to imperfect Q resolution qbroad -- quadratic peak broadening factor related to dataset @@ -56,7 +59,7 @@ class PDFDataSet(PDFComponent): refinableVars -- set (dict) of refinable variable names. """ - persistentItems = [ 'robs', 'Gobs', 'drobs', 'dGobs', 'stype', 'qmax', + persistentItems = [ 'robs', 'Gobs', 'drobs', 'dGobs', 'stype', 'pctype', 'qmax', 'qmin', 'qdamp', 'qbroad', 'dscale', 'rmin', 'rmax', 'metadata' ] refinableVars = dict.fromkeys(('qdamp', 'qbroad', 'dscale')) @@ -77,6 +80,8 @@ def clear(self): self.dGobs = [] self.stype = 'X' # user must specify qmax to get termination ripples + self.pctype = 'PC' + self.qmin = 0.0 self.qmax = 0.0 self.qdamp = 0.001 self.qbroad = 0.0 @@ -216,7 +221,7 @@ def readStr(self, datastring): if res: self.metadata['doping'] = float(res.groups()[0]) - # parsing gerneral metadata + # parsing general metadata if metadata: regexp = r"\b(\w+)\ *=\ *(%(f)s)\b" % rx while True: diff --git a/src/diffpy/pdfgui/control/pdfguicontrol.py b/src/diffpy/pdfgui/control/pdfguicontrol.py index 5f838a5f..012724da 100644 --- a/src/diffpy/pdfgui/control/pdfguicontrol.py +++ b/src/diffpy/pdfgui/control/pdfguicontrol.py @@ -529,6 +529,20 @@ def getEngineOutput(self): redirect_stdout(six.StringIO()) return txt + #long + def getCMIOutput(self): + """Get the output from the CMI engine.""" + txt = None + if self.currentFitting: + if self.currentFitting.cmiresults: + txt = self.currentFitting.cmiresults + return txt + + def resetCMIOutput(self): + self.currentFitting.cmiresults = None + return + #end Long + _pdfguicontrol = None def pdfguicontrol(*args, **kwargs): """This function will return the single instance of class PDFGuiControl""" diff --git a/src/diffpy/pdfgui/gui/mainframe.py b/src/diffpy/pdfgui/gui/mainframe.py index 24a0b578..7f338c44 100644 --- a/src/diffpy/pdfgui/gui/mainframe.py +++ b/src/diffpy/pdfgui/gui/mainframe.py @@ -352,7 +352,7 @@ def __customProperties(self): # Position other panels. Note that currently MinimizeButton does not do # anything. It is to be implemented in future versions of wx.aui self.auiManager.AddPane(self.outputPanel, wx.aui.AuiPaneInfo(). - Name("outputPanel").Caption("PDFfit2 Output"). + Name("outputPanel").Caption("DiffPy-CMI Output"). Bottom(). TopDockable(). BottomDockable(). @@ -2497,7 +2497,13 @@ def updateFittingStatus(self, job): def updateOutput(self): """Update text in outputPanel with text in stdout.""" - self.outputPanel.updateText(self.control.getEngineOutput()) + # self.outputPanel.updateText(self.control.getEngineOutput()) + #long + # TODO: append CMI result after pdffit result + if self.control.getCMIOutput(): + self.outputPanel.updateText(self.control.getCMIOutput()) + self.control.resetCMIOutput() #only output cmi results once + #end long return # end of class MainPanel diff --git a/src/diffpy/pdfgui/gui/pdfguiglobals.py b/src/diffpy/pdfgui/gui/pdfguiglobals.py index 406ed5ef..61948a22 100644 --- a/src/diffpy/pdfgui/gui/pdfguiglobals.py +++ b/src/diffpy/pdfgui/gui/pdfguiglobals.py @@ -13,13 +13,13 @@ # ############################################################################## -"""This module contains gloabal parameters needed by PDFgui.""" +"""This module contains global parameters needed by PDFgui.""" import os.path from pkg_resources import Requirement, resource_filename # Name of the program -name = "PDFgui" +name = "PDFgui 2.0" # Maximum number of files to be remembered MAXMRU = 5 # The location of the configuration file