comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-28 06:05:44 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

Merge pull request #1242 from edoapra/4basis

Browse source 
bug fixes in integral codes
This commit is contained in:
NWChem: Open Source High-Performance Computational Chemistry 2026-02-02 23:56:25 -08:00 • committed by GitHub
commit f75be727ff
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
7 changed files with 132 additions and 61 deletions
Download patch file Download diff file Expand all files Collapse all files

8
src/NWints/api/exactd_mem.F
View file

@ -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)),

38
src/NWints/deriv/hf1d.F
View file

@ -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.

7
src/NWints/dk/GNUmakefile
View file

@ -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

23
src/NWints/dk/dk_utils.F Normal file
View file

@ -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

60
src/NWints/dk/int_1edk_ga.F
View file

@ -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

5
src/util/GNUmakefile
View file

@ -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

52
src/util/ga_matpow.F
View file

@ -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