Update DFdriver

src/apps/dirac/DF.h

@@ -5,6 +5,7 @@
 #include <madness/mra/operator.h>
 #include <madness/constants.h>
 #include <madness/mra/nonlinsol.h>  // The kain solver
+#include <tuple>
 #include <vector>
 #include <math.h>
 #include <stdio.h>
@@ -110,12 +111,24 @@ class DF {
           //Creates the fermi nuclear potential from the molecule object. Also calculates the nuclear repulsion energy
           void make_fermi_potential(World& world, real_convolution_3d& op, real_function_3d& potential, double& nuclear_repulsion_energy);
 
+          //Creates the point nuclear potential from the molecule object
+          void make_point_potential(World& world, real_function_3d& potential);
+
+          //Creates the point nuclear potential from the molecule object. Also calculates the nuclear repulsion energy
+          void make_point_potential(World& world, real_function_3d& potential, double& nuclear_repulsion_energy);
+
           //Load balancing function
           void DF_load_balance(World& world, real_function_3d& Vnuc);
 
           //Does one full SCF iteration
-          bool iterate(World& world, real_function_3d& V, real_convolution_3d& op, real_function_3d& JandV, std::vector<Fcwf>& Kpsis, XNonlinearSolver<std::vector<Fcwf>, std::complex<double>, Fcwf_vector_allocator>& kainsolver, double& tolerance, int& iteration_number, double& nuclear_repulsion_energy);
-
+          std::tuple<bool, double, real_function_3d>
+          iterate(World &world, real_function_3d &V, real_convolution_3d &op,
+                  real_function_3d &JandV, std::vector<Fcwf> &Kpsis,
+                  XNonlinearSolver<std::vector<Fcwf>, std::complex<double>,
+                                   Fcwf_vector_allocator> &kainsolver,
+                  double &tolerance, int &iteration_number,
+                  double &nuclear_repulsion_energy, double &prev_energy,
+                  real_function_3d &prev_rho);
 
           //Runs the job specified in the input parameters
           void solve(World& world);

src/apps/dirac/DFParameters.h

@@ -27,6 +27,7 @@ namespace madness {
           int max_iter;                ///< Maximum number of iterations
           double small;                ///< Minimum length scale to be resolved
           double thresh;               ///< Accuracy criterion when truncating
+          double dconv;                ///< Accuracy criterion for charge density. Defaults to thresh.
           int k;                       ///< Number of legendre polynomials in scaling basis
           bool kain;                   ///< Turns on KAIN nonlinear solver 
           int maxsub;                  ///< Sets maximum subspace size for KAIN
@@ -59,6 +60,7 @@ namespace madness {
           , max_iter(20)
           , small(1e-5)
           , thresh(1e-6)
+          , dconv(1e-6)
           , k(8)
           , kain(false)
           , maxsub(10)
@@ -71,7 +73,7 @@ namespace madness {
           , nwchem(false)
           , lineplot(false)
           , no_compute(false)
-          , bohr_rad(52917.7211)
+          , bohr_rad(52917.7210544)  // bohr radius in fm from CODATA 2022
           , min_iter(2)
           , Krestricted(false)
           {}
@@ -105,6 +107,10 @@ namespace madness {
                     }
                     else if (s == "thresh"){
                          f >> thresh;
+                         dconv = thresh;
+                    }
+                    else if (s == "dconv"){
+                         f >> dconv;
                     }
                     else if (s == "k"){
                          f >> k;
@@ -175,6 +181,9 @@ namespace madness {
                if(nucleus == 1){
                     madness::print("                       Nucleus: fermi");
                }
+               else if (nucleus == 2) {
+                    madness::print("                       Nucleus: point");
+               }
                else{
                     madness::print("                       Nucleus: gaussian");
                }

src/apps/dirac/InitParameters.h

@@ -20,6 +20,7 @@ namespace madness{
           // Ground state parameters that are read in from archive
           std::string inFile;                      ///< Name of input archive to read in
           double Init_total_energy;                ///< Total energy of the nonrelativistic ground state
+          double conv_thresh;                      ///< Convergence threshold for MOs
           bool spinrestricted;                     ///< Indicates if input calc. was spin-restricted
           bool closed_shell;
           unsigned int num_occupied;               ///< Number of orbitals
@@ -82,13 +83,14 @@ namespace madness{
                     std::string localize_method;
 
                     input & version;
-                    input & Init_total_energy;              // double
+                    input & Init_total_energy;   // double
                     input & spinrestricted;      // bool
                     input & L;                   // double            box size
                     input & order;               // int               wavelet order
                     input & molecule;            // Molecule   
                     input & xc;
                     input & localize_method;
+                    input & conv_thresh;         // double
 
                     input & num_occupied;        // int
                     input & temp_energies;       // Tensor<double>    orbital energies
@@ -113,7 +115,7 @@ namespace madness{
                     complex_derivative_3d Dx(world,0);
                     complex_derivative_3d Dy(world,1);
                     complex_derivative_3d Dz(world,2);
-                    //double myc = 137.0359895; //speed of light in atomic units
+                    //double myc = 137.03599917697017; //speed of light in atomic units from CODATA 2022
                     std::complex<double> myi(0,1);
                     if(spinrestricted){ 
                          //If the calculation was spin-restricted in moldft, then we only have "spin-up" orbitals
@@ -127,7 +129,7 @@ namespace madness{
                          real_function_3d yfunc = real_factory_3d(world).f(myyfunc);
 
                          //Handle Kramers-restricted and unrestricted cases differently
-                         if(Krestricted){
+                         if(Krestricted || closed_shell){
                               //Loop over the occupied orbitals and convert
                               for(unsigned int i = 0; i < num_occupied; i++){
                                    //read in orbital

src/apps/dirac/README

@@ -38,13 +38,13 @@ List of interesting input options (from DFParameters.h)
           int maxsub;                  ///< Sets maximum subspace size for KAIN
           double maxrotn;              ///< maximum step allowed by kain
           bool restart;                ///< Indicates this is a restarted DF job
-          int nucleus;                 ///< Indicates which nucleus model to use (1 for fermi, anything else for Gaussian)
+          int nucleus;                 ///< Indicates which nucleus model to use (1 for fermi, 2 for point, anything else for Gaussian)
           bool do_save;                ///< Whether or not to save after each iteration. Defaults to true. Turn off with 'no_save'
           std::string savefile;        ///< Gives the file to save the archive each iteration Default: DFrestartdata (in working directory)
           int lb_iter;                 ///< How many iterations to load balance (after the initial load balancing)
           bool lineplot;               ///< Whether or not to make lineplots at the end of the job
           bool no_compute;             ///< If true, will skip all computation
-          double bohr_rad;             ///< bohr radius in fm (default: 52917.7211)
+          double bohr_rad;             ///< bohr radius in fm (default: 52917.7210544)
           int min_iter;                ///< minimum number of iterations (default: 2)
 
 (Kramers restricted does not seem to be working --- bug in the exchange piece?)

src/apps/dirac/fcwf.cc

@@ -186,7 +186,7 @@ Fcwf Fcwf::operator-=(const Fcwf& phi){
 //Returns the 2-norm of an initialized Fcwf
 double Fcwf::norm2(){
      MADNESS_ASSERT(m_initialized);
-     double c2 = 137.0359895*137.0359895; //speed of light in atomic units
+     double c2 = 137.03599917697017 * 137.03599917697017; //speed of light in atomic units from CODATA 2022
      std::complex<double> temp(0,0);
 
      temp += madness::inner(m_psi[0],m_psi[0]);
@@ -239,7 +239,7 @@ void Fcwf::truncate(){
 //Returns the inner product of two Fcwfs
 std::complex<double> Fcwf::inner(World& world, const Fcwf& phi) const{
      MADNESS_ASSERT(m_initialized && phi.getinitialize());
-     double c2 = 137.0359895*137.0359895; //speed of light in atomic units
+     double c2 = 137.03599917697017 * 137.03599917697017; //speed of light in atomic units from CODATA 2022
      std::complex<double> temp(0,0);
 
      temp += madness::inner(m_psi[0],phi.m_psi[0]);
@@ -294,15 +294,15 @@ Fcwf apply(World& world, complex_derivative_3d& D, const Fcwf& psi){
 //Returns the square modulus of an Fcwf, which is a real function
 real_function_3d squaremod(Fcwf& psi){
      MADNESS_ASSERT(psi.getinitialize());
-     double c2 = 137.0359895*137.0359895; //speed of light in atomic units
+     double c2 = 137.03599917697017 * 137.03599917697017; //speed of light in atomic units from CODATA 2022
      real_function_3d temp = abssq(psi[0]) + abssq(psi[1]) + abssq(psi[2]).scale(1.0/c2) + abssq(psi[3]).scale(1.0/c2);
      return temp;
 }
 
 //Returns the square modulus of the small component of an Fcwf, which is a real function
 real_function_3d squaremod_small(Fcwf& psi){
      MADNESS_ASSERT(psi.getinitialize());
-     double c2 = 137.0359895*137.0359895; //speed of light in atomic units
+     double c2 = 137.03599917697017 * 137.03599917697017; //speed of light in atomic units from CODATA 2022
      real_function_3d temp = (abssq(psi[2]) + abssq(psi[3])).scale(1.0/c2); //don't forget the factor of c^2
      return temp;
 }
@@ -317,7 +317,7 @@ real_function_3d squaremod_large(Fcwf& psi){
 //compute the function inner product between two Fcwfs. Result is a complex function.
 complex_function_3d inner_func(World& world, Fcwf& psi, Fcwf& phi){
      MADNESS_ASSERT(psi.getinitialize() && phi.getinitialize());
-     double c = 137.0359895; //speed of light in atomic units
+     double c = 137.03599917697017; //speed of light in atomic units from CODATA 2022
      std::vector<complex_function_3d> a(4);
      std::vector<complex_function_3d> b(4);
      for(unsigned int i = 0; i < 2; i++){
@@ -433,7 +433,7 @@ Tensor<std::complex<double>> matrix_inner(World& world, std::vector<Fcwf>& a, st
      unsigned int m = b.size();
      //MADNESS_ASSERT(n==m);
 
-     double c2 = 137.0359895*137.0359895; //speed of light in atomic units
+     double c2 = 137.03599917697017 * 137.03599917697017; //speed of light in atomic units from CODATA 2022
 
      //Reassign the vectors of Fcwfs to vectors of complex functions to facilitate use of vmra functions
      std::vector<complex_function_3d> a_1(n); //all first components of Fcwfs in input a

src/apps/dirac/relops.cc

@@ -7,7 +7,7 @@
 #include "../../../../mpfrc++-3/mpreal.h"
 
 static const double m = 1.0;
-static const double c = 137.0359895;
+static const double c = 137.03599917697017; // from CODATA 2022
 static const double mc2 = m*c*c;
 static const double PI = 3.14159265358979323846264338328;
 

src/apps/dirac/rk.cc

@@ -193,7 +193,7 @@ double rk_v(double Z) {
   if (NONREL) return 0.0;
   if (DK1) Z *= 0.5;
 
-    static const double c = 137.0359895;
+    static const double c = 137.03599917697017; // from CODATA 2022
     double v = -2.0*Z/(constants::pi*c);
     while (true) {
         double vnew = -2.0*std::atan(Z/(c*(v+1.0)))/constants::pi;