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

Merge pull request #506 from dmejiar/tddft

Browse source 
Density fitting in LR-TDDFT
This commit is contained in:
Edoardo Aprà 2022-01-26 14:40:26 -08:00 • committed by GitHub
commit 420116f2e8
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
34 changed files with 18203 additions and 142 deletions
Download patch file Download diff file Expand all files Collapse all files

24
QA/tests/ritddft_co/ritddft_co.nw Normal file
View file

@ -0,0 +1,24 @@
START
TITLE "TDHF/aug-cc-pVDZ CO+"
CHARGE 1
GEOMETRY
C 0.0 0.0 0.0
O 1.2 0.0 0.0
symmetry c1
END
BASIS
* library aug-cc-pVDZ
END
basis "cd basis" spherical
* library cc-pvdz-ri
* library aug-cc-pvdz-ri_diffuse
end
DFT
XC ncap
MULT 2
END
TDDFT
RPA
NROOTS 5
END
TASK TDDFT ENERGY

1825
QA/tests/ritddft_co/ritddft_co.out Normal file
View file

File diff suppressed because it is too large Load diff

21
QA/tests/ritddft_h2o/ritddft_h2o.nw Normal file
View file

@ -0,0 +1,21 @@
start
TITLE "BLYP/6-31G* H2O"
GEOMETRY
O 0.00000000 0.00000000 0.12982363
H 0.75933475 0.00000000 -0.46621158
H -0.75933475 0.00000000 -0.46621158
END
BASIS
* library 6-31G*
END
basis "cd basis" spherical
* library "Weigend Coulomb Fitting"
end
DFT
XC becke88 1.0 lyp 1.0
END
TDDFT
RPA
NROOTS 20
END
TASK TDDFT ENERGY

1396
QA/tests/ritddft_h2o/ritddft_h2o.out Normal file
View file

File diff suppressed because it is too large Load diff

46
QA/tests/ritddft_pyridine/ritddft_pyridine.nw Normal file
View file

@ -0,0 +1,46 @@
echo
start
title "RI-TDDFT/LDA geometry optimization of Pyridine with 6-31G"
geometry nocenter
N 0.00000000 0.00000000 1.41599295
C 0.00000000 -1.15372936 0.72067272
C 0.00000000 1.15372936 0.72067272
C 0.00000000 -1.20168790 -0.67391011
C 0.00000000 1.20168790 -0.67391011
C 0.00000000 0.00000000 -1.38406147
H 0.00000000 -2.07614628 1.31521089
H 0.00000000 2.07614628 1.31521089
H 0.00000000 2.16719803 -1.19243296
H 0.00000000 -2.16719803 -1.19243296
H 0.00000000 0.00000000 -2.48042299
symmetry c1
end
basis spherical
* library "6-31G"
end
basis "cd basis" spherical
* library "Weigend Coulomb Fitting"
end
driver
clear
maxiter 100
end
dft
iterations 500
xc slater 1.0 vwn_5 1.0
grid xfine
grid euler
direct
end
tddft
nroots 2
algorithm 1
notriplet
target 1
targetsym a
civecs
grad
root 1
end
end
task tddft optimize

14066
QA/tests/ritddft_pyridine/ritddft_pyridine.out Normal file
View file

File diff suppressed because it is too large Load diff

2
src/ddscf/fock_2e_file.F
View file

@ -154,6 +154,8 @@ c
j = idamax(nbf*nbf, dens(i,1), nfock)
denmax = max(denmax, abs(dens(i,j)))
end do
c return if DM is null (e.g imaginary part of RTTDFT DM at t=0)
if(denmax.lt.1d-12) return
den_tol = min(dentolmax,tol2e/denmax,tol2e/denmax**2) ! Threshold to screen integs only
c
if (ga_nodeid().eq.0 .and. util_print('fockfile',print_debug))

6
src/hessian/analytic/shell_fock_build.F
View file

@ -1267,7 +1267,7 @@ c
, ' ndens ',ndens
call util_flush(luout)
endif
#if 1
#if 0
dims(1) = nfock
dims(2) = nbf
dims(3) = nbf
@ -1282,8 +1282,8 @@ c
g_xc=ga_create3d_atom_blocked(geom, basis,'nfockxc:3dxc', nfock)
call nga_distribution(g_xc,
. ga_nodeid(), ilo,ihi)
write(6,123) ga_nodeid(),' gxc distr ilo',ilo
write(6,123) ga_nodeid(),' gxc distr ihi',ihi
c write(6,123) ga_nodeid(),' gxc distr ilo',ilo
c write(6,123) ga_nodeid(),' gxc distr ihi',ihi
123 format(i2,a,3i5)
#endif

6
src/nwdft/grid/grid_quadv0b.F
View file

