Fix issues that occured when running ppafm-generate-elff with tip's density

examples/PTCDA_Hartree_dz2/example_ptcda_hartree.py

@@ -8,7 +8,7 @@
 from ppafm.cli.relaxed_scan import main as relaxed_scan
 
 
-def example_ptcda_hartree():
+def example_ptcda_hartree_dz2():
     script_location = Path(__file__).absolute().parent
 
     # Change directory to the location of this script

examples/pyridineDensOverlap/example_pyridine_density_overlap.py

@@ -0,0 +1,42 @@
+import os
+import urllib.request
+import zipfile
+from pathlib import Path
+
+from ppafm.cli.conv_rho import main as generate_conv_rho
+from ppafm.cli.generateDFTD3 import main as generate_dftd3
+from ppafm.cli.generateElFF import main as generate_elff
+from ppafm.cli.plot_results import main as plot_results
+from ppafm.cli.relaxed_scan import main as relaxed_scan
+
+
+def example_pyridine_density_overlap():
+    script_location = Path(__file__).absolute().parent
+
+    # Change directory to the location of this script
+    os.chdir(script_location)
+
+    if not Path("LOCPOT.xsf").exists():
+        urllib.request.urlretrieve("https://zenodo.org/records/14222456/files/pyridine.zip?download=1", "sample.zip")
+        # Unzip the file using python tools only
+        with zipfile.ZipFile("sample.zip", "r") as zip_ref:
+            zip_ref.extractall(".")
+
+    if not Path("tip/density_CO.xsf").exists():
+        urllib.request.urlretrieve("https://zenodo.org/records/14222456/files/CO_tip.zip?download=1", "tip.zip")
+        # Unzip the file using python tools only
+        with zipfile.ZipFile("tip.zip", "r") as zip_ref:
+            # Make "tip" directory if it does not exist
+            Path("tip").mkdir(exist_ok=True)
+            zip_ref.extractall("tip")
+
+    generate_conv_rho("-s CHGCAR.xsf -t tip/density_CO.xsf -B 1.0 -E".split())
+    generate_elff("-i LOCPOT.xsf --tip_dens tip/density_CO.xsf --Rcore 0.7 -E --doDensity".split())
+    generate_dftd3("-i LOCPOT.xsf --df_name PBE --energy".split())
+
+    relaxed_scan("-k 0.25 -q 1.0 --noLJ --Apauli 18.0 --bDebugFFtot".split())  # Note the --noLJ for loading separate Pauli and vdW instead of LJ force field
+    plot_results("-k 0.25 -q 1.0 -a 2.0 --df".split())
+
+
+if __name__ == "__main__":
+    example_pyridine_density_overlap()

examples/pyridineDensOverlap/run.sh

@@ -1,34 +1,18 @@
-#! /bin/bash
+echo "======= STEP 1 : Download the sample and tip files."
+wget https://zenodo.org/records/14222456/files/pyridine.zip?download=1 -O sample.zip
+wget https://zenodo.org/records/14222456/files/CO_tip.zip?download=1  -O tip.zip
 
-echo " ====== STEP 0 : Download Example Data-Files "
+echo "======= STEP 2 : Unzip the sample and tip files."
+unzip sample.zip
+unzip tip.zip -d tip
 
-# You should either install this: https://megatools.megous.com/
-# Or download it manually from Mega.nz
-echo "in order to download data files we need to install megatools "
-echo "alternatively you may download it manually from mega.nz webpage and place it to proper sub-folders ./tip and ./sample "
-echo "If you have Debian based linux system with apt-get this should work : "
-echo "enter your sudo password ( I promiss not hacking you ;-) : "
-sudo apt-get install megatools
+echo "======= STEP 3 : Generate force field grid."
+ppafm-conv-rho       -s CHGCAR.xsf -t tip/density_CO.xsf -B 1.0 -E
+ppafm-generate-elff  -i LOCPOT.xsf --tip_dens tip/density_CO.xsf --Rcore 0.7 -E --doDensity
+ppafm-generate-dftd3 -i LOCPOT.xsf --df_name PBE
 
