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 #591 from edoapra/x2c-pr-pull

Browse source 
X2c code
This commit is contained in:
Edoardo Aprà 2022-04-13 23:00:38 -07:00 • committed by GitHub
commit 78f443cb72
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
38 changed files with 6303 additions and 75 deletions
Download patch file Download diff file Expand all files Collapse all files

4
QA/doafewqmtests.mpi
View file

@ -28,6 +28,7 @@ fi
./runtests.mpi.unix procs $np tddft_n2+
./runtests.mpi.unix procs $np tddft_ac_co
./runtests.mpi.unix procs $np hi_zora_sf
./runtests.mpi.unix procs $np ch3radical_unrot
#
# newproperty tests
./runtests.mpi.unix procs $np prop_ch3f
@ -54,3 +55,6 @@ fi
./runtests.mpi.unix procs $np h2o2-response-uhf-damping
./runtests.mpi.unix procs $np h2o2-prop-notrans
#
#--- relativity tests
./runtests.mpi.unix procs $np x2c-h2se
#

3
QA/doqmtests.mpi
View file

@ -328,7 +328,8 @@ np_small=$(get_limit $np 16)
# GW
./runtests.mpi.unix procs $np gw_closedshell
./runtests.mpi.unix procs $np gw_openshell
#
# x2c
./runtests.mpi.unix procs $np x2c-h2se
#--- if we only want quick tests get out now!
if [[ "$what" == "fast" ]]; then
exit

158
QA/tests/x2c-h2se/x2c-h2se.nw Normal file
View file

@ -0,0 +1,158 @@
echo
start testjob
title efgjob
geometry noautoz units bohr nocenter
symmetry c1
Se 0.0000 0.0000 0.0564
H 0.0000 1.0575 -0.9594
H 0.0000 -1.0575 -0.9594
end
BASIS "ao basis" spherical PRINT
# Sapporo-DZP uncontracted
#BASIS SET: (4s,3p) -> [4s,3p]
H S
13.0107003 1.00000000E+00
H S
1.9622571 1.00000000E+00
H S
0.4445380 1.00000000E+00
H S
0.1219496 1.00000000E+00
H P
2.8606730 1.00000000E+00
H P
0.8881900 1.00000000E+00
H P
0.3037760 1.00000000E+00
#BASIS SET: (17s,14p,9d) -> [17s,14p,9d]
Se S
241498.6382409 1.00000000E+00
Se S
36223.7408736 1.00000000E+00
Se S
8243.4547006 1.00000000E+00
Se S
2332.0080828 1.00000000E+00
Se S
760.4702669 1.00000000E+00
Se S
578.5796313 1.00000000E+00
Se S
272.0485848 1.00000000E+00
Se S
179.7486759 1.00000000E+00
Se S
100.7295011 1.00000000E+00
Se S
28.6642790 1.00000000E+00
Se S
19.9031660 1.00000000E+00
Se S
11.8693555 1.00000000E+00
Se S
3.3851284 1.00000000E+00
Se S
2.0728541 1.00000000E+00
Se S
1.4483918 1.00000000E+00
Se S
0.3352235 1.00000000E+00
Se S
0.1251479 1.00000000E+00
Se P
4987.1768459 1.00000000E+00
Se P
1182.1765262 1.00000000E+00
Se P
381.9907152 1.00000000E+00
Se P
144.3524457 1.00000000E+00
Se P
59.8682607 1.00000000E+00
Se P
54.3041566 1.00000000E+00
Se P
26.1507531 1.00000000E+00
Se P
11.1948834 1.00000000E+00
Se P
7.9587015 1.00000000E+00
Se P
3.3433343 1.00000000E+00
Se P
1.3700986 1.00000000E+00
Se P
0.5441071 1.00000000E+00
Se P
0.2107083 1.00000000E+00
Se P
0.0794330 1.00000000E+00
Se D
237.0548228 1.00000000E+00
Se D
70.8504951 1.00000000E+00
Se D
26.4912435 1.00000000E+00
Se D
10.9539115 1.00000000E+00
Se D
4.6411645 1.00000000E+00
Se D
1.9254421 1.00000000E+00
Se D
0.7322816 1.00000000E+00
Se D
0.5653420 1.00000000E+00
Se D
0.2082120 1.00000000E+00
End
charge 0
relativistic
x2c on
#zora on
#douglas-kroll on
zora:cutoff_EFG 1d-10
zora:cutoff 1d-30
#clight 1e6
clight 137.0359895
end
set lindep:n_dep 0
# set dft:no_prune T
set dft:tol_rho 1d-30
set int:cando_txs f
set dft:job_grid_acc 1d-20
dft
grid xfine
direct
#xc slater 1.0 becke88 nonlocal 1.0 lyp 1.0
xc hfexch
convergence energy 1d-10
iterations 50
#print "quadrature"
end
set lindep:n_dep 0
#set sodft:scalar true
#set sodft:spinprojection 0. 0. 1.
#set sodft:collinear true
#set sodft:bfield 0. 1e-4 0.
#set sodft:xcallinone T
#set sodft:alternativexc T
set sodft:debugenergy T
#set sodft:x2c T
set x2c:analyticpvp T
#set x2c:efg T
task dft

1981
QA/tests/x2c-h2se/x2c-h2se.out Normal file
View file

File diff suppressed because it is too large Load diff

7
src/NWints/api/apiP.fh
View file

@ -89,3 +89,10 @@ c
common /DK_log/doug_kroll, dk_even_bas
common /DK_int/dk_form_type, dk_fitbas
common /DK_real/dk_ratio, dk_clight
c
c exact decoupling settings for the API (similar to DK)
c
logical exact_2c
integer x2c_bas
common /X2C_log/exact_2c
common /X2C_int/x2c_bas

13
src/NWints/api/exact_mem.F
View file

@ -29,6 +29,7 @@ c:tex-\end{verbatim}
#include "apiP.fh"
#include "stdio.fh"
#include "util.fh"
#include "errquit.fh"
c
c:: functions
logical api_is_atomic_basis
@ -123,6 +124,18 @@ c
isz_1e = max(isz_1e, lsz_1e)
endif
c
c ... jochen:
c Memory for X2C integrals. exact_2c and x2c_bas are
c common-block variables in apiP.fh
if (exact_2c) then
if (doug_kroll) call errquit(
& 'exact_mem: x2c and DK not allowed',1,INPUT_ERR)
call emem_1e_dk(x2c_bas,x2c_bas,lsz_1e, lmem_1e)
mem_1e = max(mem_1e, lmem_1e)
isz_1e = max(isz_1e, lsz_1e)
endif
c
c.. need extra memory for H1 call externalized by int_mem_h1
c
mem_1e_h1 = mem_1e + isz_1e

44
src/NWints/api/int_init.F
View file

@ -55,6 +55,9 @@ c procedure. The wrapper sets the Douglas-Kroll stuff and adds the new
c fitting basis set to the list of active basis sets and calls
c int_init_org.
c
c ... jochen 01/22: the exact decoupling functionality needs
c the same setup as Douglas-Kroll, so we added that here
c
c::passed
c:tex-\begin{verbatim}
integer rtdb !< [Input] run time data base handle
@ -162,7 +165,46 @@ c
int_mb(mybIm_a+nbas) = dk_fitbas
nrbas = nrbas + 1
endif
if (doug_kroll .or. dyall_mod_dir) call rel_init_con (rtdb)
c ... jochen 01/22: same idea as for douglas-kroll: we initialize
c an additional uncontracted basis.
c in principle we could have used the corresponding D-K code
c and variables, but I want to avoid potentially conflicting
c variable definitions
if (.not. rtdb_get(rtdb,'x2c',MT_LOG,1,exact_2c))
& exact_2c = .false.
c write (luout,*) 'jochen: in int_init x2c = ',exact_2c
if (exact_2c) then
if (doug_kroll) call errquit(
& 'int_init: x2c and DK not allowed',1,INPUT_ERR)
c Generate the fitting basis set and add it to the list of
c used basis sets
do 300 ibas = 1, nbas
if (.not. bas_name(bases(ibas),dk_bas_name,dk_trans_name))
& call errquit('int_init: x2c failed in bas_name',0,
& INT_ERR)
if (dk_bas_name(1:8) .eq. 'ao basis') goto 350
300 continue
call errquit('int_init: x2c could not find ao basis',0,
& INT_ERR)
350 call bas_create_uncont(rtdb,bases(ibas),'auto X2C basis',
& x2c_bas)
c Add fitting basis set to list of loaded basis sets
int_mb(mybIm_a+nbas) = x2c_bas
nrbas = nrbas + 1
endif ! exact_2c
c ... jochen end
c rel_init_con initializes some relativistic constants
if (doug_kroll .or. dyall_mod_dir .or. exact_2c)
& call rel_init_con (rtdb)
call int_init_org(rtdb, nrbas, int_mb(mybIm_a))
if (.not. ma_pop_stack(mybm_a))
$ call errquit('int_init: ma_pop_stack mybm_a',0, INT_ERR)

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

@ -69,6 +69,14 @@ c
& ,0, INT_ERR)
dk_fitbas = -1
endif
c
c terminate uncontracted basis used in the X2C code
if (exact_2c) then
if (.not. bas_destroy(x2c_bas)) call errquit(
& 'int_terminate: failed to destroy x2c uncontracted basis'
& ,0, INT_ERR)
x2c_bas = -1
endif
c
end
C>

6
src/NWints/rel/GNUmakefile
View file

@ -3,7 +3,8 @@
HEADERS = rel_consts.fh \
rel_nwc.fh \
zora.fh \
modelpotential.fh
modelpotential.fh \
x2c.fh
OBJ_ENERGY = \
rel_onel.o \
@ -17,7 +18,8 @@ OBJ_ENERGY = \
rel_SSSS.o \
rel_init_con.o \
dmd_input.o \
zora_input.o
zora_input.o \
x2c_input.o
OBJ_GRAD = \
rel_oneld.o \

20
src/NWints/rel/rel_input.F
View file

@ -14,6 +14,7 @@
#include "errquit.fh"
#include "zora.fh"
#include "modelpotential.fh"
#include "x2c.fh"
c
c Read input for the relativistic block
c
@ -25,6 +26,7 @@ c
logical dk_run ! Douglas-Kroll
logical dmd_run ! Dyall modified Dirac
logical zora_run ! ZORA
logical x2c_run ! exact decoupling (X2C)
integer mpoption
c
logical status
@ -56,6 +58,11 @@ c ZORA
if (.not. rtdb_put(rtdb, 'zora', mt_log, 1, zora_run))
& call errquit('rel_input: rtdb_put failed', 0, RTDB_ERR)
c
c X2C
x2c_run = .false.
if (.not. rtdb_put(rtdb, 'x2c', mt_log, 1, x2c_run))
& call errquit('rel_input: rtdb_put failed', 0, RTDB_ERR)
c
c Set some relativistic variables we could need later on
c
1 if (.not. inp_read()) call errquit('rel_input: inp_read failed',0,
@ -75,7 +82,7 @@ c
endif
write(luout,10) clight
10 format(/,'The speed of light used in the relativistic ',
$ 'calculations is ',F16.10, 'au.'/)
$ 'calculations is ',E15.7, 'au.'/)
c
elseif (inp_compare(.false.,'douglas-kroll', test)) then
dk_run = .true.
@ -88,6 +95,12 @@ c
if (.not. rtdb_put(rtdb, 'dyall_mod_dir', mt_log, 1, dmd_run))
$ call errquit('rel_input: rtdb_put failed', 0, RTDB_ERR)
call dmd_input(rtdb)
c
elseif (inp_compare(.false.,'x2c', test)) then
x2c_run = .true.
if (.not. rtdb_put(rtdb, 'x2c', mt_log, 1, x2c_run))
& call errquit('rel_input: rtdb_put failed', 0, RTDB_ERR)
call x2c_input(rtdb)
c
elseif (inp_compare(.false.,'zora', test)) then
call zora_input(rtdb)
@ -107,6 +120,11 @@ c
if(.not. inp_f(zoracutoff)) zoracutoff = 1.0d-15 ! set a default cutoff
if (.not. rtdb_put(rtdb,'zora:cutoff',mt_dbl,1,zoracutoff))
& call errquit('rel_input: rtdb put failed',0, RTDB_ERR)
c
elseif(inp_compare(.false.,'x2c:cutoff',test)) then
if(.not. inp_f(x2ccutoff)) x2ccutoff = 1.0d-15 ! set a default cutoff
if (.not. rtdb_put(rtdb,'x2c:cutoff',mt_dbl,1,x2ccutoff))
& call errquit('rel_input: rtdb put failed',0, RTDB_ERR)
c
c zoracutoff_EFG, used in get_rhoS() from dft_zora_utils.F
elseif(inp_compare(.false.,'zora:cutoff_EFG',test)) then

15
src/NWints/rel/x2c.fh Normal file
View file

@ -0,0 +1,15 @@
c exact decoupling (X2C) related variables
logical do_x2c, scale_x2c_hartree, scale_x2c_so_hartree
common/x2c_flags/
& do_x2c,
& scale_x2c_hartree,
& scale_x2c_so_hartree
c scaling for Hartree
double precision x2chscale, x2csoscale
common/x2c_params/x2chscale, x2csoscale
c cutoff
double precision x2ccutoff
common/x2c_tols/x2ccutoff

91
src/NWints/rel/x2c_input.F Normal file
View file

@ -0,0 +1,91 @@
c Handles the X2C input
subroutine x2c_input(rtdb)
implicit none
#include "rtdb.fh"
#include "context.fh"
#include "global.fh"
#include "mafdecls.fh"
#include "inp.fh"
#include "stdio.fh"
#include "rel_consts.fh"
#include "util.fh"
#include "errquit.fh"
#include "x2c.fh"
c
c
c Read input for calculations with X2C (exact decoupling)
c
c All input is of the form:
c <keyword> <on||off>
c
integer rtdb ! [input] handle to database
c
character*255 test
double precision rtmp
c
c Preliminaries
do_x2c = .false.
scale_x2c_hartree = .false.
scale_x2c_so_hartree = .false.
x2chscale=1.0d0 ! scaling for Hartree potential
x2csoscale=1.0d0 ! scaling in the SO operator
c
c Read input data beyond x2c; store in rtdb.
10 if (inp_a(test)) then
if (inp_compare(.false.,'on',test)) then
do_x2c = .true.
else if (inp_compare(.false.,'off',test)) then
do_x2c = .false.
c
c process other x2c options
c
c Hartree scaling?
else if (inp_compare(.false.,'scalehartree',test)) then
scale_x2c_hartree = .true.
rtmp = 1.0d0
if (.not.inp_f(rtmp)) then
write (luout,*) 'rtmp = ',rtmp
call errquit('x2c: error reading Hartree scaling factor',
& 1, INPUT_ERR)
else
x2chscale = rtmp
end if
c Hartree scaling for SO only?
else if (inp_compare(.false.,'scalesoee',test)) then
scale_x2c_so_hartree = .true.
rtmp = 1.0d0
if (.not.inp_f(rtmp)) then
write (luout,*) 'rtmp = ',rtmp
call errquit('x2c: error reading SOee scaling factor',
& 1, INPUT_ERR)
else
x2csoscale = rtmp
end if
c
c keyword not recognized
else
call errquit('x2c_input: unknown directive',0, UNKNOWN_ERR)
endif
goto 10
endif
c
c Put x2c parameters in rtdb
if (.not. rtdb_put(rtdb, 'x2c:on', mt_log, 1,do_x2c))
& call errquit('x2c_input: rtdb_put failed', 0, RTDB_ERR)
c
if (scale_x2c_hartree) then
write (luout,*) 'X2C: scaling for Hartree term =', x2chscale
end if
if (scale_x2c_so_hartree) then
write (luout,*) 'X2C: scaling for SO Hartree term =',
& x2csoscale
end if
c
return
end

2
src/config/make_nwchem_config
View file

@ -318,7 +318,7 @@ STUBROUTINES := $(patsubst rism,task_rism rism_input,$(STUBROUTINES) )
STUBROUTINES := $(patsubst qmmm, qmmm_input task_qmmm_fep task_qmmm_abc qmmm_print_energy1 qmmm_forces qmmm_init qmmm_end task_qmmm_energy task_qmmm_gradient task_qmmm_optimize task_qmmm_dynamics task_qmmm_esp task_qmmm_pes task_qmmm_thermo qmmm_bq_forces,$(STUBROUTINES) )
STUBROUTINES := $(patsubst mm, mm_add_energy mm_add_egrad mm_input mmi_init mmi_end mm_task_gradient mm_task_energy,$(STUBROUTINES) )
STUBROUTINES := $(patsubst dplot,dplot dplot_input player,$(STUBROUTINES) )
STUBROUTINES := $(patsubst property,property raman_input task_raman prop_input hnd_nwhnd_tran hnd_prp_vec_read hnd_diag hnd_elfcon_symm hnd_elfcon hnd_elfcon_fde get_slctd_atoms xc_tabcd_giao_nmr_rhs int_giao_1ega fiao_b1_movecs giao_b1_movecs,$(STUBROUTINES) )
STUBROUTINES := $(patsubst property,property raman_input task_raman prop_input hnd_nwhnd_tran hnd_prp_vec_read hnd_diag hnd_elfcon_symm hnd_elfcon hnd_elfcon_fde get_slctd_atoms xc_tabcd_giao_nmr_rhs int_giao_1ega fiao_b1_movecs giao_b1_movecs int_mpole_1ega,$(STUBROUTINES) )
STUBROUTINES := $(patsubst solvation,cosmo_input cosmo_initialize cosmo_charges cosmo_charges_from_dmat cosmo_def_radii cosmo_tidy cosff cosdff grad_hnd_cos cosmo_smd_output cosmo_vem_output,$(STUBROUTINES) )
STUBROUTINES := $(patsubst vscf,task_vscf vscf_input,$(STUBROUTINES) )
STUBROUTINES := $(patsubst dntmc,task_dntmc dntmc_input,$(STUBROUTINES) )

2
src/nwdft/GNUmakefile
View file

@ -7,7 +7,7 @@
LIB_DEFINES =
SUBDIRS = include input_dft scf_dft scf_dft_cg so_dft coulomb xc util grid dftgrad lr_tddft lr_tddft_grad zora spec rt_tddft libxc
SUBDIRS = include input_dft scf_dft scf_dft_cg so_dft coulomb xc util grid dftgrad lr_tddft lr_tddft_grad zora spec rt_tddft libxc x2c
LIB_INCLUDES = -I./include -I./grid
HEADERS = ../ddscf/cosmo.fh

2
src/nwdft/scf_dft/GNUmakefile
View file

@ -40,7 +40,7 @@
dft_fant.o\
dft_printinfo.o\
dft_fock_efield.o\
dft_frac_mo.o dft_nwmolden.o dft_spinflip.o
dft_frac_mo.o dft_nwmolden.o dft_spinflip.o dft_quickguess.o
LIBRARY = libnwdft.a
#_USE_SCALAPACK = $(shell cat ${NWCHEM_TOP}/src/tools/build/config.h | awk ' /HAVE_SCALAPACK\ 1/ {print "Y"}')

29
src/nwdft/scf_dft/dft_fockbld.F
View file

