Added solid phase correction

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):

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]
     }
 }
 

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)
+

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))