Paul's July 2026 updates

include/becke.h

@@ -27,7 +27,7 @@ void becke_weight_3c(int gs, float* grid1, float* wt1, float* grid2, float* wt2,
 void becke_weight_2c(int gs, float* grid1, float* wt1, float* grid2, float* wt2,
                      int Z1, int Z2, float A2, float B2, float C2);
 void becke_weight_2d(int gs, double* grid1, double* wt1, double* grid2, double* wt2,
-                     double zeta1, double zeta2, double A2, double B2, double C2);
+                     bool is_Z, double zeta1, double zeta2, double A2, double B2, double C2);
 
 float becke_ar(float r1, float r2);
 double becke_ard(double r1, double r2);

include/cint.h

@@ -4,7 +4,8 @@
  * Parameters and function signature for libcint.
  */
 
-#define CINT_VERSION    4.4.0
+#define CINT_VERSION            "6.1.3"
+#define CINT_SOVERSION          6
 
 /* #undef I8 */
 #ifdef I8
@@ -49,9 +50,8 @@
 #define PTR_COORD       1
 #define NUC_MOD_OF      2
 #define PTR_ZETA        3
-#define PTR_FRAC_CHARGE 3
-#define RESERVE_ATMLOT1 4
-#define RESERVE_ATMLOT2 5
+#define PTR_FRAC_CHARGE 4
+#define RESERVE_ATMSLOT 5
 #define ATM_SLOTS       6
 
 
@@ -150,6 +150,78 @@ typedef struct {
     double **log_max_coeff;
     PairData **pairdata;  // NULL indicates not-initialized, NO_VALUE can be skipped
 } CINTOpt;
+
+// Add this macro def to make pyscf compatible with both v4 and v5
+#define HAVE_DEFINED_CINTENVVARS_H
+typedef struct {
+        FINT *atm;
+        FINT *bas;
+        double *env;
+        FINT *shls;
+        FINT natm;
+        FINT nbas;
+
+        FINT i_l;
+        FINT j_l;
+        FINT k_l;
+        FINT l_l;
+        FINT nfi;  // number of cartesian components
+        FINT nfj;
+        // in int1e_grids, the grids_offset and the number of grids
+        union {FINT nfk; FINT grids_offset;};
+        union {FINT nfl; FINT ngrids;};
+        FINT nf;  // = nfi*nfj*nfk*nfl;
+        FINT rys_order; // = nrys_roots for regular ERIs. can be nrys_roots/2 for SR ERIs
+        FINT x_ctr[4];
+
+        FINT gbits;
+        FINT ncomp_e1; // = 1 if spin free, = 4 when spin included, it
+        FINT ncomp_e2; // corresponds to POSX,POSY,POSZ,POS1, see cint.h
+        FINT ncomp_tensor; // e.g. = 3 for gradients
+
+        /* values may diff based on the g0_2d4d algorithm */
+        FINT li_ceil; // power of x, == i_l if nabla is involved, otherwise == i_l
+        FINT lj_ceil;
+        FINT lk_ceil;
+        FINT ll_ceil;
+        FINT g_stride_i; // nrys_roots * shift of (i++,k,l,j)
+        FINT g_stride_k; // nrys_roots * shift of (i,k++,l,j)
+        FINT g_stride_l; // nrys_roots * shift of (i,k,l++,j)
+        FINT g_stride_j; // nrys_roots * shift of (i,k,l,j++)
+        FINT nrys_roots;
+        FINT g_size;  // ref to cint2e.c g = malloc(sizeof(double)*g_size)
+
+        FINT g2d_ijmax;
+        FINT g2d_klmax;
+        double common_factor;
+        double expcutoff;
+        double rirj[3]; // diff by sign in different g0_2d4d algorithm
+        double rkrl[3];
+        double *rx_in_rijrx;
+        double *rx_in_rklrx;
+
+        double *ri;
+        double *rj;
+        double *rk;
+        // in int2e or int3c2e, the coordinates of the fourth shell
+        // in int1e_grids, the pointer for the grids coordinates
+        union {double *rl; double *grids;};
+
+        FINT (*f_g0_2e)();
+        void (*f_g0_2d4d)();
+        void (*f_gout)();
+        CINTOpt *opt;
+
+        /* values are assigned during calculation */
+        int *idx;
+        double ai[1];
+        double aj[1];
+        double ak[1];
+        double al[1];
+        double fac[1];
+        double rij[3];
+        double rkl[3];
+} CINTEnvVars;
 #endif
 
 FINT CINTlen_cart(const FINT l);