@ -9,8 +9,8 @@
, i_degen,g_rho_exact,g_tmp_exact, g_addit_exact,
, n_levels, act_levels, oprint_time,
cc AJL/Begin/FDE
c , g_zora_sf,g_frozemb)
, g_zora_sf,g_frozemb,g_vxc_fde,g_xcinv_fde,g_dens_fde,eelec )
c , g_rel_sf,g_frozemb)
, g_rel_sf,g_frozemb,g_vxc_fde,g_xcinv_fde,g_dens_fde,eelec )
cc AJL/End
c
c $Id$
@ -26,6 +26,7 @@ c
#include "dftps.fh"
#include "cosmo.fh"
#include "zora.fh"
#include "x2c.fh"
#include "case.fh"
c#include "dimqm.fh"
c
@ -65,8 +66,8 @@ c
c
double precision Ebq
c
integer g_zora_sf(*) ! spin-free zora ga
double precision Ezora_sf ! spin-free zora energy
integer g_rel_sf(*) ! spin-free zora/x2c ga
double precision Erel_sf ! spin-free zora/x2c energy
c
integer g_frozemb ! frozen embedding ga
double precision Efrozemb ! energy
@ -446,22 +447,22 @@ c == Calculate the HF exchange (exact Coulomb), DFT exchange-correlation ==
c
if (odftps) call pstat_off(ps_xc)
c
c add zora contributions
if (do_zora) then
Ezora_sf =0.d0
c add relativistic energy contributions
if (do_zora .or. do_x2c) then
Erel_sf =0.d0
if (ipol.gt.1) then
Ezora_sf = ga_ddot(g_dens(1),g_zora_sf(1))
& + ga_ddot(g_dens(2),g_zora_sf(2))
Erel_sf = ga_ddot(g_dens(1),g_rel_sf(1))
& + ga_ddot(g_dens(2),g_rel_sf(2))
else
Ezora_sf = ga_ddot(g_dens(1),g_zora_sf(1))
Erel_sf = ga_ddot(g_dens(1),g_rel_sf(1))
end if
Ecore = Ecore + Ezora_sf
Ecore = Ecore + Erel_sf
c
call ga_dadd(1.d0,g_vxc(1),1.d0,g_zora_sf(1),g_vxc(1))
call ga_dadd(1.d0,g_vxc(1),1.d0,g_rel_sf(1),g_vxc(1))
if (ipol.gt.1)
& call ga_dadd(1.d0,g_vxc(2),1.d0,g_zora_sf(2),g_vxc(2))
& call ga_dadd(1.d0,g_vxc(2),1.d0,g_rel_sf(2),g_vxc(2))
c
end if ! do_zora
end if ! do_zora .or. do_x2c
c
c == if frozen embedding ==
Efrozemb = 0.d0

12
src/nwdft/scf_dft/dft_main0d.F
View file

@ -34,6 +34,7 @@ c
c
#include "zora.fh"
#include "rel_consts.fh"
#include "x2c.fh"
#include "case.fh"
c
integer ga_create_atom_blocked
@ -212,6 +213,17 @@ c == Get the grid cutoff or set default ==
if(.not.rtdb_get(rtdb,'zora:not_zora_scale', ! FA
& mt_log,1,not_zora_scale)) ! FA
& not_zora_scale = .false. ! FA
c ... jochen:
c == X2C options ==
do_x2c = .false.
if (rtdb_get(rtdb,'x2c:on',mt_log,1,do_x2c))
& do_x2c = .true.
c we use the variable zoracutoff also for the x2c option
if(.not.rtdb_get(rtdb,'x2c:cutoff',mt_dbl,1,zoracutoff))
& zoracutoff = 1.0d-08
c
c
c == Get the light speed and other constants, if specified ==
call rel_init_con(rtdb)

41
src/nwdft/scf_dft/dft_quickguess.F Normal file
View file

@ -0,0 +1,41 @@
subroutine dft_quickguess(ipol,nbf_ao,g_dens,oskel,
N ntotel,nclosed,
D densmat)
implicit none
#include "global.fh"
#include "errquit.fh"
integer ipol
integer nbf_ao
logical oskel
integer ntotel,nclosed
integer g_dens(*)
logical densmat ![in]
c
character*256 file_mat1
c
if(densmat) then
call util_file_name('densmat', .false.,.false.,file_mat1)
if(ga_nodeid().eq.0)
W write(6,*) ' reading DM from file ',file_mat1
call ao_1prdm_read(nbf_ao,g_dens(1),file_mat1)
endif
if(ipol.eq.2) then
if(densmat) then
call errquit(' code not ready for open-shell ',0,CAPMIS_ERR)
call util_file_name('densmat', .false.,.false.,file_mat1)
if(ga_nodeid().eq.0)
W write(6,*) ' reading DM from file ',file_mat1
call ao_1prdm_read(nbf_ao,g_dens(2),file_mat1)
else
call ga_copy(g_dens(1),g_dens(2))
call ga_dscal(g_dens(1),dble(ntotel-nclosed)/(ntotel))
call ga_dscal(g_dens(2),dble(nclosed)/(ntotel))
endif
endif
if(oskel) then
call ga_symmetrize(g_dens(1))
if(ipol.eq.2) call ga_symmetrize(g_dens(2))
endif
return
end

172
src/nwdft/scf_dft/dft_scf.F
View file

