diff --git a/xpsservice/models.py b/xpsservice/models.py index bbfd6d3..188e45a 100644 --- a/xpsservice/models.py +++ b/xpsservice/models.py @@ -6,7 +6,7 @@ from ase import Atoms from pydantic import BaseModel, Field, validator from rdkit import Chem -from .settings import ALLOWED_FMAX, ALLOWED_ELEMENTS, ALLOWED_FF, ALLOWED_METHODS, transition_map +from .settings import ALLOWED_FMAX, ALLOWED_ELEMENTS, ALLOWED_FF, ALLOWED_METHODS, ALLOWED_ENERGY_REFERENCE, transition_map class TransitionValidator(BaseModel): @@ -33,9 +33,13 @@ class XPSRequest(BaseModel): description="Method for geometry optimization and XPS binding energies calculation. Allowed values: `GFNFF`, `GFN2xTB`, `GFN1xTB`.", ) fmax: Optional[float] = Field( - 0.01, + 0.001, description="Maximum force admissible during geometry optimization. Typically ranging from 0.1 to 0.0001." ) + energyReference: Optional[str] = Field( + "solid", + description="Energy referencing performed through polynomial regression, for gaz phase, solid state or none (raw data)" + ) @validator("smiles") def validate_smiles(cls, v): @@ -70,6 +74,12 @@ def validate_fmax(cls, v): if not (ALLOWED_FMAX[0] <= v <= ALLOWED_FMAX[1]): raise ValueError(f"fmax must be within the range {ALLOWED_FMAX}") return v + + @validator("energyReference") + def validate_energy_reference(cls, v): + if v not in ALLOWED_ENERGY_REFERENCE: + raise ValueError(f"Method must be in {ALLOWED_ENERGY_REFERENCE}") + return v class Position(BaseModel): diff --git a/xpsservice/settings.py b/xpsservice/settings.py index 85b7453..bd4e44b 100644 --- a/xpsservice/settings.py +++ b/xpsservice/settings.py @@ -13,7 +13,8 @@ #XTB optimization method ALLOWED_METHODS = ("GFNFF", "GFN2xTB", "GFN1xTB") MAX_ATOMS_XTB = int(os.getenv("MAX_ATOMS_XTB", 200)) -ALLOWED_FMAX = (0.000001, 1) +ALLOWED_FMAX = (0.000001, 0.1) +ALLOWED_ENERGY_REFERENCE = ("none", "solid") #timeout for the overall calculation @@ -22,18 +23,21 @@ # Define transition_map dictionary(list of orbitals for which a photoelectron emission is calculated) # Several maps possible. Adjust which one to load in xpsservice +# The regression coefficient are for polynomial a + bx + cx^2 + ..., given in an array [a,b,c,...] transition_map = { "C1s": { "element": "C", "orbital": "1s", "soap_config_filepath": os.path.abspath("SOAP_configs/soap_config_C1s.pkl"), - "model_filepath": os.path.abspath("ML_models/XPS_GPR_C1s_xtb.pkl") + "model_filepath": os.path.abspath("ML_models/XPS_GPR_C1s_xtb.pkl"), + "regression_coefficients": [-2501.002792622873, 18.12807741, -0.02938832] }, "O1s": { "element": "O", "orbital": "1s", "soap_config_filepath": os.path.abspath("SOAP_configs/soap_config_O1s.pkl"), - "model_filepath": os.path.abspath("ML_models/XPS_GPR_O1s_xtb.pkl") + "model_filepath": os.path.abspath("ML_models/XPS_GPR_O1s_xtb.pkl"), + "regression_coefficients": [-14914.578039681715, 5.67166123e+01, -5.20506217e-02] } } diff --git a/xpsservice/xps.py b/xpsservice/xps.py index 8a71164..d0c9845 100644 --- a/xpsservice/xps.py +++ b/xpsservice/xps.py @@ -7,6 +7,7 @@ import logging import pickle import hashlib +from fastapi.logger import logger from ase import Atoms from ase.build.tools import sort from quippy.descriptors import Descriptor @@ -170,7 +171,7 @@ def run_xps_calculations(ase_mol: Atoms, transition_map) -> dict: return be_predictions @wrapt_timeout_decorator.timeout(TIMEOUT, use_signals=False) -def calculate_from_molfile(molfile, method, fmax, transition_map) -> XPSResult: +def calculate_from_molfile(molfile, method, fmax, energy_reference, transition_map) -> XPSResult: ase_mol, mol = molfile2ase(molfile, get_max_atoms(method)) @@ -181,6 +182,11 @@ def calculate_from_molfile(molfile, method, fmax, transition_map) -> XPSResult: be_predictions = run_xps_calculations(opt_ase_mol, transition_map) + #reference the binding enregy if queried + if energy_reference == 'solid': + referenced_be_predictions = reference_be_predictions(be_predictions, transition_map) + be_predictions = referenced_be_predictions + if isinstance(be_predictions, dict): ordered_predictions = reorder_predictions(be_predictions, opt_ase_mol, transition_map) else: @@ -218,3 +224,48 @@ def reorder_predictions(be_predictions: dict, ase_mol: Atoms, transition_map: di return ordered_predictions +# Energy referencing of the whole predictions +def reference_be_predictions(be_predictions: dict, transition_map: dict): + + referenced_be_predictions = {} + + for transition_key, predictions in be_predictions.items(): + # get regression coeffs + transition_info = transition_map[transition_key] + regression_coefficients = transition_info['regression_coefficients'] + + referenced_predictions = [] + for prediction in predictions: + referenced_predictions.append(reference_be_prediction(prediction, regression_coefficients)) + + referenced_be_predictions[transition_key] = referenced_predictions + + return referenced_be_predictions + +''' +# Energy referencing of a single prediction based on polynomial fitting of experimental data +def reference_be_prediction(be_prediction: float, regression_coefficients: list): + + referenced_be_prediction = be_prediction + for i, coeff in enumerate(regression_coefficients): + logger.debug(f"be_prediction format {be_prediction}") + referenced_be_prediction[0] += coeff * (be_prediction[0] ** i) + + return referenced_be_prediction +''' + +# Energy referencing of a single prediction based on polynomial fitting of experimental data +def reference_be_prediction(be_prediction: tuple, regression_coefficients: list): + # Unpack the tuple + prediction, std_dev = be_prediction + + # Initialize the referenced prediction with the original prediction value + referenced_prediction = 0.0 + + # Apply the polynomial correction + for i, coeff in enumerate(regression_coefficients): + referenced_prediction += coeff * (prediction ** i) + + # Return the modified tuple with the referenced prediction and original standard deviation + return (referenced_prediction, std_dev) + diff --git a/xpsservice/xpsservice.py b/xpsservice/xpsservice.py index a930b27..610b84f 100644 --- a/xpsservice/xpsservice.py +++ b/xpsservice/xpsservice.py @@ -44,7 +44,6 @@ #Initial loading #load_soap_configs_and_models(selected_transition_map) - app = FastAPI( title="EPFL-ISIC-XRDSAP: XPS webservice", description="Offers XPS binding energy prediction tool, based on simulation, and a Gaussian process ML model. Allowed elements/orbitals to be predicted are `C1s` and `O1s`. Hydrogens are taken into account but not predicted. Allowed methods for molecular geometry optimization are `GFNFF`, `GFN2xTB`, `GFN1xTB`", @@ -54,12 +53,6 @@ ) -#ping -@app.get("/ping") -def ping(): - return {"message": "pongping"} - - #load models and descriptors @app.get("/load_soap_configs_and_models") async def load_soap_configs_and_models_endpoint() -> bool: @@ -104,6 +97,7 @@ def predict_binding_energies_endpoint(request: XPSRequest): molfile = request.molFile method = request.method fmax = request.fmax + energy_reference = request.energyReference # Perform conversions to molfile based on the provided input if smiles and not molfile: @@ -119,7 +113,7 @@ def predict_binding_energies_endpoint(request: XPSRequest): print("converted format") # Perform calculations try: - predictions = calculate_from_molfile(molfile, method, fmax, selected_transition_map) + predictions = calculate_from_molfile(molfile, method, fmax, energy_reference, selected_transition_map) except Exception as e: raise HTTPException(status_code=400, detail=str(e))