Implement PpafmParameters class based on the pydantic.BaseModel (#…

…275)

Previously, the run parameters were stored in a global dictionary called `params`. This
solution was sub-optimal, as it was very tricky to do multiple runs to avoid "contamination"
of parameters by the previous run. In this PR I introduce the parameters class `PpafmParameters`.
Now, every run will have its set of parameters. Moving forward, we can simplify the signature of
many functions relying on the parameters provided by the `PpafmParameters` object. The
original `params` dictionary is now completely removed.

Previously, the default values for parameters were scattered around. Now, all the default values are
set by the `PpafmParameters`, which allows the removal of the `params.ini` file which was part
of the package distribution.

The infrastructure built here offers an easy way to  switch input format from `ini`
to `toml`, as suggested in the issue #153. To allow easy switching `to_file` method is implemented.
It allows dumping the current parameters to the `toml` format.

Finally, tests were added to make sure the `PpafmParameters`  class behaves as expected.

ppafm/HighLevel.py

@@ -32,7 +32,7 @@ def meshgrid3d(xs, ys, zs):
     Xs, Ys = np.meshgrid(xs, ys)
 
 
-def trjByDir(n, d=[0.0, 0.0, PPU.params["scanStep"][2]], p0=[0, 0, PPU.params["scanMin"][2]]):
+def trjByDir(n, d, p0):
     trj = np.zeros((n, 3))
     trj[:, 0] = p0[0] + (np.arange(n)[::-1]) * d[0]
     trj[:, 1] = p0[1] + (np.arange(n)[::-1]) * d[1]
@@ -131,7 +131,7 @@ def relaxedScan3D_omp(xTips, yTips, zTips, trj=None, bF3d=False):
     return fzs, rs
 
 
-def perform_relaxation(lvec, FFLJ, FFel=None, FFpauli=None, FFboltz=None, FFkpfm_t0sV=None, FFkpfm_tVs0=None, tipspline=None, bPPdisp=False, bFFtotDebug=False):
+def perform_relaxation(lvec, FFLJ, FFel=None, FFpauli=None, FFboltz=None, FFkpfm_t0sV=None, FFkpfm_tVs0=None, tipspline=None, bPPdisp=False, bFFtotDebug=False, parameters=None):
     global FF  # We need FF global otherwise it is garbage collected and program crashes inside C++ e.g. in stiffnessMatrix()
     if verbose > 0:
         print(">>>BEGIN: perform_relaxation()")