include/cintwrapper.h

@@ -141,9 +141,25 @@ void gen_4c_overlap(double* ovlp4, size_t N,
                     int natm, int nbas, int nenv,
                     int* atm, int* bas, double* env);
 
-void gen_eri(double **eri, int N,
+void gen_eri_4b(double** eri, int N1, int N2, int N3, int N,
+             int natm, int nbas1, int nbas2, int nbas3, int nbas, int nenv,
+             int* atm, int* bas, double* env, int nomp);
+
+void gen_eri_2b(double** eri, int N1, int N,
+             int natm, int nbas1, int nbas, int nenv,
+             int* atm, int* bas, double* env, int nomp);
+
+void gen_eri(double** eri, int N,
+             int natm, int nbas, int nenv,
+             int* atm, int* bas, double* env);
+
+void gen_eri(double** eri, int N,
+             int natm, int nbas, int nenv,
+             int* atm, int* bas, double* env, int nomp);
+
+void gen_eri_2(double** eri, int N,
              int natm, int nbas, int nenv,
-             int *atm, int *bas, double *env);
+             int* atm, int* bas, double* env, int nomp);
 
 void gen_jMOI_gto(double **eri, int N,
                   int natm, int nbas, int nenv,

include/cpu_util.h

@@ -23,11 +23,15 @@ using namespace std;
 double randomf(double a, double b);
 
 void solve_axeb(int dim, double* A, double* b);
+void solve_axeb_ls(int m, int n, double* A, double* b);
+
+double determinant(int N, double* A);
 
 void expmat_complex_cpu(int N, double* theta, double* thetai, double* etheta);
 void expmat_cpu(int N, double *theta, double *etheta);
 
 int la_diagR(int neig, double* A, double* eigen, double* eigeni);
+void la_diag_gen(int neig, int s1, double* A, double* S, double* Ae);
 void la_diag(int neig, int s1, double* A, double* Ae);
 int invert_stable_cpu(double* A, int size, double delta);
 int invert_stable_cpu(double* A, int size, double delta, int prl);

include/cuda_util.h

@@ -22,6 +22,7 @@ void diagonalize_cusolver(int Ne, int N, double* A, double* Ae, cusolverDnHandle
 
 int mat_root_inv_cusolver(double* A, int size, cusolverDnHandle_t& cu_hdl);
 int mat_root_inv_stable_cusolver(double* A, int size, double delta, cusolverDnHandle_t& cu_hdl);
+int mat_root_inv_stable_cusolver(double* A, int size, double delta, cusolverDnHandle_t& cu_hdl, int prl);
 
 void solve_axeb_gpu(int dim, double* A, double* b, cusolverDnHandle_t& cusolverH);
 
@@ -51,4 +52,4 @@ void copy_to_all_gpu(int ngpu, int s1, float* A, int include_first);
 void copy_to_all_gpu(int ngpu, int s1, int s2, double** A, int include_first);
 
 
-#endif // CUDA_UTILH
+#endif // CUDA_UTILH

include/gauss.h

@@ -21,7 +21,7 @@ void eval_p_gh(int gs, float* grid, float* val, int n1, int l1, int m1, float no
 void eval_pd_gh(int gs, double* grid, double* val, int n1, int l1, int m1, double norm1, double zeta1);
 
 //Gaussian single-center potentials
-void eval_vghd(int gsa, double* grid, double* val, int n, int l, int m, double norm, double zt);
+void eval_vghd(int gsa, double* grid, double* val, int n, int l, double norm, double zt);
 
 void eval_hess_ghd(int gs, double* grid, double* val, int n1, int l1, int m1, double norm1, double zeta1);
 int eval_gh_full(int gs, float* grid, float** val1, int i1, int natoms, int nbas, int nenv, int N, int* atm, int* bas, double* env);
@@ -35,5 +35,7 @@ void wf_to_grid_gh(bool divide_by_rho, bool calc_rho, int natoms, int* atno, dou
 //void integrate_hole_para_gh(double* rdm, bool full_rdm, bool hfx_on, int Nc, int No, int M, int natoms, int* atno, double* coords, int gs, int gsb, vector<vector<double> >& basis,
 //                    double* Pao, double* Pmo, double* jCA, float* grid, float* gridb, float* wt, float* wtb, float* rho, float* vxch, int prl);
 
+void gen_4cj_overlap(double* ol4j, int nrad, int nang, double* ang_g, double* ang_w,
+                     int natoms, int* atno, double* coords, vector<vector<double> > basis, int prl);
 
 #endif

include/isotopes.h

@@ -0,0 +1,10 @@
+#ifndef ISOTOPESH
+#define ISOTOPESH
+
+//most abundant isotopes of elements H through Ne
+const int isos1[11] = { 0, 1, 4, 7, 9, 11, 12, 14, 16, 19, 20 };
+const int isos2[11] = { 0, 2, 3, 6, 0, 10, 13, 15, 18, 0, 22 };
+const int isos3[11] = { 0, 3, 0, 0, 0, 0, 0, 0, 17, 0, 21 };
+
+
+#endif

include/read.h

@@ -26,13 +26,16 @@ string SSTRF(float x);
 string SSTRF2(double x);
 
 vector<string> split1(const string &s, char delim);
+vector<string> split2(const string &s, char delim);
 
 int check_file(string filename);
 
 int initialize(bool gbasis, vector<vector<double> >& basis, vector<vector<double> >& basis_aux, int* atno, double* &coords, int& charge, int& unpaired, double& Enn, int prl);
 bool get_secondary_basis(string name, int natoms, int* atno, double* coords, vector<vector<double> >& basis, vector<vector<double> >& basis_aux, int prl);
 void create_basis_aux(int natoms, vector<vector<double> >& basis_std, vector<vector<double> >& basis_aux);
 string read_basis_text(string aname);
+string read_basis_text(string filename, string aname);
+int read_geoms(int& ngeom1, int natoms, vector<double*>& coords, int prl);
 void print_basis(int natoms, vector<vector<double> >& basis, vector<vector<double> >& basis_aux, int prl);
 int read_rotate(int N, double* jCA);
 double read_float(string filename);
@@ -93,6 +96,10 @@ bool read_MOI_from_file(int Mm, int M3, double** MOI, string filename, int prl);
 
 bool check_PS();
 
+double nuclear_repulsion(int natoms, double* atno, double* coords);
+double nuclear_repulsion(int natoms, int* atno, double* coords);
+
+
 string get_aname(int Z);
 
 #endif

include/write.h

@@ -51,6 +51,7 @@ void write_square_clean(int N, double* A, string fname, double thresh, int prl);
 
 void write_molden(bool gbasis, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname);
 void write_molden(bool gbasis, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, string fname);
+void write_molden_g(int btype, int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname);
 void write_molden_g(int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, double* eig, string fname);
 //void write_molden_vcf(int natoms, int* atno, double* coords, vector<vector<vcf> > &vcfs, string fname);
 void write_molden_ss(int natoms, int* atno, double* coords, vector<vector<double> > &basis, double* jCA, int No, string fname);

src/integrals/becke.cpp

@@ -954,16 +954,27 @@ void becke_weight_2c(int gs, float* grid1, float* wt1, float* grid2, float* wt2,
 }
 
 void becke_weight_2d(int gs, double* grid1, double* wt1, double* grid2, double* wt2,
-                     double zeta1, double zeta2, double A2, double B2, double C2)
+                     bool is_Z, double zeta1, double zeta2, double A2, double B2, double C2)
 {
   int mu_order = 3;
 
  //first center is at 0,0,0 and second at A2,B2,C2
   double R  = sqrt(A2*A2+B2*B2+C2*C2);
   const double oR = 1./R;
 
-  const double a1 = becke_a_zeta(zeta1,zeta2);
-  const double a2 = becke_a_zeta(zeta2,zeta1);
+  double a1 = 0.;
+  double a2 = 0.;
+  if (is_Z)
+  {
+   //not actually zeta input, but Z input to this ftn
+    a1 = becke_a(zeta1,zeta2);
+    a2 = becke_a(zeta2,zeta1);
+  }
+  else
+  {
+    a1 = becke_a_zeta(zeta1,zeta2);
+    a2 = becke_a_zeta(zeta2,zeta1);
+  }
 
   //printf(" a1/2: %8.5f %8.5f \n",a1,a2);
 
@@ -5220,9 +5231,6 @@ void compute_delt(int natoms, int* atno, double* coords, bool gbasis, vector<vec
 
   if (sgs_basis) { printf("\n ERROR: compute_rho for SGS basis does not support drho \n"); exit(-1); }
 
-  if (gbasis)
-    printf("  TESTING: compute_delt for gbasis \n");
-
   int N = basis.size();
   int N2 = N*N;
   int* n2i = new int[natoms];