diff --git a/ppafm/HighLevel.py b/ppafm/HighLevel.py index d95bf595..18b605aa 100755 --- a/ppafm/HighLevel.py +++ b/ppafm/HighLevel.py @@ -167,7 +167,7 @@ def perform_relaxation(lvec, FFLJ, FFel=None, FFpauli=None, FFboltz=None, FFkpfm 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 (np.array(parameters.stiffness) < 0.0).any(): parameters.stiffness = np.array([parameters.klat, parameters.klat, parameters.krad]) @@ -237,7 +237,7 @@ def computeLJ(geomFile, speciesFile, geometry_format=None, save_format=None, com 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]) @@ -308,7 +308,7 @@ def computeDFTD3(input_file, df_params="PBE", geometry_format=None, save_format= # 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 @@ -340,7 +340,7 @@ def computeELFF_pointCharge(geomFile, geometry_format=None, tip="s", save_format 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) + 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]) @@ -435,7 +435,9 @@ def getAtomsWhichTouchPBCcell(fname, Rcut=1.0, bSaveDebug=True, geometry_format= 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 ... ") diff --git a/ppafm/cli/generateTraining_PVE.py b/ppafm/cli/generateTraining_PVE.py index f5085839..171a805f 100644 --- a/ppafm/cli/generateTraining_PVE.py +++ b/ppafm/cli/generateTraining_PVE.py @@ -40,7 +40,7 @@ force_field, _ = prepareArrays(None, False) print("FFLJ.shape", force_field.shape) -core.setFF_shape(np.shape(force_field), lvec_t) +core.setFF_shape(np.shape(force_field), lvec_t, parameters=parameters) base_dir = os.getcwd() paths = ["out1", "out2"] diff --git a/ppafm/common.py b/ppafm/common.py index e3c7a4e0..58400545 100644 --- a/ppafm/common.py +++ b/ppafm/common.py @@ -76,58 +76,6 @@ class Config: validate_assignment = True -# fmt: off -# default parameters of simulation -params = { - "PBC": True, - "nPBC": np.array([1, 1, 1]), - "gridN": np.array([-1, -1, -1]).astype(int), - "gridO": np.array([0.0, 0.0, 0.0]), - "gridA": np.array([0.0, 0.0, 0.0]), - "gridB": np.array([0.0, 0.0, 0.0]), - "gridC": np.array([0.0, 0.0, 0.0]), - "FFgrid0": np.array([-1.0, -1.0, -1.0]), - "FFgridA": np.array([-1.0, -1.0, -1.0]), - "FFgridB": np.array([-1.0, -1.0, -1.0]), - "FFgridC": np.array([-1.0, -1.0, -1.0]), - "moleculeShift": np.array([0.0, 0.0, 0.0]), - "probeType": "O", - "charge": 0.00, - "Apauli": 18.0, - "Bpauli": 1.0, - "ffModel": "LJ", - "Rcore": 0.7, - "r0Probe": np.array([0.00, 0.00, 4.00]), - "stiffness": np.array([-1.0, -1.0, -1.0]), - "klat": 0.5, - "krad": 20.00, - "tip": "s", - "sigma": 0.7, - "scanStep": np.array([0.10, 0.10, 0.10]), - "scanMin": np.array([0.0, 0.0, 5.0]), - "scanMax": np.array([20.0, 20.0, 8.0]), - "scanTilt": np.array([0.0, 0.0, -0.1]), - "tiltedScan": False, - "kCantilever": 1800.0, - "f0Cantilever": 30300.0, - "Amplitude": 1.0, - "plotSliceFrom": 16, - "plotSliceTo": 22, - "plotSliceBy": 1, - "imageInterpolation": "bicubic", - "colorscale": "gray", - "colorscale_kpfm": "seismic", - "ddisp": 0.05, - "aMorse": -1.6, - "tip_base": np.array(["None", 0.00]), - "Rtip": 30.0, - "permit": 0.00552634959, - "vdWDampKind": 2, - "#": None, -} -# fmt: on - - class CLIParser(ArgumentParser): """ Subclass from the built-in ArgumentParser with functionality to easily add arguments commonly used in ppafm scripts. @@ -163,7 +111,7 @@ class CLIParser(ArgumentParser): "Rcore": { "action": "store", "type": float, - "default": params["Rcore"], + "default": 0.7, "help": "Width of nuclear charge density blob to achieve charge neutrality [Å]", }, "sigma": { @@ -261,7 +209,7 @@ class CLIParser(ArgumentParser): def _check_params_args(self): for arg in self._params_args: - if arg not in params: + if arg not in PpafmParameters.model_filds: raise ValueError(f"Argument name `{arg}` does not match with any parameter in global parameters dictionary.") @property @@ -555,11 +503,11 @@ def loadParams(fname, parameters): parameters.tip = parameters.tip.replace('"', "") parameters.tip = parameters.tip.replace("'", "") - ### necessary for working even with quotemarks in params.ini + # Necessary for working even with quotemarks in params.ini parameters.tip_base[0] = parameters.tip_base[0].replace('"', "") parameters.tip_base[0] = parameters.tip_base[0].replace("'", "") - ### necessary for working even with quotemarks in params.ini - print("Loaded parameters from ", fname) + if verbose > 0: + print("Loaded parameters from ", fname) def apply_options(opt, parameters): @@ -572,7 +520,7 @@ def apply_options(opt, parameters): continue if hasattr(parameters, key): setattr(parameters, key, value) - if key == "klat" and params["stiffness"][0] > 0.0: + if key == "klat" and parameters.stiffness[0] > 0.0: print("Overriding stiffness parameter with klat") parameters.stiffness = np.array([-1.0, -1.0, -1.0]) # klat overrides stiffness if verbose > 0: @@ -604,14 +552,14 @@ def loadSpecies(fname=None): # load atoms species parameters form a file ( currently used to load Lennard-Jones parameters ) def loadSpeciesLines(lines): - params = [] + parameters = [] for l in lines: l = l.split() if len(l) >= 5: # print l - params.append((float(l[0]), float(l[1]), float(l[2]), int(l[3]), l[4])) + parameters.append((float(l[0]), float(l[1]), float(l[2]), int(l[3]), l[4])) return np.array( - params, + parameters, dtype=[ ("rmin", np.float64), ("epsilon", np.float64), diff --git a/ppafm/core.py b/ppafm/core.py index 7839fee4..c64f86a7 100755 --- a/ppafm/core.py +++ b/ppafm/core.py @@ -41,13 +41,13 @@ def setGridN(n): lib.setGridCell.restype = None -def setGridCell(cell=None): +def setGridCell(cell=None, parameters=None): if cell is None: cell = np.array( [ - PPU.params["gridA"], - PPU.params["gridB"], - PPU.params["gridC"], + parameters.gridA, + parameters.gridB, + parameters.gridC, ], dtype=np.float64, ).copy() @@ -63,10 +63,10 @@ def setGridCell(cell=None): lib.setGridCell(cell) -def setFF_shape(n_, cell): +def setFF_shape(n_, cell, parameters=None): n = np.array(n_).astype(np.int32) lib.setGridN(n) - setGridCell(cell) + setGridCell(cell, parameters=parameters) # void setFF_pointer( double * gridF, double * gridE ) @@ -87,7 +87,7 @@ def setFF_Epointer(gridE): lib.setFF_Epointer(gridE) -def setFF(cell=None, gridF=None, gridE=None): +def setFF(cell=None, gridF=None, gridE=None, parameters=None): n_ = None if gridF is not None: setFF_Fpointer(gridF) @@ -98,14 +98,14 @@ def setFF(cell=None, gridF=None, gridE=None): if cell is None: cell = np.array( [ - PPU.params["gridA"], - PPU.params["gridB"], - PPU.params["gridC"], + parameters.gridA, + parameters.gridB, + parameters.gridC, ] ).copy() if n_ is not None: print("setFF() n_ : ", n_) - setFF_shape(n_, cell) + setFF_shape(n_, cell, parameters=parameters) else: "Warrning : setFF shape not set !!!" @@ -133,15 +133,15 @@ def setFIRE(finc=1.1, fdec=0.5, falpha=0.99): lib.setTip.restype = None -def setTip(lRadial=None, kRadial=None, rPP0=None, kSpring=None): +def setTip(lRadial=None, kRadial=None, rPP0=None, kSpring=None, parameters=None): if lRadial is None: - lRadial = PPU.params["r0Probe"][2] + lRadial = parameters.r0Probe[2] if kRadial is None: - kRadial = PPU.params["krad"] / -PPU.eVA_Nm + kRadial = parameters.krad / -PPU.eVA_Nm if rPP0 is None: - rPP0 = np.array((PPU.params["r0Probe"][0], PPU.params["r0Probe"][1], 0.0)) + rPP0 = np.array((parameters.r0Probe[0], parameters.r0Probe[1], 0.0)) if kSpring is None: - kSpring = np.array((PPU.params["klat"], PPU.params["klat"], 0.0)) / -PPU.eVA_Nm + kSpring = np.array((parameters.klat, parameters.klat, 0.0)) / -PPU.eVA_Nm print(" IN setTip !!!!!!!!!!!!!! ") print(" lRadial ", lRadial) print(" kRadial ", kRadial) @@ -259,10 +259,10 @@ def computeD3Coeffs(Rs, iZs, iZPP, df_params): lib.getMorseFF.restype = None -def getMorseFF(Rs, REs, alpha=None): +def getMorseFF(Rs, REs, alpha=None, parameters=None): if alpha is None: - alpha = PPU.params["aMorse"] - print("getMorseFF: alpha: %g [1/A] ", alpha) + alpha = parameters.aMorse + print(f"getMorseFF: alpha: {alpha} [1/A] ") natom = len(Rs) lib.getMorseFF(natom, Rs, REs, alpha) diff --git a/tests/human_eye/test_TipForce.py b/tests/human_eye/test_TipForce.py index 2eebeb2b..85b3df32 100755 --- a/tests/human_eye/test_TipForce.py +++ b/tests/human_eye/test_TipForce.py @@ -5,7 +5,7 @@ import matplotlib.pyplot as plt import numpy as np -import ppafm.core as PPC +from ppafm import common, core spline_ini = """ 0.0 0.4 -0.1 @@ -17,6 +17,8 @@ # S = np.genfromtxt('TipRSpline.ini') S = np.genfromtxt(StringIO(spline_ini), skip_header=1) +parameters = common.PpafmParameters() + print("TipRSpline.ini overrides harmonic tip") xs = S[:, 0].copy() print("xs: ", xs) @@ -24,9 +26,9 @@ print("ydys: ", ydys) plt.plot(xs, ydys[:, 0], "o") plt.plot(xs, ydys[:, 0] + ydys[:, 1], ".") -PPC.setTipSpline(xs, ydys) +core.setTipSpline(xs, ydys) -PPC.setTip() +core.setTip() fs = np.zeros((60, 3)) r0 = np.array([0.0, 0.0, 0.5]) @@ -35,12 +37,12 @@ xs = np.array(list(range(len(fs)))) * dr[2] + r0[2] # print "xs=",xs -print(">>> PPC.test_force( 1, r0, dr, R, fs )") -PPC.test_force(1, r0, dr, R, fs) +print(">>> core.test_force( 1, r0, dr, R, fs )") +core.test_force(1, r0, dr, R, fs) plt.plot(xs, fs[:, 2]) -print(">>> PPC.test_force( 2, r0, dr, R, fs )") -PPC.test_force(2, r0, dr, R, fs) +print(">>> core.test_force( 2, r0, dr, R, fs )") +core.test_force(2, r0, dr, R, fs) plt.plot(xs, fs[:, 2]) # print "fs:", fs diff --git a/tests/test_io.py b/tests/test_io.py index 4a8e7816..fcaba059 100755 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -1,5 +1,4 @@ #!/usr/bin/env python3 - import os import numpy as np @@ -61,9 +60,10 @@ def test_parse_comment_ase(): def test_load_aims(): - from ppafm.common import params + from ppafm import common from ppafm.io import loadGeometry + parameters = common.PpafmParameters() temp_file = "temp_geom.in" geometry_str = """ lattice_vector 4.0 0.0 0.0 @@ -76,7 +76,7 @@ def test_load_aims(): with open(temp_file, "w") as f: f.write(geometry_str) - atoms, nDim, lvec = loadGeometry(temp_file, params=params) + atoms, nDim, lvec = loadGeometry(temp_file, parameters=parameters) assert np.allclose( atoms, @@ -90,7 +90,7 @@ def test_load_aims(): ] ), ) - assert np.allclose(nDim, params["gridN"]) + assert np.allclose(nDim, parameters.gridN) assert np.allclose(lvec, np.array([[0.0, 0.0, 0.0], [4.0, 0.0, 0.0], [1.0, 5.0, 0.0], [0.0, 0.0, 6.0]])) os.remove(temp_file)