@ -40,6 +40,7 @@ c
#include "dftps.fh"
#include "cosmo.fh"
#include "zora.fh"
#include "x2c.fh"
#include "modelpotential.fh"
cold#include "dimqm.fh"
c
@ -191,7 +192,6 @@ c == zora related ==
logical dft_zora_read, dft_zora_write, dft_zora_inquire_file
external dft_zora_read, dft_zora_write, dft_zora_inquire_file
character*255 zorafilename
integer g_zora_sf(2)
integer g_zora_scale_sf(2)
double precision ener_scal
integer icalczora
@ -225,12 +225,18 @@ cc AJL/End
c
gotsmear(ss)=abs(ss).gt.1d-9
c
c ... exact decoupling and other relativistic Hamiltonians
integer g_rel(2)
integer g_x2c_u(2), g_x2c_u2c, g_x2c_c2u
c
c misc
c
logical debug_fon
logical do_zora_recalc, zora_recalc
external zora_recalc
c
external zora_recalcc
logical x2c_noatm
external x2c_noatm
c MN solvation models -->
c
double precision gcds, gspol, espol, gspoldyn, gspolneq
@ -282,6 +288,7 @@ c
c <-- MN solvation models
c
c =========================================================================
c
c
k_3cERI = 0
k_3cwhat = 0
@ -317,6 +324,38 @@ c
me = ga_nodeid()
nproc = ga_nnodes()
c ---------------------------------------------------------
c check which relativistic Hamiltonian (if any) we're using
c and write some information to the output.
c comment: Niri's zora implementation added an array g_zora_sf
c to the Fock matrix in dft_fockbld. This array contains
c T(rel) - T(nrel). The nonrelativistic kinetic energy is always
c added to the Fock matrix in dft_fockbld.
c I changed the array name for g_zora_sf to g_rel and use
c it now for both methods.
c ---------------------------------------------------------
c make sure we don't use zora, x2c, or another
c two-component relativistic method at the same time
if (do_zora .and. do_x2c) then
call errquit ('dft_scf: zora and x2c options conflict',
& 666, INPUT_ERR)
end if
c write output about the relativistic treatment
if (do_zora) then
if (me.eq.0) write (luout,*) 'dft: using ZORA'
end if
if (do_x2c) then
if (me.eq.0) write (luout,*) 'dft: using X2C'
end if
call x2c_scf_init(rtdb)
c
c !!! BGJ
c Store SCF hamiltonian type as DFT for use in BGJ routines
@ -651,6 +690,23 @@ c
c == force sleazy SCF into "direct" mode ==
call fock_force_direct(rtdb)
c
c -------------------------------------------
c calculate superposition of atomic densities
c if ZORA or X2C functionality is requested
c -------------------------------------------
if ((do_zora .or. do_x2c) .and. .not.x2c_noatm()) then
call ga_zero(g_dens_at(1))
if (ipol.gt.1) call ga_zero(g_dens_at(2))
call guess_dens(rtdb, geom, ao_bas_han, g_dens_at)
call dft_quickguess(ipol,nbf_ao,g_dens_at,oskel,
N ntotel,nclosed,
D .false.)
else if (do_x2c .and. x2c_noatm()) then
call ga_zero(g_dens_at(1))
if (ipol.gt.1) call ga_zero(g_dens_at(2))
end if ! do_zora .or. do_x2c
c
c == get atomic zora contributions ==
if (.not.rtdb_get(rtdb,'zora',mt_log,1,do_zora)) do_zora = .false.
c
@ -665,13 +721,13 @@ c == get filename for the zora data ==
call util_file_name('zora_sf',.false.,.false.,zorafilename)
c
c == zora arrays ==
if(.not.ga_duplicate(g_vxc(1),g_zora_sf(1),'sf 1'))
if(.not.ga_duplicate(g_vxc(1),g_rel(1),'sf 1'))
& call errquit('dft_scf: ga_duplicate failed',1, GA_ERR)
call ga_zero(g_zora_sf(1))
call ga_zero(g_rel(1))
if(ipol.gt.1) then
if(.not.ga_duplicate(g_vxc(2),g_zora_sf(2),'sf 2'))
if(.not.ga_duplicate(g_vxc(2),g_rel(2),'sf 2'))
& call errquit('dft_scf: ga_duplicate failed',1, GA_ERR)
call ga_zero(g_zora_sf(2))
call ga_zero(g_rel(2))
endif
c
c == zora energy scaling arrays ==
@ -691,18 +747,6 @@ c == create g_zora_Kinetic array ==
call ga_zero(g_zora_Kinetic(i))
enddo
c
c == generate an superposition of atomic densities ==
call ga_zero(g_dens_at(1))
if (ipol.gt.1) call ga_zero(g_dens_at(2))
call guess_dens(rtdb, geom, ao_bas_han, g_dens_at)
if (oskel) call ga_symmetrize(g_dens_at(1))
if(ipol.gt.1) then
call ga_copy(g_dens_at(1),g_dens_at(2))
call ga_dscal(g_dens_at(1),dble(ntotel-nclosed)/(ntotel))
call ga_dscal(g_dens_at(2),dble(nclosed)/(ntotel))
if(oskel) call ga_symmetrize(g_dens_at(2))
end if
c
c == in case fon is used together with zora ==
c pstrace is queried in the grid code ==
if (fon) then
@ -715,7 +759,7 @@ c
c == try reading the zora atomic corrections from file ==
icalczora = 0 ! initialize the flag
if (.not.dft_zora_read(zorafilename, nbf_ao, ipol, nmo,
& mult, g_zora_sf, g_zora_scale_sf)) icalczora = 1
& mult, g_rel, g_zora_scale_sf)) icalczora = 1
c
c check if zora corrections have to recalculated
do_zora_recalc = zora_recalc(rtdb)
@ -764,7 +808,7 @@ c 11 format('In dft_scf:: zetanuc_arr(',i3,')=',f35.8)
dbl_mb(k_zetanuc+iatom-1)=dsqrt(dbl_mb(k_zetanuc+iatom-1)) ! Calc sqrt(zetanuc)
enddo ! end-loop-iatom
c------- FA-code-added-09-26-11 --------------------------- END
call zora_getv_sf(rtdb, g_dens_at, g_zora_sf,g_zora_scale_sf,
call zora_getv_sf(rtdb, g_dens_at, g_rel,g_zora_scale_sf,
& ofinite,dbl_mb(k_zetanuc),
& Knucl, nexc)
if (.not. MA_free_heap(l_zetanuc)) call errquit
@ -772,7 +816,7 @@ c------- FA-code-added-09-26-11 --------------------------- END
c
c == write out the atomic zora corrections to file ==
if (.not.dft_zora_write(rtdb, ao_bas_han, zorafilename,
& nbf_ao, ipol, nmo, mult, g_zora_sf, g_zora_scale_sf))
& nbf_ao, ipol, nmo, mult, g_rel, g_zora_scale_sf))
& call errquit('dft_scf: dft_zora_write failed', 0, DISK_ERR)
end if ! icalczora
c
@ -780,10 +824,43 @@ c == NOTE.-Storing g_zora_Kinetic to be used in rohf_fock.F ==
c == for cphf routine - FA-11-08-10 ==
c == -> g_zora_Kinetic is defined in zora.fh ==
do i=1,ipol
call ga_copy(g_zora_sf(i),g_zora_Kinetic(i))
call ga_copy(g_rel(i),g_zora_Kinetic(i))
enddo
end if ! do_zora
c
c -------------------------------------------------
c set up an exact decoupling Hamiltonian (optional)
c -------------------------------------------------
if (do_x2c) then
if (me.eq.0) write (luout,*) 'allocating g_x2c',nbf_ao
if(.not.ga_duplicate(g_vxc(1),g_rel(1),'x2c 1'))
& call errquit('dft_scf: ga_duplicate failed',666, GA_ERR)
call ga_zero(g_rel(1))
if (me.eq.0) write (luout,*) 'calling calc_x2c_1e_scalar'
c note: g_vxc is used as temp array in the numerical
c integration routines.
c g_x2c_u2c is created in the x2c routine and survives
c if keepU=.T. Same for the transformation matrices g_xc_u(:).
c optionally, g_x2c_c2u is also allocated and kept.
call calc_x2c_1e_scalar_ga(rtdb,
& nexc, g_dens_at, g_vxc, g_rel(1), g_x2c_u,
& g_x2c_u2c, g_x2c_c2u, toll_s)
if (ipol.gt.1) then
if(.not.ga_duplicate(g_vxc(2),g_rel(2),'x2c 2'))
& call errquit('dft_scf: ga_duplicate failed',1, GA_ERR)
call ga_copy(g_rel(1),g_rel(2))
end if
end if ! do_x2c
c
c == determine guess ==
if(.not.rtdb_get(rtdb,'forceatguess',mt_log,1,forceatguess))
. forceatguess=.false.
@ -826,30 +903,9 @@ c == better have static ldb at high node counts ==
if (.not. rtdb_cput(rtdb, 'dft:output vectors', 1, movecs_out))
$ call errquit('dft_scf: DFT MO vectors not defined',0,
& RTDB_ERR)
if(densmat) then
call util_file_name('densmat', .false.,.false.,file_mat1)
if(ga_nodeid().eq.0)
W write(6,*) ' reading DM from file ',file_mat1
call ao_1prdm_read(nbf_ao,g_dens(1),file_mat1)
endif
if(ipol.eq.2) then
if(densmat) then
call errquit(' code not ready for open-shell ',0,CAPMIS_ERR)
call util_file_name('densmat', .false.,.false.,file_mat1)
if(ga_nodeid().eq.0)
W write(6,*) ' reading DM from file ',file_mat1
call ao_1prdm_read(nbf_ao,g_dens(2),file_mat1)
else
call ga_copy(g_dens(1),g_dens(2))
call ga_dscal(g_dens(1),dble(ntotel-nclosed)/(ntotel))
call ga_dscal(g_dens(2),dble(nclosed)/(ntotel))
endif
endif
if(oskel) then
call ga_symmetrize(g_dens(1))
if(ipol.eq.2) call ga_symmetrize(g_dens(2))
endif
call dft_quickguess(ipol,nbf_ao,g_dens,oskel,
N ntotel,nclosed,
D densmat)
goto 1789
endif
if(oskel) then
@ -1227,8 +1283,8 @@ c remove flipping to avoid multiple operation in geom opt.
, i_degen,g_rho_exact,g_tmp_exact, g_addit_exact,
, n_levels, act_levels, oprint_time,
cc AJL/Begin/FDE
c & g_zora_sf,g_frozemb)
& g_zora_sf,g_frozemb,g_vxc_fde,g_xcinv_fde,g_dens_fde,eelec)
c & g_rel,g_frozemb)
& g_rel,g_frozemb,g_vxc_fde,g_xcinv_fde,g_dens_fde,eelec)
cc AJL/End
if(util_print('dft timings', print_high))
& time2=util_wallsec() ! end fock matrix build time
@ -1889,6 +1945,12 @@ c AOR begin
endif
c AOR end
c ... jochen: calculate x2c properties if so requested:
if (do_x2c) then
call x2c_prop(rtdb,
N nexc, g_dens, g_dens_at, g_x2c_u)
end if ! do_x2c
c vdw bit
c
c activate disp if is present in rtdb
@ -2498,18 +2560,24 @@ cdft...Deallocate cdft arrays
c
c Deallocate zora related arrays
if (do_zora) then
if (.not. ga_destroy(g_zora_sf(1))) call errquit(
if (.not. ga_destroy(g_rel(1))) call errquit(
& 'dft_scf: ga_destroy failed ',0, GA_ERR)
if (.not. ga_destroy(g_zora_scale_sf(1))) call errquit(
& 'dft_scf: ga_destroy failed ',0, GA_ERR)
if(ipol.gt.1) then
if (.not. ga_destroy(g_zora_sf(2))) call errquit(
if (.not. ga_destroy(g_rel(2))) call errquit(
& 'dft_scf: ga_destroy failed ',0, GA_ERR)
if (.not. ga_destroy(g_zora_scale_sf(2))) call errquit(
& 'dft_scf: ga_destroy failed ',0, GA_ERR)
end if
end if
c
c deallocate x2c arrays if necessary
if (do_x2c) then
call x2c_scf_terminate(rtdb,ipol,g_rel,
g g_x2c_u,g_x2c_u2c,g_x2c_c2u)
end if ! do_x2c
c
c AOR begin
c deallocate xdm arrays
if (lxdm.eq.1) then

28
src/nwdft/x2c/GNUmakefile Normal file
View file

@ -0,0 +1,28 @@
#
#
OBJ = calc_x2c_w_scalar.o \
x2c_getw_scalar.o \
x2c_getw_efg_scalar.o \
x2c_getw_moment_scalar.o \
calc_x2c_1e_scalar_ga.o \
calc_x2c_efg_scalar.o \
calc_x2c_moment_scalar.o \
calc_x2c_wefg_scalar.o \
calc_x2c_wmoment_scalar.o \
write_x2c_pc_dip.o \
read_x2c_pc_dip.o \
ma_solve.o x2c_scf.o
LIBRARY = libnwdft.a
# LIB_DEFINES = -DNWCHEMX2C
LIB_DEFINES =
LIB_INCLUDES = -I../include -I../../ddscf -I../grid
USES_BLAS =
include ../../config/makefile.h
include ../../config/makelib.h

981
src/nwdft/x2c/calc_x2c_1e_scalar_ga.F Normal file
View file

@ -0,0 +1,981 @@
subroutine calc_x2c_1e_scalar_ga (rtdb,
& nexc, g_dens, g_vxc, g_x2c, g_u, g_u2c, g_c2u, toll_s)
c -----------------------------------------------------------------
c Purpose: calculate scalar X2C one-electron Hamiltonian matrix
c Integrals are supposed to be initialized by the calling
c routine
c Driver (c) 2012, 2015, 2022 J. Autschbach, jochena@buffalo.edu
c 'relham' code used for the original (2012) implementation
c (c) 2012 Daoling Peng. This code is now obsolete in NWChem.
c Contracted basis extension: J. Autschbach, 2015
c New code using GA exclusively (c) 2022 J. Autschbach
c Code is not yet (2022) well tested. Expect bugs & email me at the
c address shown above if you find one. Please provide a minimal
c working example (in-,output) for your bug report, if possible.
c The relevant equations implemented here can be found in the
c following works, among others:
c Autschbach, Peng, Reiher, J. Chem. Theory Comput. 2012, 8,
c 4239–4248. (initial NWChem implementation, field gradient calcs.)
c Liu & Peng, J. Chem. Phys. 2009, 131, 031104. (X2C formalism in
c non-orthogonal basis)
c Peng & Reiher, Theor. Chem. Acc. 2012, 131, 1081. (comparison of
c X2C approaches)
c Autschbach, J., Quantum Theory for Chemical Applications, Oxford
c University Press (2021), Chapter 24. (further details re. the
c formalism)
c The workflow for mapping onto an uncontracted basis and back is
c similar to int_1edk_ga.F. We have the problem here that anything
c that is calculated by numerical integration needs to be done in
c the contracted basis, because of the way the grid routines are
c written. I.e. unlike in int_1edk_ga.F we also need a
c back-transformation. For the time being, we use a naive
c left-inverse for this. Eventually, the numerical integrations
c will be replaced completely.
c -----------------------------------------------------------------
implicit none
#include "mafdecls.fh"
#include "rtdb.fh"
#include "apiP.fh"
#include "global.fh"
#include "tcgmsg.fh"
#include "bas.fh"
#include "util.fh"
#include "stdio.fh"
#include "errquit.fh"
#include "dftpara.fh"
#include "cdft.fh"
#include "rel_consts.fh"
#include "cint1cache.fh"
#include "numerical_constants.fh"
#include "x2c_scf.fh"
c subroutine arguments:
integer rtdb ! [inp] runtime database GA handle
integer nexc ! [inp] no. of XC terms
integer g_dens(2) ! [inp] atom densities for model pot.
integer g_vxc(4) ! scratch
integer g_x2c ! [out] exact decoupling 1e operator
integer g_u(2) ! [out] decoupling transformation
integer g_u2c ! [out] basis trafo
integer g_c2u ! [out] basis trafo
double precision toll_s ! [inp] Smat diag. tolerance
c GA handles used for temp storage in this routine:
integer g_t, g_s, g_w, g_v, g_temp, g_tmpv, g_X
integer g_smix, g_densu
integer g_tm12, g_sm12, g_h, g_m
c local MA handles
integer l_eval, k_eval
integer l_scr, k_scr
c other local variables:
integer ubas, cbas
integer nu, nu2, nunu, nc
integer me, type, i, j, ntemp
double precision sum, energy, rtemp
double precision c_2c, c_2c2, c_4c2
double precision ga_trace_diag
external ga_trace_diag
integer info
integer ga_cholesky,ga_llt_i
external ga_cholesky,ga_llt_i
double precision tol2e
logical master, analytic_pvp, modelpot_2e, use_ga, oskel_sav
integer dbg
double precision jfac(1),kfac(1)
character*13 pname
c =================================================================
c dbg = 1: debug numbers and flow. 2: debug also matrices
c dbg = 0: no debug output
dbg=1
me=ga_nodeid()
master = me.eq.0
pname = "x2c_1e_scalar"
if (dbg>0.and.master) write(luout,*) 'entering '//pname
call ga_sync()
call ga_zero(g_x2c)
tol2e=10.d0**(-itol2e)
c src/NWints/dk/int_1edk_ga.F uses similar transformations as used
c here. There, oskel is deliberately set to .false., so I'm guessing
c we should do the same here but I'm not sure yet
oskel_sav = oskel
oskel = .false.
c temporarily add alpha and beta model densities:
c TO BE REPLACED WITH 4-COMPONENT MODEL DENSITIES
if (ipol.gt.1)
& call ga_dadd(1.d0, g_dens(1), 1.d0, g_dens(2), g_dens(1))
c ---------------
c initializations
c ---------------
ubas = x2c_bas ! uncontracted basis
cbas = ao_bas_han ! contracted basis
c make sure the input GAs have the correct dimension
call ga_inquire(g_x2c,type,nc,nc)
if (nc.ne.nbf_ao) then
call errquit(pname//': g_x2c has wrong dimension',nc,GA_ERR)
end if
if (.not. bas_numbf(cbas,ntemp))
& call errquit(pname//': bas_ao error',ntemp, BASIS_ERR)
if (ntemp.ne.nc)
& call errquit(pname//': nc.ne.ntemp',ntemp, BASIS_ERR)
if (.not. bas_numbf(ubas,nu))
& call errquit(pname//': nu error',nu, BASIS_ERR)
if (master.and.dbg>0) write (luout,*)
& pname//': nu, nc =',nu,nc
if (nu.lt.nc) call errquit(
& 'x2c uncontracted basis smaller than contracted', 1,
& RTDB_ERR)
lx2c_u2c = (nu.gt.nc)
if (master) then
if (lx2c_u2c) then
write (luout,*)
& 'x2c INFO: the operator will be calculated in an'
write (luout,*)
& ' uncontracted auxiliary basis'
else
write (luout,*)
& 'x2c INFO: Assuming that the AO basis is fully uncontracted'
end if
end if
c no of elements in the AO matrices and dimension of mDE:
nunu = nu * nu
nu2 = 2* nu
c ---------------------------
c read some options from rtdb
c ---------------------------
c analytic pVp integrals are default, but it can be disabled
analytic_pvp = .true.
if (.not. rtdb_get(rtdb, 'x2c:analyticpvp',
& mt_log, 1, analytic_pvp)) continue
c model potential for 2e terms?
modelpot_2e = .false.
if (.not. rtdb_get(rtdb, 'x2c:modelpot2e',
& mt_log, 1, modelpot_2e)) continue
c if analyticpvp is .F. or if we need the model potential,
c we also need a contracted-to-uncontracted basis transformation
c unless the two bases are identical
lx2c_c2u = ((.not.analytic_pvp .or. modelpot_2e).and.lx2c_u2c)
if (master.and.lx2c_c2u) write (luout,*)
& 'x2c INFO: x2c options require constructing the c2u trafo'
c check rtdb switch to see whether we should perform decoupling
c in GA, or call relham routines with MA arrays
use_ga = .false.
if (.not.rtdb_get(rtdb, 'x2c:usega', mt_log, 1,
& use_ga)) continue
if (master.and.use_ga) write(luout,*)
& pname//': will perform decoupling in GA mode'
c -----------------------------------------
c allocate GAs for T, V, S, and for W = pVp
c -----------------------------------------
c create S matrix GA. Do we need to use that atom-block stuff
c here, like they do in the SCF calling routine dft_main0d.F ???
if (.not. ga_create(MT_DBL,nu,nc,'g_s',1,1,g_s)) call
& errquit(pname//': ga_create failed g_s',11,GA_ERR)
call ga_zero(g_s)
c the matrices for T, V, W are created as duplicates of g_s:
if (.not.ga_duplicate(g_s, g_t,'g_t')) call
& errquit(pname//': ga_create failed g_t',12,GA_ERR)
call ga_zero(g_t)
if (.not.ga_duplicate(g_s,g_v, 'g_v')) call
& errquit(pname//': ga_create failed g_v',13,GA_ERR)
call ga_zero(g_v)
if (.not.ga_duplicate(g_s, g_w, 'g_w')) call
& errquit(pname//': ga_create failed g_w',14,GA_ERR)
call ga_zero(g_w)
c ---------------------------------------
c Analytic one-electron matrices: S, T, V
c ---------------------------------------
if (ocache) call int_1e_uncache_ga()
call int_1e_ga(ubas, ubas ,g_s, 'overlap' ,oskel)
call int_1e_ga(ubas, ubas, g_t, 'kinetic' , oskel)
call int_1e_ga(ubas, ubas, g_v, 'potential', oskel)
c symmetrize, just in case ...
call ga_symmetrize(g_s)
call ga_symmetrize(g_t)
call ga_symmetrize(g_v)
call ga_sync()
if (master.and.dbg>0) write(luout,*) pname//': T,V,S done'
c ------------------------------------------------------
c later we need a transformation from the contracted (c)
c to the uncontracted (u) basis, and vice versa for the
c density matrix sooner. the sequence comes
c from routine int_1edk_ga
c ------------------------------------------------------
c u2c transforms the basis from uncontracted to contracted
c (really, these are just the contraction coefficients,
c but this is how it is done in the DKH routines)
if (lx2c_u2c) then
if (.not. ga_create(MT_DBL,nu,nc,'u2c',1,1,g_u2c))
& call errquit(pname//': create u2c failed',51,GA_ERR)
call ga_zero(g_u2c)
if (.not. ga_create(MT_DBL,nu,nc,'smix',1,1,g_smix))
& call errquit(pname//': create mixed overlap failed',52,GA_ERR)
call ga_zero(g_smix)
if (.not. ga_create(MT_DBL,nu,nu,'temp',1,1,g_temp))
& call errquit(pname//': error creating temp',53,GA_ERR)
if (master.and.dbg>0) write(luout,*) pname//': u2c mem done'
call int_1e_ga(ubas, cbas, g_smix, 'overlap' ,oskel)
call ga_copy(g_s,g_temp)
call ga_matpow(g_temp, -1.0d0, 1d-12)
call ga_dgemm('n','n',nu,nc,nu,1.0d0,
& g_temp, g_smix, 0.0d0, g_u2c)
if (.not. ga_destroy(g_temp))
& call errquit(pname//': error destroy temp',53, GA_ERR)
if (.not. ga_destroy(g_smix))
& call errquit(pname//': error destroy smix',52, GA_ERR)
if (master.and.dbg>0) write(luout,*) pname//': u2c done'
end if ! lx2c_u2c
c c2u transforms the basis from contracted to uncontracted.
c we have to do this as a generalized left-inverse of u2c.
c Use [(A^T A)**(-1) A^T] A = 1 to define the left-inverse.
if (lx2c_c2u) then
if (.not. ga_create(MT_DBL,nc,nu,'c2u',1,1,g_c2u))
& call errquit(pname//': create c2u failed',54,GA_ERR)
call ga_zero(g_c2u)
if (.not. ga_create(MT_DBL,nc,nc,'temp',1,1,g_temp))
& call errquit(pname//': failed creating temp',55,GA_ERR)
if (master.and.dbg>0) write(luout,*) pname//': c2u mem done'
call ga_dgemm('t','n',nc,nc,nu,1.0d0,g_u2c,g_u2c,
& 0.0d0,g_temp)
call ga_matpow(g_temp, -1.0d0, 1d-12)
if (master.and.dbg>0) write(luout,*) pname//': c2u matpow done'
call ga_dgemm('n','t',nc,nu,nc,1.0d0,g_temp,g_u2c,
& 0.0d0,g_c2u)
if (.not. ga_destroy(g_temp))
& call errquit(pname//': error destroy temp',55, GA_ERR)
if (master.and.dbg>0) write(luout,*) pname//': c2u done'
c check c2u * u2c
if (.not. ga_create(MT_DBL,nc,nc,'temp',1,1,g_temp))
& call errquit(pname//': error creating temp',56,GA_ERR)
call ga_dgemm('n','n',nc,nc,nu,1.0d0,g_c2u,g_u2c,
& 0.0d0,g_temp)
if (dbg>1) then
call ga_print(g_u2c)
call ga_print(g_c2u)
call ga_print(g_temp)
end if
c check the quality of the left-inverse. the product
c c2u*u2c must be a unit matrix of dimension nc
sum = ga_trace_diag(g_temp)
if (master) write (luout,'(1x,a,1x,e15.7)')
& 'x2c c2u*u2c: ',sum
if (abs(sum-real(nc)).gt.1.0d-10) then
if (master) write (luout,*)
& 'x2c WARNING: left-inverse in basis trafo is inaccurate'
end if
if (.not. ga_destroy(g_temp))
& call errquit(pname//': error destroy temp',56, GA_ERR)
if (master.and.dbg>0) write(luout,*) pname//': c2u done'
end if ! lx2c_c2u
if (lx2c_u2c) then
c transform density matrix to the uncontracted basis:
c c2u is how the basis would transform in this case.
c the density matrix is based on the MO coefficients,
c P = C n C^T, with C being a collections of MO vectors
c in the c-basis. We do a passive transformation,
c since the MOs are to remain the same. In this case,
c the transformation of the MO vectors is contra-
c gredient to the basis and we need u2c:
c P(u) = u2c P(c) u2c^T.
if (.not. ga_create(MT_DBL,nc,nu,'temp',1,1,g_temp))
& call errquit(pname//': error creating temp',57,GA_ERR)
if (.not. ga_create(MT_DBL,nu,nu,'densu',1,1,g_densu))
& call errquit(pname//': error creating densu',58,GA_ERR)
if (master.and.dbg>0) write(luout,*) pname//': dens c2u start'
call ga_zero(g_densu)
call ga_dgemm('n','t',nc,nu,nc, 1.0d0,
& g_dens(1),g_u2c,0.0d0,g_temp)
if (master.and.dbg>0) write(luout,*) pname//': dens c2u step1'
call ga_dgemm('n','n',nu,nu,nc ,1.0d0,
& g_u2c, g_temp, 0.0d0,g_densu)
call ga_symmetrize(g_densu)
if (.not. ga_destroy(g_temp))
& call errquit(pname//': error destroy temp',57, GA_ERR)
if (master.and.dbg>0) write(luout,*) pname//': dens c2u done'
else
c if there is no need for the u2c trafo, we
c just assign the GA handle of g_dens(1) to g_densu
g_densu = g_dens(1)
end if ! lx2c_u2c
if (master.and.dbg>0) write (luout,*)
& 'g_densu, g_dens(1) = ',g_densu,g_dens(1)
c ------------------------------------
c use sum of atomic densities to check
c that the AO matrices are OK
c ------------------------------------
energy = 0d0
if (.not. ga_create(MT_DBL,nu,nu,'g_temp',1,1,g_temp)) call
& errquit(pname//': ga_create g_temp failed',59,GA_ERR)
call ga_zero(g_temp)
if (master.and.dbg>0) write(luout,*) pname//': debug S:'
call ga_dgemm('n','n',nu,nu,nu, 1.0d0,
& g_densu,g_s,0.0d0,g_temp)
sum = ga_trace_diag(g_temp)
if (master) write (luout,'(1x,a,1x,e15.7)')
& pname//': tr(P * S) : ',sum
if (master.and.dbg>0) write(luout,*) pname//': debug T:'
call ga_dgemm('n','n',nu,nu,nu, 1.0d0,
& g_densu,g_t,0.0d0,g_temp)
sum = ga_trace_diag(g_temp)
if (master) write (luout,'(1x,a,1x,e15.7)')
& pname//': tr(P * T) : ',sum
energy = energy + sum
if (master.and.dbg>0) write(luout,*) pname//': debug V:'
call ga_dgemm('n','n',nu,nu,nu, 1.0d0,
& g_densu,g_v,0.0d0,g_temp)
sum = ga_trace_diag(g_temp)
if (master) write (luout,'(1x,a,1x,e15.7)')
& pname//': tr(P * V) : ',sum
energy = energy + sum
if (.not. ga_destroy(g_temp)) call errquit
& (pname//': temp ga corrupt?',59, GA_ERR)
if (master) write (luout,'(1x,a,1x,e15.7)')
& pname//': 1e energy : ',energy
c ----------------------------------------------------------
c as this routine gets embedded in the scalar DFT framework,
c we make it work similar to ZORA. I.e., for compatibility
c we need to SUBTRACT the nonrel H from the operator.
c this is done here, upfront.
c ----------------------------------------------------------
c if lx2c_u2c is .T., we now need to allocate an array for the
c x2c operator in the uncontracted basis => g_X
c otherwise, we simply use the GA handle of g_x2c.
if (lx2c_u2c) then
if (.not.ga_duplicate(g_s, g_X, 'g_X')) call
& errquit(pname//': ga_create failed g_X',61,GA_ERR)
call ga_zero(g_X)
else
g_X = g_x2c
end if ! lx2c_u2c
call ga_dadd(1.d0, g_X, -1.0d0, g_t, g_X)
call ga_dadd(1.d0, g_X, -1.0d0, g_v, g_X)
c ------------------------------------------------------
c In the case of using a model Coulomb potential in the
c transformation (V and W integrals), we need to
c subtract the corresponding Coulomb potential from the
c one-electron operator and add it to V (and later in W)
c ------------------------------------------------------
if (modelpot_2e) then
if(.not.ga_create(mt_dbl,nu,nu,'tmpv',0,0,g_tmpv))
& call errquit(pname//': error creating ga',141, GA_ERR)
jfac(1) = 1.0d0
kfac(1) = 0.0d0
call ga_zero(g_tmpv)
if (master.and.dbg>0) write (luout,*)
& 'x2c_1e: calling fock2e for model potential'
call fock_2e(geom, ubas, 1, jfac, kfac,
& tol2e, oskel, g_densu, g_tmpv, .false.)
call ga_sync()
c subtract from the 1e Hamiltonian:
call ga_add(1.0d0, g_X, -1.0d0, g_tmpv, g_X)
c add Coul. pot. to g_v:
call ga_add(1.0d0, g_v, 1.0d0, g_tmpv, g_v)
if (.not. ga_destroy(g_tmpv)) call errquit
& (pname//': ga corrupt?',141, GA_ERR)
end if ! modelpot_2e
c g_densu not needed past this point
if (lx2c_u2c) then
if (.not. ga_destroy(g_densu))
& call errquit(pname//': error destroy densu',58, GA_ERR)
else
if (g_densu.ne.g_dens(1)) call errquit(
& pname//': g_densu inconsistency 1',g_densu, UNKNOWN_ERR)
end if
call ga_zero(g_w)
c ---------------------------------------------------------------
c nuclear potential only: calculate analytic pV dot p integrals
c and add them to the array g_w(:).
c the calls to the integral routines below have been adapted
c from the douglas-kroll code in nwchem.
c ---------------------------------------------------------------
if (analytic_pvp) then
if (master) write (luout,*)
& pname//': calculating analytic p.Vp integrals'
c add scalar part p.Vp to g_w:
call int_1e_ga(ubas, ubas ,g_w, 'pvp' ,oskel)
call ga_symmetrize(g_w)
end if ! analytic_pvp
c --------------------------------------------------
c if needed: calculate all or part of
c (sig.p)V(sig.p) integrals by numerical integration,
c similar to how the ZORA spin-orbit arrays are done.
c We have the option to skip the nuclear part, if those
c integrals are done analytically
c --------------------------------------------------
if (analytic_pvp .and. .not.modelpot_2e) then
if (master) write(luout,*)
& pname//': skipping numerical W integrals'
else
if (master) write(luout,*) pname//': calculating W numerically'
if (modelpot_2e .and. master)
& write(luout,*) ' (2e model pot. only)'
if (lx2c_u2c) then
c here we need to use the contracted basis ...
c use g_tmpv to hold the array.
c g_vxc has the right dimension already
if (dbg>0.and.master) write (luout,*) 'g_w numint c2u start'
if(.not.ga_create(mt_dbl,nc,nc,'tmpw',0,0,g_tmpv))
& call errquit(pname//': error creating tmpw',147, GA_ERR)
call ga_zero(g_tmpv)
call x2c_getw_scalar(rtdb, g_dens(1), g_tmpv, g_vxc, nexc,
& analytic_pvp, .not.modelpot_2e)
if (dbg>0.and.master) write (luout,*) 'g_w numint c2u end'
c transform to ubas: tmpv = c2u^T g_tmpv c2u, add to g_w
if(.not.ga_create(mt_dbl,nc,nu,'temp',0,0,g_temp))
& call errquit(pname//': error creating temp',148, GA_ERR)
call ga_zero(g_temp)
call ga_dgemm('n','n',nc,nu,nc, 1.0d0,
& g_tmpv,g_c2u,0.0d0,g_temp)
if (master.and.dbg>0) write(luout,*) pname//': g_w c2u step1'
call ga_dgemm('t','n',nu,nu,nc ,1.0d0,
& g_c2u, g_temp, 1.0d0,g_w)
if (master.and.dbg>0) write(luout,*) pname//': g_w c2u step2'
if (.not. ga_destroy(g_tmpv)) call errquit
& (pname//': ga corrupt?',147, GA_ERR)
if (.not. ga_destroy(g_temp)) call errquit
& (pname//': ga corrupt?',148, GA_ERR)
else
c cbas = ubas. Just continue as usual
if (dbg>0.and.master) write (luout,*) 'g_w numint start'
call x2c_getw_scalar(rtdb, g_dens(1), g_w, g_vxc, nexc,
& analytic_pvp, .not.modelpot_2e)
if (dbg>0.and.master) write (luout,*) 'g_w numint end'
end if ! lx2c_u2c
end if ! need analytic integration
c W is the pVp matrix, incl model pot if requested:
if (dbg>1) call ga_print(g_w)
c sanity checks:
if (.not.lx2c_u2c) then
if (g_X.ne.g_x2c) call errquit(
& pname//': g_X/g_x2c inconsistency 1',g_X-g_x2c, UNKNOWN_ERR)
end if
c -----------------------------------------------
c Input data are all in place for decoupling step
c -----------------------------------------------
c define some constants related to speed of light:
c (define them here directly as inverses, which we need later)
c_2c = 1.0d0/(cau*2.d0) ! 1/2c
c_2c2 = 1.0d0/(cau*cau*2.d0) ! 1/2c**2
c_4c2 = 1.0d0/(cau*cau*4.d0) ! 1/4c**2
c scale W by 1/(4c^2)
call ga_scale(g_w, c_4c2)
c assemble the 'scalar 4-component' mDE Hamiltonian
if (.not. ga_create(MT_DBL,nu2,nu2,'x2c_h',1,1,g_h))
& call errquit(pname//': ga_create failed g_h',66, GA_ERR)
call ga_zero(g_h)
call ga_copy_patch('n', ! L-L block = V
& g_v, 1, nu, 1, nu,
& g_h, 1, nu, 1, nu)
call ga_copy_patch('n', ! L-S block = T
& g_t, 1, nu, 1, nu,
& g_h, 1, nu, nu+1, nu2)
call ga_copy_patch('n', ! S-L block = T
& g_t, 1, nu, 1, nu,
& g_h, nu+1, nu2, 1, nu)
call ga_add_patch( ! S-S block = W/(4c**2)-T
& 1.0d0, g_w, 1, nu, 1, nu,
& -1.0d0, g_t, 1, nu, 1, nu,
& g_h, nu+1, nu2, nu+1, nu2)
c assemble the metric for the mDE:
if (.not. ga_create(MT_DBL,nu2,nu2,'x2c_m',1,1,g_m))
& call errquit(pname//': ga_create failed g_m',67, GA_ERR)
call ga_zero(g_m)
call ga_copy_patch('n', ! L-L block = S
& g_s, 1, nu, 1, nu,
& g_m, 1, nu, 1, nu)
c call ga_scale(g_t, c_2c2)
call ga_copy_patch('n', ! S-S block = T/(2c**2)
& g_t, 1, nu, 1, nu,
& g_m, nu+1, nu2, nu+1, nu2)
call ga_scale_patch(
& g_m, nu+1, nu2, nu+1, nu2, c_2c2)
c diagonalize the mDE Hamiltonian:
if (.not. ma_push_get(mt_dbl, nu2, 'eval', l_eval, k_eval))
& call errquit(pname//': failed allocation eval', nu2,
& MA_ERR)
if (.not.ga_duplicate(g_h, g_temp,'g_temp')) call
& errquit(pname//': ga_dup failed g_temp',68,GA_ERR)
call ga_zero(g_temp)
if (dbg>1) call ga_print(g_h)
if (dbg>1) call ga_print(g_m)
#if defined(PARALLEL_DIAG)
#ifdef SCALAPACK
c
c SCALAPACK piece needs code similar to what's
c in src/NWints/dk/int_1edk_ga.F
if (master .and. dbg>0) write (luout,*)
& 'calling ga_pdsygv replacement sequence of routines'
c call ga_pdsygv(g_h, g_m, g_temp, dbl_mb(k_eval))
info= ga_cholesky('L', g_m)
info = ga_llt_i('L', g_m, -1)
info = ga_cholesky('L', g_m, -1)
call ga_dgemm('t','n',nu2, nu2, nu2, 1.0d0,
& g_m, g_h, 0.0d0, g_temp)
call ga_dgemm('n','n', nu2, nu2, nu2, 1.0d0,
& g_temp, g_m, 0.0d0, g_h)
call ga_pdsyev(g_h, g_temp, dbl_mb(k_eval), 0)
call ga_dgemm('n','n', nu2, nu2, nu2, 1.0d0,
& g_m, g_temp, 0.0d0, g_h)
call ga_copy (g_h, g_temp)
#else
if (master .and. dbg>0) write (luout,*) 'calling ga_diag'
call ga_diag(g_h, g_m, g_temp, dbl_mb(k_eval))
#endif
#else
if (master .and. dbg>0) write (luout,*) 'calling ga_diag_seq'
call ga_diag_seq(g_h, g_m, g_temp, dbl_mb(k_eval))
#endif
c we don't need the mDE eigenvalues further
if (.not.ma_chop_stack(l_eval)) call
& errquit(pname//': ma_chop_stack failed k_eval',l_eval,MA_ERR)
if (dbg>1) call ga_print(g_temp)
c destroy arrays for mDE Hamiltonian and metric,
c allocate the same handles again for scratch space
if (.not. ga_destroy(g_h))
& call errquit(pname//': ga_destroy failed g_h',66, GA_ERR)
if (.not. ga_destroy(g_m))
& call errquit(pname//': ga_destroy failed g_m',67, GA_ERR)
if (.not.ga_duplicate(g_s, g_h, 'g_h (temp)')) call
& errquit(pname//': ga_dup failed g_h(temp)',71,GA_ERR)
if (.not.ga_duplicate(g_s, g_m, 'g_m (temp)')) call
& errquit(pname//': ga_dup failed g_m(temp)',72,GA_ERR)
c set up and perform the decoupling.
call ga_copy_patch('n',
& g_temp, 1, nu, nu+1, nu2,
& g_m, 1, nu, 1, nu) ! g_m now holds array A
call ga_copy_patch('n',
& g_temp, nu+1, nu2, nu+1, nu2,
& g_h, 1, nu, 1, nu) ! g_h now holds array B
if (dbg>1) call ga_print(g_m)
if (dbg>1) call ga_print(g_h)
c destroy temp GA and reallocate with nubas**2 elements
if (.not. ga_destroy(g_temp))
& call errquit(pname//': ga_destroy failed temp',68, GA_ERR)
if (.not.ga_duplicate(g_s, g_temp,'g_temp')) call
& errquit(pname//': ga_dup failed g_temp',73,GA_ERR)
c solve matrix equation X A = B for X, store in g_h
c that is, solve A' X' = B' (' = Transpose), then transpose.
c matrix A is square, but not symmetric.
call ga_transpose(g_m, g_temp)
call ga_copy(g_temp, g_m)
call ga_transpose(g_h, g_temp)
call ga_copy(g_temp, g_h) ! keep copy of B' in g_temp, too
#ifdef SCALAPACK
if (dbg>0.and.master) write(luout,*) 'calling ga_solve'
info = ga_solve(g_m, g_h) ! g_h now has solution X
if (dbg>0.and.master) write(luout,*) 'info from ga_solve:', info
#else
if (dbg>0.and.master) write(luout,*) 'calling ma_solve'
call ma_solve(g_m, g_h)
#endif
c test quality of the solution:
call ga_dgemm('N', 'N', nu, nu, nu, 1.d0,
& g_m, g_h, -1.d0, g_temp) ! g_temp = A'X'-B', should be zero
call ga_maxelt(g_temp, rtemp)
if (master) write(luout,*)
& 'largest element in A''X''-B'': ',rtemp
call ga_transpose(g_h, g_temp)
call ga_copy(g_temp, g_h)
if (dbg>1) call ga_print(g_h)
c assemble S + X(+)TX/(2c**2) in g_s
c assemble X(+)T + TX - X(+)TX + V in g_t
call ga_dgemm('N', 'N', nu, nu, nu, 1.d0,
& g_t, g_h, 0.d0, g_temp)
call ga_dgemm('T', 'N', nu, nu, nu, 1.d0,
& g_h, g_t, 0.d0, g_m)
call ga_copy(g_m, g_t)
call ga_add(1.0d0, g_t, 1.0d0, g_temp, g_t)
call ga_dgemm('T', 'N', nu, nu, nu, 1.d0,
& g_h, g_temp, 0.d0, g_m)
call ga_add (1.0d0, g_t, -1.0d0, g_m, g_t)
if (dbg>1) call ga_print(g_t)
call ga_add (1.0d0, g_t, 1.0d0, g_v, g_t)
call ga_copy(g_s, g_temp) ! save a copy of S in g_temp
call ga_scale(g_m, c_2c2)
call ga_add (1.0d0, g_s, 1.0d0, g_m, g_s)
if (dbg>1) then
if (master) write (luout,*) 'tilde S = S + X''TX/(2c^2):'
call ga_print(g_s)
if (master) write (luout,*) 'X''T + TX - X''TX + V:'
call ga_print(g_t)
end if
c now assemble the metric correction for the FW Hamiltonian
c g_m, g_v can be used for scratch
c g_temp still holds S
c g_h still holds X
call ga_copy (g_temp, g_v)
call ga_matpow(g_v, -FP_Half, 1d-12)
call ga_dgemm('T', 'N', nu, nu, nu, 1.d0,
& g_v, g_s, 0.d0, g_m)
call ga_dgemm('N', 'N', nu, nu, nu, 1.d0,
& g_m, g_v, 0.d0, g_s)
call ga_matpow(g_s, -FP_Half, 1d-12)
call ga_dgemm('N', 'N', nu, nu, nu, 1.d0,
& g_v, g_s, 0.d0, g_m)
call ga_matpow(g_temp, FP_Half, 1d-12)
call ga_dgemm('N', 'N', nu, nu, nu, 1.d0,
& g_m, g_temp, 0.d0, g_s) ! g_s = = UL
call ga_dgemm('N', 'N', nu, nu, nu, 1.d0,
& g_h, g_s, 0.d0, g_m) ! g_m = US g_h, g_v, g_temp now free
c construct the FW Hamiltonian in g_h.
c we then need to ADD it to g_X since it might
c contain -V(coul) already.
c also, the g_X array already has T+V subtracted
c such that it ties in with the zora code
call ga_dgemm('T', 'N', nu, nu, nu, 1.d0,
& g_s, g_t, 0.d0, g_temp)
call ga_dgemm('N', 'N', nu, nu, nu, 1.d0,
& g_temp, g_s, 0.d0, g_h)
call ga_dgemm('T', 'N', nu, nu, nu, 1.d0,
& g_m, g_w, 0.d0, g_temp)
call ga_dgemm('N', 'N', nu, nu, nu, 1.d0,
& g_temp, g_m, 1.d0, g_h)
if (dbg>1) call ga_print(g_h)
call ga_add(1.0d0, g_h, 1.0d0, g_X, g_X)
c free up some memory:
if (.not. ga_destroy(g_temp))
& call errquit(pname//': ga_destroy failed temp',73, GA_ERR)
if (.not. ga_destroy(g_h))
& call errquit(pname//': ga_destroy failed g_h',71, GA_ERR)
if (.not. ga_destroy(g_t))
& call errquit(pname//': ga_destroy failed t',12, GA_ERR)
if (.not. ga_destroy(g_v))
& call errquit(pname//': ga_destroy failed v',13, GA_ERR)
if (.not. ga_destroy(g_w))
& call errquit(pname//': ga_destroy failed w',14, GA_ERR)
c -----------------------------------------------------
c save decoupling transformation, if requested by input
c this is kept in the uncontracted basis.
c if we don't keep the transformation we also don't
c need to keep the u2c [or c2u] transformation
c -----------------------------------------------------
if (lx2c_keepU) then
if (.not. ga_create(MT_DBL,nu,nu,'x2c_UL',1,1,g_u(1)))
& call errquit(pname//': ga_create failed UL',75, GA_ERR)
if (.not. ga_create(MT_DBL,nu,nu,'x2c_US',1,1,g_u(2)))
& call errquit(pname//': ga_create failed US',76, GA_ERR)
call ga_copy(g_s, g_u(1))
call ga_copy(g_m, g_u(2))
call ga_scale(g_u(2), c_2c)
end if ! lx2c_keepU
if (.not. ga_destroy(g_s))
& call errquit(pname//': ga_destroy failed s',11, GA_ERR)
if (.not. ga_destroy(g_m))
& call errquit(pname//': ga_destroy failed m',72, GA_ERR)
if (dbg>1) call ga_print(g_X)
c sanity check:
if (.not.lx2c_u2c) then
if (g_X.ne.g_x2c) call errquit(
& pname//': g_X/g_x2c inconsistency 2',g_X-g_x2c, UNKNOWN_ERR)
end if
c ----------------------------------------------
c transform x2c operator to contracted basis and
c store in g_x2c. Otherwise, we are already done
c ----------------------------------------------
if (lx2c_u2c) then
if (dbg>0.and.master) write(luout,*)
& pname//': u2c trafo of result'
if (.not. ga_create(MT_DBL,nc,nu,'temp',1,1,g_temp))
& call errquit(pname//': ga_create failed temp',77,
& GA_ERR)
if (g_X.eq.g_x2c) call errquit(
& pname//': g_X/g_x2c inconsistency 3',g_X-g_x2c, UNKNOWN_ERR)
call ga_dgemm('t','n',nc,nu,nu,1.0d0,
& g_u2c,g_X,0.0d0,g_temp)
c note: if lx2c_u2c.eq..T. then g_x2c and g_X are not the
c same. Hence, g_x2c is over-written here
call ga_dgemm('n','n',nc,nc,nu,1.0d0,
& g_temp,g_u2c, 0.0d0,g_x2c)
if (.not. ga_destroy(g_temp))
& call errquit(pname//': ga_destroy failed temp',77,GA_ERR)
end if ! lx2c_u2c
call ga_symmetrize(g_x2c)
c de-allocate g_X, if needed, and related arrays
if (lx2c_u2c) then
if (.not. ga_destroy(g_X))
& call errquit(pname//': ga_destroy failed g_X',61,GA_ERR)
if (.not.lx2c_keepU) then
if (.not. ga_destroy(g_u2c)) call errquit
& (pname//': ga_destroy failed u2c',51,GA_ERR)
end if
end if ! lx2c_u2c
if (lx2c_c2u .and. .not.lx2c_keepU) then
if (.not. ga_destroy(g_c2u)) call errquit
& (pname//': ga_destroy failed c2u',54,GA_ERR)
end if ! lx2c_c2u
c restore model alpha and beta densities:
if (ipol.gt.1)
& call ga_dadd(1.d0, g_dens(1), -1.d0, g_dens(2), g_dens(1))
c --------
c all done
c --------
call ga_sync()
if (dbg>1) then
call ga_print(g_x2c)
end if
oskel = oskel_sav ! restore oskel value
if (dbg>0.and.master) write(luout,*) 'leaving x2c_1e_scalar'
c =================================================================
return
end

606
src/nwdft/x2c/calc_x2c_efg_scalar.F Normal file
View file

@ -0,0 +1,606 @@
subroutine calc_x2c_efg_scalar (rtdb, nexc,
& g_dscf, g_dens_at, g_u)
c -----------------------------------------------------------------
c purpose: calculate one-electron property EFG using the
c two-component exact-decoupling Hamiltonian
c For the time being we MUST use an uncontracted basis !
c Integrals are supposed to be initialized by the calling
c routine.
c Driver (c) 2012, 2022 J. Autschbach, SUNY, jochena@buffalo.edu
c Approach: The transformation matrices g_u come from calc_x2c_1e
c if certain options are set.
c
c -----------------------------------------------------------------
implicit none
#include "mafdecls.fh"
#include "rtdb.fh"
#include "apiP.fh"
#include "global.fh"
#include "tcgmsg.fh"
#include "bas.fh"
#include "util.fh"
#include "stdio.fh"
#include "errquit.fh"
#include "dftpara.fh"
#include "cdft.fh"
#include "geom.fh"
c#include "rel_consts.fh"
#include "cosmo.fh"
c subroutine arguments:
integer rtdb ! [in] runtime database GA handle
integer nexc ! [in] no. of XC terms
integer g_dscf(2) ! [in] density matrix from SCF
integer g_dens_at(2) ! [in] atom densities for model pot.
integer g_u(2) ! [in] decoupling matrices
c local variables:
integer basis
integer nbas, nbas2
integer me, type, i, j
double precision sum
character*2 symbol
character*16 element, at_tag
c GA handles used for temp storage in this routine:
integer g_P(2), g_efgint, g_temp, g_efg, g_w(6)
c dimensioning arrays
integer alo(3), ahi(3), blo(3), bhi(3), clo(3), chi(3)
c other MA handles
integer l_xyzpt, k_xyzpt, l_zanpt, k_zanpt
c bq variables (MV)
logical dobq
integer bq_ncent
integer i_cbq
integer i_qbq
double precision elpotbq
double precision ga_trace_diag
external ga_trace_diag
double precision tol2e, accval
logical master, analytic_pvp, status
integer dbg
integer ndmem, ntemp
integer iat, nat
double precision xp, yp, zp, xn, yn, zn, zan
double precision rr, rr5
double precision efgel(6), efgnuc(6), etmp(6)
double precision origin(3)
double precision pifac, pi, deg, third
parameter (third=1.0d0/3.0d0)
c parameter(froth=4.0d0/3.0d0)
c parameter(pifac=froth*3.14159265358979323846264338327950288419d0)
c =================================================================
c ---------------
c initializations
c ---------------
c options / settings
dbg=0
me=ga_nodeid()
master = me.eq.0
if (dbg>0.and.master) write(luout,*) 'entering x2c_efg'
call ga_sync()
c constants:
tol2e=10.d0**(-itol2e)
call int_acc_get(accval)
if (dbg>0.and.master) write(luout,*) 'accval = ',accval
c accval = 1d-25
c call int_acc_set(accval)
pi = acos(-1.0d0)
deg = 180.0d0/pi
pifac = (4.0d0/3.0d0)*pi
basis = ao_bas_han ! AO basis handle
c make sure the input GAs have the correct dimension
call ga_inquire(g_u(1),type,nbas,nbas)
if (nbas.ne.nbf_ao) then
call errquit('x2c_efg: g_u has wrong dimension',nbas,GA_ERR)
end if
if (master.and.dbg>0) write (luout,*) 'x2c: nbas =',nbas
nbas2 = nbas * nbas
c in light of the recent updates of the x2c ham. code, it is
c important to check that we are working with an uncontracted
c basis here.
c eventually, this will be fixed
if (.not. bas_numbf(x2c_bas,ntemp))
& call errquit('x2c_efg: ntemp error',ntemp, BASIS_ERR)
if (ntemp.ne.nbas)
& call errquit(
& 'x2c_efg: cannot currently use a contracted basis',
& ntemp, BASIS_ERR)
c -------------
c Calculate EFG
c -------------
c initialize integrals that were not already initialized
c call hnd_giao_init(basis,1)
c ----- define points for calculation -----
c 1. nuclei: for now, we take ALL nuclei
c
status=geom_ncent(geom,nat)
c
if (.not. ma_push_get(mt_dbl,3*nat,'xyz pnt',l_xyzpt,k_xyzpt))
& call errquit('x2c_efg: ma failed',911,MA_ERR)
if (.not. ma_push_get(mt_dbl,nat,'zan pnt',l_zanpt,k_zanpt))
& call errquit('x2c_efg: ma failed',911,MA_ERR)
c
do iat=1,nat
status=geom_cent_get(geom,iat,at_tag,dbl_mb(k_xyzpt+3*(iat-1)),
& dbl_mb(k_zanpt+iat-1))
end do
alo(1) = nbas
alo(2) = -1
alo(3) = -1
ahi(1) = nbas
ahi(2) = nbas
ahi(3) = 6*nat
if (.not.nga_create(MT_DBL,3,ahi,'g_efgint',alo,g_efgint)) call
& errquit('x2c_efg: nga_create failed g_efgint',11,GA_ERR)
call ga_zero(g_efgint)
c calculate field gradient integrals (same as FC+SD in NMR ...)
call int_giao_1ega(basis,basis,g_efgint,'sd+fc',
& dbl_mb(k_xyzpt),nat, oskel)
c comment: The electronic integrals contain the isotropic
c `contact' piece. However, when constructing the traceless
c EFG tensor, those contributions disappear. We might as well
c take them out explicitly
c calculate FC integrals separately
ahi(3) = nat
if (.not.nga_create(MT_DBL,3,ahi,'temp',alo,g_temp)) call
& errquit('add_x2c: nga_create failed g_temp',31,GA_ERR)
call ga_zero(g_temp)
call int_giao_1ega(basis,basis,g_temp,'fc',
& dbl_mb(k_xyzpt),nat,oskel)
c
c Take out FC part from SD+FC (i.e. from xx, yy and zz components)
c
alo(1) = 1
ahi(1) = nbas
alo(2) = 1
ahi(2) = nbas
blo(1) = 1
bhi(1) = nbas
blo(2) = 1
bhi(2) = nbas
do iat = 1, nat
alo(3) = 1+(iat-1)*6
ahi(3) = 1+(iat-1)*6
blo(3) = iat
bhi(3) = iat
do i = 1, 3 ! runs over xx, yy, zz
call nga_add_patch(1.0d0,g_efgint,alo,ahi,
& pifac,g_temp,blo,bhi,
& g_efgint,alo,ahi)
alo(3) = alo(3)+1
ahi(3) = ahi(3)+1
enddo
enddo
if (.not. ga_destroy(g_temp)) call errquit
& ('x2c_efg: temp ga corrupt?',31, GA_ERR)
c the field gradient integrals carry a factor of 3
c that we don't want:
call ga_scale(g_efgint, third)
if (dbg>1 .and.master) write (luout,*) 'x2c_efg before iat loop'
c ---------------------------------------
c calculate EFGs in a loop over the atoms
c ---------------------------------------
do iat = 1,nat
xp = dbl_mb(k_xyzpt +3*(iat-1))
yp = dbl_mb(k_xyzpt+1+3*(iat-1))
zp = dbl_mb(k_xyzpt+2+3*(iat-1))
if (dbg>1 .and.master) write (luout,*) 'x2c_efg: iat=',iat
c --------------------
c nuclear contribution
c --------------------
do i = 1,6
efgnuc(i) = 0d0
end do
do i = 1,nat
xn = dbl_mb(k_xyzpt +3*(i-1)) - xp
yn = dbl_mb(k_xyzpt+1+3*(i-1)) - yp
zn = dbl_mb(k_xyzpt+2+3*(i-1)) - zp
zan = dbl_mb(k_zanpt+i-1)
rr = sqrt(xn*xn + yn*yn + zn*zn)
if (rr.ge.1.0d-3) then
c (skip if it is the same nucleus or nuclei too close)
rr5=rr*rr*rr*rr*rr
efgnuc(1) = efgnuc(1) - zan*xn*xn/rr5
efgnuc(2) = efgnuc(2) - zan*yn*yn/rr5
efgnuc(3) = efgnuc(3) - zan*zn*zn/rr5
efgnuc(4) = efgnuc(4) - zan*xn*yn/rr5
efgnuc(5) = efgnuc(5) - zan*xn*zn/rr5
efgnuc(6) = efgnuc(6) - zan*yn*zn/rr5
end if
end do ! i = loop over all other nuclei
c create traceless version of nuclear contribution to EFG:
etmp(1) = 2.0d0*efgnuc(1) - efgnuc(2) - efgnuc(3)
etmp(2) = 2.0d0*efgnuc(2) - efgnuc(1) - efgnuc(3)
etmp(3) = 2.0d0*efgnuc(3) - efgnuc(1) - efgnuc(2)
etmp(4) = 3.0d0*efgnuc(4)
etmp(5) = 3.0d0*efgnuc(5)
etmp(6) = 3.0d0*efgnuc(6)
do i = 1,6
efgnuc(i) = etmp(i)
end do
c -----------------------
c electronic contribution
c -----------------------
c First, we do the case without picture change corrections
c copy electronic EFG integrals into temp array:
c we also rearrange them such that we get the traceless
c versions of the integrals right away
alo(1) = nbas
alo(2) = -1
alo(3) = -1
ahi(1) = nbas
ahi(2) = nbas
ahi(3) = 6
if (.not.nga_create(MT_DBL,3,ahi,'efg',alo,g_efg)) call
& errquit('add_x2c: nga_create failed g_efg',32,GA_ERR)
call ga_zero(g_efg)
ahi(3) = 1
if (.not.nga_create(MT_DBL,2,ahi,'temp',alo,g_temp)) call
& errquit('add_x2c: nga_create failed g_temp',31,GA_ERR)
call ga_zero(g_temp)
alo(1) = 1
ahi(1) = nbas
alo(2) = 1
ahi(2) = nbas
blo(1) = 1
bhi(1) = nbas
blo(2) = 1
bhi(2) = nbas
clo(1) = 1
chi(1) = nbas
clo(2) = 1
chi(2) = nbas
clo(3) = 1
chi(3) = 1
c create traceless version of EFG matrices in g_efg
c for the selected reference point
c step 1: accumulate 3xx, 3yy, 3zz,
c and accumulate xx+yy+zz in g_temp
do i = 1,3
alo(3) = (iat-1)*6 + i
ahi(3) = alo(3)
blo(3) = i
bhi(3) = blo(3)
call nga_add_patch(3.0d0,g_efgint,alo,ahi,
& 1d0,g_efg,blo,bhi,
& g_efg,blo,bhi)
call nga_add_patch(1.0d0,g_efgint,alo,ahi,
& 1d0,g_temp,clo,chi,
& g_temp,clo,chi)
end do
c step 2: subtract xx+yy+zz from diagonal
do i = 1,3
blo(3) = i
bhi(3) = blo(3)
call nga_add_patch(1.0d0,g_efg,blo,bhi,
& -1d0,g_temp,clo,chi,
& g_efg,blo,bhi)
end do
c indices 1, 2, 3 of g_efg are now the components
c 3xx-rr, 3yy-rr, 3zz-rr
c step 3: assign 3xy, 3xz, 3yz to the remaining elements 4, 5, 6
do i = 4,6
alo(3) = (iat-1)*6 + i
ahi(3) = alo(3)
blo(3) = i
bhi(3) = blo(3)
call nga_add_patch(3.0d0,g_efgint,alo,ahi,
& 1d0,g_efg,blo,bhi,
& g_efg,blo,bhi)
end do
c calculate electronic EFG components:
alo(3) = 1 ! for density matrix
ahi(3) = 1
do i = 1,6 ! loop over EFG components
if (dbg>1 .and.master) write (luout,*) 'x2c_efg: i=',i
efgel(i) = 0d0
c if we use array g_efg, use this:
blo(3) = i ! select this section from array g_efg
bhi(3) = blo(3)
c if we use g_efgint, use the boundaries below:
c$$$ blo(3) = i + (iat-1)*6
c$$$ bhi(3) = blo(3)
c tr[P*EFG]
alo(1) = 1
ahi(1) = nbas
alo(2) = 1
ahi(2) = nbas
call ga_zero(g_temp)
call nga_matmul_patch('n','n',1d0,0d0,
& g_dscf(1),alo,ahi,
& g_efg, blo,bhi,
& g_temp,clo,chi)
sum = ga_trace_diag(g_temp)
efgel(i) = efgel(i) + sum
end do ! i = terms in EFG
c form the traceless EFG tensor elements
c 3x**2 - r**2, etc, for diagonal elements
c$$$ etmp(1) = 2.0d0*efgel(1) - efgel(2) - efgel(3)
c$$$ etmp(2) = 2.0d0*efgel(2) - efgel(1) - efgel(3)
c$$$ etmp(3) = 2.0d0*efgel(3) - efgel(1) - efgel(2)
c$$$ etmp(4) = 3.0d0*efgel(4)
c$$$ etmp(5) = 3.0d0*efgel(5)
c$$$ etmp(6) = 3.0d0*efgel(6)
c$$$
c$$$ do i = 1,6
c$$$ efgel(i) = etmp(i)
c$$$ end do
c print results:
if (master) then
write (luout,*)
write (luout,*) 'X2C EFG: Before picture-change corr.'
write (luout,'(1x,a,i4,6e15.7)') 'x2cefg: inuc, efgnuc =',
& iat, (efgnuc(i), i=1,6)
write (luout,'(1x,a,i4,6e15.7)') 'x2cefg: inuc, efgel =',
& iat, (efgel(i), i = 1,6)
write (luout,'(1x,a,i4,6e15.7)') 'x2cefg: inuc, efgtot =',
& iat, (efgel(i)+efgnuc(i), i = 1,6)
write (luout,*)
end if
c ---------------------------------------------------
c calculate picture-change corrected EFG operator
c for the selected reference point (nuclear position)
c ---------------------------------------------------
c allocate arrays to hold the AO matrices of the EFG
c operator in the small-component basis
c (we'll calculate all 6 matrices in one go)
alo(1) = nbas
alo(2) = -1
alo(3) = -1
ahi(1) = nbas
ahi(2) = nbas
ahi(3) = 1
do j = 1,6
if (.not.nga_create(MT_DBL,2,ahi,'g_w',alo,g_w(j)))
& call errquit(
& 'x2cefg: nga_create failed g_w',1+10*j,GA_ERR)
call ga_zero(g_w(j))
end do
origin(1) = xp
origin(2) = yp
origin(3) = zp
if (dbg>0.and.master) write (luout,*) 'calling getw'
call x2c_getw_efg_scalar(
& rtdb, g_dens_at, origin,
& g_w,
& nexc)
if (dbg>0.and.master) write (luout,*) 'return from getw'
c --------------------
c allocate operator GA
c --------------------
if(.not.ga_create(mt_dbl, nbas, nbas,'g_P',0,0,
& g_P))
& call errquit(
& 'x2c_efg: error creating g_P',666, GA_ERR)
c
c ------------------------------------------------------------
c in a loop over the EFG components, transform the
c operator to two-component form, and calculate the electronic
c EFG component again
c ------------------------------------------------------------
do i = 1,6 ! loop over EFG components
call ga_zero(g_P)
c array bounds needed to extract EFG component from GA:
alo(1) = 1
ahi(1) = nbas
alo(2) = 1
ahi(2) = nbas
alo(3) = i
ahi(3) = i
c g_temp always uses the same dimensions clo, chi
c we use the same ones for g_u, g_P
clo(1) = 1
chi(1) = nbas
clo(2) = 1
chi(2) = nbas
c transform EFG operator
c UL' EFG UL, store in g_P
call nga_matmul_patch('t','n',1d0,0d0,
& g_u(1),clo,chi,
& g_efg,alo,ahi,
& g_temp,clo,chi)
call nga_matmul_patch('n','n',1d0,0d0,
& g_temp,clo,chi,
& g_u(1),clo,chi,
& g_P,clo,chi)
c US' pEFGp US, add to g_P
call nga_matmul_patch('t','n',1d0,0d0,
& g_u(2),clo,chi,
& g_w(i),clo,chi,
& g_temp,clo,chi)
call nga_matmul_patch('n','n',1d0,1d0,
& g_temp,clo,chi,
& g_u(2),clo,chi,
& g_P,clo,chi)
c call ga_print(g_P)
c -----------------------------------------
c calculate electronic EFG again, this time
c using the transformed operator
c -----------------------------------------
efgel(i) = 0.0d0
alo(1) = 1
ahi(1) = nbas
alo(2) = 1
ahi(2) = nbas
blo(1) = 1
bhi(1) = nbas
blo(2) = 1
bhi(2) = nbas
c tr[P*EFG]
call ga_zero(g_temp)
call nga_matmul_patch('n','n',1d0,0d0,
& g_dscf(1),alo,ahi,
& g_P,blo,bhi,
& g_temp,clo,chi)
sum = ga_trace_diag(g_temp)
efgel(i) = efgel(i) + sum
end do ! i = loop over EFG components
c print results:
if (master) then
write (luout,*)
write (luout,*) 'X2C EFG: After picture-change corr.'
write (luout,'(1x,a,i4,6e15.7)') 'x2cefg: inuc, efgnuc =',
& iat, (efgnuc(i), i=1,6)
write (luout,'(1x,a,i4,6e15.7)') 'x2cefg: inuc, efgel =',
& iat, (efgel(i), i = 1,6)
write (luout,'(1x,a,i4,6e15.7)') 'x2cefg: inuc, efgtot =',
& iat, (efgel(i)+efgnuc(i), i = 1,6)
write (luout,*)
end if
c --------------------
c clean up some GA memory
c --------------------
if (.not. ga_destroy(g_efg)) call errquit
& ('x2cefg: ga corrupt?',32, GA_ERR)
if (.not. ga_destroy(g_temp)) call errquit
& ('x2cefg: ga corrupt?',31, GA_ERR)
if (.not. ga_destroy(g_P)) call errquit
& ('x2cefg: ga P corrupt?',666, GA_ERR)
do j = 1,6
if (.not. ga_destroy(g_w(j))) call errquit
& ('x2cefg: failed destroying g_w',1+10*j, GA_ERR)
end do
c --------------------------------------------------------
c end loop over refernce points at which EFG is calculated
c (typically the list of nuclei)
c --------------------------------------------------------
end do ! iat, loop over nuclei
if (.not. ga_destroy(g_efgint)) call errquit
& ('x2cefg: ga corrupt?',11, GA_ERR)
c --------
c all done
c --------
call ga_sync()
if (dbg>0.and.master) write(luout,*) 'leaving x2c_efg'
c =================================================================
return
end

410
src/nwdft/x2c/calc_x2c_moment_scalar.F Normal file
View file

@ -0,0 +1,410 @@
subroutine calc_x2c_moment_scalar (rtdb, nexc,
& g_dscf, g_dens_at, g_u, type, ncomp)
c -----------------------------------------------------------------
c purpose: calculate some multipole moments using the
c two-component exact-decoupling Hamiltonian
c For the time being we MUST use an uncontracted basis !
c Integrals are supposed to be initialized by the calling
c routine.
c Driver (c) 2012,2022 J. Autschbach, SUNY, jochena@buffalo.edu
c Approach: The transformation matrices g_u come from calc_x2c_1e
c if certain options are set.
c
c -----------------------------------------------------------------
implicit none
c
#include "mafdecls.fh"
#include "rtdb.fh"
#include "apiP.fh"
#include "global.fh"
#include "tcgmsg.fh"
#include "bas.fh"
#include "util.fh"
#include "stdio.fh"
#include "errquit.fh"
#include "dftpara.fh"
#include "cdft.fh"
#include "geom.fh"
c#include "rel_consts.fh"
#include "cosmo.fh"
#include "msgids.fh"
#include "dra.fh"
#include "inp.fh"
c
c subroutine arguments:
integer rtdb ! [in] runtime database GA handle
integer nexc ! [in] no. of XC terms
integer g_dscf(2) ! [in] density matrix from SCF
integer g_dens_at(2) ! [in] atom densities for model pot.
integer g_u(2) ! [in] decoupling matrices
integer ncomp ! [in] number of components for moment
character*(*) type ! type of moment
c local variables:
integer basis
integer nbas, nbas2
integer me, idum, i, j
double precision sum
character*2 symbol
character*16 element, at_tag
c GA handles used for temp storage in this routine:
integer g_P, g_int, g_temp, g_moment, g_w(ncomp)
c dimensioning arrays
integer alo(3), ahi(3), blo(3), bhi(3), clo(3), chi(3)
double precision ga_trace_diag
external ga_trace_diag
double precision tol2e, accval
logical master, analytic_pvp, status
integer dbg
integer ndmem, ntemp
double precision moment(ncomp)
double precision origin(3)
c
integer i_comp
character*3 ch_comp
logical dosavepcdip
character*(nw_max_path_len) fn_pcdip
c
logical dmat_to_file
external dmat_to_file
c
c =================================================================
c ---------------
c initializations
c ---------------
c options / settings
dbg=1
me=ga_nodeid()
master = me.eq.0
if (dbg>0.and.master) write(luout,*) 'entering x2c_moment'
call ga_sync()
c
dosavepcdip = .false.
if(.not.(rtdb_get(rtdb,'x2c:savepcdip',mt_log,1,dosavepcdip)))
& dosavepcdip = .false.
c sanity check: type of integral, and dimensioning
if (type.eq.'dipole') then
if (ncomp.ne.3) call errquit
& ('x2c_moment: ncomp has wrong value',666,INPUT_ERR)
if(master) write(luout,*) 'x2c: dipole moment calculation'
else
if (master) write(luout,*)
& 'x2c_moment: integral type not supported. Aborting'
return
end if
tol2e=10.d0**(-itol2e)
call int_acc_get(accval)
if (dbg>0.and.master) write(luout,*) 'accval = ',accval
c accval = 1d-25
c call int_acc_set(accval)
basis = ao_bas_han ! AO basis handle
c make sure the input GAs have the correct dimension
call ga_inquire(g_u(1),idum,nbas,nbas)
if (nbas.ne.nbf_ao) then
call errquit('x2c_moment: g_u has wrong dimension',nbas,GA_ERR)
end if
if (master.and.dbg>0) write (luout,*) 'nbas =',nbas
nbas2 = nbas * nbas
c in light of the recent updates of the x2c ham. code, it is
c important to check that we are working with an uncontracted
c basis here.
c eventually, this will be fixed
if (.not. bas_numbf(x2c_bas,ntemp))
& call errquit('x2c_moment: ntemp error',ntemp, BASIS_ERR)
if (ntemp.ne.nbas)
& call errquit(
& 'x2c_moment: cannot currently use a contracted basis',
& ntemp, BASIS_ERR)
c coordinate origin
origin(1) = 0d0
origin(2) = 0d0
origin(3) = 0d0
c --------------------------
c Calculate electric moments
c --------------------------
c initialize integrals that were not already initialized
c call hnd_giao_init(basis,1)
c allocate integral array and calculate AO matrix.
c the integral type in variable 'type' must match
c the available types and no. of components in
c routine int_mpole_1ega
alo(1) = nbas
alo(2) = -1
alo(3) = -1
ahi(1) = nbas
ahi(2) = nbas
ahi(3) = ncomp
if (.not.nga_create(MT_DBL,3,ahi,'g_int',alo,g_int)) call
& errquit('x2c_moment: nga_create failed g_int',11,GA_ERR)
call ga_zero(g_int)
call int_mpole_1ega(basis,basis,g_int,type,origin,
& oskel)
if(.not.ga_create(mt_dbl, nbas, nbas,'g_temp',0,0,g_temp))
& call errquit(
& 'x2c_moment: error creating g_temp',12, GA_ERR)
c First, we do the case without picture change corrections
alo(1) = 1
ahi(1) = nbas
alo(2) = 1
ahi(2) = nbas
blo(1) = 1
bhi(1) = nbas
blo(2) = 1
bhi(2) = nbas
clo(1) = 1
chi(1) = nbas
clo(2) = 1
chi(2) = nbas
clo(3) = 1
chi(3) = 1
c calculate electronic moment components:
alo(3) = 1 ! for density matrix
ahi(3) = 1
do i = 1,ncomp ! loop over moment components
moment(i) = 0.0d0
c
blo(3) = i ! select this section from array g_int
bhi(3) = blo(3)
c tr[P*Moment]
moment(i) = 0.0d0
call ga_zero(g_temp)
call nga_matmul_patch('n','n',1d0,0d0,
& g_dscf(1),alo,ahi,
& g_int, blo,bhi,
& g_temp,clo,chi)
sum = ga_trace_diag(g_temp)
moment(i) = moment(i) + sum
end do ! i = terms in moment
c print results:
if (master) then
write (luout,*)
write (luout,*) 'X2C Moment: Before picture-change corr.'
write (luout,'(1x,a,50e15.7)') 'moment = ',
& (moment(i), i=1,ncomp)
write (luout,*)
end if
c allocate arrays to hold the AO matrices of the EFG
c operator in the small-component basis
c (we'll calculate all matrices in one go)
alo(1) = nbas
alo(2) = -1
alo(3) = -1
ahi(1) = nbas
ahi(2) = nbas
ahi(3) = 1
do j = 1,ncomp
if (.not.nga_create(MT_DBL,2,ahi,'g_w',alo,g_w(j)))
& call errquit(
& 'x2cefg: nga_create failed g_w',1+10*j,GA_ERR)
call ga_zero(g_w(j))
end do
if (dbg>0.and.master) write (luout,*) 'calling getw'
call x2c_getw_moment_scalar(
& rtdb, g_dens_at, origin,
& g_w,
& nexc, type, ncomp)
if (dbg>0.and.master) write (luout,*) 'return from getw'
c --------------------
c allocate operator GA
c --------------------
if(.not.ga_create(mt_dbl, nbas, nbas,'g_P',0,0,
& g_P))
& call errquit(
& 'x2c_efg: error creating g_P',666, GA_ERR)
c ------------------------------------------------------------
c in a loop over the moment components, transform the
c operator to two-component form, and recalculate moment
c ------------------------------------------------------------
do i = 1,ncomp ! loop over components
call ga_zero(g_P)
c array bounds needed to extract EFG component from GA:
alo(1) = 1
ahi(1) = nbas
alo(2) = 1
ahi(2) = nbas
alo(3) = i
ahi(3) = i
c g_temp always uses the same dimensions clo, chi
c we use the same ones for g_u, g_P
clo(1) = 1
chi(1) = nbas
clo(2) = 1
chi(2) = nbas
c transform moment operator
c UL' MomenT UL, store in g_P
call nga_matmul_patch('t','n',1d0,0d0,
& g_u(1),clo,chi,
& g_int,alo,ahi,
& g_temp,clo,chi)
call nga_matmul_patch('n','n',1d0,0d0,
& g_temp,clo,chi,
& g_u(1),clo,chi,
& g_P,clo,chi)
c US' pMomenTp US, add to g_P
call nga_matmul_patch('t','n',1d0,0d0,
& g_u(2),clo,chi,
& g_w(i),clo,chi,
& g_temp,clo,chi)
call nga_matmul_patch('n','n',1d0,1d0,
& g_temp,clo,chi,
& g_u(2),clo,chi,
& g_P,clo,chi)
c call ga_print(g_P)
c ---------------------------------
c calculate moment again, this time
c using the transformed operator
c ---------------------------------
moment(i) = 0.0d0
alo(1) = 1
ahi(1) = nbas
alo(2) = 1
ahi(2) = nbas
blo(1) = 1
bhi(1) = nbas
blo(2) = 1
bhi(2) = nbas
c tr[P*Moment]
call ga_zero(g_temp)
call nga_matmul_patch('n','n',1d0,0d0,
& g_dscf(1),alo,ahi,
& g_P,blo,bhi,
& g_temp,clo,chi)
sum = ga_trace_diag(g_temp)
moment(i) = moment(i) + sum
c The GA g_P now has the picture-change transformed
c electric moment AO matrix for component no. i
c (as returned from int_mpole_1ega).
cng check if we need to save the g_P matrices for use in TDDFT
c
if (dosavepcdip) call write_x2c_pc_dip(i,g_P) ! save components
c
end do ! i = loop over moment components
c print results:
if (master) then
write (luout,*)
write (luout,*) 'X2C Moment: After picture-change corr.'
write (luout,'(1x,a,50e15.7)') 'moment = ',
& (moment(i), i=1,ncomp)
write (luout,*)
end if
c --------------------
c clean up some GA memory
c --------------------
if (.not. ga_destroy(g_P)) call errquit
& ('x2c_moment: ga P corrupt?',666, GA_ERR)
do j = 1,ncomp
if (.not. ga_destroy(g_w(j))) call errquit
& ('x2c_moment: failed destroying g_w',1+10*j, GA_ERR)
end do
if (.not. ga_destroy(g_temp)) call errquit
& ('x2c_moment: ga corrupt?',12, GA_ERR)
if (.not. ga_destroy(g_int)) call errquit ('x2cefg: ga corrupt?'
& ,11, GA_ERR)
c --------
c all done
c --------
call ga_sync()
if (dbg>0.and.master) write(luout,*) 'leaving x2c_moment'
c =================================================================
return
end

189
src/nwdft/x2c/calc_x2c_w_scalar.F Normal file
View file

@ -0,0 +1,189 @@
c
c == calculate scalar x2c contribution ==
subroutine calc_x2c_w_scalar (
& ao_bas_han, ! in: AO basis handle
& geom, ! in: geometry handle
& ipol, ! in: nr. of spin components in g_dens
& g_dens, ! in: atomic superposition density
& skipnuc, ! in: whether to exclude nuclear potential
& skipcoul, ! in: whether to skip atomic Hartree pot.
& delchi_ao, ! in: derivatives of basis fcts.
& qxyz, ! in: grid points
& qwght, ! in: weight coefficients
& nbf, ! in: nr. basis functions
& npts, ! in: nr. grid points
& natoms, ! in: nr. atoms
& x2c0) ! out: W0
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "stdio.fh"
cinclude "x2c.fh"
#include "global.fh"
#include "bas.fh"
c
integer nbf,npts,ao_bas_han,natoms,geom
logical skipnuc, skipcoul
integer g_dens(2),ipol
double precision qwght(npts)
double precision pot(npts)
double precision qxyz(3,npts)
double precision delchi_ao(npts,3,nbf)
double precision x2c0(nbf,nbf)
integer i,j,k
double precision amat_coul(npts,ipol)
double precision amat_nucl(npts)
integer ipt,closegridpts(npts)
double precision clight_au2,tol
double precision amat_tot
double precision denom
double precision fac1_arr(npts)
double precision ac_sf
external get_ints_x2c_scalar
integer dbg
c =================================================================
dbg = 0
c debug:
if (ga_nodeid().eq.0 .and. dbg>0) then
write (luout,*) 'hello from calc_x2c_w_scalar'
write (luout,*) 'skipnuc, skipcoul = ',skipnuc, skipcoul
end if
c == preliminaries ==
do ipt = 1,npts
amat_coul(ipt,1) = 0.d0
if (ipol.gt.1) amat_coul(ipt,2) = 0.d0
amat_nucl(ipt) = 0.d0
closegridpts(ipt) = 0
end do
tol = 1d-8
c calculate sum of atomic Hartree potentials, if
c the skipcoul option is not set:
if (.not.skipcoul) then
do i=1,ipol
call potential_list(ao_bas_han, g_dens(i), npts, qxyz,
& amat_coul(1,i), tol)
enddo
if (ipol.gt.1) then ! enter-ONLY-if ipol>1
do k = 1,npts
amat_coul(k,1)=amat_coul(k,1)+amat_coul(k,2)
enddo
endif
c write (luout,*) 'x2chscale =', x2chscale
c scale Hartree potential globally if requested by input
c$$$ if (scale_x2c_hartree) then
c$$$ do k = 1,npts
c$$$ amat_coul(k,1)=amat_coul(k,1) * x2chscale
c$$$ end do
c$$$ end if
else
if (ga_nodeid().eq.0 .and. dbg>0)
& write (luout,*) 'calc_x2c_w_sc: skipped Coul. pot'
end if ! .not.skipcoul
c calculate the total point nuclear potential on the grid
c if the option skipnuc is not set
if (.not.skipnuc) then
call gridNuclearPotentialPoint(geom,natoms,npts,qxyz,qwght,
& closegridpts,amat_nucl)
do k = 1,npts
if (k.eq.closegridpts(k)) qwght(k) = 0.d0
end do
else
if (ga_nodeid().eq.0 .and. dbg>0)
& write (luout,*) 'calc_x2c_w_sc: skipped nuc. pot'
end if ! .not.skipnuc
c === define fac1_arr
do k = 1,npts
c == assemble hartree and nuclear contributions ==
c == and pre-multiply with grid point weight factor ==
amat_tot =-amat_coul(k,1)+amat_nucl(k)
fac1_arr(k)=amat_tot * qwght(k)
end do
c scale Hartree potential if requested by input
c (caution, it may have been scaled already once, see above)
c$$$ if (scale_x2c_so_hartree) then
c$$$ write (luout,*) 'x2csoscale =', x2csoscale
c$$$ do k = 1,npts
c$$$ amat_coul(k,1)=amat_coul(k,1) * x2csoscale
c$$$ end do
c$$$ end if
c == assemble x2c correction ==
c ----- main diagonal --- START
do i = 1, nbf
j=i
call get_ints_x2c_scalar(nbf,npts,delchi_ao,i,j,
& fac1_arr,
& ac_sf) ! out
x2c0(i,j) = x2c0(i,j) + ac_sf
enddo ! end-loop-i
c ----- main diagonal --- END
c ----- off diagonal --- START
do i = 1, nbf
do j = i+1, nbf
call get_ints_x2c_scalar(nbf,npts,delchi_ao,i,j,
& fac1_arr,
& ac_sf) ! out
x2c0(i,j) = x2c0(i,j) + 2.0d0*ac_sf
enddo ! end-loop-j
enddo ! end-loop-i
c ----- off diagonal --- END
if (ga_nodeid().eq.0 .and. dbg>0) then
write (luout,*) 'leaving calc_x2c_w_scalar'
end if
return
end
subroutine get_ints_x2c_scalar(
& nbf, ! in: # basis fncts.
& npts, ! in: grid points
& delchi_ao, ! in: deriv. of basis fcts.
& i,j, ! in: (i,j) indices for delchi_ao
& fac1_arr, ! in
& ac_sf) ! out
c
implicit none
#include "errquit.fh"
#include "stdio.fh"
#include "global.fh"
integer nbf,npts,i,j,k
double precision delchi_ao(npts,3,nbf)
double precision fac1_arr(npts)
double precision ac_sf
double precision prod0
ac_sf = 0.0d0
do k = 1, npts
prod0 = delchi_ao(k,1,i)*delchi_ao(k,1,j)
& +delchi_ao(k,2,i)*delchi_ao(k,2,j)
& +delchi_ao(k,3,i)*delchi_ao(k,3,j)
ac_sf = ac_sf + fac1_arr(k)*prod0
end do ! k
return
end

152
src/nwdft/x2c/calc_x2c_wefg_scalar.F Normal file
View file

@ -0,0 +1,152 @@
c
c calculate EFG operator matrix in small component basis, for x2c
subroutine calc_x2c_wefg_scalar(
& ao_bas_han, ! in: AO basis handle
& geom, ! in: geometry handle
& ipol, ! in: nr. of polarizations
& g_dens, ! in: superposit. atomic
& delchi_ao, ! in: deriv. basis fcts.
& qxyz, ! in: grid points
& qwght, ! in: weighting coeffs.
& nbf, ! in: nr. basis functions
& npts, ! in: nr. grid points
& natoms, ! in: nr. atoms
& origin, ! in: EFG center
& x2c0) ! out
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "stdio.fh"
cinclude "x2c.fh"
#include "global.fh"
#include "bas.fh"
c
integer nbf,npts,ao_bas_han,natoms,geom
integer g_dens(2),ipol
double precision qwght(npts)
double precision pot(npts)
double precision qxyz(3,npts)
double precision delchi_ao(npts,3,nbf)
double precision origin(3)
double precision x2c0(nbf,nbf,6)
integer i,j,k, iefg
integer closegridpts(npts)
double precision clight_au2
double precision amat_tot
double precision denom
double precision fac1_arr(npts,6)
double precision ac_sf(6)
external get_ints_x2c_so_efg
c == preliminaries ==
do k = 1,npts
closegridpts(k) = 0
end do
do iefg = 1,6
do k = 1,npts
fac1_arr(k,iefg) = 0.0d0
end do
end do
c == calculate the EFG operator on the grid ==
call gridQpqPotential_all(origin,
& npts,qxyz,
& fac1_arr, ! out: EFG potential
& closegridpts)
c black out grid points too close to a nucleus;
c the numerical integration may get noisy otherwise:
do k = 1,npts
if (k.eq.closegridpts(k)) qwght(k) = 0.d0
end do
c === premultiply fac1_arr with integration weights
do iefg = 1,6
do k = 1,npts
fac1_arr(k,iefg)= fac1_arr(k,iefg) * qwght(k)
end do
end do
c == assemble integrals ==
c ----- main diagonal --- START
do i = 1, nbf
j=i
call get_ints_x2c_sc_efg(
& nbf,npts,delchi_ao,i,j,
& fac1_arr,
& ac_sf) ! out
do iefg = 1,6
x2c0(i,j,iefg) = x2c0(i,j,iefg) + ac_sf(iefg)
end do
enddo ! end-loop-i
c ----- main diagonal --- END
c ----- off diagonal --- START
do i = 1, nbf
do j = i+1, nbf
call get_ints_x2c_sc_efg(
& nbf,npts,delchi_ao,i,j,
& fac1_arr,
& ac_sf) ! out
do iefg = 1,6
x2c0(i,j,iefg) = x2c0(i,j,iefg) + 2.0d0*ac_sf(iefg)
end do
enddo ! end-loop-j
enddo ! end-loop-i
c ----- off diagonal --- END
c all done
return
end
subroutine get_ints_x2c_sc_efg(
& nbf, ! in: # basis fncts.
& npts, ! in: grid points
& delchi_ao, ! in: deriv. of basis fcts.
& i,j, ! in: (i,j) indices for delchi_ao
& fac1_arr, ! operator on grid points * weights
& ac_sf) ! out
implicit none
#include "errquit.fh"
#include "stdio.fh"
#include "global.fh"
integer nbf,npts,i,j,k, iefg
double precision delchi_ao(npts,3,nbf)
double precision fac1_arr(npts,6)
double precision
& ac_sf(6)
double precision prod0
do iefg = 1,6
ac_sf(iefg) = 0.0d0
end do
do k = 1, npts
prod0 = delchi_ao(k,1,i)*delchi_ao(k,1,j)
& +delchi_ao(k,2,i)*delchi_ao(k,2,j)
& +delchi_ao(k,3,i)*delchi_ao(k,3,j)
do iefg = 1,6
ac_sf(iefg) = ac_sf(iefg) + fac1_arr(k,iefg)*prod0
end do ! iefg
end do ! k
return
end

157
src/nwdft/x2c/calc_x2c_wmoment_scalar.F Normal file
View file

@ -0,0 +1,157 @@
c
c calculate moment operator matrix in small component basis, for x2c
subroutine calc_x2c_wmoment_scalar(
& ao_bas_han, ! in: AO basis handle
& geom, ! in: geometry handle
& ipol, ! in: nr. of polarizations
& g_dens, ! in: superposit. atomic
& delchi_ao, ! in: deriv. basis fcts.
& qxyz, ! in: grid points
& qwght, ! in: weighting coeffs.
& nbf, ! in: nr. basis functions
& npts, ! in: nr. grid points
& natoms, ! in: nr. atoms
& origin, ! in: center
& type, ! in: moment type
& ncomp, ! in: no. of components
& x2c0) ! out
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "stdio.fh"
cinclude "x2c.fh"
#include "global.fh"
#include "bas.fh"
c
integer nbf,npts,ao_bas_han,natoms,ncomp,geom
integer g_dens(2),ipol
double precision qwght(npts)
double precision pot(npts)
double precision qxyz(3,npts)
double precision delchi_ao(npts,3,nbf)
double precision origin(3)
character*(*) type
double precision x2c0(nbf,nbf,ncomp)
integer i,j,k, icomp
integer closegridpts(npts)
double precision clight_au2
double precision amat_tot
double precision denom
double precision fac1_arr(npts,ncomp)
double precision ac_sf(6)
external get_ints_x2c_so_efg
c == preliminaries ==
do k = 1,npts
closegridpts(k) = 0
end do
do icomp = 1,ncomp
do k = 1,npts
fac1_arr(k,icomp) = 0.0d0
end do
end do
c == calculate the EFG operator on the grid ==
call grid_moments(origin,type,ncomp,
& npts,qxyz,
& fac1_arr, ! out: moments on grid
& closegridpts)
c black out grid points too close to a nucleus;
c the numerical integration may get noisy otherwise:
c$$$ do k = 1,npts
c$$$ if (k.eq.closegridpts(k)) qwght(k) = 0.d0
c$$$ end do
c === premultiply fac1_arr with integration weights
do icomp = 1,ncomp
do k = 1,npts
fac1_arr(k,icomp)= fac1_arr(k,icomp) * qwght(k)
end do
end do
c == assemble integrals ==
c ----- main diagonal --- START
do i = 1, nbf
j=i
call get_ints_x2c_sc_moment(
& nbf,npts,delchi_ao,i,j,ncomp,
& fac1_arr,
& ac_sf) ! out
do icomp = 1,ncomp
x2c0(i,j,icomp) = x2c0(i,j,icomp) + ac_sf(icomp)
end do
enddo ! end-loop-i
c ----- main diagonal --- END
c ----- off diagonal --- START
do i = 1, nbf
do j = i+1, nbf
call get_ints_x2c_sc_moment(
& nbf,npts,delchi_ao,i,j,ncomp,
& fac1_arr,
& ac_sf) ! out
do icomp = 1,ncomp
x2c0(i,j,icomp) = x2c0(i,j,icomp) + 2.0d0*ac_sf(icomp)
end do
enddo ! end-loop-j
enddo ! end-loop-i
c ----- off diagonal --- END
c all done
return
end
subroutine get_ints_x2c_sc_moment(
& nbf, ! in: # basis fncts.
& npts, ! in: grid points
& delchi_ao, ! in: deriv. of basis fcts.
& i,j, ! in: (i,j) indices for delchi_ao
& ncomp, ! in: number of components,
& fac1_arr, ! operator on grid points * weights
& ac_sf) ! out
implicit none
#include "errquit.fh"
#include "stdio.fh"
#include "global.fh"
integer nbf,npts,i,j,ncomp
integer k, icomp
double precision delchi_ao(npts,3,nbf)
double precision fac1_arr(npts,ncomp)
double precision
& ac_sf(6)
double precision prod0
do icomp = 1,ncomp
ac_sf(icomp) = 0.0d0
end do
do k = 1, npts
prod0 = delchi_ao(k,1,i)*delchi_ao(k,1,j)
& +delchi_ao(k,2,i)*delchi_ao(k,2,j)
& +delchi_ao(k,3,i)*delchi_ao(k,3,j)
do icomp = 1,ncomp
ac_sf(icomp) = ac_sf(icomp) + fac1_arr(k,icomp)*prod0
end do ! icomp
end do ! k
return
end

128
src/nwdft/x2c/ma_solve.F Normal file
View file

@ -0,0 +1,128 @@
subroutine ma_solve(g_a, g_bx)
#include "mafdecls.fh"
#include "global.fh"
#include "rtdb.fh"
#include "util.fh"
#include "stdio.fh"
#include "errquit.fh"
c purpose: solve linear equation system AX = B with square
c non-symmetric matrix. ga_solve relies on SCALAPACK for reliable
c results, which we may not have available in a given build
c routine uses 2n**2 words of extra floating point storage, with
c input matrices of dimension n x n
integer g_a ! inout: coefficient matrix A
integer g_bx ! in: rhs. B. out: solution matrix X
c local variables
integer g_x, g_y
integer l_amat, k_amat
integer l_bmat, k_bmat
integer l_xmat, k_xmat
integer l_ipiv, k_ipiv
integer l_tmp, k_tmp
integer n, n2, ntmp, inf1, inf2, dbg, type
logical master
character*8 pname
c ================================================================
dbg = 0 ! debug setting during code development
pname = "ma_solve"
master = ga_nodeid().eq.0
if (dbg>0.and.master) write(luout,*) 'entering '//pname
n = 0
call ga_inquire(g_a,type,n,n)
call ga_inquire(g_bx,type,ntmp,ntmp)
if (ntmp.ne.n) then
call errquit(
& pname//': input ga dimensions do not match',n,GA_ERR)
end if
if (dbg>0.and.master) write(luout,*) pname//'n = ',n
n2 = n*n
info1 = 0
info2 = 0
c allocate MA arrays and retrieve data from GAs
c amat
if (.not. ma_push_get(mt_dbl, n2, 'amat', l_amat, k_amat))
& call errquit(pname//': failed allocation amat', n2, MA_ERR)
c tmp
if (.not. ma_push_get(mt_dbl, n, 'tmp', l_tmp, k_tmp))
& call errquit(pname//': failed allocation tmp', n, MA_ERR)
c ipiv
if (.not. ma_push_get(mt_int, n, 'ipiv', l_ipiv, k_ipiv))
& call errquit(pname//': failed allocation ipiv', n, MA_ERR)
call dfill(n2, 0.0d0, dbl_mb(k_amat), 1)
call dfill(n, 0.0d0, dbl_mb(k_tmp), 1)
call ga_get(g_a, 1, n, 1, n, dbl_mb(k_amat), n)
c invert matrix A and store in temp GA
call dgetrf(n, n, dbl_mb(k_amat), n, int_mb(k_ipiv), inf1)
call dgetri(n, dbl_mb(k_amat), n, int_mb(k_ipiv),
& dbl_mb(k_tmp), n, inf2)
if(dbg>0.and.master)
& write(luout,*) pname//': inf1, inf2 = ',inf1, inf2
c allocate temp GA and store inverse of A
if (.not.ga_duplicate(g_a, g_x, 'g_x')) call
& errquit(pname//': ga_dup failed g_x ',101, GA_ERR)
call ga_put(g_x, 1,n, 1,n, dbl_mb(k_amat), n)
c MA arrays no longer needed
if (.not.ma_chop_stack(l_amat)) call
& errquit(pname//': ma_chop_stack failed k_amat',l_tmp,MA_ERR)
c allocate one more GA, perform X = A**(-1) B, and
c store the result in g_bx
if (.not.ga_duplicate(g_a, g_y, 'g_y')) call
& errquit(pname//': ga_dup failed g_y ',102, GA_ERR)
call ga_dgemm('N', 'N', n, n, n, 1.d0,
& g_x, g_bx, 0.0d0, g_y)
call ga_copy(g_y, g_bx)
c deallocate GAs
if (.not. ga_destroy(g_y))
& call errquit(pname//': error destroy g_y',102, GA_ERR)
if (.not. ga_destroy(g_x))
& call errquit(pname//': error destroy g_x',101, GA_ERR)
c all done. return
if (dbg>0.and.master) write(luout,*) 'leaving '//pname
return
end

50
src/nwdft/x2c/read_x2c_pc_dip.F Normal file
View file

@ -0,0 +1,50 @@
c
subroutine read_x2c_pc_dip(icomp,g_pcdip)
c
implicit none
c
#include "mafdecls.fh"
#include "global.fh"
#include "tcgmsg.fh"
#include "util.fh"
#include "stdio.fh"
#include "errquit.fh"
#include "msgids.fh"
#include "dra.fh"
#include "inp.fh"
#include "msgtypesf.h"
c
integer icomp
integer g_pcdip
c
character*(nw_max_path_len) fn_pcdip
character*3 ch_comp
c
logical dmat_from_file
external dmat_from_file
c
integer inntsize,ok
c
c prepare file name
call util_file_name('pcdip',.false.,.false.,fn_pcdip)
if (icomp.eq.1) ch_comp='.1'
if (icomp.eq.2) ch_comp='.2'
if (icomp.eq.3) ch_comp='.3'
fn_pcdip = fn_pcdip(1:inp_strlen(fn_pcdip))//ch_comp ! append component
c
c resolve path
call util_file_name_resolve(fn_pcdip, .false.)
c
c read ga from file
if (.not. dmat_from_file(g_pcdip,fn_pcdip))
& call errquit('read_x2c_pc_dip: dmat_from_file',0,
& UNKNOWN_ERR)
c
c propagate status
ok = 1
inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0) ! Propagate status
call ga_sync()
c
return
end

41
src/nwdft/x2c/write_x2c_pc_dip.F Normal file
View file

@ -0,0 +1,41 @@
c
subroutine write_x2c_pc_dip(icomp,g_pcdip)
c
implicit none
c
#include "mafdecls.fh"
#include "global.fh"
#include "tcgmsg.fh"
#include "util.fh"
#include "stdio.fh"
#include "errquit.fh"
#include "msgids.fh"
#include "dra.fh"
#include "inp.fh"
c
integer icomp
integer g_pcdip
character*(nw_max_path_len) fn_pcdip
character*3 ch_comp
c
logical dmat_to_file
external dmat_to_file
c
c prepare file name
call util_file_name('pcdip',.false.,.false.,fn_pcdip)
if (icomp.eq.1) ch_comp='.1'
if (icomp.eq.2) ch_comp='.2'
if (icomp.eq.3) ch_comp='.3'
fn_pcdip = fn_pcdip(1:inp_strlen(fn_pcdip))//ch_comp ! append component
c
c resolve path name
call util_file_name_resolve(fn_pcdip, .false.)
c
c write ga to file
if (.not. dmat_to_file(g_pcdip, fn_pcdip))
& call errquit('write_x2c_pc_dip: dmat_to_file error',0,
& UNKNOWN_ERR)
c
return
end

266
src/nwdft/x2c/x2c_getw_efg_scalar.F Normal file
View file

@ -0,0 +1,266 @@
c
c == compute the matrix elements for EFG operator in
c == small component basis, for x2c
subroutine x2c_getw_efg_scalar(rtdb, g_dens, origin,
& g_w,
& nexc)
c
c modified from zora_getv_so
implicit none
#include "rtdb.fh"
#include "bas.fh"
#include "cdft.fh"
#include "errquit.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "geom.fh"
#include "msgtypesf.h"
#include "msgids.fh"
#include "stdio.fh"
#include "cgridfile.fh"
#include "grid_cube.fh"
c
c == arguments ==
integer rtdb
integer g_dens(2)
double precision origin(3) ! [in] EFG center
integer g_w(6) ! [out]
integer nexc
c
c == local variables ==
integer i,j,k,iefg,ind,nij
double precision rho_n
double precision tmat
double precision dummy(2)
integer iqsh, istep, nxyz, ncontrset
integer ixyz, lxyz, icharge, lcharge, itags, ltags
integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube,
& ictr_buf,iqpts
double precision rad,ke
integer lbas_cset_info, ibas_cset_info,
& lbas_cent_info, ibas_cent_info,
& ldocset, idocset,
& l_rchi_atom,i_rchi_atom,
& l_rq,i_rq,lniz, iniz,
& lchi_ao, ichi_ao,
& ldelchi_ao, idelchi_ao
integer lx2c ,ix2c
integer dbg
integer inntsize,ddblsize,ok
logical grid_file_rewind
external grid_file_rewind,ga_antisymmetrize
logical me
c
me = (ga_nodeid().eq.0)
dbg=0
if (dbg>0 .and. me) write (luout,*) 'entering getw_efg'
c == allocate memory ==
c$$$ do i=1,4
c$$$ do iefg = 1,6
c$$$ if (.not.MA_Push_Get(mt_dbl,nbf_ao*nbf_ao,
c$$$ & 'x2ci',lx2c(i,iefg),ix2c(i,iefg)))
c$$$ & call errquit('x2c_getv_so: x2ci',i + iefg*10, MA_ERR)
c$$$ end do
c$$$ enddo
if (.not.MA_Push_Get(mt_dbl, 6*(nbf_ao**2),
& 'x2ci',lx2c,ix2c))
& call errquit('x2c_getv_so: x2ci', 6*(nbf_ao**2), MA_ERR)
c == set output MA arrays to zero
c$$$ do i = 1,4
c$$$ do iefg = 1,6
c$$$ call dfill(nbf_ao**2, 0.0d0, dbl_mb(ix2c(i,iefg)), 1)
c$$$ end do
c$$$ end do
call dfill(6*(nbf_ao**2), 0.0d0, dbl_mb(ix2c), 1)
c == generate the grid ==
dummy(1) = 0.d0
dummy(2) = 0.d0
c g_w is a dummy array in the following call
if (dbg>0 .and.me) write (luout,*) 'calling grid_quadv0'
call grid_quadv0(rtdb,g_dens,g_w,nexc,rho_n,dummy,tmat)
if (me) write (luout,*) 'rho_n =',rho_n
c
c == ao basis set info used by xc_eval_basis ==
if (.not.bas_numcont(AO_bas_han, ncontrset))
& call errquit('x2c_getv_so:bas_numcont',0, BASIS_ERR)
if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info',
& lbas_cent_info, ibas_cent_info))
& call errquit('x2c_getv_so: cannot allocate bas_cent_info',0,
& MA_ERR)
if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info',
& lbas_cset_info, ibas_cset_info))
& call errquit('x2c_getv_so: cannot allocate bas_cset_info',0,
& MA_ERR)
call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info),
& int_mb(ibas_cset_info), ncenters)
c
if (.not.MA_Push_Get(mt_log, ncontrset, 'docset',
& ldocset, idocset))
& call errquit('x2c_getv_so: cannot allocate ccdocset',
. ncontrset, MA_ERR)
do i=1,ncontrset
log_mb(idocset+i-1)=.true.
enddo
if(.not.MA_push_get(MT_int, ncenters, 'iniz',
& lniz, iniz))
& call errquit("x2c_getv_so:iniz",0, MA_ERR)
do i= 1, ncenters
int_mb(iniz+i-1)=1
enddo
c
nxyz = 3*ncenters
if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz))
& call errquit('x2c_getv_so: cannot allocate xyz',0, MA_ERR)
if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge))
& call errquit('x2c_getv_so: cannot allocate charge',0, MA_ERR)
if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags))
& call errquit('x2c_getv_so: cannot allocate tags',0, MA_ERR)
if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags),
& Dbl_MB(ixyz), Dbl_MB(icharge)))
& call errquit('x2c_getv_so: geom_cart_get failed',74, GEOM_ERR)
if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz))
& call errquit('x2c_getv_so: cannot allocate qxyz',0, MA_ERR)
if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght))
& call errquit('x2c_getv_so: cannot allocate qwght',0, MA_ERR)
if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4,
& 'quad pts buffer', lrqbuf, irqbuf))
& call errquit('x2c_getv_so: quad buffer', 3, MA_ERR)
if (.not. grid_file_rewind())
$ call errquit('x2c_getv_so: rewinding gridpts?', 0,
& UNKNOWN_ERR)
c
c == loop over records in the grid file ==
do iqsh = 1, n_rec_in_file
c
c == define current range of rad. shells and integration center ==
call grid_file_read(n_per_rec, nqpts, ictr_buf,
& rad,dbl_mb(irqbuf),nsubb)
if(nqpts.gt.buffer_size)
& call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR)
c
c == loop over a subset of the grid ==
istep=0
do ncube=1,nsubb
c
c put buf into currently used arrays qxyz and qwght
call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz),
& dbl_mb(iqwght), nqpts, rad,istep)
if(nqpts.ne.0) then
c
c == compute the basis functions over the grid ==
if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom',
& l_rchi_atom,i_rchi_atom))
& call errquit("x2c_getv:rchi_atom",0, MA_ERR)
c
if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq',
& l_rq,i_rq))
& call errquit("x2c_getv_so:rq",0, MA_ERR)
c
c == delchi ==
if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao,
& 'delchi_ao', ldelchi_ao, idelchi_ao))
& call errquit('x2c_getv: delchi_ao',0, MA_ERR)
c
c == chi ==
if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao,
& 'chi_ao', lchi_ao, ichi_ao))
& call errquit('x2c_getv: chi_ao',0, MA_ERR)
call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq),
& dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters)
call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao),
& dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq),
& dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters,
& int_mb(iniz), log_mb(idocset),
& int_mb(ibas_cent_info), int_mb(ibas_cset_info))
c
c calculate EFG integrals
if(dbg>0 .and. me) write (luout,*) 'calling wefg_sc'
call calc_x2c_wefg_scalar(ao_bas_han, geom, ipol, g_dens,
& dbl_mb(idelchi_ao),
& dbl_mb(iqxyz),dbl_mb(iqwght),
& nbf_ao, nqpts, ncenters,
& origin,
& dbl_mb(ix2c)) ! out
if(dbg>0 .and. me) write (luout,*) 'return from wefg_sc'
c == delete memory ==
if(.not.MA_pop_stack(lchi_ao))
& call errquit("x2c_getv: pop chi_ao", 100, MA_ERR)
if(.not.MA_pop_stack(ldelchi_ao))
& call errquit("x2c_getv: pop delchi_ao", 100, MA_ERR)
if(.not.MA_pop_stack(l_rq))
& call errquit("x2c_getv: pop rq", 100, MA_ERR)
if(.not.MA_pop_stack(l_rchi_atom))
& call errquit("x2c_getv: pop rchi_atom",100,MA_ERR)
endif ! nqpts
enddo ! ncube
end do ! iqsh
c
c == delete memory ==
if(.not.MA_pop_stack(lrqbuf))
& call errquit("x2c_getv_so: pop rqbuf", 100, MA_ERR)
if(.not.MA_pop_stack(lqwght))
& call errquit("x2c_getv_so: pop qwght", 100, MA_ERR)
if(.not.MA_pop_stack(lqxyz))
& call errquit("x2c_getv_so: pop qxyz", 100, MA_ERR)
if(.not.MA_pop_stack(ltags))
& call errquit("x2c_getv_so: pop tags", 100, MA_ERR)
if(.not.MA_pop_stack(lcharge))
& call errquit("x2c_getv_so: pop charge", 100, MA_ERR)
if(.not.MA_pop_stack(lxyz))
& call errquit("x2c_getv_so: pop xyz", 100, MA_ERR)
if(.not.MA_pop_stack(lniz))
& call errquit("x2c_getv_so: pop niz", 100, MA_ERR)
if(.not.MA_pop_stack(ldocset))
& call errquit("x2c_getv_so: pop docset", 100, MA_ERR)
if(.not.MA_pop_stack(lbas_cset_info))
& call errquit("x2c_getv_so: pop bas_cset_info", 100, MA_ERR)
if(.not.MA_pop_stack(lbas_cent_info))
& call errquit("x2c_getv_so: pop bas_cent_info", 100, MA_ERR)
c
c == tally up over all the nodes ==
nij = 0
do i = 1, nbf_ao
do j = 1, nbf_ao
do iefg = 1,6
nij = nij + 1
call ga_dgop(msg_excrho,dbl_mb(ix2c +nij-1),1,'+')
end do
enddo
enddo
c == pack into GA and symmetrize ==
do iefg = 1,6
call ga_zero(g_w(iefg))
c call ga_put(g_w(i,iefg),1,nbf_ao,1,nbf_ao,
c & dbl_mb(ix2c(i,iefg)),nbf_ao)
call ga_put(g_w(iefg),1,nbf_ao,1,nbf_ao,
& dbl_mb(ix2c + ((iefg-1)*(nbf_ao**2)) ),nbf_ao)
call ga_symmetrize(g_w(iefg))
end do
call ga_sync()
c
if(.not.MA_chop_stack(lx2c))
& call errquit("x2c_getv_so: pop lx2c", 100, MA_ERR)
if (dbg>0 .and. me) write (luout,*) 'leaving getw_efg'
c
return
end

249
src/nwdft/x2c/x2c_getw_moment_scalar.F Normal file
View file

@ -0,0 +1,249 @@
c
c == compute the matrix elements for electric moment operator in
c == small component basis, for x2c
subroutine x2c_getw_moment_scalar(rtdb, g_dens, origin,
& g_w,
& nexc, type, ncomp)
c
c modified from zora_getv_so
implicit none
#include "rtdb.fh"
#include "bas.fh"
#include "cdft.fh"
#include "errquit.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "geom.fh"
#include "msgtypesf.h"
#include "msgids.fh"
#include "stdio.fh"
#include "cgridfile.fh"
#include "grid_cube.fh"
c
c == arguments ==
integer rtdb
integer g_dens(2)
double precision origin(3) ! [in] EFG center
integer nexc, ncomp
character*(*) type
integer g_w(ncomp) ! [out]
c
c == local variables ==
integer i,j,k,icomp,ind,nij
double precision rho_n
double precision tmat
double precision dummy(2)
integer iqsh, istep, nxyz, ncontrset
integer ixyz, lxyz, icharge, lcharge, itags, ltags
integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube,
& ictr_buf,iqpts
double precision rad,ke
integer lbas_cset_info, ibas_cset_info,
& lbas_cent_info, ibas_cent_info,
& ldocset, idocset,
& l_rchi_atom,i_rchi_atom,
& l_rq,i_rq,lniz, iniz,
& lchi_ao, ichi_ao,
& ldelchi_ao, idelchi_ao
integer lx2c ,ix2c
integer dbg
integer inntsize,ddblsize,ok
logical grid_file_rewind
external grid_file_rewind,ga_antisymmetrize
logical me
c
me = (ga_nodeid().eq.0)
dbg=0
if (dbg>0 .and. me) write (luout,*) 'entering getw_moment'
if (.not.MA_Push_Get(mt_dbl, ncomp*(nbf_ao**2),
& 'x2ci',lx2c,ix2c))
& call errquit('x2c_getv_so: x2ci', ncomp*(nbf_ao**2), MA_ERR)
call dfill(ncomp*(nbf_ao**2), 0.0d0, dbl_mb(ix2c), 1)
c == generate the grid ==
dummy(1) = 0.d0
dummy(2) = 0.d0
c g_w is a dummy array in the following call
if (dbg>0 .and.me) write (luout,*) 'calling grid_quadv0'
call grid_quadv0(rtdb,g_dens,g_w,nexc,rho_n,dummy,tmat)
if (me) write (luout,*) 'rho_n =',rho_n
c
c == ao basis set info used by xc_eval_basis ==
if (.not.bas_numcont(AO_bas_han, ncontrset))
& call errquit('x2c_getv_so:bas_numcont',0, BASIS_ERR)
if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info',
& lbas_cent_info, ibas_cent_info))
& call errquit('x2c_getv_so: cannot allocate bas_cent_info',0,
& MA_ERR)
if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info',
& lbas_cset_info, ibas_cset_info))
& call errquit('x2c_getv_so: cannot allocate bas_cset_info',0,
& MA_ERR)
call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info),
& int_mb(ibas_cset_info), ncenters)
c
if (.not.MA_Push_Get(mt_log, ncontrset, 'docset',
& ldocset, idocset))
& call errquit('x2c_getv_so: cannot allocate ccdocset',
. ncontrset, MA_ERR)
do i=1,ncontrset
log_mb(idocset+i-1)=.true.
enddo
if(.not.MA_push_get(MT_int, ncenters, 'iniz',
& lniz, iniz))
& call errquit("x2c_getv_so:iniz",0, MA_ERR)
do i= 1, ncenters
int_mb(iniz+i-1)=1
enddo
c
nxyz = 3*ncenters
if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz))
& call errquit('x2c_getv_so: cannot allocate xyz',0, MA_ERR)
if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge))
& call errquit('x2c_getv_so: cannot allocate charge',0, MA_ERR)
if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags))
& call errquit('x2c_getv_so: cannot allocate tags',0, MA_ERR)
if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags),
& Dbl_MB(ixyz), Dbl_MB(icharge)))
& call errquit('x2c_getv_so: geom_cart_get failed',74, GEOM_ERR)
if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz))
& call errquit('x2c_getv_so: cannot allocate qxyz',0, MA_ERR)
if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght))
& call errquit('x2c_getv_so: cannot allocate qwght',0, MA_ERR)
if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4,
& 'quad pts buffer', lrqbuf, irqbuf))
& call errquit('x2c_getv_so: quad buffer', 3, MA_ERR)
if (.not. grid_file_rewind())
$ call errquit('x2c_getv_so: rewinding gridpts?', 0,
& UNKNOWN_ERR)
c
c == loop over records in the grid file ==
do iqsh = 1, n_rec_in_file
c
c == define current range of rad. shells and integration center ==
call grid_file_read(n_per_rec, nqpts, ictr_buf,
& rad,dbl_mb(irqbuf),nsubb)
if(nqpts.gt.buffer_size)
& call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR)
c
c == loop over a subset of the grid ==
istep=0
do ncube=1,nsubb
c
c put buf into currently used arrays qxyz and qwght
call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz),
& dbl_mb(iqwght), nqpts, rad,istep)
if(nqpts.ne.0) then
c
c == compute the basis functions over the grid ==
if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom',
& l_rchi_atom,i_rchi_atom))
& call errquit("x2c_getv:rchi_atom",0, MA_ERR)
c
if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq',
& l_rq,i_rq))
& call errquit("x2c_getv_so:rq",0, MA_ERR)
c
c == delchi ==
if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao,
& 'delchi_ao', ldelchi_ao, idelchi_ao))
& call errquit('x2c_getv: delchi_ao',0, MA_ERR)
c
c == chi ==
if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao,
& 'chi_ao', lchi_ao, ichi_ao))
& call errquit('x2c_getv: chi_ao',0, MA_ERR)
call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq),
& dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters)
call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao),
& dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq),
& dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters,
& int_mb(iniz), log_mb(idocset),
& int_mb(ibas_cent_info), int_mb(ibas_cset_info))
c
c calculate integrals
if(dbg>0 .and. me) write (luout,*) 'calling wmoment_sc'
call calc_x2c_wmoment_scalar(
& ao_bas_han, geom, ipol, g_dens,
& dbl_mb(idelchi_ao),
& dbl_mb(iqxyz),dbl_mb(iqwght),
& nbf_ao, nqpts, ncenters,
& origin,type,ncomp,
& dbl_mb(ix2c)) ! out
if(dbg>0 .and. me) write (luout,*) 'return wmoment_sc'
c == delete memory ==
if(.not.MA_pop_stack(lchi_ao))
& call errquit("x2c_getv: pop chi_ao", 100, MA_ERR)
if(.not.MA_pop_stack(ldelchi_ao))
& call errquit("x2c_getv: pop delchi_ao", 100, MA_ERR)
if(.not.MA_pop_stack(l_rq))
& call errquit("x2c_getv: pop rq", 100, MA_ERR)
if(.not.MA_pop_stack(l_rchi_atom))
& call errquit("x2c_getv: pop rchi_atom",100,MA_ERR)
endif ! nqpts
enddo ! ncube
end do ! iqsh
c
c == delete memory ==
if(.not.MA_pop_stack(lrqbuf))
& call errquit("x2c_getv_so: pop rqbuf", 100, MA_ERR)
if(.not.MA_pop_stack(lqwght))
& call errquit("x2c_getv_so: pop qwght", 100, MA_ERR)
if(.not.MA_pop_stack(lqxyz))
& call errquit("x2c_getv_so: pop qxyz", 100, MA_ERR)
if(.not.MA_pop_stack(ltags))
& call errquit("x2c_getv_so: pop tags", 100, MA_ERR)
if(.not.MA_pop_stack(lcharge))
& call errquit("x2c_getv_so: pop charge", 100, MA_ERR)
if(.not.MA_pop_stack(lxyz))
& call errquit("x2c_getv_so: pop xyz", 100, MA_ERR)
if(.not.MA_pop_stack(lniz))
& call errquit("x2c_getv_so: pop niz", 100, MA_ERR)
if(.not.MA_pop_stack(ldocset))
& call errquit("x2c_getv_so: pop docset", 100, MA_ERR)
if(.not.MA_pop_stack(lbas_cset_info))
& call errquit("x2c_getv_so: pop bas_cset_info", 100, MA_ERR)
if(.not.MA_pop_stack(lbas_cent_info))
& call errquit("x2c_getv_so: pop bas_cent_info", 100, MA_ERR)
c
c == tally up over all the nodes ==
nij = 0
do i = 1, nbf_ao
do j = 1, nbf_ao
do icomp = 1,ncomp
nij = nij + 1
call ga_dgop(msg_excrho,dbl_mb(ix2c +nij-1),1,'+')
end do
enddo
enddo
c == pack into GA and symmetrize ==
do icomp = 1,ncomp
call ga_zero(g_w(icomp))
call ga_put(g_w(icomp),1,nbf_ao,1,nbf_ao,
& dbl_mb(ix2c + ((icomp-1)*(nbf_ao**2)) ),nbf_ao)
call ga_symmetrize(g_w(icomp))
end do
call ga_sync()
c
if(.not.MA_chop_stack(lx2c))
& call errquit("x2c_getv_so: pop lx2c", 100, MA_ERR)
if (dbg>0 .and. me) write (luout,*) 'leaving getw_moment'
c
return
end

262
src/nwdft/x2c/x2c_getw_scalar.F Normal file
View file

@ -0,0 +1,262 @@
c
c == compute the matrix elements for scalar x2c ==
subroutine x2c_getw_scalar(rtdb, g_dens,
& g_w, g_v,
& nexc, skipnuc, skipcoul)
c
c modified from zora_getv type routines
implicit none
#include "rtdb.fh"
#include "bas.fh"
#include "cdft.fh"
#include "errquit.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "geom.fh"
#include "msgtypesf.h"
#include "msgids.fh"
#include "stdio.fh"
#include "cgridfile.fh"
#include "grid_cube.fh"
c
c == arguments ==
integer rtdb ! [inp] runtime database GA handle
integer g_dens(2) ! [inp] model pot. density matrix
integer g_w ! [out] W-matrix (=pV.p)
integer g_v(4) ! scratch
integer nexc ! dummy needed for numerical integration code
logical skipnuc ! [inp] whether to skip nuclear potential
logical skipcoul ! [inp] whether to skip atomic Hartree pot
c
c == local variables ==
integer i,j,k,ind,nij
double precision rho_n
double precision tmat
double precision dummy(2)
integer iqsh, istep, nxyz, ncontrset
integer ixyz, lxyz, icharge, lcharge, itags, ltags
integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube,
& ictr_buf,iqpts
double precision rad,ke
integer lbas_cset_info, ibas_cset_info,
& lbas_cent_info, ibas_cent_info,
& ldocset, idocset,
& l_rchi_atom,i_rchi_atom,
& l_rq,i_rq,lniz, iniz,
& lchi_ao, ichi_ao,
& ldelchi_ao, idelchi_ao
integer lx2c,ix2c
integer inntsize,ddblsize,ok
integer g_tmp
logical grid_file_rewind
external grid_file_rewind,ga_antisymmetrize
c ================================================================
if (skipnuc .and. skipcoul) then
if (ga_nodeid().eq.0) write (luout,*)
& 'x2c_getw_scalar: noting to do. Input error?'
call errquit ('X2C_GETW_SC skipnuc=T, skipcoul=T', 1,
& INPUT_ERR)
end if
c debug:
if (ga_nodeid().eq.0) then
write (luout,*) 'hello from x2c_getw_scalar'
write (luout,*) 'skipnuc, skipcoul = ',skipnuc, skipcoul
end if
c == allocate memory ==
if (.not.MA_Push_Get(mt_dbl,nbf_ao*nbf_ao,
& 'x2ci',lx2c,ix2c))
& call errquit('x2c_getv_scalar: x2ci',666, MA_ERR)
c == preliminaries ==
do i= 1, nbf_ao*nbf_ao
dbl_mb(ix2c+i-1)=0.d0
enddo
c == generate the grid ==
dummy(1) = 0.d0
dummy(2) = 0.d0
c g_v is a scratch array in the following call, must
c have dimension (4)
call grid_quadv0(rtdb,g_dens,g_v,nexc,rho_n,dummy,tmat)
if (ga_nodeid().eq.0) write (luout,*) 'rho_n =',rho_n
c
c == ao basis set info used by xc_eval_basis ==
if (.not.bas_numcont(AO_bas_han, ncontrset))
& call errquit('x2c_getv_sc:bas_numcont',0, BASIS_ERR)
if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info',
& lbas_cent_info, ibas_cent_info))
& call errquit('x2c_getv_sc: cannot allocate bas_cent_info',0,
& MA_ERR)
if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info',
& lbas_cset_info, ibas_cset_info))
& call errquit('x2c_getv_sc: cannot allocate bas_cset_info',0,
& MA_ERR)
call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info),
& int_mb(ibas_cset_info), ncenters)
c
if (.not.MA_Push_Get(mt_log, ncontrset, 'docset',
& ldocset, idocset))
& call errquit('x2c_getv_sc: cannot allocate ccdocset',
. ncontrset, MA_ERR)
do i=1,ncontrset
log_mb(idocset+i-1)=.true.
enddo
if(.not.MA_push_get(MT_int, ncenters, 'iniz',
& lniz, iniz))
& call errquit("x2c_getv_sc:iniz",0, MA_ERR)
do i= 1, ncenters
int_mb(iniz+i-1)=1
enddo
c
nxyz = 3*ncenters
if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz))
& call errquit('x2c_getv_sc: cannot allocate xyz',0, MA_ERR)
if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge))
& call errquit('x2c_getv_sc: cannot allocate charge',0, MA_ERR)
if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags))
& call errquit('x2c_getv_sc: cannot allocate tags',0, MA_ERR)
if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags),
& Dbl_MB(ixyz), Dbl_MB(icharge)))
& call errquit('x2c_getv_sc: geom_cart_get failed',74, GEOM_ERR)
if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz))
& call errquit('x2c_getv_sc: cannot allocate qxyz',0, MA_ERR)
if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght))
& call errquit('x2c_getv_sc: cannot allocate qwght',0, MA_ERR)
if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4,
& 'quad pts buffer', lrqbuf, irqbuf))
& call errquit('x2c_getv_sc: quad buffer', 3, MA_ERR)
if (.not. grid_file_rewind())
$ call errquit('x2c_getv_sc: rewinding gridpts?', 0,
& UNKNOWN_ERR)
c
c == loop over records in the grid file ==
do iqsh = 1, n_rec_in_file
c
c == define the current range of radial shells and integration center ==
call grid_file_read(n_per_rec, nqpts, ictr_buf,
& rad,dbl_mb(irqbuf),nsubb)
if(nqpts.gt.buffer_size)
& call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR)
c
c == loop over a subset of the grid ==
istep=0
do ncube=1,nsubb
c
c put buf into currently used arrays qxyz and qwght
call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz),
& dbl_mb(iqwght), nqpts, rad,istep)
if(nqpts.ne.0) then
c
c == compute the basis functions over the grid ==
if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom',
& l_rchi_atom,i_rchi_atom))
& call errquit("x2c_getv:rchi_atom",0, MA_ERR)
c
if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq',
& l_rq,i_rq))
& call errquit("x2c_getv_sc:rq",0, MA_ERR)
c
c == delchi ==
if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao,
& 'delchi_ao', ldelchi_ao, idelchi_ao))
& call errquit('x2c_getv: delchi_ao',0, MA_ERR)
c
c == chi ==
if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao,
& 'chi_ao', lchi_ao, ichi_ao))
& call errquit('x2c_getv: chi_ao',0, MA_ERR)
call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq),
& dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters)
call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao),
& dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq),
& dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters,
& int_mb(iniz), log_mb(idocset),
& int_mb(ibas_cent_info), int_mb(ibas_cset_info))
c
c == calculate spin-orbit x2c integrals (pVp) ==
call calc_x2c_w_scalar(
& ao_bas_han, geom, ipol, g_dens, skipnuc, skipcoul,
& dbl_mb(idelchi_ao),
& dbl_mb(iqxyz),dbl_mb(iqwght),
& nbf_ao, nqpts, ncenters,
& dbl_mb(ix2c)) ! out
c == delete memory ==
if(.not.MA_pop_stack(lchi_ao))
& call errquit("x2c_getv: pop chi_ao", 100, MA_ERR)
if(.not.MA_pop_stack(ldelchi_ao))
& call errquit("x2c_getv: pop delchi_ao", 100, MA_ERR)
if(.not.MA_pop_stack(l_rq))
& call errquit("x2c_getv: pop rq", 100, MA_ERR)
if(.not.MA_pop_stack(l_rchi_atom))
& call errquit("x2c_getv: pop rchi_atom",100,MA_ERR)
endif ! nqpts
enddo ! ncube
end do ! iqsh
c
c == delete memory ==
if(.not.MA_pop_stack(lrqbuf))
& call errquit("x2c_getv_sc: pop rqbuf", 100, MA_ERR)
if(.not.MA_pop_stack(lqwght))
& call errquit("x2c_getv_sc: pop qwght", 100, MA_ERR)
if(.not.MA_pop_stack(lqxyz))
& call errquit("x2c_getv_sc: pop qxyz", 100, MA_ERR)
if(.not.MA_pop_stack(ltags))
& call errquit("x2c_getv_sc: pop tags", 100, MA_ERR)
if(.not.MA_pop_stack(lcharge))
& call errquit("x2c_getv_sc: pop charge", 100, MA_ERR)
if(.not.MA_pop_stack(lxyz))
& call errquit("x2c_getv_sc: pop xyz", 100, MA_ERR)
if(.not.MA_pop_stack(lniz))
& call errquit("x2c_getv_sc: pop niz", 100, MA_ERR)
if(.not.MA_pop_stack(ldocset))
& call errquit("x2c_getv_sc: pop docset", 100, MA_ERR)
if(.not.MA_pop_stack(lbas_cset_info))
& call errquit("x2c_getv_sc: pop bas_cset_info", 100, MA_ERR)
if(.not.MA_pop_stack(lbas_cent_info))
& call errquit("x2c_getv_sc: pop bas_cent_info", 100, MA_ERR)
c
c == tally up over all the nodes ==
nij = 0
do i = 1, nbf_ao
do j = 1, nbf_ao
nij = nij + 1
call ga_dgop(msg_excrho, dbl_mb(ix2c+nij-1), 1, '+')
enddo
enddo
c == pack into GA and symmetrize ==
c Note: we need to ADD to the g_w GA because it may alreday hold
c analytic W integrals with the nuclear potential
if(.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'getw_tmp',0,0, g_tmp))
& call errquit('x2c_getw_sc: error cdreating tmp GA',2, GA_ERR)
call ga_zero(g_tmp)
call ga_put(g_tmp,1,nbf_ao,1,nbf_ao,
& dbl_mb(ix2c),nbf_ao)
call ga_symmetrize(g_tmp)
call ga_add(1.0d0, g_w, 1.0d0, g_tmp, g_w)
if (.not. ga_destroy(g_tmp)) call errquit
& ('x2c_getw_sc: could not destroy g_tmp', 2, GA_ERR)
call ga_sync()
c
if(.not.MA_chop_stack(lx2c))
& call errquit("x2c_getv_sc: pop lx2c", 100, MA_ERR)
c
return
end

102
src/nwdft/x2c/x2c_scf.F Normal file
View file

@ -0,0 +1,102 @@
subroutine x2c_scf_init(rtdb)
implicit none
integer rtdb
#include "rtdb.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "stdio.fh"
#include "x2c_scf.fh"
c
c check for additional x2c keywords and set
c lx2c_keepU to T if we need it for any of the properties:
lx2c_keepU = .false.
lx2c_efg = .false.
if (rtdb_get(rtdb, 'x2c:efg', mt_log, 1,
& lx2c_efg)) lx2c_keepU = .true.
lx2c_hfc = .false.
if (rtdb_get(rtdb, 'x2c:hfc', mt_log, 1,
& lx2c_hfc)) lx2c_keepU = .true.
lx2c_moment = .false.
if (rtdb_get(rtdb, 'x2c:moment', mt_log, 1,
& lx2c_moment)) lx2c_keepU = .true.
lx2c_noatm = .false.
if (rtdb_get(rtdb, 'x2c:noatm', mt_log, 1,
& lx2c_noatm)) continue
if (lx2c_efg .and. ga_nodeid().eq.0) write (luout,*)
& 'x2c: will calculate EFGs, too'
return
end
subroutine x2c_scf_terminate(rtdb,
i ipol,
g g_rel,g_x2c_u,g_x2c_u2c,g_x2c_c2u)
implicit none
integer rtdb
integer ipol
integer g_rel(*)
integer g_x2c_u(*),g_x2c_u2c,g_x2c_c2u
#include "rtdb.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "stdio.fh"
#include "errquit.fh"
#include "x2c_scf.fh"
c
c deallocate x2c arrays if necessary
if (.not. ga_destroy(g_rel(1))) call errquit(
& 'dft_scf: ga_destroy failed g_rel',0, GA_ERR)
if (ipol.gt.1) then
if (.not. ga_destroy(g_rel(2))) call errquit(
& 'dft_scf: ga_destroy failed g_rel',0, GA_ERR)
end if
if (lx2c_keepU) then
if(.not.ga_destroy(g_x2c_u(1))) call errquit(
& 'dft_scf: ga_destroy failed x2c_u(1)',1, GA_ERR)
if(.not.ga_destroy(g_x2c_u(2))) call errquit(
& 'dft_scf: ga_destroy failed x2c_u(2)',2, GA_ERR)
if (lx2c_u2c) then
if(.not.ga_destroy(g_x2c_u2c)) call errquit(
& 'dft_scf: ga_destroy failed x2c_u2c',51, GA_ERR)
end if
if (lx2c_c2u) then
if(.not.ga_destroy(g_x2c_c2u)) call errquit(
& 'dft_scf: ga_destroy failed x2c_c2u',54, GA_ERR)
end if
endif
return
end
subroutine x2c_prop(rtdb,
N nexc, g_dens, g_dens_at, g_x2c_u)
implicit none
integer rtdb
integer nexc, g_dens, g_dens_at, g_x2c_u
#include "rtdb.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "stdio.fh"
#include "x2c_scf.fh"
c
if (lx2c_efg) then
c x2c EFG
if (ga_nodeid().eq.0)
& write (luout,*) 'calling calc_x2c_efg'
call calc_x2c_efg_scalar(rtdb,
& nexc, g_dens, g_dens_at, g_x2c_u)
end if
if (lx2c_moment) then
c x2c dipole moment
if (ga_nodeid().eq.0)
& write (luout,*) 'calling x2c_moment'
call calc_x2c_moment_scalar(rtdb,
& nexc, g_dens, g_dens_at, g_x2c_u, 'dipole', 3)
end if
return
end
logical function x2c_noatm()
implicit none
#include "x2c_scf.fh"
x2c_noatm=lx2c_noatm
return
end

8
src/nwdft/x2c/x2c_scf.fh Normal file
View file

@ -0,0 +1,8 @@
logical lx2c_efg, lx2c_moment, lx2c_hfc, lx2c_keepU
logical lx2c_c2u, lx2c_u2c, lx2c_noatm
c
common /x2c_scf/
& lx2c_efg, lx2c_moment, lx2c_hfc, lx2c_keepU,
& lx2c_c2u, lx2c_u2c, lx2c_noatm

1
src/nwdft/zora/GNUmakefile
View file

@ -8,6 +8,7 @@
gridNuclearPotential.o\
gridHartreePotential.o\
gridQpqPotential.o\
grid_moments.o\
zora_getv_EFGZ4_SR.o\
calc_zora_EFGZ4_SR.o\
zora_getv_EFGZ4_SO.o\

79
src/nwdft/zora/gridQpqPotential.F
View file

@ -122,4 +122,81 @@ c == distance from the grid points to given xyz_EFGcoords() ==
endif
return
end
c $Id$
subroutine gridQpqPotential_all(
& xyz_EFGcoords, ! in : EFG-nuclear coord.
& nqpts, ! in : nr. grid points
& qxyz, ! in : grid points
& amat_Qnucl, ! out: EFG potential
& closegridpts)
c Purpose: Evaluates Q_pq(\vec{r},\vec{r}')
c for all 6 unique pq simultaneously
c index = 1 -> Evaluates Qxx
c = 2 -> Evaluates Qyy
c = 3 -> Evaluates Qzz
c = 4 -> Evaluates Qxy
c = 5 -> Evaluates Qxz
c = 6 -> Evaluates Qyz
c nqpts , number of grid points
c qxyz , grid points
c xyz_EFGcoords, Quadrupole potential is evaluated
c in this point (\vec{r})
c Output:
c amat_Qnucl, Quadrupole potential ev. in the grid
c for integration purpose (\vec{r}')
c
c Note: these are the traceless versions where xx is actually
c 3xx - rr, etc
implicit none
#include "global.fh"
#include "stdio.fh"
#include "zora.fh"
integer igrid,nqpts
integer closegridpts(*)
double precision xyz_EFGcoords(3)
double precision qxyz(3,nqpts)
double precision rx,ry,rz,dist,dist5
double precision amat_Qnucl(nqpts,6) ! output
integer dbg
logical me
c =================================================================
dbg = 0
me = (ga_nodeid().eq.0)
if (dbg>0 .and. me) write (luout,*) 'gridQpq: cutoff =',
& zoracutoff_EFG
c == loop over the grid points ==
do igrid = 1,nqpts
c == distance from the grid points to given xyz_EFGcoords() ==
rx = xyz_EFGcoords(1) - qxyz(1,igrid)
ry = xyz_EFGcoords(2) - qxyz(2,igrid)
rz = xyz_EFGcoords(3) - qxyz(3,igrid)
dist = dsqrt(rx*rx + ry*ry + rz*rz)
if (dist.gt.zoracutoff_EFG) then ! check-cutoff
dist5=dist*dist*dist*dist*dist
amat_Qnucl(igrid,1) = (2.d0*rx*rx-(ry*ry+rz*rz))/dist5
amat_Qnucl(igrid,2) = (2.d0*ry*ry-(rx*rx+rz*rz))/dist5
amat_Qnucl(igrid,3) = (2.d0*rz*rz-(rx*rx+ry*ry))/dist5
amat_Qnucl(igrid,4) = (3.d0*rx*ry)/dist5
amat_Qnucl(igrid,5) = (3.d0*rx*rz)/dist5
amat_Qnucl(igrid,6) = (3.d0*ry*rz)/dist5
else
closegridpts(igrid) = igrid
end if
end do ! end-grid
return
end
c $Id: gridQpqPotential.F 22517 2012-05-24 18:37:27Z jochen $

58
src/nwdft/zora/grid_moments.F Normal file
View file

@ -0,0 +1,58 @@
subroutine grid_moments(
& origin, ! in : origin
& type, ! in: integral type
& ncomp, ! in: no of components
& nqpts, ! in : nr. grid points
& qxyz, ! in : grid points
& amat, ! out: moments on grid
& closegridpts)
c nqpts , number of grid points
c qxyz , grid points
implicit none
#include "global.fh"
#include "stdio.fh"
#include "zora.fh"
#include "errquit.fh"
integer igrid,nqpts,ncomp,dbg
integer closegridpts(*)
double precision origin(3)
double precision qxyz(3,nqpts)
character*(*) type
double precision rx,ry,rz,dist
double precision amat(nqpts,ncomp) ! output
logical master
c =================================================================
dbg = 0
master = (ga_nodeid().eq.0)
if (dbg>0 .and. master) write (luout,*) 'grid_moments: cutoff =',
& zoracutoff_EFG
c == loop over the grid points ==
do igrid = 1,nqpts
c == distance from the grid points to given origin ==
rx = qxyz(1,igrid)-origin(1)
ry = qxyz(2,igrid)-origin(2)
rz = qxyz(3,igrid)-origin(3)
c dist = dsqrt(rx*rx + ry*ry + rz*rz)
if (ncomp.eq.3) then
c dipole
amat(igrid,1) = rx
amat(igrid,2) = ry
amat(igrid,3) = rz
else
call errquit("grid_moments: ncomp unsupported", 666,
& INPUT_ERR)
end if
end do ! end-grid
return
end

1
travis/run_qas.sh
View file

@ -157,6 +157,7 @@ fi
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs dft_he2p_wb97
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs ritddft_pyridine
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs au2-sarc-zora-mp
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs x2c-h2se
if [[ ! -z "$USE_LIBXC" ]]; then
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs libxc_he2+
fi