@ -389,6 +389,7 @@ c AJL/End
& ttau, kske)
endif
else if (do_cpks_l.or.do_gxc) then
if (odftps) call pstat_on(ps_xcrho)
call xc_rhogen(1,
T tol_rho, AO_bas_han, g_dens, nbf_ao_mxnbf_ce,
N natoms, 1,1, nmat+1,
@ -399,6 +400,7 @@ c AJL/End
& rchi_atom, rdelchi_atom,
& rdens_atom, cetobfr,wmax,
& ttau, kske, .false.)
if (odftps) call pstat_off(ps_xcrho)
else if (do_cpks_r) then
call errquit(pname//'what do do for cpks rhs?',0, CALC_ERR)
endif
@ -744,11 +746,14 @@ c
c Combine derivs of functional with deriv densities to get coeffs
c of bfn pairs in CPKS LHS --> these are written back into rho, delrho, ttau
c
if (odftps) call pstat_on(ps_cpkscoef)
call xc_cpks_coeff(rho, delrho, ttau,
& Amat2, Cmat2, Cmat, Mmat2,
& delrho(iscf_delrho), nmat, ipol, nqpts, grad, kske,
& triplet)
if (odftps) call pstat_off(ps_cpkscoef)
c
if (odftps) call pstat_on(ps_tabcd)
call xc_tabcd(1,l3d,
T tol_rho,Fmat,Pmat, rho, Bmat,delrho, Dbl_MB(iDmat),
N 0d0,0d0,0d0,0d0,
@ -760,6 +765,7 @@ c
& rdens_atom, cetobfr,
& kske, ttau,
& .false.,0d0,0d0,0d0)
if (odftps) call pstat_off(ps_tabcd)
c
c Daniel (1-11-13): Need specific control over contraction for the
c XC-third derivatives

10
src/nwdft/lr_tddft/GNUmakefile
View file

@ -6,13 +6,17 @@ tddft_subspace.o tddft_diagon.o tddft_residual.o tddft_analysis.o tddft_dftinit.
tddft_finalize.o tddft_schmidt.o tddft_fold.o \
tddft_bdtensor.o tddft_patch.o tddft_civecs_phase.o tddft_oscstr.o \
tddft_checkpoint.o tddft_restart.o tddft_s2.o tddft_simpleso.o \
tddft_lanczos.o tddft_lanczos_tda.o tddft_dipole.o tddft_prod_v.o
tddft_lanczos.o tddft_lanczos_tda.o tddft_dipole.o tddft_prod_v.o \
tddft_basis.o
OBJ_OPTIMIZE =
OBJ_OPTIMIZE = tddft_fock_fit.o
LIBRARY = libnwdft.a
USES_BLAS = tddft_diagon.F tddft_transfm.F tddft_analysis.F tddft_restart.F tddft_civecs_phase.F
LIB_INCLUDES = -I../include
USES_BLAS = tddft_diagon.F tddft_transfm.F tddft_analysis.F tddft_restart.F tddft_civecs_phase.F \
tddft_fock_fit.F
LIB_DEFINES = -DDEBUG_PRINT
#

158
src/nwdft/lr_tddft/tddft_basis.F Normal file
View file

@ -0,0 +1,158 @@
subroutine tddft_basis(rtdb,geom_in,aobasis,lhascd)
implicit none
#include "rtdb.fh"
#include "cdft.fh"
#include "errquit.fh"
#include "bas.fh"
#include "geom.fh"
#include "global.fh"
#include "mafdecls.fh"
#include "stdio.fh"
logical lhascd
integer rtdb, geom_in
integer aobasis
character(len=13),parameter :: pname = 'tddft_basis: '
logical ok, iolgc
logical, external :: int_normalize, int_norm_2c
integer g_tmpb,lmiss,me,info
lhascd = .false.
me = ga_nodeid()
iolgc = noio.ne.1
geom = geom_in
c
c Define the contracted Gaussian functions for the AO basis set.
c
ok = bas_destroy(AO_bas_han)
ok = bas_destroy(CD_bas_han)
if (.not. bas_create(AO_bas_han, 'ao basis'))
& call errquit(pname//'bas_create?', 86, BASIS_ERR)
if (.not. bas_rtdb_load(rtdb, geom, AO_bas_han, 'ao basis'))
& call errquit(pname//'no ao basis ', 86, BASIS_ERR)
if (.not. int_normalize(rtdb,AO_bas_han))
& call errquit(pname//'int_normalize failed', 86, INT_ERR)
if (.not. bas_numbf(AO_bas_han, nbf))
& call errquit(pname//'basis set error:', 86, BASIS_ERR)
nbf_ao = nbf
if (.not. bas_nprim_cn_max(AO_bas_han,nbf_ao_mxprim))
& call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (.not. bas_high_angular(AO_bas_han,nbf_ao_mxang))
& call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (.not. bas_ncontr_cn_max(AO_bas_han,nbf_ao_mxcont))
& call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (.not. bas_nbf_cn_max(AO_bas_han,nbf_ao_mxnbf_cn))
& call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (.not. bas_nbf_ce_max(AO_bas_han,nbf_ao_mxnbf_ce))
& call errquit(pname//'basis set error:', 86, BASIS_ERR)
if (.not. bas_numcont(AO_bas_han,nshells_ao))
& call errquit(pname//'basis set error:', 86, BASIS_ERR)
aobasis = ao_bas_han
if (.not. bas_create(CD_bas_han, 'cd basis'))
& call errquit(pname//'bas_create?',102, BASIS_ERR)
CDFIT=bas_rtdb_load(rtdb, geom, CD_bas_han, 'cd basis')
if (CDFIT)then
call int_init(rtdb, 1, cd_bas_han)
if (.not. int_norm_2c(rtdb,CD_bas_han))
& call errquit('rdinput: int_norm_2c failed', 208, BASIS_ERR)
call int_terminate()
if (.not. bas_numbf(CD_bas_han, nbf_cd))
& call errquit('rdinput: basis set error:', 86, BASIS_ERR)
if (.not. bas_nprim_cn_max(CD_bas_han,nbf_cd_mxprim))
& call errquit('rdinput: basis set error:', 86, BASIS_ERR)
if (.not. bas_high_angular(CD_bas_han,nbf_cd_mxang))
& call errquit('rdinput: basis set error:', 86, BASIS_ERR)
if (.not. bas_ncontr_cn_max(CD_bas_han,nbf_cd_mxcont))
& call errquit('rdinput: basis set error:', 86, BASIS_ERR)
if (.not. bas_nbf_cn_max(CD_bas_han,nbf_cd_mxnbf_cn))
& call errquit('rdinput: basis set error:', 86, BASIS_ERR)
if (.not. bas_nbf_ce_max(CD_bas_han,nbf_cd_mxnbf_ce))
& call errquit('rdinput: basis set error:', 86, BASIS_ERR)
if (.not. bas_numcont(CD_bas_han,nshells_cd))
& call errquit('rdinput: basis set error:', 86, BASIS_ERR)
if (.not. rtdb_put(rtdb, 'dft:cdfit', mt_log, 1, .true.))
& call errquit('rdinput: failed put_rtdb:', 99, RTDB_ERR)
c
c texas breaks when AO_spherical and CD_cartesian
c
if(bas_is_spherical(AO_bas_han).and.
N (.not.bas_is_spherical(CD_bas_han)))
T call int_app_set_no_texas(rtdb)
else
if (.not. bas_destroy(CD_bas_han))
& call errquit('rdinput:not able to destroy CD_bas:',86,
& BASIS_ERR)
CD_bas_han=-99999
nbf_cd = 0
nbf_cd_mxprim = 0
nbf_cd_mxang = 0
nbf_cd_mxcont = 0
nbf_cd_mxnbf_cn = 0
nbf_cd_mxnbf_ce = 0
nshells_cd = 0
endif
if (cd_bas_han.eq.-99999) then
call int_init(rtdb,1,ao_bas_han)
else
call int_init(rtdb,2,(/ao_bas_han,cd_bas_han/))
endif
if (CDFIT)then
if (.not. ga_create(mt_dbl, nbf_cd, nbf_cd, 'CD 2cERI',
& 0, nbf_cd, g_2ceri))
& call errquit(pname//'Error creating g_2ceri',0,
& GA_ERR)
call ga_zero(g_2ceri)
call dft_get2eri(CD_bas_han, g_2ceri,.false.)
call ga_sync()
if (.not. ga_duplicate(g_2ceri, g_cdinv, 'CD 2cERInv'))
& call errquit(pname//'Error creating g_cdinv',0, GA_ERR)
call ga_zero(g_cdinv)
lmiss = 1
info = 0
call ga_chol_seq(g_2ceri, g_cdinv, info)
if (info.ne.0)then
if (me.eq.0)then
write(LuOut,*)' Problem in performing a Cholesky '
write(LuOut,*)' decomposition of the 2-ctr ERI '
write(LuOut,*)' matrix using CD fitting basis. '
write(LuOut,*)' Attempting a diag/inverse. '
endif
endif
if (info.eq.0) then
g_tmpb = g_2ceri
call ga_copy(g_cdinv, g_tmpb)
call ga_inv_seq(g_tmpb, g_cdinv)
else
call dft_invdiag(g_2ceri, g_cdinv,nbf_cd)
endif
if (IOLGC.and.(me.eq.0)) then
lmiss = 0
call dft_invio('CDI', g_cdinv, nbf_cd, 'WRITE', lmiss)
if (lmiss.ne.0)call errquit
& (pname//'dft_invio - abnormal write of CDI ', 0,DISK_ERR)
lmiss = 0
call dft_invio('CD', g_2ceri, nbf_cd, 'WRITE', lmiss)
if (lmiss.ne.0)call errquit
& (pname//'dft_invio - abnormal write of CD ', 0,DISK_ERR)
endif
if (IOLGC) then
if (.not. ga_destroy(g_cdinv)) call errquit
& (pname//'Could not destroy g_xcinv', 0, GA_ERR)
if (.not. ga_destroy(g_2ceri)) call errquit
& (pname//'Could not destroy g_xcinv', 0, GA_ERR)
endif
lhascd = .true.
endif
return
end

33
src/nwdft/lr_tddft/tddft_contract.F
View file

@ -1,6 +1,6 @@
subroutine tddft_contract(rtdb,geom,ao_bas_han,
1 g_td,g_pd,g_qd,nbf_ao,ipol,tol2e,tda,oskel,kfac,
2 lhashf,triplet)
2 lhashf,triplet,lhascd)
c
c $Id$
c
@ -46,6 +46,7 @@ c
double precision kfac ! Prefactor for HF exchange
logical lhashf ! True if XC has HF contribution
logical triplet ! True if restricted triplet
logical lhascd ! True if CD basis is available
c
integer idim(3) ! Dimension array for NGA create
integer ichunk(3) ! Chunk array for NGA create
@ -313,7 +314,7 @@ c
c
subroutine tddft_nga_cont(rtdb,geom,ao_bas_han,
1 g_ntrl,g_nprd_p,g_nprd_m,nbf_ao,ipol,tol2e,tda,oskel,kfac,
2 lhashf,triplet,nvectors,iwhich)
2 lhashf,triplet,nvectors,iwhich,lhascd)
c
c $Id: tddft_contract.F 26045 2014-08-26 04:09:46Z niri $
c
@ -352,6 +353,7 @@ c
double precision jfac ! Prefactor for Coulomb exchange
logical lhashf ! True if XC has HF contribution
logical triplet ! True if restricted triplet
logical lhascd ! True if CD basis available
integer nvectors ! Number of vectors to be processed
c
integer idim(3) ! Dimension array for NGA create
@ -458,8 +460,13 @@ c
enddo
if(.not.cam_exch) then ! normal calculations
c write(6,*)'HERE'
call shell_fock_build(geom,ao_bas_han,0,nvectors,
1 dbl_mb(k_jf),dbl_mb(k_kf),tol2e,g_nsym,g_nprd_p,.true.)
if (lhascd.and..not.lhashf) then
call tddft_fock_fit(0,nvectors,dbl_mb(k_jf),tol2e,g_nsym,
1 g_nprd_p)
else
call shell_fock_build(geom,ao_bas_han,0,nvectors,
1 dbl_mb(k_jf),dbl_mb(k_kf),tol2e,g_nsym,g_nprd_p,.true.)
endif
else ! attenuated calculations
c write(6,*)'NOT HERE'
@ -472,8 +479,13 @@ c write(6,*)'NOT HERE'
dbl_mb(k_kf+n-1)=-kfac*0.5d0*scale_k
enddo
if (.not.cam_exch) then ! normal calculations
call shell_fock_build(geom,ao_bas_han,0,nvectors,
1 dbl_mb(k_jf),dbl_mb(k_kf),tol2e,g_nsym,g_nprd_p,.true.)
if (lhascd.and..not.lhashf) then
call tddft_fock_fit(0,nvectors,dbl_mb(k_jf),tol2e,g_nsym,
1 g_nprd_p)
else
call shell_fock_build(geom,ao_bas_han,0,nvectors,
1 dbl_mb(k_jf),dbl_mb(k_kf),tol2e,g_nsym,g_nprd_p,.true.)
endif
else ! attenuated calculations
call shell_fock_build_cam(geom,ao_bas_han,0,nvectors,
1 dbl_mb(k_jf),dbl_mb(k_kf),tol2e,g_nsym,g_nprd_p,.true.)
@ -488,8 +500,13 @@ c write(6,*)'NOT HERE'
call uhf_fock_setup(g_nsym,g_nprd_p,dbl_mb(k_jf),dbl_mb(k_kf),
1 ndu,nfu,nbf_ao,.false.)
if (.not.cam_exch) then ! normal calculations
call shell_fock_build(geom,ao_bas_han,0,ndu,
1 dbl_mb(k_jf),dbl_mb(k_kf),tol2e,g_nsym,g_nprd_p,.true.)
if (lhascd.and..not.lhashf) then
call tddft_fock_fit(0,ndu,dbl_mb(k_jf),tol2e,g_nsym,
1 g_nprd_p)
else
call shell_fock_build(geom,ao_bas_han,0,ndu,
1 dbl_mb(k_jf),dbl_mb(k_kf),tol2e,g_nsym,g_nprd_p,.true.)
endif
else
call shell_fock_build_cam(geom,ao_bas_han,0,ndu,
1 dbl_mb(k_jf),dbl_mb(k_kf),tol2e,g_nsym,g_nprd_p,.true.)

23
src/nwdft/lr_tddft/tddft_davidson.F
View file

@ -3,7 +3,8 @@
2 g_buffer1,g_buffer2,d_trials,d_prod_p,d_prod_m,
3 g_movecs,evl,ntrials,nproducts,
4 tda,oskel,singlet,triplet,tol2e,kfac,lhashf,thresh,
5 maxiter,target,targetsym,algorithm,mult,l_irs,k_irs)
5 maxiter,target,targetsym,algorithm,mult,l_irs,k_irs,
6 lhascd)
c
c $Id$
c
@ -130,6 +131,7 @@ c
logical singlet ! True if singlet excited states
logical triplet ! True if triplet excited states
logical oskel ! True if use symmetry in fock_2e
logical lhascd ! True if a CD basis is available
double precision tol2e ! 2-e integral tolerance
double precision kfac ! Prefactor for HF exchange contribution
logical lhashf ! True if XC has HF contribution
@ -188,8 +190,8 @@ c
integer ispintag
integer iwhich
c
integer ga_create_atom_blocked
external ga_create_atom_blocked
integer ga_create_atom_blocked, ga_create3d_atom_blocked
external ga_create_atom_blocked, ga_create3d_atom_blocked
c
nodezero=(ga_nodeid().eq.0)
iwhich = 0 ! call to tddft_nga_cont()
@ -349,6 +351,7 @@ c
idim(1)=ipol*nvectors
idim(2)=nbf_ao
idim(3)=nbf_ao
#if 0
ichunk(1)=ipol*nvectors
ichunk(2)=-1
ichunk(3)=-1
@ -363,6 +366,16 @@ c
1 ichunk,g_nprd_m)) call errquit
2 ('tddft_davidson: failed to create g_nprd_m',0, GA_ERR)
endif
#else
g_ntrl=ga_create3d_atom_blocked(geom,ao_bas_han,'Work trial',
1 ipol*nvectors)
g_nprd_p=ga_create3d_atom_blocked(geom,ao_bas_han,
1 'Work product p',ipol*nvectors)
if (.not.tda) then
g_nprd_m=ga_create3d_atom_blocked(geom,ao_bas_han,
1 'Work product m',ipol*nvectors)
endif
#endif
c
c -----------------------
c MO to AO transformation
@ -414,7 +427,7 @@ c ---------------------------------------------------------
c
call tddft_nga_cont(rtdb,geom,ao_bas_han,
1 g_ntrl,g_nprd_p,g_nprd_m,nbf_ao,ipol,tol2e,tda,oskel,kfac,
2 lhashf,triplet,nvectors,iwhich)
2 lhashf,triplet,nvectors,iwhich,lhascd)
ccpu=ccpu+util_cpusec()
cwall=cwall+util_wallsec()
if (nodezero.and.util_print('timings',print_high))
@ -617,7 +630,7 @@ c ---------------------------------------------------------
c
call tddft_contract(rtdb,geom,ao_bas_han,
1 g_td,g_pd,g_qd,nbf_ao,ipol,tol2e,tda,oskel,kfac,
2 lhashf,triplet)
2 lhashf,triplet,lhascd)
ccpu=ccpu+util_cpusec()
cwall=cwall+util_wallsec()
if (nodezero.and.util_print('timings',print_high))

13
src/nwdft/lr_tddft/tddft_energy.F
View file

@ -66,6 +66,7 @@ c
logical tda ! True if Tamm-Dancoff approximation
logical singlet ! True if singlet excited states
logical triplet ! True if triplet excited states
logical lhascd ! True if a CD basis was specified
integer ispintag ! 1 = restricted singlet, 2 = restricted triplet, 3 = unrestricted
integer nroots ! Number of roots sought
integer maxvecs ! Maximum number of trial vectors
@ -194,7 +195,8 @@ c
4 tda,singlet,triplet,nroots,maxvecs,tol2e,thresh,
5 maxiter,target,targetsym,symmetry,algorithm,g_movecs,
6 l_evl,k_evl,l_occ,k_occ,l_irs,k_irs,mult,lecut,ecut,diff_max,
7 lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
7 lowin,owstart,owend,lewin,ewinl,ewinh,vshift,
8 lhascd)
c
c ----------------------------------
c Cosmo initialization, if necessary
@ -386,7 +388,8 @@ c tddft_analysis is called at the end of tddft_davidson
2 g_buffer1,g_buffer2,d_trials,d_prod_p,d_prod_m,
3 g_movecs,dbl_mb(k_evl),ntrials,nproducts,
4 tda,oskel,.true.,.false.,tol2e,xfac(1),lhashf,thresh,
5 maxiter,target,targetsym,algorithm,mult,l_irs,k_irs)
5 maxiter,target,targetsym,algorithm,mult,l_irs,k_irs,
6 lhascd)
c
endif ! singlet
c
@ -427,7 +430,8 @@ c tddft_analysis is called at the end of tddft_davidson
2 g_buffer1,g_buffer2,d_trials,d_prod_p,d_prod_m,
3 g_movecs,dbl_mb(k_evl),ntrials,nproducts,
4 tda,oskel,.false.,.true.,tol2e,xfac(1),lhashf,thresh,
5 maxiter,target,targetsym,algorithm,mult,l_irs,k_irs)
5 maxiter,target,targetsym,algorithm,mult,l_irs,k_irs,
6 lhascd)
c
endif ! triplet
c
@ -463,7 +467,8 @@ c
2 g_buffer1,g_buffer2,d_trials,d_prod_p,d_prod_m,
3 g_movecs,dbl_mb(k_evl),ntrials,nproducts,
4 tda,oskel,.false.,.false.,tol2e,xfac(1),lhashf,thresh,
5 maxiter,target,targetsym,algorithm,mult,l_irs,k_irs)
5 maxiter,target,targetsym,algorithm,mult,l_irs,k_irs,
6 lhascd)
endif
c
c Clean up cosmo stuff, if necessary ==

