{"id":"23b06ea8-55a7-49c6-9b85-525858bd11ee","arxiv_id":"2608.08474","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"A single variational electrostatic framework, implemented in SPARC, gives exact open-boundary Dirichlet conditions and analytic energy, forces, and stress for isolated, 1D periodic, and 2D periodic systems under a uniform electric field.","lead":"This paper builds a unified way to handle electrostatics in real-space DFT for molecules, wires, and slabs, including under a uniform electric field, using boundary values computed from Green's functions. It adds a stress tensor expression and implements the method in the SPARC code, showing faster convergence with less vacuum than plane-wave codes.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Eq. (24) asserts an isotropic half-dipole electric-field stress with no derivation; the strain derivative of the field energy in Eqs. (5)/(22) suggests an anisotropic E_i μ_j term, so the claimed exact stress is unsupported.","rationale":"The reader's verdict is CONDITIONAL, and I agree with that verdict, but for a sharper reason. The reader's identified weakest assumption, that the electron density vanishes on and beyond the open boundaries, is standard, explicitly stated in Sec. 2, and supported by the exponential convergence tests; it is a practical limitation rather than an internal inconsistency. The more load-bearing gap is the stress tensor: Eq. (24) is asserted without derivation, and it appears structurally inconsistent with the field-energy term from which it should follow. The field-dipole energy ∫(ρ+b)x·E under an affine strain yields an anisotropic contribution E_i μ_j/|Ω|, with no 1/2 and no off-diagonal cancellation, whereas Eq. (24) gives a purely isotropic half-dipole term. If this is correct, the paper's claim to provide the first stress for such systems fails. The numerical validation against energy derivatives is the right idea, but the reported magnitudes and percentages do not constrain the term sharply enough, so a targeted test is needed. I am not asserting that Eq. (24) is false; the authors may have a derivation that was omitted from the preprint, and the surface terms in Eq. (3) could alter the naive affine-scaling argument. That is exactly why a concrete numerical check on a classical dipole and on the production code should settle it. The rest of the work is strong: the variational formulation is internally consistent, the Dirichlet expressions for 0D/1D/2D are carefully derived, the dipole-energy consistency check passes, and the agreement with Quantum ESPRESSO is convincing. Because the stress issue is testable and localized, the appropriate disposition remains CONDITIONAL rather than REJECT. If the test confirms the anisotropic prediction, the verdict should be moved to REJECT or a major revision; if Eq. (24) survives the test, the paper can be accepted after adding the missing derivation.","tokens_in":27425,"tokens_out":36534,"duration_ms":413303,"concrete_test":"Build a classical model of two opposite Gaussian pseudocharges ±q at ±d/2 along z in an orthorhombic cell with vacuum L, with no Kohn-Sham orbitals, and evaluate the functional Eq. (5) exactly. Compute the stress by numerical differentiation of this energy with respect to ε_xx, ε_yy, ε_zz at fixed E_z using central differences with step 1e-4. Compare with Eq. (24) and with the analytic affine-scaling result σ_zz = E_z q d/V, σ_xx=σ_yy=0. If the numerical derivative matches the affine result and not Eq. (24), the stress formula in the paper is wrong. Repeat in SPARC on the polycarbonitrile wire with a transverse field E_z = 0.01 Ha/(e·Bohr), comparing the analytical axial stress (Eq. 24) to finite-difference derivatives of Eq. (22) with respect to L_y; a mismatch of order E_z μ_z/V at this field would confirm the issue in the production code.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central exactness claim includes the stress tensor, and Eq. (24) is the only statement of the field-dependent stress: σ_αβ(R,E)=σ_αβ(R,0)+δ_αβ (1/(2|Ω|))∫(ρ+b)x·E. No derivation is provided. This is not a cosmetic gap: the proposed term has the wrong tensor structure for the energy from which it is supposed to follow. In Eq. (5)/(22), the field-dependent part of the electrostatic energy is the full dipole coupling ∫_Ω (ρ+b)x·E, not half. Under an affine strain ε, the coordinates entering this term transform as x→(1+ε)x, so the first-order change is ∫(ρ+b) ε_ij E_i x_j = E_i μ_j ε_ij. The corresponding stress contribution is (1/|Ω|) E_i μ_j, which is anisotropic and carries no factor 1/2. Eq. (24) instead gives an isotropic half-dipole term. For a 1D wire with a transverse field E_z, Eq. (24) predicts an axial stress change E_z μ_z/(2V), whereas the stated energy has no axial strain dependence unless μ_z changes with ε_yy; any such change is a piezoelectric response E_z ∂μ_z/∂ε_yy, not the dipole density itself. The reported stress validation against numerical energy derivatives is too weak to distinguish these: the field-induced stress in the test cases is small relative to the zero-field stress, and the percentage errors are quoted against unspecified stress magnitudes. The burden is on the paper to derive Eq. (24) or replace it; as written, the claimed exactness of the stress is unsupported and likely incorrect.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript develops a unified real-space electrostatics framework for Kohn–Sham density functional theory covering isolated, 1D-periodic, and 2D-periodic systems with open boundaries and a uniform applied electric field along the open directions. The central construction is the local variational functional of Eq. (3), whose stationarity is claimed to yield the Poisson equation for the total (electron plus pseudocharge) density, with Dirichlet values supplied by the analytical Green's-function/multipole expressions in Eqs. (9), (13), and (17). The paper further derives ground-state energy, atomic forces, and stress, implements them in SPARC, validates exponential vacuum convergence and agreement with Quantum ESPRESSO, and applies the framework to static polarizabilities and piezoelectric coefficients.","tokens_in":27697,"tokens_out":11414,"duration_ms":136002,"significance":"If the central claim holds, the framework provides a unified, parameter-free treatment of open-boundary electrostatics with applied fields, eliminating multipole-image errors and substantially reducing the vacuum required in real-space DFT calculations, while also supplying a stress tensor not currently available in other implementations. The strengths of the paper are its explicit variational derivation, the absence of fitted parameters in the boundary values, the exponential-convergence validation, the direct comparisons with an established plane-wave code, and the open data repository. These features make the work a potentially important contribution to real-space electronic structure methods. The main caveat is the stress expression in Eq. (24), which is asserted without derivation and appears inconsistent with the energy functional from which it is supposed to follow; because stress is part of the paper's central claim, this issue is load-bearing.","major_comments":[{"comment":"The field-dependent stress contribution is stated without derivation and its tensor structure is inconsistent with the energy functional from which it is supposed to follow. In Eqs. (5) and (22), the applied-field part of the electrostatic energy is the dipole coupling ∫_Ω (ρ+b) x·E. Under an affine strain ε, the coordinates transform as x → (1+ε)x, so the first-order change of this term is E_i μ_j ε_ij, with μ_j = ∫ (ρ+b) x_j. The corresponding stress contribution is the symmetric part of (1/|Ω|) E_i μ_j, which is anisotropic and carries no factor of 1/2. Equation (24) instead asserts an isotropic half-dipole term, δ_αβ (1/(2|Ω|)) ∫ (ρ+b) x·E. The numerical validation in Sec. 5.2 is too weak to distinguish the two: the errors are quoted as percentages of unspecified stress magnitudes, and the field-induced change is small relative to the zero-field stress. The authors should either derive Eq. (24) from Eq. (22), including all surface, boundary-value, and pseudocharge terms, or replace it; as written, the claimed exactness of the stress is unsupported.","section":"Sec. 3, Eq. (24)"},{"comment":"The force expression is obtained by differentiating the ground-state energy and invoking the Hellmann–Feynman theorem, but the admissible space of potentials in Eq. (3) has Dirichlet values φ0(R) + x·E that depend on atomic positions through φ0. For a variational problem with R-dependent constraints, the stationarity condition alone does not eliminate the derivative of the constraint. The manuscript does not show that the boundary terms involving ∂φ0/∂R_I vanish, and the asserted expression in Eq. (23) has no such terms. Please provide this step or state the conditions under which the contribution of the R-dependent Dirichlet values cancels. The agreement with Quantum ESPRESSO in Table 1 is reassuring, but it does not replace the missing derivation.","section":"Sec. 3, Eq. (23)"}],"minor_comments":[{"comment":"The phrase 'the the' should be corrected to 'the'.","section":"Sec. 2, after Eq. (6)"},{"comment":"The symbols x2D and θ are used in Eq. (13) and defined only afterward; please move the definitions before the equation or restructure the sentence for readability.","section":"Sec. 2.2, Eq. (13)"},{"comment":"The quantities b̃ and V_c appear in Eq. (22) but are defined only in the following text; please define them before first use.","section":"Sec. 3, Eq. (22)"},{"comment":"The stress validation reports percentage errors without giving the reference stress magnitudes. Please report the actual zero-field stresses and the field-induced changes so that the reader can assess the significance of the 0.1% level agreement.","section":"Sec. 5.2"},{"comment":"The convergence plots show multiple truncation-parameter curves, but the caption and legends do not identify which curve corresponds to which truncation parameter. Please add this information for reproducibility.","section":"Fig. 4"}],"recommendation":"major_revision","confidential_remarks":"The manuscript is a solid methods contribution, and the variational core, energy, and force validation appear sound. The main blocker is the stress expression in Eq. (24), which needs a proper derivation or correction, along with a targeted numerical test of the field-dependent stress. If the authors can resolve this, I would support publication in this journal. The literature review and credit to prior real-space electrostatics work are appropriate."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Quick take: the variational electrostatics is real and the numerics are convincing; the stress part needs serious work before I'd trust it.\n\nWhat's new: one local variational functional that handles 0D, 1D periodic, and 2D periodic systems with non-vanishing Dirichlet boundary conditions, includes a uniform field, and gives energy, forces, and stress. The Dirichlet expressions (spherical multipoles, cylindrical multipoles, 2D Fourier/dipole) are standard and carefully assembled. The implementation in SPARC with O(N) boundary evaluation is a genuine engineering achievement. Validation is strong: exponential convergence with vacuum size, 2x–8x less vacuum than plane-wave, agreement with Quantum ESPRESSO for energy/forces/polarization, and a clean energy–dipole consistency check.\n\nSoft spots. First, Eq. (24) is not derived. The claim that the field contribution to the stress is isotropic with a factor of 1/2 is not what you get from a simple strain derivative of the field-coupling term in the energy. Under strain, the coordinate in the integral ∫(ρ+b)x·E transforms as x→(1+ε)x, which gives a contribution symmetrizing to (1/2V)(E_i μ_j + E_j μ_i), not δ_ij (1/2V)E·μ. The numerical stress validation cannot distinguish these: the field-induced stress in the test cases is small compared to the zero-field stress, and the quoted percent errors do not reveal which tensor is being tested. So either Eq. (24) needs a real derivation, or the field-dependent stress claim should be withdrawn or qualified. As written, the asserted exactness of the stress is not established. Second minor issue: the text says e311=e322=0.0012 e/Bohr while the Fig. 7 caption says 0.0011 e/Bohr; one is a typo. Neither issue breaks the central electrostatics.\n\nThe weaker assumption—that the charge density vanishes at the open boundary—is standard, clearly stated, and not a concern. The variational derivation in Sec. 2 is internally consistent.\n\nBottom line: this is a useful methods paper for real-space DFT, and the electrostatics part is solid. The stress tensor needs a fix before final. I'd send it to peer review, with an explicit request that a referee check the derivation of Eq. (24). The central contribution is real and the flaw is localized.","headline":"Unified open-boundary electrostatics in real-space DFT is a genuine, well-tested advance, but the field-dependent stress tensor (Eq. 24) is asserted without derivation and looks wrong on tensorial grounds.","tokens_in":28315,"tokens_out":11291,"would_cite":true,"duration_ms":140151,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":["71.15.Mb"],"model":"deepseek-v4-flash","headline":"One local variational functional with Green's-function Dirichlet values yields the exact electrostatic potential, energy, forces, and stress for isolated, 1D, and 2D periodic systems in DFT, with or without an applied uniform field.","keywords":["real-space density functional theory","open-boundary electrostatics","Dirichlet boundary conditions","Poisson equation","uniform electric field","multipole expansion","static polarizability","piezoelectric coefficients"],"falsifier":"Run the method on a neutral molecule at progressively smaller vacuum sizes until the electron density at the boundary is no longer negligible, and compare the energy and forces against a well-converged plane-wave calculation: the exponential-convergence curve should flatten or diverge at the point where the density support reaches the boundary. Equivalently, apply the formulation to a charged molecule, for which the zeroth multipole term is dropped: the computed energy will drift from the exact value as the cell grows, exposing the neutrality condition as the binding assumption.","tokens_in":27154,"feed_emoji":"⚡","tokens_out":10146,"duration_ms":94754,"temperature":0.7,"pith_summary":"This paper claims that the open-boundary electrostatics of isolated molecules, 1D-periodic wires, and 2D-periodic slabs can be handled in real-space density functional theory by a single local variational functional, Eq. (3), whose maximizer is the exact electrostatic potential of the total charge density plus any applied uniform electric field along the open directions. The required Dirichlet boundary values are written analytically from Green's functions: spherical multipoles for isolated systems, cylindrical multipoles with modified Bessel functions for wires, and a dipole step with in-plane Fourier terms for slabs. If the claim is right, the vacuum region only needs to be large enough for the electron density to decay, not large enough to suppress multipole-image errors, so real-space DFT would need substantially less vacuum than plane-wave supercell calculations at the same accuracy. The paper verifies exponential convergence with vacuum size, agreement with a plane-wave reference, stresses that match numerical derivatives of the energy, and polarizabilities and piezoelectric coefficients that agree with published values.","feed_headline":"One functional unifies open-boundary electrostatics in DFT","feed_subtitle":"Real-space DFT can now drop excess vacuum around molecules, wires, and slabs without losing accuracy.","key_machinery":"The load-bearing object is the local variational functional of Eq. (3), a maximum over $\\phi \\in \\mathcal{U} \\subset H^1(\\Omega)$, the affine space of functions periodic along the periodic directions and fixed to $\\phi_0(\\mathbf{x}) + \\mathbf{x}\\cdot\\mathbf{E}$ on the open boundaries, with boundary surface integrals that absorb the slow or growing decay of the potential along the open directions. This functional replaces the nonlocal Coulomb kernel with a differential operator whose inversion reproduces the long-range interaction exactly, so stationarity of the functional guarantees consistency between the energy and its derivatives. The Dirichlet values $\\phi_0$ come from the Green's functions of each geometry: spherical harmonics for the 0D case, cylindrical multipole moments together with zeroth-order modified Bessel functions $K_0$ for the 1D case, and a dipole step term with exponentially screened in-plane Fourier components for the 2D case, each valid where the charge density has vanished and convergent for a charge-neutral system.","core_discovery":"The central discovery is a maximization formulation of electrostatics that makes no assumption that the potential or its gradient vanishes on the open boundaries. The paper defines an affine space of functions $\\phi$ that are periodic along the periodic directions and take prescribed values $\\phi_0(\\mathbf{x}) + \\mathbf{x}\\cdot\\mathbf{E}$ on the open faces, where $\\phi_0$ is the potential of the total charge density alone. Eq. (3) expresses the electrostatic energy as the maximum over this space of a local functional containing surface integrals on the open boundaries; the maximizer solves the Poisson equation $-\\frac{1}{4\\pi}\\nabla^2\\phi = \\rho + b$ for the total charge density, and substituting it back yields the closed-form energy, Eq. (5). The Dirichlet data $\\phi_0$ are derived by a Green's-function ansatz for each dimensionality: the spherical multipole series of Eq. (9) for isolated systems, the cylindrical multipole and Bessel-function series of Eq. (13) for 1D periodic systems, and the dipole step plus exponentially screened in-plane Fourier series of Eq. (17) for 2D periodic systems. The force and stress follow by differentiation, with the electric field entering only through $\\phi$; the stress, Eq. (24), is the zero-field stress plus an isotropic diagonal term, and because no existing code provides stresses for such systems they are validated against numerical derivatives of the energy.","pith_inferences":["Charged systems are the natural next stress test: the authors explicitly defer them to future work, and the omission of the $\\ell = 0$ term in Eq. (9) suggests a neutralizing-background term would need to be added to $\\phi_0$ before the formulation extends to net-charged species.","The built-in consistency check between the energy derivative and the dipole moment could be turned into an automated diagnostic for the truncation parameters, flagging when $\\ell_{\\max}$, $m_{\\max}$, or $Q^{\\max}_{mn}$ is too small for the system's multipole content.","Because the electric field enters only through $\\phi$, the first-order variation of the boundary values should slot naturally into real-space density functional perturbation theory, a path the authors name as future work.","The Green's-function ansatz is the only part that changes with geometry, so the machinery could likely be re-derived for cyclic and helical symmetry, giving open-boundary electrostatics for bent and twisted nanostructures."],"forward_implications":["One framework covers isolated, 1D-periodic, and 2D-periodic systems, with an applied uniform electric field entering entirely through the electrostatic potential rather than through an explicit change to the Kohn–Sham Hamiltonian.","Energy, forces, and stresses converge exponentially with vacuum size, reaching roughly $10^{-6}$ Ha/atom and $10^{-5}$ Ha/Bohr at about 10 Bohr of vacuum, several times less vacuum than the plane-wave reference calculations needed in the test cases.","The stress tensor is available for systems with a non-vanishing potential on the open boundaries, so cell relaxation and equation-of-state studies need no finite-difference energy derivatives for such systems.","Static polarizabilities and piezoelectric coefficients of polar low-dimensional systems follow from the same machinery without separate correction schemes, matching published values.","The additional cost of evaluating the Dirichlet boundary values is negligible relative to the Poisson solve, so the parallel scalability of the real-space solver is retained.","The boundary conditions also accommodate fixed-potential electrodes and could be matched to bulk boundary conditions for semi-infinite surface calculations."],"supporting_citations":[{"why":"Supplies the Green's functions for a periodic line of charges and a doubly periodic charge distribution that seed the 1D and 2D Dirichlet boundary-value derivations.","marker":"[67]"},{"why":"Provides the spherical multipole expansion for isolated systems that the 0D boundary condition of Eq. (9) builds on.","marker":"[61]"},{"why":"Derives cylindrical multipole Dirichlet values for 1D periodic real-space systems, the direct antecedent of Eq. (13).","marker":"[62]"},{"why":"Gives the asymptotic potential for 2D slabs used previously as the open-boundary condition, the antecedent of Eq. (17).","marker":"[63]"},{"why":"Prior real-space slab formulation with a uniform out-of-plane electric field under bias, which the present work extends to multipole effects and all dimensionalities.","marker":"[64]"},{"why":"The prior real-space electrostatic formulation in the code, recovered exactly when the electric field and boundary values vanish.","marker":"[48]"},{"why":"Provides the zero-field stress tensor that the electric-field stress expression of Eq. (24) extends.","marker":"[75]"},{"why":"The plane-wave code used as the accuracy reference for energy, forces, and polarization comparisons.","marker":"[10]"},{"why":"The real-space code in which the formulation and the Dirichlet boundary conditions are implemented.","marker":"[65]"}],"fun_headline_variants":["DFT electrostatics: less vacuum, same accuracy","Maximize electrostatics, minimize vacuum in DFT","One functional, all boundaries: DFT electrostatics unified","No excess vacuum: unified open-boundary DFT electrostatics","SPARC cuts vacuum with unified open-boundary electrostatics"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the total charge density (electrons plus ionic pseudocharge) is exactly zero on and beyond the open boundary faces and that the system is charge neutral, so the multipole and Fourier series for the Dirichlet values converge; if the vacuum is so small that the density has not decayed, an atom crosses the boundary, or the system carries a net charge, the boundary values are no longer exact and the claimed exactness collapses.","fun_headline_variants_meta":{"raw":{"variants":["DFT electrostatics: less vacuum, same accuracy","Maximize electrostatics, minimize vacuum in DFT","One functional, all boundaries: DFT electrostatics unified","No excess vacuum: unified open-boundary DFT electrostatics","SPARC cuts vacuum with unified open-boundary electrostatics"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001103,"raw_usage":{"total_tokens":4654,"prompt_tokens":1055,"completion_tokens":3599,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":671,"completion_tokens_details":{"reasoning_tokens":3517}},"tokens_in":671,"tokens_out":3599,"duration_ms":24867,"temperature":1.0,"reasoning_tokens":3517,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T04:36:46.032841+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the method on a neutral molecule at progressively smaller vacuum sizes until the electron density at the boundary is no longer negligible, and compare the energy and forces against a well-converged plane-wave calculation: the exponential-convergence curve should flatten or diverge at the point where the density support reaches the boundary. Equivalently, apply the formulation to a charged molecule, for which the zeroth multipole term is dropped: the computed energy will drift from the exact value as the cell grows, exposing the neutrality condition as the binding assumption.","supporting_citations":[{"cited_title":"Lennard-Jones, B","cited_arxiv_id":null,"evidence_quote":"Supplies the Green's functions for a periodic line of charges and a doubly periodic charge distribution that seed the 1D and 2D Dirichlet boundary-value derivations."},{"cited_title":"Hirose, T","cited_arxiv_id":null,"evidence_quote":"Provides the spherical multipole expansion for isolated systems that the 0D boundary condition of Eq. (9) builds on."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Derives cylindrical multipole Dirichlet values for 1D periodic real-space systems, the direct antecedent of Eq. (13)."},{"cited_title":"Natan, A","cited_arxiv_id":null,"evidence_quote":"Gives the asymptotic potential for 2D slabs used previously as the open-boundary condition, the antecedent of Eq. (17)."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Prior real-space slab formulation with a uniform out-of-plane electric field under bias, which the present work extends to multipole effects and all dimensionalities."},{"cited_title":"Ghosh, P","cited_arxiv_id":null,"evidence_quote":"The prior real-space electrostatic formulation in the code, recovered exactly when the electric field and boundary values vanish."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the zero-field stress tensor that the electric-field stress expression of Eq. (24) extends."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"The real-space code in which the formulation and the Dirichlet boundary conditions are implemented."}],"review_version":1}