Feature: periodic self-consistent continuum solvation (SCCS) Step2: SCCS with PCC 0D/2D

source/Makefile.Objects

@@ -284,6 +284,8 @@ OBJS_ELECSTAT=elecstate.o\
     pcc_2d.o\
     pot_pcc.o\
     pot_sccs.o\
+    sccs_pcc_0d.o\
+    sccs_pcc_2d.o\
     elecstate_op.o\
     efield.o\
     gatefield.o\
@@ -631,10 +633,14 @@ OBJS_SURCHEM=surchem.o\
     sol_force.o\
 
 OBJS_SCCS=sccs_cavity.o\
+    sccs_cav_deriv.o\
     sccs_functional.o\
     sccs_parameters.o\
     sccs_ion_charge.o\
     sccs_ion_force.o\
+    sccs_coul_0d.o\
+    sccs_coul_2d.o\
+    sccs_moments.o\
     sccs_pw_coulomb.o\
     sccs_response.o
 

source/source_estate/CMakeLists.txt

@@ -22,6 +22,8 @@ list(APPEND objects
     module_pot/pot_local.cpp
     module_pot/pot_pcc.cpp
     module_pot/pot_sccs.cpp
+    module_pot/sccs_pcc_0d.cpp
+    module_pot/sccs_pcc_2d.cpp
     module_pot/potential_new.cpp
     module_pot/potential_types.cpp
     module_pot/pot_surchem.cpp

source/source_estate/module_pot/pot_sccs.cpp

@@ -1,4 +1,6 @@
 #include "pot_sccs.h"
+#include "sccs_pcc_0d.h"
+#include "sccs_pcc_2d.h"
 
 #include "source_base/parallel_reduce.h"
 #include "source_base/timer.h"
@@ -8,6 +10,7 @@
 #include "source_hamilt/module_sccs/sccs_ion_charge.h"
 #include "source_hamilt/module_sccs/sccs_ion_force.h"
 #include "source_hamilt/module_sccs/sccs_response.h"
+#include "source_hamilt/module_sccs/sccs_pw_coulomb.h"
 #include "source_io/module_parameter/input_parameter.h"
 
 #include <cmath>
@@ -16,6 +19,13 @@
 
 namespace elecstate
 {
+namespace
+{
+// Charge tolerance in e, the same as the PotPcc electron-count check: far
+// above the grid round-off (~1e-12 e) and far below any physical net charge.
+const double charge_tolerance = 1e-6;
+} // namespace
+
 void sccs_cfg_input(const Input_para& input,
                     ModuleSccs::SccsConfig& config,
                     ModuleSccs::PolarizationSolverParameters& solver)
@@ -35,6 +45,19 @@ void sccs_cfg_input(const Input_para& input,
         config = ModuleSccs::preset_config(preset);
     }
     config.surface_regularization = input.sccs_surface_eta;
+    if (input.assume_isolated == "pcc_0d")
+    {
+        config.boundary = ModuleSccs::Boundary::Pcc0d;
+    }
+    else if (input.assume_isolated == "pcc_2d")
+    {
+        config.boundary = ModuleSccs::Boundary::Pcc2d;
+        config.pcc_2d_axis = input.pcc_2d_axis;
+    }
+    else
+    {
+        config.boundary = ModuleSccs::Boundary::Periodic;
+    }
     solver.max_iterations = input.sccs_maxiter;
     solver.tolerance_rms = input.sccs_tol_rms;
     solver.tolerance_max = input.sccs_tol_max;
@@ -48,9 +71,14 @@ std::string check_charge(const ModuleSccs::SccsConfig& config, const UnitCell& c
         ionic_charge += cell.atoms[it].ncpp.zv * cell.atoms[it].na;
     }
     const double net_charge = ionic_charge - electron_count;
-    if (config.cavity.epsilon_bulk <= 1.0 || std::abs(net_charge) <= 1e-6) { return std::string(); }
+    const bool periodic = config.boundary == ModuleSccs::Boundary::Periodic;
+    if (!periodic || config.cavity.epsilon_bulk <= 1.0 || std::abs(net_charge) <= charge_tolerance)
+    {
+        return std::string();
+    }
     return "charged SCCS with periodic boundaries: the periodic Poisson solver drops the G = 0 component "
-           "of the net charge, so the energy depends on the cell size.";
+           "of the net charge, so the energy depends on the cell size. Use assume_isolated pcc_0d or pcc_2d "
+           "for converged charged energies.";
 }
 
 PotSccs::PotSccs(const ModulePW::PW_Basis* basis,
@@ -87,22 +115,36 @@ void PotSccs::cal_v_eff(const Charge* charge, const UnitCell* cell, ModuleBase::
     double valence_charge = 0.0;
     for (const unitcell::AtomData& atom : atoms) { valence_charge += atom.valence_charge; }
     const double ionic_error = grid_charges[1] * dv - valence_charge;
-    if (std::abs(ionic_error) > 1e-6)
+    if (std::abs(ionic_error) > charge_tolerance)
     {
         ModuleBase::WARNING_QUIT("PotSccs::cal_v_eff", "SCCS ionic density normalization does not match the valence charge");
     }
-    if (std::abs(electron_error) > 1e-6)
+    // With PCC the PCC component reports the same mismatch.
+    if (config_.boundary == ModuleSccs::Boundary::Periodic && std::abs(electron_error) > charge_tolerance)
     {
         std::ostringstream message;
         message << "SCCS grid electron count differs from the expected value by " << electron_error
                 << " e; the solute charge uses the grid density";
         ModuleBase::WARNING("PotSccs::cal_v_eff", message.str());
     }
+    std::unique_ptr<ModuleSccs::CoulombOperator> coulomb;
+    if (config_.boundary == ModuleSccs::Boundary::Pcc0d)
+    {
+        pcc_0d_operator(*cell, basis, atoms, coulomb);
+    }
+    else if (config_.boundary == ModuleSccs::Boundary::Pcc2d)
+    {
+        pcc_2d_operator(*cell, basis, atoms, config_.pcc_2d_axis, coulomb);
+    }
+    else
+    {
+        coulomb.reset(new ModuleSccs::PeriodicCoulombOperator(basis, cell->tpiba));
+    }
     ModuleSccs::SccsResponse response;
     ModuleSccs::solve_response(density, solute_charge, config_.cavity, solver_, restart_potential_, basis,
-                               cell->tpiba, response);
+                               cell->tpiba, *coulomb, response);
     ModuleSccs::FunctionalResult functional;
-    ModuleSccs::eval_functional(solute_charge, response, config_, basis, cell->tpiba, functional);
+    ModuleSccs::eval_functional(solute_charge, response, config_, basis, cell->tpiba, *coulomb, functional);
     electrostatic_rydberg_ = 2.0 * functional.reaction_energy;
     non_electrostatic_rydberg_ = 2.0 * (functional.surface_energy + functional.volume_energy);
     electrostatic_potential_.resize(basis.nrxx);

source/source_estate/module_pot/sccs_pcc_0d.cpp

@@ -0,0 +1,58 @@
+#include "sccs_pcc_0d.h"
+
+#include "source_base/tool_quit.h"
+#include "source_basis/module_pw/pw_basis.h"
+#include "source_basis/module_pw/pw_grid_geometry.h"
+#include "source_cell/cell_geometry.h"
+#include "source_cell/cell_tools.h"
+#include "source_cell/unitcell.h"
+#include "source_hamilt/module_sccs/sccs_coul_0d.h"
+
+namespace elecstate
+{
+void pcc_0d_geometry(const UnitCell& cell,
+                     const std::vector<unitcell::AtomData>& atoms,
+                     unitcell::OrthogonalCell& geometry,
+                     Pcc0dParameters& parameters)
+{
+    // Same orthogonality and equal-edge tolerance as PotPcc.
+    const double cell_tolerance = 1e-10;
+    const bool orthogonal = unitcell::make_orthogonal_cell(cell.latvec, cell.lat0, cell_tolerance, geometry);
+    const bool cubic = orthogonal && make_pcc_0d_parameters(geometry, cell_tolerance, parameters);
+    if (!cubic)
+    {
+        ModuleBase::WARNING_QUIT("pcc_0d_geometry", "PCC 0D requires an equal-edge cubic cell");
+    }
+    std::vector<ModuleBase::Vector3<double>> positions;
+    std::vector<double> masses;
+    for (const unitcell::AtomData& atom : atoms)
+    {
+        positions.push_back(atom.position);
+        masses.push_back(atom.mass);
+    }
+    geometry.origin = unitcell::weighted_center(positions, masses, geometry);
+}
+
+void pcc_0d_relative(const unitcell::OrthogonalCell& geometry,
+                     std::vector<ModuleBase::Vector3<double>>& positions)
+{
+    for (ModuleBase::Vector3<double>& position : positions)
+    {
+        position = unitcell::relative_position(position, geometry);
+    }
+}
+
+void pcc_0d_operator(const UnitCell& cell,
+                     const ModulePW::PW_Basis& basis,
+                     const std::vector<unitcell::AtomData>& atoms,
+                     std::unique_ptr<ModuleSccs::CoulombOperator>& coulomb)
+{
+    unitcell::OrthogonalCell geometry;
+    Pcc0dParameters parameters;
+    pcc_0d_geometry(cell, atoms, geometry, parameters);
+    std::vector<ModuleBase::Vector3<double>> relative_positions;
+    ModulePW::grid_positions(basis, cell.latvec, cell.lat0, relative_positions);
+    pcc_0d_relative(geometry, relative_positions);
+    coulomb.reset(new ModuleSccs::Pcc0dCoulombOperator(basis, cell.tpiba, relative_positions, parameters));
+}
+} // namespace elecstate

source/source_estate/module_pot/sccs_pcc_0d.h

@@ -0,0 +1,48 @@
+#ifndef SCCS_PCC_0D_ADAPTER_H
+#define SCCS_PCC_0D_ADAPTER_H
+
+#include <memory>
+#include <vector>
+
+class UnitCell;
+namespace ModuleBase
+{
+template <typename T> class Vector3;
+}
+namespace ModulePW
+{
+class PW_Basis;
+}
+namespace ModuleSccs
+{
+class CoulombOperator;
+}
+namespace unitcell
+{
+struct AtomData;
+struct OrthogonalCell;
+}
+namespace elecstate
+{
+struct Pcc0dParameters;
+
+// Cubic PCC cell with the mass-weighted ionic origin, as PotPcc; a cell that
+// is not an equal-edge cube stops the run.
+void pcc_0d_geometry(const UnitCell& cell,
+                     const std::vector<unitcell::AtomData>& atoms,
+                     unitcell::OrthogonalCell& geometry,
+                     Pcc0dParameters& parameters);
+
+// Minimum-image displacements of Cartesian positions from the origin, in place.
+void pcc_0d_relative(const unitcell::OrthogonalCell& geometry,
+                     std::vector<ModuleBase::Vector3<double>>& positions);
+
+// Assemble an SCCS grid-charge operator using the existing cell geometry and
+// mass-weighted ionic origin. An unsupported cell stops the run.
+void pcc_0d_operator(const UnitCell& cell,
+                     const ModulePW::PW_Basis& basis,
+                     const std::vector<unitcell::AtomData>& atoms,
+                     std::unique_ptr<ModuleSccs::CoulombOperator>& coulomb);
+} // namespace elecstate
+
+#endif

source/source_estate/module_pot/sccs_pcc_2d.cpp

@@ -0,0 +1,65 @@
+#include "sccs_pcc_2d.h"
+
+#include "source_base/tool_quit.h"
+#include "source_basis/module_pw/pw_basis.h"
+#include "source_basis/module_pw/pw_grid_geometry.h"
+#include "source_cell/cell_geometry.h"
+#include "source_cell/cell_tools.h"
+#include "source_cell/unitcell.h"
+#include "source_hamilt/module_sccs/sccs_coul_2d.h"
+
+namespace elecstate
+{
+void pcc_2d_geometry(const UnitCell& cell,
+                     const std::vector<unitcell::AtomData>& atoms,
+                     int open_axis,
+                     unitcell::SlabCell& geometry,
+                     Pcc2dParameters& parameters)
+{
+    // Same alignment tolerance as PotPcc, loose enough for lattice vectors
+    // typed with rounded digits.
+    const double alignment_tolerance = 1e-6;
+    const bool perpendicular = unitcell::make_slab_cell(cell.latvec, cell.lat0, open_axis, alignment_tolerance,
+                                                        geometry);
+    if (!perpendicular)
+    {
+        ModuleBase::WARNING_QUIT("pcc_2d_geometry",
+                                 "PCC 2D requires the open lattice vector to be perpendicular to the periodic plane");
+    }
+    parameters = make_pcc_2d_parameters(geometry);
+    std::vector<ModuleBase::Vector3<double>> positions;
+    std::vector<double> masses;
+    for (const unitcell::AtomData& atom : atoms)
+    {
+        positions.push_back(atom.position);
+        masses.push_back(atom.mass);
+    }
+    geometry.origin = unitcell::weighted_center(positions, masses, geometry);
+}
+
+void pcc_2d_project(const unitcell::SlabCell& geometry,
+                    std::vector<ModuleBase::Vector3<double>>& positions)
+{
+    for (ModuleBase::Vector3<double>& position : positions)
+    {
+        const double coordinate = unitcell::relative_coordinate(position, geometry);
+        position = geometry.normal * coordinate;
+    }
+}
+
+void pcc_2d_operator(const UnitCell& cell,
+                     const ModulePW::PW_Basis& basis,
+                     const std::vector<unitcell::AtomData>& atoms,
+                     int open_axis,
+                     std::unique_ptr<ModuleSccs::CoulombOperator>& coulomb)
+{
+    unitcell::SlabCell geometry;
+    Pcc2dParameters parameters;
+    pcc_2d_geometry(cell, atoms, open_axis, geometry, parameters);
+    std::vector<ModuleBase::Vector3<double>> projected_positions;
+    ModulePW::grid_positions(basis, cell.latvec, cell.lat0, projected_positions);
+    pcc_2d_project(geometry, projected_positions);
+    coulomb.reset(new ModuleSccs::Pcc2dCoulombOperator(basis, cell.tpiba, projected_positions, geometry.normal,
+                                                       parameters));
+}
+} // namespace elecstate

source/source_estate/module_pot/sccs_pcc_2d.h

@@ -0,0 +1,51 @@
+#ifndef SCCS_PCC_2D_ADAPTER_H
+#define SCCS_PCC_2D_ADAPTER_H
+
+#include <memory>
+#include <vector>
+
+class UnitCell;
+namespace ModuleBase
+{
+template <typename T> class Vector3;
+}
+namespace ModulePW
+{
+class PW_Basis;
+}
+namespace ModuleSccs
+{
+class CoulombOperator;
+}
+namespace unitcell
+{
+struct AtomData;
+struct SlabCell;
+}
+namespace elecstate
+{
+struct Pcc2dParameters;
+
+// Slab PCC cell with the mass-weighted ionic origin along the normal, with
+// the PotPcc alignment tolerance; a tilted open axis stops the run.
+void pcc_2d_geometry(const UnitCell& cell,
+                     const std::vector<unitcell::AtomData>& atoms,
+                     int open_axis,
+                     unitcell::SlabCell& geometry,
+                     Pcc2dParameters& parameters);
+
+// Replace Cartesian positions, in place, by the normal times their
+// minimum-image coordinate from the slab origin.
+void pcc_2d_project(const unitcell::SlabCell& geometry,
+                    std::vector<ModuleBase::Vector3<double>>& positions);
+
+// Assemble an SCCS grid-charge operator using the existing cell geometry and
+// mass-weighted ionic origin. An unsupported cell stops the run.
+void pcc_2d_operator(const UnitCell& cell,
+                     const ModulePW::PW_Basis& basis,
+                     const std::vector<unitcell::AtomData>& atoms,
+                     int open_axis,
+                     std::unique_ptr<ModuleSccs::CoulombOperator>& coulomb);
+} // namespace elecstate
+
+#endif