2
src/nwdft/lr_tddft/tddft_finalize.F
View file

@ -37,6 +37,8 @@ c Basis set
c ---------
if (.not.bas_destroy(ao_bas_han))
1 call errquit('tddft_finalize: bas_destroy failed',0, BASIS_ERR)
call tddft_finalize_fit()
c ----------
c MO vectors
c ----------

463
src/nwdft/lr_tddft/tddft_fock_fit.F Normal file
View file

@ -0,0 +1,463 @@
subroutine tddft_fock_fit(nder,ndens,jfac,tol2e,g_dens,g_fock)
implicit none
#include "errquit.fh"
c
#include "cdft.fh"
#include "bas.fh"
#include "mafdecls.fh"
#include "stdio.fh"
#include "global.fh"
#include "schwarz.fh"
#include "sym.fh"
#include "util.fh"
integer ndens
integer nfocks
integer nder
integer g_dens
integer g_fock
double precision tol2e
double precision jfac(*)
integer lo(3), hi(3), ld(2)
integer ilo, ihi, jlo, jhi, adrc, ldc
integer nproc, me, n1, i1
logical dowork, v_nonzero, doshell, iolgc
integer nao2
integer lTvec, iTvec
integer l_p, k_p
integer lgvec, igvec
integer leri, ieri
integer lscr, iscr, nscr
integer lcdc, icdc
integer lcd_shell, icd_shell
integer maxg, mscratch_2e3c
integer atom_c,ishc,sh_lo_c,sh_hi_c,nbf_at_c,nshbfc,nshbf_cd
integer bf_lo_c,bf_hi_c
integer bf_lo_d,bf_hi_d
integer atom_d,ishd,sh_lo_d,sh_hi_d,nbf_at_d,nshbfd
integer ishd_0, ishd_1, ifirstc, ifirstd, ilastc, ilastd
integer atom_p,ishp,ifirstp,ilastp,nshbfp
integer size_atblock
integer lmiss,iptr,j
integer ishbf,lshbf
integer icount, icount_c, icount_d
integer istart, nintegrals
integer idens,i,nnii,nfunz,ninteg_api
integer g_tmp
integer ityp, ndim, dims(3)
double precision threshold, az, shmax, cd_coef_max
double precision pmat_max,eri_est,fac,vcoul_est,q1,qscale,qsym
double precision maxf, t0
logical,external :: xc_gotxc
integer,external :: nxtask_static, nxtask
double precision, external :: ddot, dabsmax
! GA information
me = ga_nodeid()
nproc = ga_nnodes()
call nga_inquire(g_fock, ityp, ndim, dims)
nfocks = dims(1)
! Whether we read or not CDINV from disk
iolgc = noio.ne.1
! Integral screening threshold
THRESHOLD = 10.d0**(-itol2e-2)
! Initialize Fock matrices
call ga_zero(g_fock)
! Skip ERIs and go straight to XC part
if (sum(abs(jfac(1:ndens))).eq.0d0) goto 100
! Memory allocation
if(.not.MA_Push_Get(MT_Dbl,nbf_cd*ndens,'Tvec',lTvec,iTvec))
& call errquit('tddft_fock_fit: cannot allocate Tvec',0, MA_ERR)
call dcopy(nbf_cd*ndens, 0.D0, 0, Dbl_MB(iTvec), 1)
if(.not.MA_Push_Get(MT_Dbl,nbf_cd*ndens,'CDCvec',lcdc,icdc))
& call errquit('tddft_fock_fit: cannot allocate CDCvec',0, MA_ERR)
call dcopy(nbf_cd*ndens, 0.D0, 0, Dbl_MB(icdc), 1)
Nao2 = nbf_ao_mxnbf_ce*nbf_ao_mxnbf_ce
if(.not.MA_Push_Get(MT_Dbl,Nao2*ndens,'gather vector',l_p,k_p))
& call errquit('tddft_fock_fit: cannot allocate gather vector',0,
& MA_ERR)
Nao2 = nbf_ao_mxnbf_cn*nbf_ao_mxnbf_cn
if(.not.MA_Push_Get(MT_Dbl,Nao2*ndens,'gathervector',
& lgvec,igvec))
& call errquit('tddft_fock_fit: cannot allocate gather vector',0,
& MA_ERR)
call int_mem_2e3c(maxg, mscratch_2e3c)
if(.not.MA_Push_Get(MT_Dbl,maxg*nshells_cd,'ERI',lERI,iERI))
& call errquit('dft_fitcd: cannot allocate ERI',0, MA_ERR)
nscr = mscratch_2e3c
if(.not.MA_Push_Get(MT_Dbl,nscr,'scr',lscr,iscr))
& call errquit('dft_fitcd: cannot allocate scr',0, MA_ERR)
! Limits for NGA_put and NGA_get
lo(1) = 1
hi(1) = ndens
ld(1) = ndens
! Part 1. Contract 3c-ERIs with density matrices
icount = 0
n1 = nxtask_static(nproc,1)
do atom_c=1,ncenters
if (.not. bas_ce2cnr( AO_bas_han, atom_c, sh_lo_c, sh_hi_c))
& call errquit('Exiting in tddft_fock_fit',110, BASIS_ERR)
if (.not. bas_ce2bfr( AO_bas_han, atom_c, bf_lo_c, bf_hi_c))
& call errquit('Exiting in tddft_fock_fit',110, BASIS_ERR)
nbf_at_c=bf_hi_c - bf_lo_c + 1
do atom_d=1,atom_c
dowork=n1.eq.icount
! Skip until assigned task
if(.not.dowork) then
icount = icount + 1
cycle
endif
if (.not.bas_ce2bfr(AO_bas_han,atom_d,bf_lo_d,bf_hi_d))
& call errquit('Exiting in tddft_fock_fit',110, BASIS_ERR)
if (.not.bas_ce2cnr(AO_bas_han,atom_d,sh_lo_d,sh_hi_d))
& call errquit('Exiting in tddft_fock_fit',110, BASIS_ERR)
nbf_at_d=bf_hi_d - bf_lo_d + 1
! Get DM block for current atom pair
size_atblock= nbf_at_d*nbf_at_c
lo(2:3) = (/ bf_lo_d, bf_lo_c /)
hi(2:3) = (/ bf_hi_d, bf_hi_c /)
ld(2) = nbf_at_d
call nga_get(g_dens,lo,hi,dbl_mb(k_p),ld)
icount_c=0
do ishc = sh_lo_c,sh_hi_c
if (.not. bas_cn2bfr(AO_bas_han,ishc,ifirstc,ilastc))
& call errquit('Exiting in dft_fitcd',111, BASIS_ERR)
nshbfc = ilastc-ifirstc+1
icount_d=0
do ishd = sh_lo_d, sh_hi_d
c
if (.not. bas_cn2bfr(AO_bas_han,ishd,ifirstd,ilastd))
& call errquit('Exiting in dft_fitcd',1111, BASIS_ERR)
nshbfd = ilastd-ifirstd+1
c
! First schwarz screening
ERI_est = schwarz_shell(ishc,ishd)
if(ERI_est.lt.THRESHOLD) then
icount_d = icount_d + nshbfd
cycle
endif
! Copy DM block for current shell pair
call tddft_copyfrom(ndens,nbf_at_d,nbf_at_c,nshbfd,
+ nshbfc,dbl_mb(k_p),dbl_mb(igvec),icount_d,icount_c,
+ jfac)
c
Nao2 = nshbfd*nshbfc
! Second Schwarz screening
Pmat_max = dabsmax(Nao2*ndens,DBL_MB(igvec))
if( Pmat_max*ERI_est.lt.THRESHOLD )then
icount_d = icount_d + nshbfd
cycle
endif
fac = 2.0d0
if(atom_c.eq.atom_d) fac = 1.d0
istart = iERI
do ishp = 1,nshells_cd
c
if (.not.bas_cn2bfr(CD_bas_han,ishp,ifirstp,ilastp))
& call errquit('fitcd.',5, BASIS_ERR)
nshbfp = ilastp-ifirstp+1
c
nintegrals = nshbfp*Nao2
Ninteg_api = maxg
c
call int_2e3c(CD_bas_han,ishp,AO_bas_han,
I ishc,ishd,nscr,Dbl_MB(iscr),
N Ninteg_api,Dbl_MB(istart))
istart = istart+nintegrals
enddo
! Contract DM blocks with 3c-ERIs
call dgemm('T','N',nbf_cd,ndens,nao2,fac,
. dbl_mb(ieri),nao2,dbl_mb(igvec),nao2,1d0,
. dbl_mb(itvec),nbf_cd)
icount_d = icount_d + nshbfd
enddo
icount_c = icount_c + nshbfc
enddo
n1 = nxtask_static(nproc,1)
icount = icount + 1
enddo
enddo
n1 = nxtask_static(-nproc,1)
c
! Collect results from all ranks
call ga_mask_sync(.false.,.true.)
call ga_dgop(1,dbl_mb(iTvec),nbf_cd*ndens,'+')
call ga_sync()
! Part 2. Multiply by inverse metric
! Read CDINV if necessary
if (IOLGC)then
if (.not. ga_create(mt_dbl, nbf_cd, nbf_cd, 'CD 2cERInv',
& nbf_cd, 0, g_cdinv))
& call errquit('dft_fitcd: error creating g_cdinv',0, GA_ERR)
if (me.eq.0)then
lmiss = 0
call dft_invio('CDI',g_cdinv,nbf_cd,'READ',lmiss)
endif
call ga_brdcst(2000, lmiss,MA_sizeof(MT_INT,1,MT_BYTE),0)
if(lmiss.ne.0)call errquit
& (' dft_fitcd; dft_invio - abnormal read of CDI ',0, GA_ERR)
endif
!Multiply T and by CDINV
call ga_distribution(g_cdinv,me,ilo,ihi,jlo,jhi)
if (ilo.gt.0) then
call ga_access(g_cdinv, ilo, ihi, jlo, jhi, adrc, ldc)
nnii=ihi-ilo+1
do j = jlo,jhi
iptr=adrc+(j-jlo)*ldc
do idens = 1, ndens
AZ = ddot(nnii,DBL_mb(iptr),1,
& DBL_MB(iTvec+(idens-1)*nbf_cd+ilo-1),1)
dbl_mb(icdc+(idens-1)*nbf_cd + j - 1) = az
enddo
enddo
call ga_release(g_cdinv, ilo, ihi, jlo, jhi)
endif
call ga_mask_sync(.true.,.false.)
call ga_dgop(1,dbl_mb(icdc),nbf_cd*ndens,'+')
! Destroy GA if necessary
if (iolgc) then
if (.not.ga_destroy(g_cdinv))
$ call errquit("tddft_fock_fit: destroy GA",g_cdinv,GA_ERR)
endif
! Obtain maximum values for all Shells
if(.not.MA_Push_Get(MT_Dbl,nshells_cd,'CD shell max',
& lcd_shell,icd_shell))
& call errquit('tddft_fock_fit: cannot allocate CD shell max',0,0)
c
CD_coef_max = 0.D0
i1 = 1
do ishp = 1,nshells_cd
if( .not. bas_cn2bfr(CD_bas_han,ishp,ishbf,lshbf) )then
call errquit('tddft_fock_fit',4,0)
end if
shmax = 0.D0
do i = ishbf,lshbf
do idens = 1, ndens
shmax = max(abs(dbl_mb(icdc+(idens-1)*nbf_cd+i-1)), shmax)
enddo
enddo
Dbl_MB((icd_shell-1)+ishp) = shmax
CD_coef_max = max(shmax, CD_coef_max)
enddo
! Part 3. Obtain Fock matrices
icount = 0
n1 = nxtask(nproc,1)
do atom_c = 1,ncenters
if (.not. bas_ce2cnr( AO_bas_han, atom_c, sh_lo_c, sh_hi_c))
& call errquit('Exiting in dft_fitcd',110, BASIS_ERR)
if (.not. bas_ce2bfr( AO_bas_han, atom_c, bf_lo_c, bf_hi_c))
& call errquit('Exiting in dft_fitcd',110, BASIS_ERR)
nbf_at_c=bf_hi_c - bf_lo_c + 1
c
do atom_d= 1,atom_c
dowork=n1.eq.icount
if(dowork) then
v_nonzero=.false.
if (.not.bas_ce2bfr(AO_bas_han,atom_d,bf_lo_d,bf_hi_d))
& call errquit('Exiting in dft_fitcd',110, BASIS_ERR)
nbf_at_d=bf_hi_d - bf_lo_d + 1
if (.not.bas_ce2cnr(AO_bas_han,atom_d,sh_lo_d,sh_hi_d))
& call errquit('Exiting in dft_fitcd',110, BASIS_ERR)
icount_c=0
czero k_p
call dcopy(nbf_at_d*nbf_at_c*ndens,0d0,0,dbl_mb(k_p),1)
maxf = 0d0
do ishc=sh_lo_c,sh_hi_c
if( .not. bas_cn2bfr( AO_bas_han,ishc,ifirstc,ilastc))
& call errquit('Exiting in fitvcoul.',3,0)
nshbfc=ilastc-ifirstc+1
ishd_0 = sh_lo_d
ishd_1 = sh_hi_d
if (atom_c.eq.atom_d) ishd_1=ishc
icount_d=0
do ishd = ishd_0, ishd_1
if (.not. bas_cn2bfr( AO_bas_han, ishd, ifirstd,
I ilastd))
& call errquit('Exiting in fitvcoul.',4,0)
nshbfd=ilastd-ifirstd+1
Nao2 = nshbfc*nshbfd
c
ERI_est=schwarz_shell(ishc,ishd)
if( ERI_est*CD_coef_max.lt.THRESHOLD )then
icount_d = icount_d + nshbfd
cycle
endif
v_nonzero=.true.
qscale=1d0
if(ishc.ne.ishd) qscale=2d0
istart=ieri
nfunz=0
call dfill(nao2*ndens,0d0,dbl_mb(igvec),1)
do ishp = 1,nshells_cd
c
if( .not.bas_cn2bfr(CD_bas_han, ishp,
I ifirstp, ilastp))
& call errquit('fitvc',0,BASIS_ERR)
nshbf_cd=ilastp-ifirstp+1
Vcoul_est=dbl_mb(icd_shell+ishp-1)*ERI_est
nintegrals = nao2*nshbf_cd
if( Vcoul_est.gt.THRESHOLD )then
Ninteg_api=maxg
call int_2e3c(CD_bas_han,ishp,ao_bas_han,ishc,
A ishd,nscr,dbl_mb(iscr),ninteg_api,
& Dbl_MB(istart))
call dgemm('N','N',nao2,ndens,
. nshbf_cd,qscale,dbl_mb(istart),nao2,
. dbl_mb(icdc+nfunz),nbf_cd,1d0,dbl_mb(igvec),
. nao2)
endif ! schw on VCoul_est
c
istart = istart + nintegrals
nfunz=nfunz+nshbf_cd
c
enddo ! ishp
c
call tddft_copyto(ndens,nbf_at_d,nbf_at_c,
. nshbfd,
. nshbfc,dbl_mb(k_p),dbl_mb(igvec),icount_d,
. icount_c)
icount_d=icount_d+nshbfd
enddo ! ishd
icount_c=icount_c+nshbfc
enddo ! ishc
if(v_nonzero) then
lo(2:3) = (/ bf_lo_d, bf_lo_c /)
hi(2:3) = (/ bf_hi_d, bf_hi_c /)
ld(2) = nbf_at_d
call nga_put(g_fock,lo,hi,dbl_mb(k_p),ld)
endif
n1 = nxtask(nproc,1)
endif
icount = icount + 1
enddo ! atom_d
enddo ! atom_c
c
! Chop stack
if(.not.ma_chop_stack(ltvec))
& call errquit('dft_fitcd: cannot pop stack',12, MA_ERR)
n1 = nxtask(-nproc,1)
call ga_sync()
100 continue
if (xc_gotxc().and.nder.eq.0) then
call xc_newfock(geom,ao_bas_han,jfac,jfac,nbf_ao,
& nder, ndens, nfocks, g_fock,g_dens)
endif
call ga_mask_sync(.false.,.false.)
call ga_symmetrize(g_fock)
end
subroutine tddft_finalize_fit()
implicit none
#include "cdft.fh"
#include "errquit.fh"
#include "bas.fh"
#include "global.fh"
if (cd_bas_han.eq.-99999) return
if (.not.bas_destroy(cd_bas_han)) then
call errquit("tddft: cannot destroy cdba",cd_bas_han,BASIS_ERR)
endif
if (noio.eq.1) then
if (.not.ga_destroy(g_2ceri) .or. .not.ga_destroy(g_cdinv))
$ call errquit("tddft: cannot destroy GAs",0,GA_ERR)
endif
end
subroutine tddft_copyfrom(nmat,ld1,ld2,row_sub,col_sub,mat,matsub,
$ row_ptr,col_ptr,fac)
implicit none
integer nmat, ld1, ld2, row_sub, col_sub
double precision mat(nmat,ld1,ld2), matsub(row_sub,col_sub,nmat)
integer row_ptr, col_ptr
double precision fac(*)
integer imat,j,i
do imat=1,nmat
if (fac(imat).eq.0d0) then
matsub(:,:,imat) = 0d0
cycle
endif
do j=1,col_sub
do i=1,row_sub
matsub(i,j,imat) = mat(imat,row_ptr+i,col_ptr+j)
enddo
enddo
enddo
end
subroutine tddft_copyto(nmat,ld1,ld2,row_sub,col_sub,mat,matsub,
$ row_ptr,col_ptr)
implicit none
integer nmat, ld1, ld2, row_sub, col_sub
c double precision mat(ld1,ld2,nmat), matsub(row_sub,col_sub,nmat)
double precision mat(nmat,ld1,ld2), matsub(row_sub,col_sub,nmat)
integer row_ptr, col_ptr
integer imat,j,i
c do imat=1,nmat
do j=1,col_sub
do i=1,row_sub
do imat=1,nmat
c mat(row_ptr+i,col_ptr+j,imat) = matsub(i,j,imat)
mat(imat,row_ptr+i,col_ptr+j) = matsub(i,j,imat)
enddo
enddo
enddo
c enddo
end

