diff --git a/src/bsemol/GNUmakefile b/src/bsemol/GNUmakefile new file mode 100644 index 0000000000..116ecb3dd4 --- /dev/null +++ b/src/bsemol/GNUmakefile @@ -0,0 +1,32 @@ +OBJ = bsemol.o bse_input.o + +USES_BLAS = +OBJ_OPTIMIZE = + +LIBRARY = libbsemol.a + +HEADERS = ../nwdft/include/cdft.fh ../nwdft/grid/cgridfile.fh + +LIB_DEFINES = +LIB_INCLUDES = -I./ + +include ../config/makefile.h + + +ifeq ($(_USE_SCALAPACK),Y) +_USE_ELPA = $(shell ${GA_PATH}/bin/ga-config --use_elpa| awk ' /1/ {print "Y"}') + +ifeq ($(_USE_ELPA),Y) + LIB_DEFINES += -DUSE_ELPA +endif +#assuming we have GAPZHEEVD +GOTGAPZHEEVD= Y +ifeq ($(GOTGAPZHEEVD),Y) + LIB_DEFINES += -DUSE_GAPZHEEVD +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 +endif +include ../config/makelib.h diff --git a/src/bsemol/bse.fh b/src/bsemol/bse.fh new file mode 100644 index 0000000000..d89e5f3ba6 --- /dev/null +++ b/src/bsemol/bse.fh @@ -0,0 +1,35 @@ + + double precision, parameter :: ha2ev = 27.211386245988d0 + + type :: bse_params_t + sequence + integer :: rtdb + integer :: me, nprocs, iMaxThreads + integer :: ipol + logical :: main, odft +c + integer :: natoms + integer :: nmo + integer :: nocc(2),nvir(2),npoles(2) + integer :: noqp(2),nvqp(2),nstates(2),mynpoles(2) + integer :: oolo(2),oohi(2) + integer :: ovlo(2),ovhi(2) + integer :: vvlo(2),vvhi(2) + integer :: llmo(2),nqp(2) + integer :: ngl + double precision :: tol2e, eta +c + integer :: l_evals, k_evals + integer :: l_occ, k_occ +c + integer :: fd +c + integer :: g_movecs(2), g_moeris, g_sols + integer :: g_erioo(2), g_eriov(2), g_erivv(2), g_erim + logical :: noio +c + logical :: ncap, davidson, lanczos + double precision :: thresh, fdstep, vxddp(2), vxddm(2) + + end type bse_params_t + diff --git a/src/bsemol/bse_input.F b/src/bsemol/bse_input.F new file mode 100644 index 0000000000..34429f731a --- /dev/null +++ b/src/bsemol/bse_input.F @@ -0,0 +1,76 @@ + subroutine bse_input(rtdb) +c +c Input parser for molecular BSE module +c +c Daniel Mejia-Rodriguez +c November 2022 +c + implicit none +c +#include "errquit.fh" +#include "util.fh" +#include "inp.fh" +#include "rtdb.fh" +#include "mafdecls.fh" +#include "stdio.fh" +c + integer rtdb + character*255 test + character(*), parameter :: pname = "bse_input: " +c + integer, parameter :: num_dirs = 2 + character*11 dirs(num_dirs) + data dirs / 'method', 'end' / +c + integer ind, qpiter,eviter,noqp(2),nvqp(2), ngl, first + double precision threshold, eta + double precision, parameter :: ha2ev = 27.211386245988d0 +c +c Set default values and put them in RTDB +c + call bse_defaults(rtdb) +c +c Parse the input file +c + 10 if (.not. inp_read()) then + call errquit(pname//'failed to read input',0,INPUT_ERR) + endif +c + if (.not. inp_a(test)) then + call errquit(pname//'failed to read keyword',0,INPUT_ERR) + endif +c + if (.not.inp_match(num_dirs, .false., test, dirs, ind)) then + call errquit(pname//'unknown directive', 0, INPUT_ERR) + end if +c + goto (100, 9999) ind + + call errquit(pname//'unknown directive',ind, INPUT_ERR) +c +c ----------------------------------------------------------------- +c +c method +c + 100 if (.not. inp_a(test)) + & call errquit(pname//'could not read method', 100, INPUT_ERR) + + if (inp_compare(.false.,'analytic',test)) then + if (.not.rtdb_put(rtdb, 'bse:analytic', mt_log, 1, .true.)) + & call errquit(pname//'rtdb_put failed', 101, RTDB_ERR) + elseif (inp_compare(.false.,'bse:davidson',test)) then + if (.not.rtdb_put(rtdb, 'bse:davidson', mt_log, 1, .true.)) + & call errquit(pname//'rtdb_put failed', 102, RTDB_ERR) + elseif (inp_compare(.false.,'bse:lanczos',test)) then + if (.not.rtdb_put(rtdb, 'bse:davidson', mt_log, 1, .true.)) + & call errquit(pname//'rtdb_put failed', 103, RTDB_ERR) + else + call errquit(pname//'unknown BSE method',104,INPUT_ERR) + endif + goto 10 +c +c normal termination +c + 9999 return +c + end subroutine bse_input diff --git a/src/bsemol/bsemol.F b/src/bsemol/bsemol.F new file mode 100644 index 0000000000..2ed16f1ff7 --- /dev/null +++ b/src/bsemol/bsemol.F @@ -0,0 +1,73 @@ + logical function bsemol(rtdb) + implicit none +#include "rtdb.fh" + integer rtdb +c + logical bse_doit + external bse doit +c + logical status +c + status = rtdb_cput(rtdb,'dft:theory', 1, 'dft') + if (status) status = bse_doit(rtdb) + bsemol = status + call grid_cleanup(.false.) +c + end function +c +c ----------------------------------------------------------------- +c + logical function bse_doit(rtdb) + implicit none +#include "errquit.fh" +#include "global.fh" +#include "rtdb.fh" +#include "stdio.fh" +#include "util.fh" +#include "bse.fh" + + integer, intent(in) :: rtdb + character(*), parameter :: pname = "bse_doit: " + + type(bse_params_t) :: params + character*80 theory + integer me +c +c ----------------------------------------------------------------- +c +c Initialize +c + bse_doit = .false. +c +c Determine the theory we are using +c + if (.not. rtdb_cget(rtdb, "dft:theory", 1, theory)) + $ call errquit(pname//"theory not specified",0,RTDB_ERR) +c +c Initialize GW calculation +c + params%rtdb = rtdb + params%me = ga_nodeid() +c + call bse_init(params) +c + if (theory.eq."dft") then + + if (params%davidson) then + call bse_davidson(params) + else + call bse_analytic(params) + endif +c + elseif (theory.eq.'sodft') then + call errquit(pname//"spin-orbit not implemented",0,0) + else + call errquit(pname//"invalid theory: "//trim(theory),0,0) + endif +c + call bse_finalize(params) + bse_doit = .true. +c +c ----------------------------------------------------------------- +c + end function diff --git a/src/config/make_nwchem_config b/src/config/make_nwchem_config index 81ac74ce91..03427ae572 100644 --- a/src/config/make_nwchem_config +++ b/src/config/make_nwchem_config @@ -21,7 +21,7 @@ KNOWN_MODULE_SUBDIRS = NWints atomscf ddscf develop gradients moints nwdft \ property solvation nwpw fft rimp2_grad python argos analyz diana nwmd cafe space \ drdy vscf uccsdt qmmm rism qmd etrans tce geninterface bq mm cons\ perfm dntmc smd nbo dangchang leps ccca lucia dimqm rdmft fcidump\ - gwmol xtb nwc_columbus + gwmol xtb nwc_columbus bsemol # These are the libraries for the high-level modules. They should be # specified in an order that will link correctly, but that shouldn't @@ -36,7 +36,7 @@ KNOWN_MODULE_LIBS = -ltest -lccsd -llucia -lmcscf -lselci -lmp2 -lmoints \ -lcafe -lspace -lanalyze -ldiana -lqhop \ -lpfft -ldplot -lnwpython -ldrdy -lvscf -lqmmm -lrism -lqmd -letrans\ -lpspw -ltce -lgeninterface -lbq -lmm -lcons -lperfm -lnbo -ldntmc -lccca\ - -ldimqm -lrdmft -lfcidump -lgwmol -lnwxtb + -ldimqm -lrdmft -lfcidump -lgwmol -lnwxtb -lbsemol # This include file handles configuration of the NW_MODULE_SUBDIRS and # NW_MODULE_LIBS macros for what we actually want to build. It @@ -77,22 +77,22 @@ EXP_alldev = develop $(EXP_all) EXP_pnnldev = develop $(EXP_pnnl) EXP_qm = NWints atomscf ddscf gradients moints nwdft rimp2 hessian bq cons \ optim stepper driver ccsd vib mcscf property solvation dplot mp2_grad esp etrans tce vscf dimqm \ - gwmol + gwmol bsemol EXP_smallqm = NWints atomscf ddscf gradients moints nwdft rimp2 hessian bq cons \ optim stepper driver vib qmd mcscf property solvation dplot \ - mp2_grad esp etrans vscf nwmd dimqm gwmol + mp2_grad esp etrans vscf nwmd dimqm gwmol bsemol EXP_smallqmmm = NWints atomscf ddscf gradients moints nwdft rimp2 hessian bq cons \ optim stepper driver vib qmd mcscf property solvation dplot \ - mp2_grad esp etrans vscf nwmd qmmm prepar analyz dimqm gwmol + mp2_grad esp etrans vscf nwmd qmmm prepar analyz dimqm gwmol bsemol EXP_notce = NWints atomscf ddscf gradients moints nwdft rimp2 hessian bq cons \ optim stepper driver ccsd vib mcscf property solvation dplot mp2_grad esp etrans vscf dimqm \ - gwmol + gwmol bsemol EXP_tinyqmpw = nwdft property solvation driver ccsd mp2_grad hessian vib nwpw EXP_ver2 = prepar gradients mp2_grad stepper vib mcscf \ optim driver analyz diana #plane_wave ... broken? RJH EXP_qmandpw = nwpw NWints atomscf ddscf gradients moints nwdft rimp2 hessian bq cons \ - optim stepper driver ccsd vib mcscf property solvation dplot mp2_grad esp etrans vscf dimqm gwmol + optim stepper driver ccsd vib mcscf property solvation dplot mp2_grad esp etrans vscf dimqm gwmol bsemol EXP_tps = prepar nwmd analyz diana EXP_md = prepar nwmd analyz ddscf nwdft gradients esp EXP_classical = prepar nwmd analyz @@ -187,6 +187,7 @@ MODULES := $(patsubst nwpw,nwpw fft,$(MODULES) ) MODULES := $(patsubst ddscf,ddscf atomscf nwdft NWints gradients,$(MODULES) ) MODULES := $(patsubst gwmol, gwmol nwdft ddscf gradients hessian atomscf NWints,$(MODULES) ) +MODULES := $(patsubst bsemol, bsemol gwmol nwdft ddscf gradients hessian atomscf NWints,$(MODULES) ) MODULES := $(patsubst lucia, lucia mp2_grad ddscf nwdft property cphf gradients hessian atomscf NWints moints,$(MODULES) ) @@ -274,6 +275,7 @@ MOD_LIBS := $(patsubst dimqm,-ldimqm,$(MOD_LIBS) ) MOD_LIBS := $(patsubst rdmft,-lrdmft,$(MOD_LIBS) ) MOD_LIBS := $(patsubst fcidump,-lfcidump,$(MOD_LIBS) ) MOD_LIBS := $(patsubst gwmol,-lgwmol,$(MOD_LIBS) ) +MOD_LIBS := $(patsubst bsemol,-lbsemol,$(MOD_LIBS) ) MOD_LIBS := $(patsubst xtb,-lnwxtb,$(MOD_LIBS) ) MOD_LIBS := $(patsubst nwc_columbus,-lnwc_columbus,$(MOD_LIBS) ) MOD_LIBS := $(strip $(MOD_LIBS)) @@ -339,6 +341,7 @@ STUBROUTINES := $(patsubst solvation, solvation,$(STUBROUTINES)) STUBROUTINES := $(patsubst dimqm,dimqm_addop dimqm_addop_uhf_damp dimqm_addop_uhf dimqm_tidy dimqm_input dimpar_input dimqm_main fock_dim dimqm_prep dim_grid_quadv0 dimqm_lclfld dimqm_polarizability dimqm_broadcast dimqm_Dens2Dipoles dimqm_rohf_calcDIM dimqm_rohf_wrap dimqm_used dimqm_setgo dimqm_getgo dimqm_getenergy dimqm_setlfirst dimqm_getlrsp dimqm_lrsp dimqm_getlclfld dimqm_getgdipeli dimqm_energy dimqm_init,$(STUBROUTINES) ) #STUBROUTINES := $(patsubst cons, cons_end cons_input cons_add_egrad cons_add_energy cons_init cons_print_summary cons_active,$(STUBROUTINES) ) STUBROUTINES := $(patsubst gwmol, gw_input gwmol, $(STUBROUTINES) ) +STUBROUTINES := $(patsubst bsemol, bse_input bsemol, $(STUBROUTINES) ) STUBROUTINES := $(patsubst xtb, xtb_input xtb_energy xtb_energy_gradient, $(STUBROUTINES) ) STUBROUTINES := $(patsubst nwc_columbus, sifcfg sifwh sodft_aoints, $(STUBROUTINES) ) STUBROUTINES := $(strip $(STUBROUTINES)) diff --git a/src/task/GNUmakefile b/src/task/GNUmakefile index 14fe35f204..0b1b8b9041 100644 --- a/src/task/GNUmakefile +++ b/src/task/GNUmakefile @@ -9,7 +9,7 @@ task_save_state.o task.o task_input.o task_shell_input.o task_num_grad.o \ task_vscf.o task_property.o task_dntmc.o task_bsse.o \ task_jefftce.o task_ncc.o task_rfcf.o \ - task_gwmol.o \ + task_gwmol.o task_bsemol.o task_aoints.o task_nadtrc.o USES_BLAS = task_num_grad.F task_bsse.F diff --git a/src/task/task.F b/src/task/task.F index 5f112136fc..b98b94fc26 100644 --- a/src/task/task.F +++ b/src/task/task.F @@ -163,6 +163,10 @@ c logical task_nadtrc external task_nadtrc + logical task_bsemol + external task_bsemol + + c logical embed_optimize c external embed_optimize c @@ -468,6 +472,8 @@ c status=embed_prepare(rtdb) status=task_dntmc(rtdb) else if (inp_compare(.false., 'gw', operation)) then status = task_gwmol(rtdb) + else if (inp_compare(.false., 'bse', operation)) then + status = task_bsemol(rtdb) c c **** check for theories that have specific operations **** c diff --git a/src/task/task_bsemol.F b/src/task/task_bsemol.F new file mode 100644 index 0000000000..3a3d1fb9c4 --- /dev/null +++ b/src/task/task_bsemol.F @@ -0,0 +1,140 @@ +c +c Task handler for Molecular BSE calcualtions +c Daniel Mejia-Rodriguez, 2022 +c + logical function task_bsemol(rtdb) +c + implicit none +c +#include "errquit.fh" +#include "rtdb.fh" +#include "mafdecls.fh" +#include "inp.fh" +#include "stdio.fh" +#include "util.fh" +#include "global.fh" +c + integer rtdb +c + character*32 theory + character*255 title + integer iproc, ipol + logical oprint, rodft, dftmp2, status, converged, ignore + double precision cputime, walltime +c + logical bsemol, task_energy + external bsemol, task_energy +c +c call ecce_print_module_entry('task gwmol') +c +c Initialization +c + task_bsemol = .false. +c +c Time at the start +c + cputime = util_cpusec() + walltime = util_wallsec() +c +c Invalidate status +c + if (.not. rtdb_put( rtdb, 'task:status', mt_log, 1, .false.)) + & call errquit('task_bse: failed to invalidate status',0,RTDB_ERR) +c +c Print BSE header +c + oprint = util_print('information', print_low) + if (oprint) then + if (.not. rtdb_cget(rtdb, 'title', 1, title)) title = ' ' + if (ga_nodeid().eq.0) then + write(LuOut,*) + call util_print_centered( + & LuOut, 'NWChem BSE Module', 40, .true.) + write(LuOut,*) + write(LuOut,*) + if (title.ne.' ') then + call util_print_centered(LuOut, title, 40, .false.) + write(LuOut,*) + write(LuOut,*) + endif + endif + endif +c +c Check if theory has been specified +c + if (.not. rtdb_cget(rtdb, 'task:theory', 1, theory)) + & call errquit('task_bsemol: theory not specified',0, RTDB_ERR) +c +c Restrict to DFT for now (SCF could also be used) +c + if (.not.inp_compare(.false.,'dft',theory)) + & call errquit('task_bsemol: BSE is only implemented for DFT', 0, + & INPUT_ERR) +c + if (rtdb_get(rtdb, 'dft:rodft', mt_log, 1, rodft)) + & call errquit('task_bsemol: BSE is not implemented for RODFT', 0, + & INPUT_ERR) +c + if (rtdb_get(rtdb, 'dft:dftmp2', mt_log, 1, dftmp2)) + & call errquit('task_bsemol: BSE is not implemented for double'// + & ' hybrid functionals',0,INPUT_ERR) +c + if (.not.rtdb_get(rtdb, 'dft:ipol', mt_int, 1, ipol)) + & call errquit('task_bsemol: failed to read ipol', 0, RTDB_ERR) + + if (.not. rtdb_get(rtdb, 'task:ignore', mt_log, 1, ignore)) + & ignore = .true. +c +c Make sure a DFT solution is available +c + status = rtdb_get(rtdb, 'dft:converged', mt_log, 1, converged) + if (.not.status) then + status = task_energy(rtdb) + if (.not.rtdb_get(rtdb, 'dft:converged', mt_log, 1, converged)) + & call errquit('task_gwmol: failed to read convergence info',0, + & RTDB_ERR) + endif + + status = converged + + + if (ga_nodeid().eq.0) then + write(luout,*) " " + write(luout,'(29X,A)') "***********************" + write(luout,'(29X,A)') "*** BSE ***" + write(luout,'(29X,A)') "***********************" + write(luout,*) " " + + if (.not.converged) then + write(LuOut,*) + write(LuOut,'(27X,A)') '!!! WARNING: !!!' + write(LuOut,'(27X,A)') '!!! DFT DID NOT CONVERGE !!!' + write(LuOut,*) + endif + + call util_flush(luout) + endif + + if ((.not.status) .and. (.not.ignore)) + & call errquit('task_bsemol: Energy failure',0,CALC_ERR) +c +c Actually compute GW +c + status = bsemol(rtdb) +c + cputime = util_cpusec() - cputime + walltime = util_wallsec() - walltime +c + if (.not. rtdb_put(rtdb, 'task:status', mt_log, 1, .true.)) + & call errquit('task_bse: failed to set status',0, + & RTDB_ERR) + if (.not. rtdb_put(rtdb, 'task:cputime', mt_dbl, 1, cputime)) + & call errquit('task_bse: failed storing cputime',0, + & RTDB_ERR) + if (.not. rtdb_put(rtdb, 'task:walltime', mt_dbl, 1, walltime)) + & call errquit('task_bse: failed storing walltime',0, + & RTDB_ERR) +c + task_bsemol = status +c + end function