@@ -151,36 +151,36 @@ def perform_relaxation(lvec, FFLJ, FFel=None, FFpauli=None, FFboltz=None, FFkpfm
             if verbose > 0:
                 print("cannot load tip spline from " + tipspline)
             sys.exit()
-    xTips, yTips, zTips, lvecScan = PPU.prepareScanGrids()
+    xTips, yTips, zTips, lvecScan = PPU.prepareScanGrids(parameters=parameters)
     FF = FFLJ.copy()
     if FFel is not None:
-        FF += FFel * PPU.params["charge"]
+        FF += FFel * parameters.charge
         if verbose > 0:
-            print("adding charge:", PPU.params["charge"])
+            print("adding charge:", parameters.charge)
     if FFkpfm_t0sV is not None and FFkpfm_tVs0 is not None:
-        FF += (PPU.params["charge"] * FFkpfm_t0sV - FFkpfm_tVs0) * PPU.params["Vbias"]
+        FF += (parameters.charge * FFkpfm_t0sV - FFkpfm_tVs0) * parameters.Vbias
         if verbose > 0:
-            print("adding charge:", PPU.params["charge"], "and bias:", PPU.params["Vbias"], "V")
+            print("adding charge:", parameters.charge, "and bias:", parameters.Vbias, "V")
     if FFpauli is not None:
-        FF += FFpauli * PPU.params["Apauli"]
+        FF += FFpauli * parameters.Apauli
     if FFboltz != None:
         FF += FFboltz
     if bFFtotDebug:
         io.save_vec_field("FFtotDebug", FF, lvec)
-    core.setFF_shape(np.shape(FF), lvec)
+    core.setFF_shape(np.shape(FF), lvec, parameters=parameters)
     core.setFF_Fpointer(FF)
-    if (PPU.params["stiffness"] < 0.0).any():
-        PPU.params["stiffness"] = np.array([PPU.params["klat"], PPU.params["klat"], PPU.params["krad"]])
+    if (np.array(parameters.stiffness) < 0.0).any():
+        parameters.stiffness = np.array([parameters.klat, parameters.klat, parameters.krad])
     if verbose > 0:
-        print("stiffness:", PPU.params["stiffness"])
-    core.setTip(kSpring=np.array((PPU.params["stiffness"][0], PPU.params["stiffness"][1], 0.0)) / -PPU.eVA_Nm, kRadial=PPU.params["stiffness"][2] / -PPU.eVA_Nm)
+        print("stiffness:", parameters.stiffness)
+    core.setTip(kSpring=np.array((parameters.stiffness[0], parameters.stiffness[1], 0.0)) / -PPU.eVA_Nm, kRadial=parameters.stiffness[2] / -PPU.eVA_Nm, parameters=parameters)
 
     # grid origin has to be moved to zero, hence the subtraction of lvec[0,:] from trj and xTip, yTips, zTips
     trj = None
-    if PPU.params["tiltedScan"]:
-        trj = trjByDir(len(zTips), d=PPU.params["scanTilt"], p0=[0.0, 0.0, PPU.params["scanMin"][2] - lvec[0, 2]])
-    # fzs, PPpos = relaxedScan3D(xTips - lvec[0, 0], yTips - lvec[0, 1], zTips - lvec[0, 2], trj=trj, bF3d=PPU.params["tiltedScan"])
-    fzs, PPpos = relaxedScan3D_omp(xTips - lvec[0, 0], yTips - lvec[0, 1], zTips - lvec[0, 2], trj=trj, bF3d=PPU.params["tiltedScan"])
+    if parameters.tiltedScan:
+        trj = trjByDir(len(zTips), d=parameters.scanTilt, p0=[0.0, 0.0, parameters.scanMin[2] - lvec[0, 2]])
+    # fzs, PPpos = relaxedScan3D(xTips - lvec[0, 0], yTips - lvec[0, 1], zTips - lvec[0, 2], trj=trj, bF3d=parameters.tiltedScan)
+    fzs, PPpos = relaxedScan3D_omp(xTips - lvec[0, 0], yTips - lvec[0, 1], zTips - lvec[0, 2], trj=trj, bF3d=parameters.tiltedScan)
 
     # transform probe-particle positions back to the original coordinates
     PPpos[:, :, :, 0] += lvec[0, 0]
@@ -189,7 +189,7 @@ def perform_relaxation(lvec, FFLJ, FFel=None, FFpauli=None, FFboltz=None, FFkpfm
 
     if bPPdisp:
         PPdisp = PPpos.copy()
-        init_pos = np.array(np.meshgrid(xTips, yTips, zTips)).transpose(3, 1, 2, 0) + np.array([PPU.params["r0Probe"][0], PPU.params["r0Probe"][1], -PPU.params["r0Probe"][2]])
+        init_pos = np.array(np.meshgrid(xTips, yTips, zTips)).transpose(3, 1, 2, 0) + np.array([parameters.r0Probe[0], parameters.r0Probe[1], -parameters.r0Probe[2]])
         PPdisp -= init_pos
     else:
         PPdisp = None
@@ -201,14 +201,14 @@ def perform_relaxation(lvec, FFLJ, FFel=None, FFpauli=None, FFboltz=None, FFkpfm
 # ==== Forcefield grid generation
 
 
-def prepareArrays(FF, Vpot):
-    if PPU.params["gridN"][0] <= 0:
+def prepareArrays(FF, Vpot, parameters):
+    if parameters.gridN[0] <= 0:
         PPU.autoGridN()
     if FF is None:
-        gridN = PPU.params["gridN"]
+        gridN = parameters.gridN
         FF = np.zeros((gridN[2], gridN[1], gridN[0], 3))
     else:
-        PPU.params["gridN"] = np.shape(FF)
+        parameters.gridN = np.shape(FF)
     core.setFF_Fpointer(FF)
     if Vpot:
         V = np.zeros((gridN[2], gridN[1], gridN[0]))
@@ -218,26 +218,26 @@ def prepareArrays(FF, Vpot):
     return FF, V
 
 
-def computeLJ(geomFile, speciesFile, geometry_format=None, save_format=None, computeVpot=False, Fmax=Fmax_DEFAULT, Vmax=Vmax_DEFAULT, ffModel="LJ"):
+def computeLJ(geomFile, speciesFile, geometry_format=None, save_format=None, computeVpot=False, Fmax=Fmax_DEFAULT, Vmax=Vmax_DEFAULT, ffModel="LJ", parameters=None):
     if verbose > 0:
         print(">>>BEGIN: computeLJ()")
     # --- load species (LJ potential)
     FFparams = PPU.loadSpecies(speciesFile)
     elem_dict = PPU.getFFdict(FFparams)
     # print elem_dict
     # --- load atomic geometry
-    atoms, nDim, lvec = io.loadGeometry(geomFile, format=geometry_format, params=PPU.params)
+    atoms, nDim, lvec = io.loadGeometry(geomFile, format=geometry_format, parameters=parameters)
     atomstring = io.primcoords2Xsf(PPU.atoms2iZs(atoms[0], elem_dict), [atoms[1], atoms[2], atoms[3]], lvec)
     if verbose > 0:
-        print(PPU.params["gridN"], PPU.params["gridO"], PPU.params["gridA"], PPU.params["gridB"], PPU.params["gridC"])
-    iZs, Rs, Qs = PPU.parseAtoms(atoms, elem_dict, autogeom=False, PBC=PPU.params["PBC"], lvec=lvec)
+        print(parameters.gridN, parameters.gridO, parameters.gridA, parameters.gridB, parameters.gridC)
+    iZs, Rs, Qs = PPU.parseAtoms(atoms, elem_dict, autogeom=False, PBC=parameters.PBC, lvec=lvec, parameters=parameters)
     # --- prepare LJ parameters
-    iPP = PPU.atom2iZ(PPU.params["probeType"], elem_dict)
+    iPP = PPU.atom2iZ(parameters.probeType, elem_dict)
     # --- prepare arrays and compute
-    FF, V = prepareArrays(None, computeVpot)
+    FF, V = prepareArrays(None, computeVpot, parameters=parameters)
     if verbose > 0:
         print("FFLJ.shape", FF.shape)
-    core.setFF_shape(np.shape(FF), lvec)
+    core.setFF_shape(np.shape(FF), lvec, parameters=parameters)
 
     # shift atoms to the coordinate system in which the grid origin is zero
     Rs0 = shift_positions(Rs, -lvec[0])
@@ -246,7 +246,7 @@ def computeLJ(geomFile, speciesFile, geometry_format=None, save_format=None, com
         REs = PPU.getAtomsRE(iPP, iZs, FFparams)
         core.getMorseFF(Rs0, REs)  # THE MAIN STUFF HERE
     elif ffModel == "vdW":
-        vdWDampKind = PPU.params["vdWDampKind"]
+        vdWDampKind = parameters.vdWDampKind
         if vdWDampKind == 0:
             cLJs = PPU.getAtomsLJ(iPP, iZs, FFparams)
             core.getVdWFF(Rs0, cLJs)  # THE MAIN STUFF HERE
@@ -277,7 +277,7 @@ def computeLJ(geomFile, speciesFile, geometry_format=None, save_format=None, com
     return FF, V, nDim, lvec
 
 
-def computeDFTD3(input_file, df_params="PBE", geometry_format=None, save_format=None, compute_energy=False):
+def computeDFTD3(input_file, df_params="PBE", geometry_format=None, save_format=None, compute_energy=False, parameters=None):
     """
     Compute the Grimme DFT-D3 force field and optionally save to a file. See also :meth:`.add_dftd3`.
 
@@ -297,18 +297,18 @@ 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, params=PPU.params)
+    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=PPU.params["PBC"], lvec=lvec)
-    iPP = PPU.atom2iZ(PPU.params["probeType"], elem_dict)
+    iZs, Rs, _ = PPU.parseAtoms(atoms, elem_dict, autogeom=False, PBC=parameters.PBC, lvec=lvec)
+    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)
-    core.setFF_shape(np.shape(FF), lvec)
+    core.setFF_shape(np.shape(FF), lvec, parameters=parameters)
     core.getDFTD3FF(shift_positions(Rs, -lvec[0]), coeffs)
 
     # Save to file
@@ -321,7 +321,7 @@ def computeDFTD3(input_file, df_params="PBE", geometry_format=None, save_format=
     return FF, V, lvec
 
 
-def computeELFF_pointCharge(geomFile, geometry_format=None, tip="s", save_format=None, computeVpot=False, Fmax=Fmax_DEFAULT, Vmax=Vmax_DEFAULT):
+def computeELFF_pointCharge(geomFile, geometry_format=None, tip="s", save_format=None, computeVpot=False, Fmax=Fmax_DEFAULT, Vmax=Vmax_DEFAULT, parameters=None):
     if verbose > 0:
         print(">>>BEGIN: computeELFF_pointCharge()")
     tipKinds = {"s": 0, "pz": 1, "dz2": 2}
@@ -333,14 +333,14 @@ def computeELFF_pointCharge(geomFile, geometry_format=None, tip="s", save_format
     elem_dict = PPU.getFFdict(FFparams)
     # print elem_dict
 
-    atoms, nDim, lvec = io.loadGeometry(geomFile, format=geometry_format, params=PPU.params)
+    atoms, nDim, lvec = io.loadGeometry(geomFile, format=geometry_format, parameters=parameters)
     atomstring = io.primcoords2Xsf(PPU.atoms2iZs(atoms[0], elem_dict), [atoms[1], atoms[2], atoms[3]], lvec)
     # --- prepare arrays and compute
     if verbose > 0:
-        print(PPU.params["gridN"], PPU.params["gridA"], PPU.params["gridB"], PPU.params["gridC"])
-    iZs, Rs, Qs = PPU.parseAtoms(atoms, elem_dict=elem_dict, autogeom=False, PBC=PPU.params["PBC"], lvec=lvec)
-    FF, V = prepareArrays(None, computeVpot)
-    core.setFF_shape(np.shape(FF), lvec)
+        print(parameters.gridN, parameters.gridA, parameters.gridB, parameters.gridC)
+    _, Rs, Qs = PPU.parseAtoms(atoms, elem_dict=elem_dict, autogeom=False, PBC=parameters.PBC, lvec=lvec, parameters=parameters)
+    FF, V = prepareArrays(None, computeVpot, parameters=parameters)
+    core.setFF_shape(np.shape(FF), lvec, parameters=parameters)
 
     # shift atoms to the coordinate system in which the grid origin is zero
     Rs0 = shift_positions(Rs, -lvec[0])
@@ -367,13 +367,13 @@ def computeELFF_pointCharge(geomFile, geometry_format=None, tip="s", save_format
     return FF, V, nDim, lvec
 
 
-def computeElFF(V, lvec, nDim, tip, computeVpot=False, tilt=0.0, sigma=None, deleteV=True):
+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 = PPU.params["sigma"]
+        sigma = parameters.sigma
     if type(tip) is np.ndarray:
         rho = tip
     elif type(tip) is dict:
@@ -428,14 +428,16 @@ def _getAtomsWhichTouchPBCcell(Rs, elems, nDim, lvec, Rcut, bSaveDebug, fname=No
 
 
 def getAtomsWhichTouchPBCcell(fname, Rcut=1.0, bSaveDebug=True, geometry_format=None):
-    atoms, nDim, lvec = io.loadGeometry(fname, format=geometry_format, params=PPU.params)
+    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])
     Rs, elems = _getAtomsWhichTouchPBCcell(Rs, elems, nDim, lvec, Rcut, bSaveDebug, fname)
     return Rs, elems
 
 
-def subtractCoreDensities(rho, lvec, elems=None, Rs=None, fname=None, valElDict=None, Rcore=0.7, bSaveDebugDens=False, bSaveDebugGeom=True, head=io.XSF_HEAD_DEFAULT):
+def subtractCoreDensities(
+    rho, lvec, elems=None, Rs=None, fname=None, valElDict=None, Rcore=0.7, bSaveDebugDens=False, bSaveDebugGeom=True, head=io.XSF_HEAD_DEFAULT, parameters=None
+):
     nDim = rho.shape
     if fname is not None:
         elems, Rs = getAtomsWhichTouchPBCcell(fname, Rcut=Rcore, bSaveDebug=bSaveDebugDens)
@@ -451,7 +453,7 @@ def subtractCoreDensities(rho, lvec, elems=None, Rs=None, fname=None, valElDict=
         print("V : ", V, " N: ", N, " dV: ", dV)
     if verbose > 0:
         print("sum(RHO): ", rho.sum(), " Nelec: ", rho.sum() * dV, " voxel volume: ", dV)  # check sum
-    core.setFF_shape(rho.shape, lvec)  # set grid sampling dimension and shape
+    core.setFF_shape(rho.shape, lvec, parameters=parameters)  # set grid sampling dimension and shape
     core.setFF_Epointer(rho)  # set pointer to array with density data (to write into)
     if verbose > 0:
         print(">>> Projecting Core Densities ... ")

ppafm/cli/fitting/fitting.py

@@ -1,15 +1,14 @@
 #!/usr/bin/python
 import os
+from collections import OrderedDict
 from optparse import OptionParser
 
-import __main__ as main
 import numpy as np
 from scipy.interpolate import RegularGridInterpolator
-from scipy.optimize import basinhopping, minimize
+from scipy.optimize import minimize
 
-import ppafm as PPU
 import ppafm.HighLevel as PPH
-from ppafm import io
+from ppafm import common, io
 
 iteration = 0
 
@@ -53,31 +52,31 @@ def getFzlist(BIGarray, MIN, MAX, points):
 FFparams = None
 if os.path.isfile("atomtypes.ini"):
     print(">> LOADING LOCAL atomtypes.ini")
-    FFparams = PPU.loadSpecies("atomtypes.ini")
+    FFparams = common.loadSpecies("atomtypes.ini")
     print(FFparams)
     elem_dict = {}
     for i, ff in enumerate(FFparams):
         elem_dict[ff[3]] = i
 else:
     raise ValueError('Please provide the file "atomtypes.ini"')
 
+
+parameters = common.PpafmParameters.from_file("params.ini")
 print(" >> OVEWRITING SETTINGS by params.ini  ")
-PPU.loadParams("params.ini", FFparams=FFparams)
-scan_min = PPU.params["scanMin"]
-scan_max = PPU.params["scanMax"]
-atoms, nDim, lvec = io.loadGeometry("p_eq.xyz", params=PPU.params)
+scan_min = parameters.scanMin
+scan_max = parameters.scanMax
+atoms, nDim, lvec = io.loadGeometry("p_eq.xyz", parameters=parameters)
 # The function automatically loads the geometry from the file of any
 # supported format. The decision about the file format is based on the
 # filename extension or it can be supplied by an extra "format" argument
-PPU.params["gridN"] = nDim
-PPU.params["gridA"] = lvec[1]
-PPU.params["gridB"] = lvec[2]
-PPU.params["gridC"] = lvec[3]
+parameters.gridN = nDim
+parameters.gridA = lvec[1]
+parameters.gridB = lvec[2]
+parameters.gridC = lvec[3]
 V, lvec_bak, nDim_bak, head = io.loadCUBE("hartree.cube")
 loaded_forces = np.loadtxt("frc_tip.txt", converters={0: pm2a, 1: pm2a, 2: pm2a}, skiprows=2, usecols=(0, 1, 2, 5))
 points = loaded_forces[:, :3]
-iZs, Rs, Qs = PPH.parseAtoms(atoms, autogeom=False, PBC=PPU.params["PBC"], FFparams=FFparams)
-from collections import OrderedDict
+iZs, Rs, Qs = PPH.parseAtoms(atoms, autogeom=False, PBC=parameters.PBC, FFparams=FFparams)
 
 fit_dict = OrderedDict()
 
@@ -160,12 +159,12 @@ def comp_msd(x=[]):
     global iteration
     iteration += 1
     update_fit_dict(x)  # updating the array with the optimized values
-    PPU.apply_options(fit_dict)  # setting up all the options according to their
+    parameters.apply_options(fit_dict)  # setting up all the options according to their
     # current values
     update_atoms(atms=fit_dict["atom"])
     print(FFparams)
     FFLJ, VLJ = PPH.computeLJFF(iZs, Rs, FFparams)
-    FFel = PPH.computeElFF(V, lvec_bak, nDim_bak, PPU.params["tip"])
+    FFel = PPH.computeElFF(V, lvec_bak, nDim_bak, parameters.tip)
     FFboltz = None
     fzs, PPpos, PPdisp, lvecScan = PPH.perform_relaxation(lvec, FFLJ, FFel, FFboltz, tipspline=None)
     Fzlist = getFzlist(BIGarray=fzs, MIN=scan_min, MAX=scan_max, points=points)
@@ -191,7 +190,7 @@ def comp_msd(x=[]):
 
     (options, args) = parser.parse_args()
     opt_dict = vars(options)
-    PPU.apply_options(opt_dict)  # setting up all the options according to their
+    parameters.apply_options(opt_dict)  # setting up all the options according to their
     x_new, bounds = set_fit_dict(opt=opt_dict)
     print("params", x_new)
     print("bounds", bounds)