19
src/nwdft/lr_tddft/tddft_init.F
View file

@ -5,7 +5,8 @@
4 tda,singlet,triplet,nroots,maxvecs,tol2e,thresh,
5 maxiter,target,targetsym,symmetry,algorithm,g_movecs,
6 l_evl,k_evl,l_occ,k_occ,l_irs,k_irs,mult,lecut,ecut,diff_max,
7 lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
7 lowin,owstart,owend,lewin,ewinl,ewinh,vshift,
8 lhascd)
c
c $Id$
c
@ -81,6 +82,7 @@ c
logical lb94 ! LB94 asymptotic correction
logical cs00 ! CS00 asymptotic correction
logical ncap ! NCAP derivative discontinuity
logical lhascd ! Has CD basis available
double precision shift ! Shift for CS00 or NCAP
C character*80 vector ! CI vector filename
c
@ -156,14 +158,16 @@ c
c ---------
c Basis set
c ---------
if (.not.bas_create(ao_bas_han,'ao basis'))
1 call errquit('tddft_init: bas_create failed',0, BASIS_ERR)
if (.not.bas_rtdb_load(rtdb,geom,ao_bas_han,'ao basis'))
1 call errquit('tddft_init: no ao basis found',0, BASIS_ERR)
if (.not.int_normalize(rtdb,ao_bas_han))
1 call errquit('tddft_init: int_normalize failed',0, BASIS_ERR)
c if (.not.bas_create(ao_bas_han,'ao basis'))
c 1 call errquit('tddft_init: bas_create failed',0, BASIS_ERR)
c if (.not.bas_rtdb_load(rtdb,geom,ao_bas_han,'ao basis'))
c 1 call errquit('tddft_init: no ao basis found',0, BASIS_ERR)
c if (.not.int_normalize(rtdb,ao_bas_han))
c 1 call errquit('tddft_init: int_normalize failed',0, BASIS_ERR)
call tddft_basis(rtdb,geom,ao_bas_han,lhascd)
if (.not.bas_numbf(ao_bas_han,nbf_ao))
1 call errquit('tddft_init: failed to read nbf_ao',0, BASIS_ERR)
c ----------
c MO vectors
c ----------
@ -390,7 +394,6 @@ c ----------------------
c ----------------------
c Initialize integrators
c ----------------------
call int_init(rtdb,1,ao_bas_han)
call schwarz_init(geom,ao_bas_han)
call scf_get_fock_param(rtdb,tol2e)
c --------------------------