-echo "Downloading Sample data-files ... "
-megadl 'https://mega.nz/file/bWhTUQDQ#7mS9E-wArUzHqOevCepckHezzO8uLLC0S1PbRDsiQfs'
-megadl 'https://mega.nz/file/rPgFxYAB#kQ6J90i4qQ4LlDFKq-k8PCX0rBie75_zdReNgNYbaKY'
-mkdir sample
-mv CHGCAR.xsf sample
-mv LOCPOT.xsf sample
-
-echo "Downloading Tip data-files ... "
-megadl 'https://mega.nz/#!2CgjyCZR!b0E-vPUp6TmkV0_uTwAYHDrMXv8mJcShYOVBZVSUYs0'
-mkdir tip
-mv CHGCAR.xsf tip
-
-echo "======= STEP 1 : Generate force field grid."
-ppafm-conv-rho       -s sample/CHGCAR.xsf -t tip/CHGCAR.xsf -B 1.0 -E
-ppafm-generate-elff  -i sample/LOCPOT.xsf --tip_dens tip/CHGCAR.xsf --Rcore 0.7 -E --doDensity
-ppafm-generate-dftd3 -i sample/LOCPOT.xsf --df_name PBE
-
-echo "======= STEP 2 : Relax Probe Particle using that force field grid."
+echo "======= STEP 4 : Relax Probe Particle using that force field grid."
 ppafm-relaxed-scan -k 0.25 -q 1.0 --noLJ --Apauli 18.0 --bDebugFFtot # Note the --noLJ for loading separate Pauli and vdW instead of LJ force field
 
-echo "======= STEP 3 : Plot the results."
+echo "======= STEP 5 : Plot the results."
 ppafm-plot-results -k 0.25 -q 1.0 -a 2.0 --df

ppafm/HighLevel.py

