Add density cutoff option to fix FDBM for all-electron densities.

doc/sphinx/source/tutorials/afmulator-tutorial.md

@@ -222,6 +222,7 @@ The full-density-based model (FDBM) can provide a better approximation of the Pa
 In order to use the FDBM, more input files and parameters are required:
 ```python
 from ppafm.ocl.AFMulator import AFMulator, HartreePotential, ElectronDensity, TipDensity
+
 pot, xyzs, Zs = HartreePotential.from_file("LOCPOT.xsf", scale=-1) # Sample Hartree potential
 rho_sample, _, _ = ElectronDensity.from_file("CHGCAR.xsf") # Sample electron density
 rho_tip, xyzs_tip, Zs_tip = TipDensity.from_file("density_CO.xsf") # Tip electron density
@@ -242,11 +243,22 @@ Above, we first load the input files.
 In addition to the Hartree potential we also load the sample electron density `rho_sample`, and two densities for the tip, the full electron density `rho_tip`, and a delta density `rho_tip_delta`.
 When constructing the `AFMulator`, we use both of the densities.
 In this case the `rho=rho_tip` is used for the Pauli density overlap intergral, and the `rho_delta=rho_tip_delta` is used in the electrostatics calculation as in the previous section.
+
+````{warning}
+All-electron densities from codes such as FHI-aims can sometimes produce severe artifacts when used with the FDBM.
+A simple work around is to simply cutoff extremely high values of the electron density around the nuclei, which can be done using the {meth}`.clamp` method for the electron density grid:
+```python
+rho_sample = rho_sample.clamp(maximum=100)
+```
+A value of 100 for the cutoff is usually safe to use, although sometimes when both the prefactor ``'A_pauli'`` and exponent ``'B_pauli'`` are large, a lower cutoff may be required.
+````
+
 The delta density can be separately calculated and loaded from a file as above, or it can be estimated by subtracting the valence electrons from the full electron density:
 ```python
 rho_tip_delta = rho_tip.subCores(xyzs_tip, Zs_tip, Rcore=0.7)
 ```
 
+
 We also specify the two parameters of the Pauli overlap integral, the prefactor `A_pauli` and the exponent `B_pauli`, as well as the DFT-D3 parameters based on the functional name via `d3_params`.
 In this case we suppose that the electron densities were calculated using the `PBE` functional, so we use the corresponding parameters.
 Other predefined parameter sets can be found in the description of the method {meth}`.add_dftd3`.
@@ -257,7 +269,3 @@ afm_images = afmulator(xyzs, Zs, qs=pot, rho_sample=rho_sample)
 ```
 
 For examples of using the FDBM, see the [pyridine example](https://github.com/Probe-Particle/ppafm/blob/main/examples/pyridineDensOverlap/run_gpu.py), as well as another [more advanced example](https://github.com/Probe-Particle/ppafm/blob/main/examples/paper_figure/run_simulation.py) where multiple simulations are performed using all of the different force field models described here.
-
-```{warning}
-Electron densities from FHI-aims are known to not work well with the FDBM. This simulation model is generally less explored compared to the other ones, so expect more problems. Known working configurations have used densities calculated with VASP.
-```

examples/BrClPyridine_FDBM_aims_density/params.ini

@@ -0,0 +1,8 @@
+probeType   O
+r0Probe     0.00 0.00 3.00
+krad        30.00
+PBC         True
+nPBC        1       1       0
+scanStep    0.10    0.10    0.10
+scanMin      2.0     2.0     6.0
+scanMax     18.0    18.0    10.0

examples/BrClPyridine_FDBM_aims_density/run.sh

@@ -0,0 +1,33 @@
+#! /bin/bash
+
+# An example of using the full-density based model (FDBM) with an sample electron density obtained with FHI-aims.
+# The elecron density in aims is all-electron, so the density at the nuclei positions can be extremely high.
+# This sometimes causes artifacts to appear in the resuting simulation images. This can be mitigated by cutting
+# out the high electron density values, which is demonstrated in this example.
+
+# Download input data
+tip_dir="./tip"
+sample_dir="./sample"
+if [ ! -d "$tip_dir" ]; then
+    wget https://zenodo.org/records/10563098/files/CO_tip_densities.tar.gz
+    mkdir $tip_dir
+    tar -xvf CO_tip_densities.tar.gz --directory $tip_dir
+fi
+if [ ! -d "$sample_dir" ]; then
+    wget "https://zenodo.org/records/14222456/files/pyridineBrCl.zip?download=1" -O hartree.zip
+    mkdir $sample_dir
+    unzip hartree.zip -d $sample_dir
+fi
+
+# Generate the force field
+# Notice here the `--density_cutoff` option. Try removing it and observe how artifacts appear in the image at the top-right
+# at the position of the bromine atom.
+ppafm-conv-rho       -s sample/density.xsf -t tip/density_CO.xsf -B 1.1 --density_cutoff 100
+ppafm-generate-elff  -i sample/hartree.xsf --tip dz2
+ppafm-generate-dftd3 -i sample/hartree.xsf --df_name PBE
+
+# Relax the probe particle in the force field
+ppafm-relaxed-scan -k 0.24 -q -0.1 --noLJ --Apauli 25.0
+
+# Plot the results
+ppafm-plot-results -k 0.24 -q -0.1 -a 1.0 --df

examples/BrClPyridine_FDBM_aims_density/run_gpu.py

@@ -0,0 +1,63 @@
+#!/usr/bin/env python3
+
+"""
+An example of using the full-density based model (FDBM) with an sample electron density obtained with FHI-aims.
+The elecron density in aims is all-electron, so the density at the nuclei positions can be extremely high.
+This sometimes causes artifacts to appear in the resuting simulation images. This can be mitigated by cutting
+out the high electron density values, which is demonstrated in this example.
+
+Simulations are printed both with and without the cutoff into two separate directories, `sims_cutoff100` and `sims_no_cutoff`.
+Compare the simulation image and observe the unphysical artifacts around the bromine atom at the top-right without the cutoff
+and how the artifacts disappear with the cutoff.
+"""
+
+from pathlib import Path
+
+import ppafm.ocl.field as FFcl
+import ppafm.ocl.oclUtils as oclu
+from ppafm.data import download_dataset
+from ppafm.ocl.AFMulator import AFMulator
+
+if __name__ == "__main__":
+
+    # Initialize an OpenCL environment. You can change i_platform to select the device to use
+    oclu.init_env(i_platform=0)
+
+    # Download input files (if not already downloaded)
+    tip_dir = Path("tip")
+    sample_dir = Path("sample")
+    download_dataset("CO-tip-densities", target_dir=tip_dir)
+    download_dataset("BrClPyridine", target_dir=sample_dir)
+
+    # Load all input files
+    rho_tip, xyzs_tip, Zs_tip = FFcl.TipDensity.from_file(tip_dir / "density_CO.xsf")
+    rho_tip_delta, _, _ = FFcl.TipDensity.from_file(tip_dir / "CO_delta_density_aims.xsf")
+    pot, xyzs, Zs = FFcl.HartreePotential.from_file(sample_dir / "hartree.xsf", scale=-1.0)  # Scale=-1.0 for correct units of potential (V) instead of energy (eV)
+    rho_sample, _, _ = FFcl.ElectronDensity.from_file(sample_dir / "density.xsf")
+
+    # Construct the simulator
+    afmulator = AFMulator(
+        pixPerAngstrome=10,
+        scan_dim=(160, 160, 40),
+        scan_window=((2, 2, 6), (18, 18, 10)),
+        rho=rho_tip,
+        rho_delta=rho_tip_delta,
+        A_pauli=25.0,  # Arbitrarily chosen values, not fitted to anything.
+        B_pauli=1.1,  # Higher values for the A and B parameters result in more severe artifacts.
+        fdbm_vdw_type="D3",
+        d3_params="PBE",
+        tipR0=[0, 0, 3],
+        tipStiffness=[0.25, 0.25, 0.0, 20.0],
+        npbc=(1, 1, 0),
+        df_steps=10,
+    )
+
+    # Let's run the simulation both with and without the cutoff
+    for cutoff in [None, 100]:
+
+        # Apply cutoff to sample electron density
+        rho_sample_cutoff = rho_sample.clamp(maximum=cutoff, in_place=False)
+        plot_dir = f"sims_cutoff{cutoff}" if cutoff else "sims_no_cutoff"
+
+        # Run simulation and plot images
+        X = afmulator(xyzs, Zs, pot, rho_sample_cutoff, plot_to_dir=plot_dir)

ppafm/cli/conv_rho.py

@@ -27,6 +27,7 @@ def main(argv=None):
     parser.add_argument('-o', '--output', action='store', default='pauli', help='Name of output energy/force files.')
     parser.add_argument('--saveDebugXsfs', action='store_true', help='Save auxiliary xsf files for debugging.')
     parser.add_argument('--no_negative_check', action='store_true', help='Input density files may contain negative voxels. This is handled by default by setting negative values to zero and rescaling the density so that the total charge is conserved. Setting this option disables the check.' )
+    parser.add_argument("--density_cutoff", action="store", default=None, type=float, help="Apply a cutoff to the electron densities to cut out high values. Helpful when using all-electron densities where extremely high values at the nuclei positions can cause artifacts in the resulting simulations. In these cases, 100 is usually a safe value to use, although sometimes when both the prefactor 'Apauli' and exponent 'Bpauli' are large, a lower cutoff may be required.")
     parser.add_arguments(['output_format', 'energy', 'Apauli', 'Bpauli'])
     # fmt: on
 
@@ -57,6 +58,11 @@ def main(argv=None):
             io.save_scal_field("sample_density_pow_%03.3f.xsf" % args.Bpauli, rho_sample, lvec_sample, data_format=args.output_format, head=head_sample)
             io.save_scal_field("tip_density_pow_%03.3f.xsf" % args.Bpauli, rho_tip, lvec_tip, data_format=args.output_format, head=head_tip)
 
+    if args.density_cutoff:
+        print(f">>> Applying a density cutoff of {args.density_cutoff} to sample and tip electron densities.")
+        rho_sample[rho_sample > args.density_cutoff] = args.density_cutoff
+        rho_tip[rho_tip > args.density_cutoff] = args.density_cutoff
+
     print(">>> Evaluating convolution E(R) = A*Integral_r ( rho_tip^B(r-R) * rho_sample^B(r) ) using FFT ... ")
     f_x, f_y, f_z, energy = fieldFFT.potential2forces_mem(rho_sample, lvec_sample, n_dim_sample, rho=rho_tip, doForce=True, doPot=True, deleteV=True)
 

ppafm/cli/gui/ppafm_gui.py

@@ -74,6 +74,7 @@
     "Rotation": "Rotation: Set sample counter-clockwise rotation angle around center of atom coordinates.",
     "fdbm_V0": "FDBM V0: Prefactor in Pauli interaction integral in the full-density based model.",
     "fdbm_alpha": "FDBM alpha: Exponent in Pauli interaction integral in the full-density based model.",
+    "fdbm_cutoff": "Density cutoff: cut out high values of the electron density. Useful when working with all-electron densities where extremely high values can arise near nuclei, causing artifacts in the images.",
     "PBCz": "z periodicity: When checked, the lattice is also periodic in z direction. This is usually not required, since the scan is aligned with the xy direction.",
     "PBC": "Periodic Boundaries: Lattice vectors for periodic images of atoms. Does not affect electrostatics calculated from a Hartree potential file, which is always assumed to be periodic.",
     "k": "k: Cantilever spring constant. Only appears as a scaling constant.",
@@ -123,11 +124,13 @@ def __init__(self, input_files=None, device=0, verbose=0):
         self.Zs = None
         self.qs = None
         self.rho_sample = None
+        self.rho_sample_original = None  # Stores the original density when a cutoff is applied
         self.rho_tip = None
         self.rho_tip_delta = None
         self.sample_lvec = None
         self.rot = np.eye(3)
         self.df_points = []
+        self.previous_cutoff = None
 
         # Initialize OpenCL environment on chosen device and create an afmulator instance to use for simulations
         oclu.init_env(device)
@@ -141,6 +144,7 @@ def __init__(self, input_files=None, device=0, verbose=0):
         self.afmulator.forcefield.verbose = verbose
         self.afmulator.scanner.verbose = verbose
         FFcl.bRuntime = verbose > 1
+        self.afmulator.bRuntime = verbose > 1
 
         # --- init QtMain
         QtWidgets.QMainWindow.__init__(self)
@@ -302,17 +306,42 @@ def updateParams(self, preset_none=True):
         if self.verbose > 0:
             print("updateParams", Q, sigma, multipole, tipStiffness, tipR0, use_point_charge, A_pauli, B_pauli, type(self.qs))
 
-        self.afmulator.iZPP = int(self.bxZPP.value())
         if self.rho_sample is not None:  # Use FDBM
-            if self.afmulator.B_pauli != B_pauli:
-                self.afmulator.setRho(self.rho_tip, B_pauli=B_pauli)
+
+            if self.bxCutoffEnable.isChecked():
+                self.bxCutoffValue.setDisabled(False)
+                cutoff = float(self.bxCutoffValue.value())
+            else:
+                self.bxCutoffValue.setDisabled(True)
+                cutoff = None
+
+            if cutoff is None or self.previous_cutoff is None:
+                cutoff_changed = cutoff != self.previous_cutoff
+            else:
+                cutoff_changed = not np.isclose(cutoff, self.previous_cutoff)
+            self.previous_cutoff = cutoff
+
+            if cutoff_changed:
+                if cutoff is not None:
+                    self.rho_sample = self.rho_sample_original.clamp(maximum=cutoff, in_place=False)
+                else:
+                    self.rho_sample = self.rho_sample_original
+
+            if self.afmulator.B_pauli != B_pauli or cutoff_changed:
+                if cutoff is not None:
+                    rho_tip = self.rho_tip.clamp(maximum=cutoff, in_place=False)
+                else:
+                    rho_tip = self.rho_tip
+                self.afmulator.setRho(rho_tip, B_pauli=B_pauli)
         else:
             if self.rho_tip is None:  # else we just keep the tip density loaded from file
                 if use_point_charge:
                     self.afmulator.setQs(Qs, QZs)
                     self.afmulator.setRho(None, sigma)
                 elif isinstance(self.qs, HartreePotential):
                     self.afmulator.setRho({multipole: Q}, sigma)
+
+        self.afmulator.iZPP = int(self.bxZPP.value())
         self.afmulator.scanner.stiffness = np.array(tipStiffness, dtype=np.float32) / -PPU.eVA_Nm
         self.afmulator.tipR0 = tipR0
         self.afmulator.A_pauli = A_pauli
@@ -471,6 +500,8 @@ def loadInput(self, main_input, rho_sample=None, rho_tip=None, rho_tip_delta=Non
             self.qs.release()
         if self.rho_sample is not None:
             self.rho_sample.release()
+        if self.rho_sample_original is not None:
+            self.rho_sample_original.release()
         if self.rho_tip is not None:
             self.rho_tip.release()
         if self.rho_tip_delta is not None:
@@ -502,6 +533,7 @@ def loadInput(self, main_input, rho_sample=None, rho_tip=None, rho_tip_delta=Non
         # Load auxiliary files (if any)
         if rho_sample:
             self.rho_sample, _, _ = ElectronDensity.from_file(rho_sample)
+            self.rho_sample_original = self.rho_sample
         else:
             self.rho_sample = None
         if rho_tip:
@@ -577,9 +609,14 @@ def loadInput(self, main_input, rho_sample=None, rho_tip=None, rho_tip_delta=Non
         if self.rho_sample is None:
             self.bxV0.setDisabled(True)
             self.bxAlpha.setDisabled(True)
+            self.bxCutoffEnable.setDisabled(True)
+            self.bxCutoffEnable.setChecked(False)
+            self.bxCutoffValue.setDisabled(True)
         else:
             self.bxV0.setDisabled(False)
             self.bxAlpha.setDisabled(False)
+            self.bxCutoffEnable.setDisabled(False)
+            self.previous_cutoff = None
 
         # Create geometry editor widget
         self.createGeomEditor()
@@ -1000,6 +1037,18 @@ def _create_probe_settings_ui(self):
         l.addWidget(lb)
         self.bxAlpha = _spin_box((0.1, 9.9), 1.0, 0.02, self.updateParams, TTips["fdbm_alpha"], l)
 
+        vb = QtWidgets.QHBoxLayout()
+        self.l0.addLayout(vb)
+        self.bxCutoffEnable = QtWidgets.QCheckBox()
+        self.bxCutoffEnable.setChecked(False)
+        self.bxCutoffEnable.toggled.connect(self.updateParams)
+        vb.addWidget(self.bxCutoffEnable)
+        lb = QtWidgets.QLabel("Density cutoff")
+        lb.setToolTip(TTips["fdbm_cutoff"])
+        vb.addWidget(lb)
+        self.bxCutoffValue = _spin_box((0.1, 10000.0), 100.0, 10.0, self.updateParams, TTips["fdbm_cutoff"], vb)
+        self.bxCutoffValue.setDisabled(True)
+
     def _create_scan_settings_ui(self):
         # Title
         vb = QtWidgets.QHBoxLayout()