11
src/nwdft/lr_tddft/tddft_lanczos.F
View file

@ -48,6 +48,7 @@ c
logical tda ! True if Tamm-Dancoff approximation
logical singlet ! True if singlet excited states
logical triplet ! True if triplet excited states
logical lhascd ! True if CD basis available
integer nroots ! Number of roots sought
integer maxvecs ! Maximum number of trial vectors
integer maxiter ! Maximum number of Davidson iterations
@ -134,7 +135,7 @@ c
4 tda,singlet,triplet,nroots,maxvecs,tol2e,thresh,
5 maxiter,target,targetsym,symmetry,algorithm,g_movecs,
6 l_evl,k_evl,l_occ,k_occ,l_irs,k_irs,mult,lecut,ecut,diff_max,
7 lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
7 lowin,owstart,owend,lewin,ewinl,ewinh,vshift,lhascd)
if (lhasdft) then
if (.not.rtdb_get(rtdb,'bgj:xc_active',mt_log,1,
@ -366,7 +367,7 @@ c********************************
1 nocc,nmo,nfc,nfv,nroots,maxvecs,g_trials,g_prod_p,g_prod_m,
3 g_movecs,dbl_mb(k_evl),1,nproducts,
4 tda,oskel,singlet,triplet,tol2e,xfac(1),lhashf,thresh,
5 target,targetsym,algorithm,mult,l_irs,k_irs,1)
5 target,targetsym,algorithm,mult,l_irs,k_irs,1,lhascd)
c
c keep result in memory
c
@ -400,7 +401,7 @@ c
1 nocc,nmo,nfc,nfv,nroots,maxvecs,g_prod_p,g_prod_p,g_prod_m,
3 g_movecs,dbl_mb(k_evl),1,nproducts,
4 tda,oskel,singlet,triplet,tol2e,xfac(1),lhashf,thresh,
5 target,targetsym,algorithm,mult,l_irs,k_irs,2)
5 target,targetsym,algorithm,mult,l_irs,k_irs,2,lhascd)
cjb spectral transform
@ -418,13 +419,13 @@ cjb spectral transform
1 nocc,nmo,nfc,nfv,nroots,maxvecs,g_prod_m,g_prod_p,g_prod_m,
3 g_movecs,dbl_mb(k_evl),1,nproducts,
4 tda,oskel,singlet,triplet,tol2e,xfac(1),lhashf,thresh,
5 target,targetsym,algorithm,mult,l_irs,k_irs,1)
5 target,targetsym,algorithm,mult,l_irs,k_irs,1,lhascd)
call tddft_prod_vec(rtdb,geom,ao_bas_han,nbf_ao,ipol,
1 nocc,nmo,nfc,nfv,nroots,maxvecs,g_prod_p,g_prod_p,g_prod_m,
3 g_movecs,dbl_mb(k_evl),1,nproducts,
4 tda,oskel,singlet,triplet,tol2e,xfac(1),lhashf,thresh,
5 target,targetsym,algorithm,mult,l_irs,k_irs,2)
5 target,targetsym,algorithm,mult,l_irs,k_irs,2,lhascd)
do i = 1,ipol