@@ -295,15 +295,15 @@ def computeDFTD3(input_file, df_params="PBE", geometry_format=None, save_format=
     # Load atomic geometry
     atoms, nDim, lvec = io.loadGeometry(input_file, format=geometry_format, parameters=parameters)
     elem_dict = PPU.getFFdict(PPU.loadSpecies())
-    iZs, Rs, _ = PPU.parseAtoms(atoms, elem_dict, autogeom=False, PBC=parameters.PBC, lvec=lvec)
+    iZs, Rs, _ = PPU.parseAtoms(atoms, elem_dict, autogeom=False, PBC=parameters.PBC, lvec=lvec, parameters=parameters)
     iPP = PPU.atom2iZ(parameters.probeType, elem_dict)
 
     # Compute coefficients for each atom
     df_params = d3.get_df_params(df_params)
     coeffs = core.computeD3Coeffs(Rs, iZs, iPP, df_params)
 
     # Compute the force field
-    FF, V = prepareArrays(None, compute_energy)
+    FF, V = prepareArrays(None, compute_energy, parameters=parameters)
     core.setFF_shape(np.shape(FF), lvec, parameters=parameters)
     core.getDFTD3FF(shift_positions(Rs, -lvec[0]), coeffs)
 
@@ -364,15 +364,13 @@ def computeELFF_pointCharge(geomFile, geometry_format=None, tip="s", save_format
 
 
 def computeElFF(V, lvec, nDim, tip, computeVpot=False, tilt=0.0, sigma=None, deleteV=True, parameters=None):
-    if verbose > 0:
-        print(" ========= get electrostatic forcefiled from hartree ")
     rho = None
     multipole = None
     if sigma is None:
         sigma = parameters.sigma
-    if type(tip) is np.ndarray:
+    if isinstance(tip, (list, np.ndarray)):
         rho = tip
-    elif type(tip) is dict:
+    elif isinstance(tip, dict):
         multipole = tip
     else:
         if tip in {"s", "px", "py", "pz", "dx2", "dy2", "dz2", "dxy", "dxz", "dyz"}:
@@ -383,8 +381,6 @@ def computeElFF(V, lvec, nDim, tip, computeVpot=False, tilt=0.0, sigma=None, del
             if any(nDim_tip != nDim):
                 sys.exit("Error: Input file for tip charge density has been specified, but the dimensions are incompatible with the Hartree potential file!")
             rho *= -1  # Negative charge density from positive electron density
-    if verbose > 0:
-        print(" computing convolution with tip by FFT ")
     Fel_x, Fel_y, Fel_z, Vout = fFFT.potential2forces_mem(V, lvec, nDim, rho=rho, sigma=sigma, multipole=multipole, doPot=computeVpot, tilt=tilt, deleteV=deleteV)
     FFel = io.packVecGrid(Fel_x, Fel_y, Fel_z)
     del Fel_x, Fel_y, Fel_z
@@ -423,7 +419,7 @@ def _getAtomsWhichTouchPBCcell(Rs, elems, nDim, lvec, Rcut, bSaveDebug, fname=No
     return Rs_, elems
 
 
-def getAtomsWhichTouchPBCcell(fname, Rcut=1.0, bSaveDebug=True, geometry_format=None):
+def getAtomsWhichTouchPBCcell(fname, Rcut=1.0, bSaveDebug=True, geometry_format=None, parameters=None):
     atoms, nDim, lvec = io.loadGeometry(fname, format=geometry_format, parameters=parameters)
     Rs = np.array(atoms[1:4])  # get just positions x,y,z
     elems = np.array(atoms[0])

ppafm/cli/conv_rho.py

@@ -15,7 +15,7 @@ def handle_negative_density(rho):
     rho *= q / rho.sum()
 
 
-def main():
+def main(argv=None):
     parser = common.CLIParser(
         description="Calculate the density overlap integral for Pauli force field in the full-density based model. "
         "The integral has two parameters A and B, and is of form A*Integral( rho_tip^B * rho_sample^B )"
@@ -30,7 +30,7 @@ def main():
     parser.add_arguments(['output_format', 'energy', 'Apauli', 'Bpauli'])
     # fmt: on
 
-    args = parser.parse_args()
+    args = parser.parse_args(argv)
 
     print(">>> Loading sample from ", args.sample, " ... ")
     rho_sample, lvec_sample, n_dim_sample, head_sample = io.loadXSF(args.sample)

ppafm/cli/generateDFTD3.py

@@ -8,7 +8,7 @@
 from ..HighLevel import computeDFTD3
 
 
-def main():
+def main(argv=None):
     parser = common.CLIParser(
         description="Generate Grimme DFT-D3 vdW force field using the Becke-Johnson damping function. The generated force field is saved to FFvdW_{x,y,z}.[ext]."
     )
@@ -32,7 +32,7 @@ def main():
         metavar=("s6", "s8", "a1", "a2"),
         help="Manually specify scaling parameters s6, s8, a1, a2. Overwrites --df_name.",
     )
-    args = parser.parse_args()
+    args = parser.parse_args(argv)
 
     parameters = common.PpafmParameters.from_file("params.ini")
 
@@ -48,7 +48,7 @@ def main():
             sys.exit(1)
         df_params = args.df_name
 
-    computeDFTD3(args.input, df_params=df_params, geometry_format=args.input_format, save_format=args.output_format, compute_energy=args.energy)
+    computeDFTD3(args.input, df_params=df_params, geometry_format=args.input_format, save_format=args.output_format, compute_energy=args.energy, parameters=parameters)
 
 
 if __name__ == "__main__":

ppafm/cli/generateElFF.py

@@ -44,7 +44,7 @@ def main(argv=None):
         if args.tip_dens is None:
             raise Exception("Rcore>0 but no tip density provided!")
         valence_electrons_dictionary = loadValenceElectronDict()
-        rs_tip, elems_tip = getAtomsWhichTouchPBCcell(args.tip_dens, Rcut=args.Rcore)
+        rs_tip, elems_tip = getAtomsWhichTouchPBCcell(args.tip_dens, Rcut=args.Rcore, parameters=parameters)
 
     atoms_samp, _, lvec_samp = io.loadGeometry(args.input, format=args.input_format, parameters=parameters)
     head_samp = io.primcoords2Xsf(atoms_samp[0], [atoms_samp[1], atoms_samp[2], atoms_samp[3]], lvec_samp)

ppafm/common.py

@@ -1,5 +1,6 @@
 #!/usr/bin/python
 
+import copy
 import os
 import typing
 from argparse import ArgumentParser
@@ -320,7 +321,7 @@ def add_arguments(self, arg_names):
         for name in arg_names:
             if name not in self.cli_args:
                 raise ValueError(f"Invalid argument name `{name}`")
-            arg_dict = self.cli_args[name]
+            arg_dict = copy.deepcopy(self.cli_args[name])
             if "help" not in arg_dict:
                 raise ValueError(f"No help message defined for `{name}`")
             if "short_name" in arg_dict: