Skip to content

Commit

Permalink
Fix issues that occured when running ppafm-generate-elff with tip's d…
Browse files Browse the repository at this point in the history
…ensity.
  • Loading branch information
yakutovicha committed Nov 5, 2024
1 parent 12a9034 commit b0aaf15
Show file tree
Hide file tree
Showing 3 changed files with 7 additions and 11 deletions.
14 changes: 5 additions & 9 deletions ppafm/HighLevel.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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"}:
Expand All @@ -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
Expand Down Expand Up @@ -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])
Expand Down
2 changes: 1 addition & 1 deletion ppafm/cli/generateDFTD3.py
Original file line number Diff line number Diff line change
Expand Up @@ -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__":
Expand Down
2 changes: 1 addition & 1 deletion ppafm/cli/generateElFF.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit b0aaf15

Please sign in to comment.