5
src/nwdft/lr_tddft/tddft_lanczos_tda.F
View file

@ -48,6 +48,7 @@ c
logical tda ! True if Tamm-Dancoff approximation
logical singlet ! True if singlet excited states
logical triplet ! True if triplet excited states
logical lhascd
integer nroots ! Number of roots sought
integer maxvecs ! Maximum number of trial vectors
integer maxiter ! Maximum number of Davidson iterations
@ -134,7 +135,7 @@ c
4 tda,singlet,triplet,nroots,maxvecs,tol2e,thresh,
5 maxiter,target,targetsym,symmetry,algorithm,g_movecs,
6 l_evl,k_evl,l_occ,k_occ,l_irs,k_irs,mult,lecut,ecut,diff_max,
7 lowin,owstart,owend,lewin,ewinl,ewinh,vshift)
7 lowin,owstart,owend,lewin,ewinl,ewinh,vshift,lhascd)
if (lhasdft) then
if (.not.rtdb_get(rtdb,'bgj:xc_active',mt_log,1,
@ -400,7 +401,7 @@ c 1 g_bigvp(i),1,nov(i),ideg,ideg)
1 nocc,nmo,nfc,nfv,nroots,maxvecs,g_trials,g_prod_p,g_prod_m,
3 g_movecs,dbl_mb(k_evl),1,nproducts,
4 tda,oskel,singlet,triplet,tol2e,xfac(1),lhashf,thresh,
5 target,targetsym,algorithm,mult,l_irs,k_irs,1)
5 target,targetsym,algorithm,mult,l_irs,k_irs,1,lhascd)
cc call ga_print(g_prod_p)

5
src/nwdft/lr_tddft/tddft_prod_v.F
View file

@ -2,7 +2,7 @@
1 nocc,nmo,nfc,nfv,nroots,maxvecs,g_trials,g_prod_p,g_prod_m,
3 g_movecs,evl,ntrials,nproducts,
4 tda,oskel,singlet,triplet,tol2e,kfac,lhashf,thresh,
5 target,targetsym,algorithm,mult,l_irs,k_irs,iwhich)
5 target,targetsym,algorithm,mult,l_irs,k_irs,iwhich,lhascd)
c
c
c
@ -47,6 +47,7 @@ c
logical tda ! True if Tamm-Dancoff approximation
logical singlet ! True if singlet excited states
logical triplet ! True if triplet excited states
logical lhascd
logical oskel ! True if use symmetry in fock_2e
double precision tol2e ! 2-e integral tolerance
double precision kfac ! Prefactor for HF exchange contribution
@ -283,7 +284,7 @@ c if(nodezero)write(6,*)'Before contractions'
call tddft_nga_cont(rtdb,geom,ao_bas_han,
1 g_ntrl,g_nprd_p,g_nprd_m,nbf_ao,ipol,tol2e,tda,oskel,kfac,
2 lhashf,triplet,nvectors,iwhich)
2 lhashf,triplet,nvectors,iwhich,lhascd)
ccpu=ccpu+util_cpusec()
cwall=cwall+util_wallsec()

30
src/nwdft/lr_tddft/tddft_transfm.F
View file

@ -60,15 +60,15 @@ c I choose to carry out full matrix multiplication in the space
c of nbf_ao*nbf_ao, while leaving the meaningless part of the
c space zeroed.
c
if (.not.ma_push_get(mt_dbl,nbf_ao**2,'Work array',
1 l_w1,k_w1)) call errquit
2 ('tddft_transfm: failed to allocate w1',0, MA_ERR)
if (.not.ma_push_get(mt_dbl,nbf_ao**2,'Work array',
1 l_w2,k_w2)) call errquit
2 ('tddft_transfm: failed to allocate w2',0, MA_ERR)
if (.not.ma_push_get(mt_dbl,nbf_ao**2,'Work MO vectors',
1 l_w3,k_w3)) call errquit
2 ('tddft_transfm: failed to allocate w2',0, MA_ERR)
c if (.not.ma_push_get(mt_dbl,nbf_ao**2,'Work array',
c 1 l_w1,k_w1)) call errquit
c 2 ('tddft_transfm: failed to allocate w1',0, MA_ERR)
c if (.not.ma_push_get(mt_dbl,nbf_ao**2,'Work array',
c 1 l_w2,k_w2)) call errquit
c 2 ('tddft_transfm: failed to allocate w2',0, MA_ERR)
c if (.not.ma_push_get(mt_dbl,nbf_ao**2,'Work MO vectors',
c 1 l_w3,k_w3)) call errquit
c 2 ('tddft_transfm: failed to allocate w2',0, MA_ERR)
c
c ------
c Type 1
@ -255,12 +255,12 @@ c ------
c Return
c ------
c
if (.not.ma_pop_stack(l_w3)) call errquit
1 ('tddft_transfm: failed to deallocate w3',0, MA_ERR)
if (.not.ma_pop_stack(l_w2)) call errquit
1 ('tddft_transfm: failed to deallocate w2',0, MA_ERR)
if (.not.ma_pop_stack(l_w1)) call errquit
1 ('tddft_transfm: failed to deallocate w1',0, MA_ERR)
c if (.not.ma_pop_stack(l_w3)) call errquit
c 1 ('tddft_transfm: failed to deallocate w3',0, MA_ERR)
c if (.not.ma_pop_stack(l_w2)) call errquit
c 1 ('tddft_transfm: failed to deallocate w2',0, MA_ERR)
c if (.not.ma_pop_stack(l_w1)) call errquit
c 1 ('tddft_transfm: failed to deallocate w1',0, MA_ERR)
return
9000 format(' Vector: ',i5,' Spin: ',i5)
end

9
src/nwdft/lr_tddft_grad/tddft_grad_compute_energy.F
View file

@ -4,7 +4,7 @@ c HvD 9/2006, NG 11/2012, DS 1/2013
c
subroutine tddft_grad_compute_energy(rtdb,ihdl_geom,ihdl_bfao,
+ tol2e,tda,ipol,nao,nfc,naoc,nocc,nav,nfv,
+ nroot,o_eig,g_mo,g_xpy,g_xmy,kfac,lhashf,otriplet)
+ nroot,o_eig,g_mo,g_xpy,g_xmy,kfac,lhashf,otriplet,lhascd)
c
c This subroutine simply evaluates the TDDFT energy using the
c information read in. The energy expression evaluated is
@ -59,6 +59,7 @@ c
logical otriplet ! =.true. triplet excited states
! =.false. singlet excited states
! value does not matter for TDUDFT
logical lhascd
c
double precision kfac ! the weight of the Hartree-Fock
! exchange contributions
@ -174,9 +175,11 @@ c
endif
c
call tddft_nga_cont(rtdb,ihdl_geom,ihdl_bfao,g_x,g_apbx,g_ambx,
+ nao,ipol,tol2e,tda,oskel,kfac,lhashf,otriplet,nroot,iwhich)
+ nao,ipol,tol2e,tda,oskel,kfac,lhashf,otriplet,nroot,iwhich,
+ lhascd)
call tddft_nga_cont(rtdb,ihdl_geom,ihdl_bfao,g_y,g_apby,g_amby,
+ nao,ipol,tol2e,tda,oskel,kfac,lhashf,otriplet,nroot,iwhich)
+ nao,ipol,tol2e,tda,oskel,kfac,lhashf,otriplet,nroot,iwhich,
+ lhascd)
c
c Dispose of X and Y
c

9
src/nwdft/lr_tddft_grad/tddft_grad_compute_r.F
View file

@ -65,6 +65,7 @@ c
logical otriplet ! =.true. triplet excited states
! =.false. singlet excited states
! value does not matter for TDUDFT
logical lhascd
c
double precision kfac ! the weight of the Hartree-Fock exchange
! contributions
@ -178,7 +179,7 @@ c fock_xc:triplet here for restricted triplet calculations to work.
endif
c
call tddft_nga_cont(rtdb,ihdl_geom,ihdl_bfao,g_tuv,g_apbt,g_ambt,
+ nao,ipol,tol2e,tda,oskel,kfac,lhashf,.false.,nroot,iwhich)
+nao,ipol,tol2e,tda,oskel,kfac,lhashf,.false.,nroot,iwhich,lhascd)
c Daniel (2-26-13): Reset fock_xc:triplet here for restricted triplet
c calculations to work.
if (otriplet) then
@ -304,9 +305,11 @@ c the global one.
endif
endif
call tddft_nga_cont(rtdb,ihdl_geom,ihdl_bfao,g_x,g_apbx,g_ambx,
+ nao,ipol,tol2e,tdaloc,oskel,kfac,lhashf,otriplet,nroot,iwhich)
+nao,ipol,tol2e,tdaloc,oskel,kfac,lhashf,otriplet,nroot,iwhich,
+lhascd)
call tddft_nga_cont(rtdb,ihdl_geom,ihdl_bfao,g_y,g_apby,g_amby,
+ nao,ipol,tol2e,tdaloc,oskel,kfac,lhashf,otriplet,nroot,iwhich)
+ nao,ipol,tol2e,tdaloc,oskel,kfac,lhashf,otriplet,nroot,iwhich,
+ lhascd)
c
c 2.6 Dispose of X and Y
c

10
src/nwdft/lr_tddft_grad/tddft_grad_compute_w.F
View file

