Electrostatic potential calculation

Dear Shankar,

I plan to use jdftx to calculate the distribution of electrostatic potential on the surface of two-dimensional material TM@BN under neutral implicit solvation and applied potential U, in order to expect to obtain information about interfacial double layers deviating from PZC. But the problem is that I according to the tutorial: http://jdftx.org/OxideSurfaces.html, the electrostatic potential distribution of TM@BN is calculated and reported in the literature (Potential Effects on the Catalytic Mechanisms of OER and ORR, J. Phys. Chem. C 2023, 127, 16346−16356) of similar electrostatic potential distribution of the two-dimensional material has obvious different, and there is a positive electrostatic potential at the plane of TM@BN. Electrostatic potentials

I was puzzled by this result, and I ended by attaching details of the calculation and result processing. Can you point out any errors in my calculation and processing? Attached is the calculation command I submitted according to the tutorial,

#jdftx command
ion-species  SG15/$ID_ONCV_PBE-1.0.upf  
include lattice.lattice 
include in.ionpos                                       
coords-type lattice                     

elec-ex-corr gga-PBEsol
fluid LinearPCM
pcm-variant CANDLE
fluid-solvent H2O
fluid-cation Na+ 1.
fluid-anion F- 1.

coulomb-interaction Slab  001                                               
coulomb-truncation-embed 0  0  1                                            
elec-cutoff  30                                                            
kpoint-folding 2 2 1                                                        
elec-smearing Fermi 0.01
spintype z-spin                                       
elec-initial-magnetization  3  no                                               
electronic-minimize  \                                                      
    energyDiffThreshold  1e-06

dump-name surface.$VAR                                              
initial-state surface.$VAR                                              

dump End ElecDensity                                       
dump End EigStats                              
dump End Dtot BoundCharge

and plot the data with a script and the resulting picture, where TM@BN is located at z=0 angstroms.

#python script
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator    
x_major_locator=MultipleLocator(1.5)              

ha = 27.211396641308
A=0.52918       #https://www.convertunits.com/from/bohr/to/angstrom
S = [   120  ,120  ,140 ]            #Sample count (fftbox from log file), 
#Lattice vector lengths (bohrs),Respectively sqrt(2)a, c and chosen height with spacing
L = [ 23.73549643, 23.73549643, 28.3457424695 ]    
z = np.arange(S[2])*L[2]/S[2]  #z-coordinates of grid points

def readPlanarlyAveraged(fileName, S):
    out = np.fromfile(fileName, dtype=np.float64)
    out = np.reshape(out, S)   #Reshape data
    out = np.mean(out, axis=1) #y average
    out = np.mean(out, axis=0) #x average
    return out

CCM = readPlanarlyAveraged('surfaceCCM.d_tot',S)*ha
dShift = CCM #-(-0.135926)*ha

print ("VBMshift =", dShift[0])  #Report VBM shift in eV
plt.plot(z*A, dShift);     #Plot with unit conversions
plt.xlabel('z [Angstroms]')
plt.ylabel('Potential [eV]')
ax=plt.gca()                                     
ax.xaxis.set_major_locator(x_major_locator)     
plt.xlim([0, 15])     
plt.show()

Best, lyy

shankar1729 / jdftx

Electrostatic potential calculation #315