From b0aaf1522183f10f9c1924d3eaeef2cf533b94cc Mon Sep 17 00:00:00 2001 From: Aliaksandr Yakutovich Date: Tue, 5 Nov 2024 13:46:52 +0100 Subject: [PATCH] Fix issues that occured when running ppafm-generate-elff with tip's density. --- ppafm/HighLevel.py | 14 +++++--------- ppafm/cli/generateDFTD3.py | 2 +- ppafm/cli/generateElFF.py | 2 +- 3 files changed, 7 insertions(+), 11 deletions(-) diff --git a/ppafm/HighLevel.py b/ppafm/HighLevel.py index 2969bd90..70012c3f 100755 --- a/ppafm/HighLevel.py +++ b/ppafm/HighLevel.py @@ -295,7 +295,7 @@ 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 @@ -303,7 +303,7 @@ def computeDFTD3(input_file, df_params="PBE", geometry_format=None, save_format= 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]) diff --git a/ppafm/cli/generateDFTD3.py b/ppafm/cli/generateDFTD3.py index ece16bbe..8b59c500 100755 --- a/ppafm/cli/generateDFTD3.py +++ b/ppafm/cli/generateDFTD3.py @@ -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__": diff --git a/ppafm/cli/generateElFF.py b/ppafm/cli/generateElFF.py index f0f5fb25..e3253e5b 100755 --- a/ppafm/cli/generateElFF.py +++ b/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)