@ -43,6 +43,7 @@ c
logical otriplet ! =.true. triplet excited states
! =.false. singlet excited states
! value does not matter for TDUDFT
logical lhascd
c
double precision kfac ! the weight of the Hartree-Fock
! exchange contributions
@ -472,7 +473,8 @@ c fock_xc:triplet here for restricted triplet calculations to work.
1 call errquit(pname//'failed to set triplet',0,RTDB_ERR)
endif
call tddft_nga_cont(rtdb,ihdl_geom,ihdl_bfao,g_puv,g_apbp,g_ambp,
+ nao,ipol,tol2e,tda,oskel,kfac,lhashf,.false.,nroots,iwhich)
+ nao,ipol,tol2e,tda,oskel,kfac,lhashf,.false.,nroots,iwhich,
+ lhascd)
c Daniel (2-26-13): Reset fock_xc:triplet here for restricted triplet
c calculations to work.
if (otriplet) then
@ -768,9 +770,11 @@ c the global one.
endif
endif
call tddft_nga_cont(rtdb,ihdl_geom,ihdl_bfao,g_x,g_apbx,g_ambx,
+ nao,ipol,tol2e,tdaloc,oskel,kfac,lhashf,otriplet,nroots,iwhich)
+ nao,ipol,tol2e,tdaloc,oskel,kfac,lhashf,otriplet,nroots,iwhich,
+ lhascd)
call tddft_nga_cont(rtdb,ihdl_geom,ihdl_bfao,g_y,g_apby,g_amby,
+ nao,ipol,tol2e,tdaloc,oskel,kfac,lhashf,otriplet,nroots,iwhich)
+ nao,ipol,tol2e,tdaloc,oskel,kfac,lhashf,otriplet,nroots,iwhich,
+ lhascd)
c
c - Dispose of X and Y
c

2
src/nwdft/lr_tddft_grad/tddft_grad_finalize.F
View file

@ -42,6 +42,8 @@ c
+ call errquit(pname//'failed to destroy geometry')
if (.not.bas_destroy(ihdl_bfao))
+ call errquit(pname//'failed to destroy basis set')
c
call tddft_finalize_fit()
c
end
c $Id$

16
src/nwdft/lr_tddft_grad/tddft_grad_init.F
View file

@ -11,7 +11,7 @@ c
+ ntroots,iptr_troots,ihdl_troots,fn_civecs_triplet,
+ osinglet,otriplet,ntotroots,ihdl_bfao,ihdl_geom,tol2e,
+ solve_thresh,solve_maxit,solve_maxsub,
+ fn_solve,solve_when,kfac,title)
+ fn_solve,solve_when,kfac,title,lhascd)
c
implicit none
c
@ -71,6 +71,7 @@ c
logical lhybrid ! True if hybrid DFT and HF
logical lhasdft ! True if it has DFT
logical lhashf ! True if it has HF
logical lhascd
logical use_nwxc ! True if using NWXC library
c
logical cs00
@ -160,12 +161,13 @@ c get the geometry
+ call errquit(pname//'geom_ncent failed',0,GEOM_ERR)
c
c get the basis set
if (.not.bas_create(ihdl_bfao,'ao basis'))
+ call errquit(pname//'bas_create failed',0, BASIS_ERR)
if (.not.bas_rtdb_load(rtdb,ihdl_geom,ihdl_bfao,'ao basis'))
+ call errquit(pname//'bas_rtdb_load failed',0, BASIS_ERR)
if (.not.int_normalize(rtdb,ihdl_bfao))
+ call errquit(pname//'int_normalize failed',0, BASIS_ERR)
call tddft_basis(rtdb, ihdl_geom, ihdl_bfao, lhascd)
c if (.not.bas_create(ihdl_bfao,'ao basis'))
c + call errquit(pname//'bas_create failed',0, BASIS_ERR)
c if (.not.bas_rtdb_load(rtdb,ihdl_geom,ihdl_bfao,'ao basis'))
c + call errquit(pname//'bas_rtdb_load failed',0, BASIS_ERR)
c if (.not.int_normalize(rtdb,ihdl_bfao))
c + call errquit(pname//'int_normalize failed',0, BASIS_ERR)
if (.not.bas_numbf(ihdl_bfao,nao))
+ call errquit(pname//'bas_numbf failed',0, BASIS_ERR)
c

2
src/nwdft/lr_tddft_grad/tddft_grad_ints.F
View file

@ -33,7 +33,7 @@ c Code
+ calc_type))
+ call errquit(pname//'failed to set calc_type_save',0,RTDB_ERR)
c
call int_init(rtdb,1,ihdl_bfao)
c call int_init(rtdb,1,ihdl_bfao)
call schwarz_init(ihdl_geom,ihdl_bfao)
call scf_get_fock_param(rtdb,tol2e)
c

4
src/nwdft/lr_tddft_grad/tddft_grad_solve.F
View file

@ -5,7 +5,7 @@ c
subroutine tddft_grad_solve(rtdb,ihdl_geom,ihdl_bfao,ipol,nroots,
+ nao,nocc,naoc,nav,nfc,nfv,g_mo,iptr_mo_e,tda,oskel,
+ kfac,lhashf,otriplet,solve_thresh,solve_maxsub,
+ solve_maxit,fn_solve,solve_when,tol2e,g_r,g_z)
+ solve_maxit,fn_solve,solve_when,tol2e,g_r,g_z,lhascd)
c
c At this point Ria has just been computed. The quantities
c are stored in 3D global arrays as R(1)(1:nroot,1:nocca,1:nvirta),
@ -70,6 +70,7 @@ c
logical oskel ! Use of symmetry in Fock build
logical lhashf ! has the functional got HF exchange
logical otriplet ! are we doing triplet excitations
logical lhascd
integer nroots ! the number of roots
integer nao ! the number of AOs
integer nfc(2) ! the number of frozen core orbitals
@ -182,6 +183,7 @@ c tdg_tda = tda
tdg_lhashf = lhashf
cold tdg_triplet = otriplet
tdg_triplet = .false.
tdg_lhascd = lhascd
c Daniel (2-26-13): It was not obvious that we need to unset
c fock_xc:triplet here for restricted triplet calculations to work.
if (otriplet) then

3
src/nwdft/lr_tddft_grad/tddft_grad_solve_common.fh
View file

@ -28,8 +28,9 @@ c
logical tdg_oskel
logical tdg_lhashf
logical tdg_triplet
logical tdg_lhascd
common/tddft_grad_solve_common/tdg_tol2e,tdg_kfac,tdg_rtdb,
+ tdg_geom,tdg_bfao,tdg_ipol,tdg_nroots,tdg_nao,tdg_nocc,
+ tdg_naoc,tdg_nav,tdg_nfc,tdg_nfv,tdg_nmo,
+ tdg_g_mo,tdg_e_mo,tdg_tda,
+ tdg_oskel,tdg_lhashf,tdg_triplet
+ tdg_oskel,tdg_lhashf,tdg_triplet,tdg_lhascd

3
src/nwdft/lr_tddft_grad/tddft_grad_solve_product.F
View file

@ -106,7 +106,8 @@ c
call ga_zero(g_Ax_3d_ao)
call tddft_nga_cont(tdg_rtdb,tdg_geom,tdg_bfao,g_x_3d_ao,
+ g_Ax_3d_ao,g_ambt,tdg_nao,tdg_ipol,tdg_tol2e,tdg_tda,
+ tdg_oskel,tdg_kfac,tdg_lhashf,tdg_triplet,tdg_nroots,iwhich)
+ tdg_oskel,tdg_kfac,tdg_lhashf,tdg_triplet,tdg_nroots,iwhich,
+ tdg_lhascd)
if (.not.tdg_tda) then
if (.not.ga_destroy(g_ambt))
+ call errquit(pname//'failed to destroy g_ambt',0,GA_ERR)

11
src/nwdft/lr_tddft_grad/tddft_gradient.F
View file

@ -41,6 +41,7 @@ c
logical oskel ! Use of symmetry in Fock build
logical osinglet ! Compute singlet excited states?
logical otriplet ! Compute triplet excited states?
logical lhascd
integer nsroots ! the number of singlet roots
integer ihdl_sroots ! the singlet roots handle
integer iptr_sroots ! the singlet roots list
@ -181,7 +182,7 @@ c Parameter initialization
+ ntroots,iptr_troots,ihdl_troots,fn_civecs_triplet,
+ osinglet,otriplet, ntotroots,ihdl_bfao,ihdl_geom,tol2e,
+ solve_thresh,solve_maxit,solve_maxsub,
+ fn_solve,solve_when,kfac,title)
+ fn_solve,solve_when,kfac,title,lhascd)
c
c Cosmo initialization, if necessary
cosmo_on = .false.
@ -277,7 +278,7 @@ c Compute the TDDFT energy for each root using the solution vectors
& write(luout,*)"call tddft_grad_compute_energy"
call tddft_grad_compute_energy(rtdb,ihdl_geom,ihdl_bfao,tol2e,tda,
+ ipol,nao,nfc,naoc,nocc,nav,nfv,ntotroots,
+ dbl_mb(iptr_mo_e),g_mo,g_xpy,g_xmy,kfac,lhashf,otriplet)
+ dbl_mb(iptr_mo_e),g_mo,g_xpy,g_xmy,kfac,lhashf,otriplet,lhascd)
if (oroot.and.odebug)
& write(luout,*)"done tddft_grad_compute_energy"
c
@ -297,7 +298,7 @@ c Assemble right hand side: R
if (oroot.and.odebug) write(luout,*) "call tddft_grad_compute_r"
call tddft_grad_compute_r(rtdb,ihdl_geom,ihdl_bfao,tol2e,tda,ipol,
+ nao,nfc,naoc,nocc,nav,nfv,ntotroots,g_mo,g_xpy,g_xmy,g_tp,
+ g_r,kfac,lhashf,otriplet)
+ g_r,kfac,lhashf,otriplet,lhascd)
if (oroot.and.odebug) write(luout,*) "done tddft_grad_compute_r"
c
c Solve the Z equation: (A+B)Z = -R
@ -306,7 +307,7 @@ c Note: g_z has the solutions
call tddft_grad_solve(rtdb,ihdl_geom,ihdl_bfao,ipol,ntotroots,
+ nao,nocc,naoc,nav,nfc,nfv,g_mo,iptr_mo_e,tda,oskel,kfac,
+ lhashf,otriplet,solve_thresh,solve_maxsub,solve_maxit,
+ fn_solve,solve_when,tol2e,g_r,g_z)
+ fn_solve,solve_when,tol2e,g_r,g_z,lhascd)
if (oroot.and.odebug) write(luout,*) "done tddft_grad_solve"
c
c Deallocate right hand side: R
@ -330,7 +331,7 @@ c W is the energy-weighted difference density matrix
call tddft_grad_compute_w(rtdb,ihdl_geom,ihdl_bfao,tol2e,tda,ipol,
+ ntotroots,nfc,naoc,nocc,nav,nfv,nao,g_mo,g_tp,g_z,
+ g_xpy,g_xmy,dbl_mb(iptr_mo_e),dbl_mb(iptr_root_e),g_w,
+ kfac,lhashf,otriplet)
+ kfac,lhashf,otriplet,lhascd)
if (oroot.and.odebug) write(luout,*) "done tddft_grad_compute_w"
c
c Finalize integrals

