From 1562986ddae01711ef484ca5c19d0a087b351a63 Mon Sep 17 00:00:00 2001 From: edoapra Date: Mon, 2 Feb 2026 15:30:44 -0800 Subject: [PATCH 1/4] fix for issue https://github.com/nwchemgit/nwchem/issues/1241 --- src/NWints/api/exactd_mem.F | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/NWints/api/exactd_mem.F b/src/NWints/api/exactd_mem.F index bc2aedc5de..54525685b4 100644 --- a/src/NWints/api/exactd_mem.F +++ b/src/NWints/api/exactd_mem.F @@ -172,7 +172,7 @@ c::local integer bra_prim, ket_prim, ket_lim integer bra_gen, ket_gen integer max_nbf_bra,max_nbf_ket, nat - integer bra_fact, ket_fact + integer bra_fact, ket_fact, lp double precision factor c #include "bas_exndcf_sfn.fh" @@ -238,7 +238,11 @@ c ket_gen = 1 endif c - memdum = 60 000 000 +c must be larger than NPP*Lp3*ncenters (line 340 of hf1d.F) + lp=bra_type+ket_type+1 + memdum= (((lp+1)*(lp+2)*(lp+3))/6* + * nat*bra_prim*ket_prim)*125/100 + memdum=max(memdum,60000000) nint = int_unint(brain,ibc,ketin,ikc,0,0,0,0) call hf1d( & dbl_mb(mb_exndcf(1,1)), dbl_mb(mb_exndcf(1,1)), From cefe6e2878f75f269437949248a7c911672fd101 Mon Sep 17 00:00:00 2001 From: edoapra Date: Mon, 2 Feb 2026 15:31:04 -0800 Subject: [PATCH 2/4] implicit none --- src/NWints/deriv/hf1d.F | 38 ++++++++++++++++++++++---------------- 1 file changed, 22 insertions(+), 16 deletions(-) diff --git a/src/NWints/deriv/hf1d.F b/src/NWints/deriv/hf1d.F index fcc20818e3..5d37f1f1ff 100644 --- a/src/NWints/deriv/hf1d.F +++ b/src/NWints/deriv/hf1d.F @@ -16,29 +16,30 @@ C> Lachlan Belcher 14 Feb 2023 & O2I,KEI,NAI,canAB,dryrun,W0,maxW0) c $Id$ - Implicit real*8 (a-h,o-z) - Implicit integer (i-n) + implicit none logical hlfdrv ! [input] if true, calculate same atom as well ! as different atom integrals. This is necessary ! when calculating the half-derivative overlap ! matrix for NADCTs - integer ictrA, ictrB !NB already used implicitly + integer ictrA, ictrB !NB already used implicitly + integer NPA,NCA,La,NPB,NCB,Lb,ncenters,nint,maxw0 Logical dryrun Logical O2I,KEI,NAI,canAB,GenCon c--> Cartesian Coordinates, Primitives & Contraction Coefficients - Dimension Axyz(3),Aprims(NPA),Acoefs(NPA,NCA) - Dimension Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB) + double precision Axyz(3),Aprims(NPA),Acoefs(NPA,NCA) + double precision Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB) c--> Nuclear Cartesian Coordinates, Charges & Inverse Exponents - Dimension Cxyz(3,ncenters),zan(ncenters),exinv(ncenters) + double precision Cxyz(3,ncenters),zan(ncenters),exinv(ncenters) c--> Blocks of Overlap, Kinetic Energy & Nuclear Attraction Integrals - Dimension bO2I(Nint,*),bKEI(Nint,*),bNAI(Nint,*) + Double precision bO2I(Nint,*),bKEI(Nint,*),bNAI(Nint,*) + Double Precision W0(maxW0) hlfdrv=.false. @@ -66,9 +67,7 @@ C------------------------------------------------------------------------- & hlfdrv) c $Id$ - Implicit real*8 (a-h,o-z) - Implicit integer (i-n) - + implicit none #include "errquit.fh" #include "stdio.fh" LOGICAL smatm ! if true, both integrals use the same atom @@ -80,30 +79,31 @@ c $Id$ ! when calculating the half-derivative overlap ! matrix for NADCTs integer ictrA, ictrB !NB already used implicitly + integer NPA,NCA,La,NPB,NCB,Lb,ncenters,nint,maxw0 Logical dryrun Logical O2I,KEI,NAI,canAB,GenCon c--> Cartesian Coordinates, Primitives & Contraction Coefficients - Dimension Axyz(3),Aprims(NPA),Acoefs(NPA,NCA) - Dimension Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB) + Double Precision Axyz(3),Aprims(NPA),Acoefs(NPA,NCA) + Double Precision Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB) c--> Nuclear Cartesian Coordinates, Charges & Inverse Exponents - Dimension Cxyz(3,ncenters),zan(ncenters),exinv(ncenters) + Double Precision Cxyz(3,ncenters),zan(ncenters),exinv(ncenters) c--> Blocks of Overlap, Kinetic Energy & Nuclear Attraction Integrals - Dimension bO2I(Nint,*),bKEI(Nint,*),bNAI(Nint,*) + Double Precision bO2I(Nint,*),bKEI(Nint,*),bNAI(Nint,*) c--> Derivative Indices - Dimension inder1(3,3),inder2(6,6) + integer inder1(3,3),inder2(6,6) c--> Scratch Space. - Dimension W0(maxW0) + Double Precision W0(maxW0) Data inder1/ 1,0,0, & 0,1,0, @@ -115,6 +115,12 @@ c--> Scratch Space. & 0,0,0,1,0,0, & 0,0,0,0,1,0, & 0,0,0,0,0,1 / +c + integer i_ALPHAp,i_Ep,i_ESp,i_ff,i_IJK,i_IPAIRp, + I i_left,i_P,i_PC,i_pf,i_R0,i_R0C,i_right, + I i_Rj,i_RS,i_Ti,i_top,i_VA,i_VB,i_VP,i_VR, + I Li, Lp,Lp3, lprod, MaxMem, MXD, + I n, NCP, nd, nder, npp c c Compute gradient of the overlap, kinetic energy, and nuclear attraction c integrals for 2 shells of contracted Gaussians functions. From b4b4cd4db1c210fc8e0ceda83f0034b89f717e59 Mon Sep 17 00:00:00 2001 From: edoapra Date: Mon, 2 Feb 2026 17:46:28 -0800 Subject: [PATCH 3/4] added serial version of ga_matpow and use ga_pdsyevr if available --- src/util/GNUmakefile | 5 +++++ src/util/ga_matpow.F | 52 +++++++++++++------------------------------- 2 files changed, 20 insertions(+), 37 deletions(-) diff --git a/src/util/GNUmakefile b/src/util/GNUmakefile index e5ea3fed11..76b6e09217 100644 --- a/src/util/GNUmakefile +++ b/src/util/GNUmakefile @@ -346,6 +346,11 @@ ifdef USE_NO_WALLTIME LIB_DEFINES += -DNO_WALLTIME endif +_PDSYEVR := $(shell [ -e ${NWCHEM_TOP}/src/tools/install/lib/libga.a ] && nm ${NWCHEM_TOP}/src/tools/install/lib/libga.a|egrep T |egrep ga_pdsyevr |awk '/ga_pdsyevr/ {print "Y"}') +ifeq ($(_PDSYEVR),Y) + LIB_DEFINES += -DUSE_PDSYEVR +endif + include ../config/makelib.h diff --git a/src/util/ga_matpow.F b/src/util/ga_matpow.F index cd54d4bc6a..0c9fbea104 100644 --- a/src/util/ga_matpow.F +++ b/src/util/ga_matpow.F @@ -6,7 +6,12 @@ C REVISION C $Id$ C C SYNOPSIS +#if !defined SERIAL_DIAG Subroutine GA_MatPow(G_V, Pow, MinEval) +#elif defined(SERIAL_DIAG) +#warning serial sub + Subroutine GA_MatPow_serial(G_V, Pow, MinEval) +#endif Implicit NONE #include "errquit.fh" Integer G_V @@ -57,9 +62,6 @@ C LOCAL VARIABLES Integer D, DiagLo, DiagHi Integer G_EV, G_Tmp Integer I, J, VType, Rows, Cols, ILo, IHi, JLo, JHi -#if defined(USE_GEN_DIAG) - Integer G_Junk -#endif Double precision Elem C C *************************** @@ -109,44 +111,13 @@ C C C Now we can solve the stupid thing... C -#if defined(USE_GEN_DIAG) -C Create an array for the metric and turn it into the unit matrix -C - If ( .NOT. GA_Create( VType, Rows, Cols, 'Metric', - $ -MinChunk, -MinChunk, G_Junk) ) - $ Call ErrQuit('GA_MatPow: Unable to allocate metric', - $ Rows*Cols, GA_ERR) -C -C This is clearly serial, but hopefully soon we should be able to -C go back to the standard eigensolver instead of the generalized one. -C - Call GA_Zero( G_Junk) - If ( GA_NodeID() .eq. 0) then - Do I = 1, Rows - Call GA_Put(G_Junk, I, I, I, I, FP_One, 1) - EndDo - EndIf -C -#ifdef SCALAPACK - call ga_pdsygv(g_v, g_junk, g_ev, dbl_mb(eval_idx),0) -#elif defined(PEIGS) - Call GA_Diag( G_V, G_Junk, G_EV, Dbl_MB(Eval_Idx)) -#else - Call GA_Diag_Seq( G_V, G_Junk, G_EV, Dbl_MB(Eval_Idx)) -#endif -C - If ( .NOT. GA_Destroy( G_Junk) ) - $ Call ErrQuit('GA_MatPow: Unable to free metric', 0, GA_ERR) -C -#else /* defined(USE_GEN_DIAG) */ -#ifdef SCALAPACK - call ga_pdsyev(g_v, g_ev, dbl_mb(eval_idx),0) -#elif defined(PEIGS) +#if defined(SCALAPACK) && defined(USE_PDSYEVR) && !defined(SERIAL_DIAG) + call ga_pdsyevr(g_v, g_ev, dbl_mb(eval_idx),0,rows) +#elif defined(PEIGS) && !defined(SERIAL_DIAG) Call GA_Diag_Std( G_V, G_EV, Dbl_MB(Eval_Idx)) #else Call GA_Diag_Std_Seq( G_V, G_EV, Dbl_MB(Eval_Idx)) #endif -#endif /* defined(USE_GEN_DIAG) */ C C Sometimes people may just want a look at the eigenvalues. C @@ -301,3 +272,10 @@ C C Return End +#ifndef SERIAL_DIAG +#define SERIAL_DIAG +c +c Compile source again for the serial diag case +c +#include "ga_matpow.F" +#endif From cbc8c25a18647bf5161426051055f6f447aee759 Mon Sep 17 00:00:00 2001 From: edoapra Date: Mon, 2 Feb 2026 17:48:25 -0800 Subject: [PATCH 4/4] use np<=3 to fix numerical issues old code still available by uncommen ting the PGROUP define --- src/NWints/dk/GNUmakefile | 7 ++++- src/NWints/dk/dk_utils.F | 23 ++++++++++++++ src/NWints/dk/int_1edk_ga.F | 60 +++++++++++++++++++++++++++++++++---- 3 files changed, 84 insertions(+), 6 deletions(-) create mode 100644 src/NWints/dk/dk_utils.F diff --git a/src/NWints/dk/GNUmakefile b/src/NWints/dk/GNUmakefile index 9fc69209e0..8255380961 100644 --- a/src/NWints/dk/GNUmakefile +++ b/src/NWints/dk/GNUmakefile @@ -4,7 +4,7 @@ include ../../config/makefile.h OBJ_OPTIMIZE = int_1edk_ga.o \ dk_input.o dkdiv_energy.o int_pxvp_ga.o \ - int_1edkso_ga.o + int_1edkso_ga.o dk_utils.o ifdef XLF10 OBJ = bas_create_uncont.o @@ -16,5 +16,10 @@ endif USES_BLAS = int_1edk_ga.F int_1edkso_ga.F +_PDSYEVR := $(shell [ -e ${NWCHEM_TOP}/src/tools/install/lib/libga.a ] && nm ${NWCHEM_TOP}/src/tools/install/lib/libga.a|egrep T |egrep ga_pdsyevr |awk '/ga_pdsyevr/ {print "Y"}') +ifeq ($(_PDSYEVR),Y) + LIB_DEFINES += -DUSE_PDSYEVR +endif + include ../../config/makelib.h diff --git a/src/NWints/dk/dk_utils.F b/src/NWints/dk/dk_utils.F new file mode 100644 index 0000000000..62e3e63489 --- /dev/null +++ b/src/NWints/dk/dk_utils.F @@ -0,0 +1,23 @@ + subroutine ga_serialcopy(g_inpg,g_inworld) + implicit none +#include "mafdecls.fh" +#include "errquit.fh" +#include "global.fh" + integer g_inpg,g_inworld +c + integer itype,dim + integer k_a,l_a,i,j +c + if(ga_nodeid().eq.0) then + call ga_inquire(g_inpg, itype, dim, dim) + if (.not. ma_push_get(MT_DBL,dim*dim,'MA temp',l_a,k_a)) + $ call errquit('serialc: ma_push_get failed MA temp', + Z 0 , MA_ERR) + call ga_get(g_inpg,1,dim,1,dim,dbl_mb(k_a),dim) + call ga_put(g_inworld,1,dim,1,dim,dbl_mb(k_a),dim) + if (.not. ma_pop_stack(l_a)) + $ call errquit('serialc: ma_pop_stack l_a',0, MA_ERR) + endif + call ga_sync() + return + end diff --git a/src/NWints/dk/int_1edk_ga.F b/src/NWints/dk/int_1edk_ga.F index f6b126f214..05bc877ab1 100644 --- a/src/NWints/dk/int_1edk_ga.F +++ b/src/NWints/dk/int_1edk_ga.F @@ -39,6 +39,8 @@ c:tex-\noindent Author: W. A. de Jong c:tex- c:tex-{\it Syntax:} c:tex-\begin{verbatim} +#define PGROUP +#define PGROUP_MAXPROC 3 subroutine int_1edk_ga(basis_in, g_in, type, kskel) c:tex-\end{verbatim} implicit none @@ -122,9 +124,13 @@ c integer ZpxVp_ga ! Unmodified pxVp for DK3 integer tmp3g_a ! Temporary g_a matrix #3 handle for ! full DK3 + integer max2,list2(PGROUP_MAXPROC),p_handle,g_inpg,g_inworld + integer comm_world + logical status integer info integer ga_cholesky,ga_llt_i external ga_cholesky,ga_llt_i + integer itype,dim c c cau = 137.0360000d0 ! Speed of light by Bernd Hess c @@ -135,6 +141,18 @@ c c c Define the dimensions of the arrays c + call ga_sync() +#ifdef PGROUP + max2=min(ga_nnodes(),PGROUP_MAXPROC) + do iloop=1,PGROUP_MAXPROC + list2(iloop)=iloop-1 + enddo + p_handle=ga_pgroup_create(list2,max2) + call ga_inquire(g_in, itype, dim, dim) + comm_world=ga_pgroup_get_world() + if(ga_nodeid().lt.max2) then + call ga_pgroup_set_default(p_handle) +#endif if (.not. bas_numbf(dk_fitbas,nbf)) $ call errquit('int_1edk_ga: bas_numbf failed',0, BASIS_ERR) c @@ -182,7 +200,7 @@ c if (.not. ma_push_get(MT_DBL,nbf,'psquared',p2m_a,p2Im_a)) $ call errquit('int_1edk_ga: ma_push_get failed psquared',0, & MA_ERR) -#ifdef SCALAPACK +#if defined(SCALAPACK) && !defined(PGROUP) c pdsygv is unstable since it uses the inaccurate pdsyevx cold call ga_pdsygv(Tg_a,tempg_a,Ug_a,dbl_mb(p2Im_a)) c S=chol*chol_t @@ -190,7 +208,7 @@ c S=chol*chol_t c inv(a) info = ga_llt_i('L',tempg_a,-1) c U=chol(inv(a))_t - info = ga_cholesky('L',tempg_a,-1) + info = ga_cholesky('L',tempg_a) c f_prime=U_t*f*inv(chol)*U call ga_dgemm('t','n',nbf,nbf,nbf,1.0d0,tempg_a,Tg_a,0.0d0,tmpg_a) @@ -201,7 +219,7 @@ c diag f_prime c_prime c c=inv(chol)_t*c_prime call ga_dgemm('n','n',nbf,nbf,nbf,1.0d0,tempg_a,tmpg_a, Z 0.0d0,ug_a) -#elif defined(PEIGS) +#elif defined(PEIGS) && !defined(PGROUP) call ga_diag(Tg_a,tempg_a,Ug_a,dbl_mb(p2Im_a)) #else call ga_diag_seq(Tg_a,tempg_a,Ug_a,dbl_mb(p2Im_a)) @@ -1744,6 +1762,7 @@ c c S_-1 is inverse overlap in the uncontracted basis c Scontr_uncon is the mixed basis overlap c + if (.not. bas_numbf(basis_in,nbf_in)) $ call errquit('int_1edk_ga: bas_numbf failed',0, BASIS_ERR) if (.not. ga_create(MT_DBL,nbf,nbf_in,'Mixed overlap',1,1, @@ -1765,7 +1784,11 @@ c Create overlap matrix dk_fitbas, invert to S^-1 and multiply with Mixg_a c call ga_zero(tempg_a) call int_1e_oldga(dk_fitbas,dk_fitbas,tempg_a,'overlap',.false.) +#ifdef PGROUP + call ga_matpow_serial(tempg_a, -1.0d0, 0.0d0) +#else call ga_matpow(tempg_a, -1.0d0, 0.0d0) +#endif call ga_dgemm('n','n',nbf,nbf_in,nbf,1.0d0,tempg_a,Mixg_a, $ 0.0d0,tmpg_a) call ga_copy(tmpg_a,Mixg_a) @@ -1782,9 +1805,18 @@ c call ga_zero(tempg_a) call ga_dgemm('t','n',nbf_in,nbf,nbf,1.0d0,Mixg_a, $ DKg_a,0.0d0,tempg_a) + call ga_zero(DKg_a) +#ifdef PGROUP + if (.not. ga_create(MT_DBL,dim,dim,'GA tmp',1,1,g_inpg)) + $ call errquit('int_1edk_ga: ga_create failed GA temp',0, + & GA_ERR) + call ga_zero(g_inpg) +#else + g_inpg=g_in +#endif call ga_dgemm('n','n',nbf_in,nbf_in,nbf,1.0d0,tempg_a,Mixg_a, - $ 1.0d0,g_in) - call ga_symmetrize(g_in) + $ 1.0d0,g_inpg) + call ga_symmetrize(g_inpg) c c Destroy transformation matrices Mixg_a, tempg_a and DKg_a c @@ -1795,6 +1827,24 @@ c if (.not. ga_destroy(DKg_a)) $ call errquit('int_1edk_ga: ga_destroy DKg_a',0, GA_ERR) c +#ifdef PGROUP + endif + call ga_pgroup_set_default(comm_world) + call ga_sync() + if (.not. ga_create(MT_DBL,dim,dim,'GA tmp',1,1,g_inworld)) + $ call errquit('int_1edk_ga: ga_create failed GA temp',0, + & GA_ERR) + call ga_serialcopy(g_inpg,g_inworld) + if(ga_nodeid().lt.max2) then + if (.not. ga_destroy(g_inpg)) + $ call errquit('int_1edk_ga: ga_destroy g_inpg',0, GA_ERR) + endif + call ga_add(1d0,g_inworld,1d0,g_in, g_in) + if (.not. ga_destroy(g_inworld)) + $ call errquit('int_1edk_ga: ga_destroy g_inworld',0, GA_ERR) + if(.not.ga_pgroup_destroy(p_handle)) + $ call errquit('int_1edk_ga: pgroupdestroy',0, GA_ERR) +#endif return c end