2
src/nwdft/scf_dft/dft_pstat.F
View file

@ -66,6 +66,8 @@ c
$ ps_scfend)) call errquit('dft: dft_pstat_init', 0,0)
if (.not. pstat_allocate('dft:scf', pstat_qstat, 0, junk,
$ ps_scfdft)) call errquit('dft: dft_pstat_init', 0,0)
if (.not. pstat_allocate('dft:cpks_coef', pstat_qstat, 0, junk,
$ ps_cpkscoef)) call errquit('dft: dft_pstat_init', 0,0)
if (.not. pstat_allocate('dft:total', pstat_qstat, 0, junk,
$ ps_dft)) call errquit('dft: dft_pstat_init', 0,0)
c

6
src/nwdft/scf_dft/dftps.fh
View file

@ -7,7 +7,8 @@ c
, ps_scfend,ps_incore,ps_bld12,ps_ssw,ps_xcrho,
, ps_tabcd,ps_ebf,ps_excf,ps_bdgemm,ps_xmrge,
, ps_xcopy,ps_inv,ps_2cer,ps_qv0b,ps_f2e,
, ps_chol,ps_getvxc
, ps_chol,ps_getvxc,
, ps_cpkscoef
logical odftps ! True if gathering stats
common /cdftps/ odftps,
@ -16,6 +17,7 @@ c
, ps_scfend,ps_incore,ps_bld12,ps_ssw,ps_xcrho,
, ps_tabcd,ps_ebf,ps_excf,ps_bdgemm,ps_xmrge,
, ps_xcopy,ps_inv,ps_2cer,ps_qv0b,ps_f2e,
, ps_chol,ps_getvxc
, ps_chol,ps_getvxc,
, ps_cpkscoef
c

104
src/nwdft/xc/xc_tabcd.F
View file

@ -180,6 +180,7 @@ c
double precision chi1
double precision scr(nq) !< [Scratch] Temporary stuff in
!< quadrature
double precision scr2(nq)
double precision dabsmax
double precision tolrho15
external dabsmax
@ -205,8 +206,8 @@ c
data indh / 1, 2, 3,
& 2, 4, 5,
& 3, 5, 6 /
data nbhandl1 /0./
data nbhandl2 /0./
data nbhandl1 /0/
data nbhandl2 /0/
save nbhandl1
save nbhandl2
double precision ddot
@ -295,53 +296,40 @@ c
nnia = ifinia - inizia + 1
nbfia = ilast - ifirst + 1
do ii = 1, ipol
call starttimer(monitor_mult1)
do mu = 1, nnia
call starttimer(monitor_mult1)
mu1 = mu+inizia-1
Bmat(1:nq,mu) = Amat(1:nq,ii,imat)*chi(1:nq,mu1)
if (GRAD) then
do n = 1, nq
chi1 = chi(n,mu1)
Bmat(n,mu) = Amat(n,ii,imat)*chi1 +
& delchi(n,1,mu1)*Cmat(n,1,ii,imat) +
& delchi(n,2,mu1)*Cmat(n,2,ii,imat) +
& delchi(n,3,mu1)*Cmat(n,3,ii,imat)
Dmat(n,1,mu) = Cmat(n,1,ii,imat)*chi1
Dmat(n,2,mu) = Cmat(n,2,ii,imat)*chi1
Dmat(n,3,mu) = Cmat(n,3,ii,imat)*chi1
if(kske) then
Dmat(n,1,mu) = Dmat(n,1,mu) +
& Mmat(n,ii,imat)*delchi(n,1,mu1)
Dmat(n,2,mu) = Dmat(n,2,mu) +
& Mmat(n,ii,imat)*delchi(n,2,mu1)
Dmat(n,3,mu) = Dmat(n,3,mu) +
& Mmat(n,ii,imat)*delchi(n,3,mu1)
endif
if(kslap) then
Dmat(n,1,mu) = Dmat(n,1,mu) +
& 2d0*Lmat(n,ii,imat)*delchi(n,1,mu1)
Dmat(n,2,mu) = Dmat(n,2,mu) +
& 2d0*Lmat(n,ii,imat)*delchi(n,2,mu1)
Dmat(n,3,mu) = Dmat(n,3,mu) +
& 2d0*Lmat(n,ii,imat)*delchi(n,3,mu1)
Bmat(n,mu) = Bmat(n,mu) + Lmat(n,ii,imat)*
& (heschi(n,1,mu1)+heschi(n,4,mu1)+heschi(n,6,mu1))
Gmat(n,mu) = Lmat(n,ii,imat)*chi1
endif
enddo
else
do n = 1, nq
Bmat(n,mu) = chi(n,mu1)*Amat(n,ii,imat)
enddo
endif ! GRAD
enddo ! mu
Bmat(1:nq,mu) = Bmat(1:nq,mu) +
& delchi(1:nq,1,mu1)*Cmat(1:nq,1,ii,imat) +
& delchi(1:nq,2,mu1)*Cmat(1:nq,2,ii,imat) +
& delchi(1:nq,3,mu1)*Cmat(1:nq,3,ii,imat)
Dmat(1:nq,1,mu) = chi(1:nq,mu1)*Cmat(1:nq,1,ii,imat)
Dmat(1:nq,2,mu) = chi(1:nq,mu1)*Cmat(1:nq,2,ii,imat)
Dmat(1:nq,3,mu) = chi(1:nq,mu1)*Cmat(1:nq,3,ii,imat)
endif
if (kske.and.grad) then
Dmat(1:nq,1,mu) = Dmat(1:nq,1,mu) +
& Mmat(1:nq,ii,imat)*delchi(1:nq,1,mu1)
Dmat(1:nq,2,mu) = Dmat(1:nq,2,mu) +
& Mmat(1:nq,ii,imat)*delchi(1:nq,2,mu1)
Dmat(1:nq,3,mu) = Dmat(1:nq,3,mu) +
& Mmat(1:nq,ii,imat)*delchi(1:nq,3,mu1)
endif
if (kslap.and.grad) then
Dmat(1:nq,1,mu) = Dmat(1:nq,1,mu) +
& 2d0*Lmat(1:nq,ii,imat)*delchi(1:nq,1,mu1)
Dmat(1:nq,2,mu) = Dmat(1:nq,2,mu) +
& 2d0*Lmat(1:nq,ii,imat)*delchi(1:nq,2,mu1)
Dmat(1:nq,3,mu) = Dmat(1:nq,3,mu) +
& 2d0*Lmat(1:nq,ii,imat)*delchi(1:nq,3,mu1)
Bmat(1:nq,mu) = Bmat(1:nq,mu) +
& Lmat(1:nq,ii,imat)*(heschi(1:nq,1,mu1)+
& heschi(1:nq,4,mu1)+heschi(1:nq,6,mu1))
Gmat(1:nq,mu) = Lmat(1:nq,ii,imat)*chi(1:nq,mu1)
endif
enddo
c Monitoring
call endtimer(monitor_mult1)
@ -477,20 +465,16 @@ c leading dimension of delchi must be set correctly
enddo
endif
endif ! GRAD
if(kslap) then
if (w01) then
indT=0
do nu = 1, nnja
nu1 = nu+inizja-1
do mu = 1, nnia
indT = indT + 1
Tmat(indT) = Tmat(indT) +
& ddot(nq,Gmat(1,mu),1,heschi(1,1,nu1),1) +
& ddot(nq,Gmat(1,mu),1,heschi(1,4,nu1),1) +
& ddot(nq,Gmat(1,mu),1,heschi(1,6,nu1),1)
enddo
enddo
endif
if(kslap.and.w01) then
indT=1
do nu = 1, nnja
nu1 = nu+inizja-1
scr2(1:nq) = heschi(1:nq,1,nu1) +
& heschi(1:nq,4,nu1) + heschi(1:nq,6,nu1)
call dgemv('T',nq,nnia,1d0,Gmat,nq,scr2,1,1d0,
& Tmat(indT),1)
indT = indT + nnia
enddo
endif
c Daniel (2-7-13): For the CPKS RHS stuff, the first call has what=1,
c so n3d=1. Also, dofull is true in that case.