mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-28 22:25:48 -04:00
Implementation of spin-polarised ECPs, including 1st and 2nd derivatives wrt coordinates. I'll put together an example and add it to the QA tests before merging my frozen-density embedding implementation
This commit is contained in:
Andrew Logsdail
parent
d2f7f3186d
commit
3b054e9e56
24 changed files with 2180 additions and 192 deletions
|
|
@ -68,7 +68,9 @@ c:tex-
|
|||
c:tex-{\it Syntax:}
|
||||
c:tex-\begin{verbatim}
|
||||
subroutine int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lstv,S,T,V,
|
||||
& doS,doT,doV)
|
||||
cc AJL/Begin
|
||||
& doS,doT,doV,msg)
|
||||
cc AJL/End
|
||||
c:tex-\end{verbatim}
|
||||
implicit none
|
||||
#include "nwc_const.fh"
|
||||
|
|
@ -130,6 +132,10 @@ c
|
|||
integer WarnP
|
||||
save WarnP
|
||||
data WarnP /0/
|
||||
cc AJL/Begin
|
||||
cc Soft-coded the value of msg
|
||||
character*(*) msg
|
||||
cc AJL/End
|
||||
c
|
||||
#include "bas_exndcf_sfn.fh"
|
||||
#include "bas_ibs_sfn.fh"
|
||||
|
|
@ -279,7 +285,10 @@ c
|
|||
& coords(1,1,i_geom),charge(1,i_geom),
|
||||
& geom_invnucexp(1,i_geom),ncenter(i_geom),
|
||||
& S,T,V,lstv,doS,doT,doV,canAB,.false.,
|
||||
& scr,lscr,'int_1eov')
|
||||
cc AJL/Begin
|
||||
c & scr,lscr,'int_1eov')
|
||||
& scr,lscr,msg)
|
||||
cc AJL/End
|
||||
elseif (onw_ok) then
|
||||
call hf1(
|
||||
& coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)),
|
||||
|
|
@ -414,8 +423,11 @@ c:tex-\begin{verbatim}
|
|||
c:tex-\end{verbatim}
|
||||
c
|
||||
call int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lov,Ov,scr,scr,
|
||||
& .true.,.false.,.false.)
|
||||
c doS doT doV
|
||||
cc AJL/Begin
|
||||
c & .true.,.false.,.false.)
|
||||
& .true.,.false.,.false.,'int_1eov')
|
||||
c doS doT doV msg
|
||||
cc AJL/End
|
||||
return
|
||||
end
|
||||
c----------------------------------------------------------------------
|
||||
|
|
@ -455,8 +467,11 @@ c:tex-\begin{verbatim}
|
|||
c:tex-\end{verbatim}
|
||||
c
|
||||
call int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lke,scr,Ke,scr,
|
||||
& .false.,.true.,.false.)
|
||||
c doS doT doV
|
||||
cc AJL/Begin
|
||||
c & .false.,.true.,.false.)
|
||||
& .false.,.true.,.false.,'int_1eke')
|
||||
c doS doT doV msg
|
||||
cc AJL/End
|
||||
return
|
||||
end
|
||||
c----------------------------------------------------------------------
|
||||
|
|
@ -503,11 +518,67 @@ c:tex-\begin{verbatim}
|
|||
c:tex-\end{verbatim}
|
||||
c
|
||||
call int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lpe,scr,scr,Pe,
|
||||
& .false.,.false.,.true.)
|
||||
c doS doT doV
|
||||
cc AJL/Begin
|
||||
c & .false.,.false.,.true.)
|
||||
& .false.,.false.,.true.,'int_1epe')
|
||||
c doS doT doV msg
|
||||
cc AJL/End
|
||||
return
|
||||
end
|
||||
C> @}
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c----------------------------------------------------------------------
|
||||
*
|
||||
C> \brief Computes the nuclear attraction integrals
|
||||
C>
|
||||
C> This routine calculates the 1-electron nuclear attraction integrals
|
||||
C> for a pair of shells. If any ECPs are present then the ECP terms
|
||||
C> are directly summed into these integrals.
|
||||
C>
|
||||
C> The integrals are defined as:
|
||||
C> \f{eqnarray*}{
|
||||
C> V & = &
|
||||
C({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\\\
|
||||
C> & = & \int_{-\infty}^{\infty}
|
||||
Cg_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac
|
||||
C> {-Z_{\alpha}}{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1}
|
||||
C> \f}
|
||||
c:tex-% this is part of the API Standard Integral routines.
|
||||
c:tex-\subsection{int\_1epe}
|
||||
c:tex-This routine computes the 1 electron potential integrals, ($V$):
|
||||
c:tex-\begin{eqnarray*}
|
||||
c:tex-V & = &({\mu}|\sum_{\alpha}\frac{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}|{\nu}) \\
|
||||
c:tex- & = & \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac
|
||||
c:tex-{-Z_{\alpha}}{|r_{1}-R_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1}
|
||||
c:tex-\end{eqnarray*}
|
||||
c:tex-If a *BETA* ECP is defined then the ECP integral contributions are summed
|
||||
c:tex-directly into the V integrals.
|
||||
c:tex-
|
||||
c:tex-{\it Syntax:}
|
||||
c:tex-\begin{verbatim}
|
||||
subroutine int_1epe_beta(i_basis,ish,j_basis,jsh,lscr,scr,lpe,Pe)
|
||||
#include "util.fh"
|
||||
c:tex-\end{verbatim}
|
||||
c
|
||||
c::passed
|
||||
c:tex-\begin{verbatim}
|
||||
integer i_basis !< [Input] basis set handle for ish
|
||||
integer ish !< [Input] i shell/contraction
|
||||
integer j_basis !< [Input] basis set handle for jsh
|
||||
integer jsh !< [Input] j shell/contraction
|
||||
integer lscr !< [Input] length of scratch array
|
||||
double precision scr(lscr) !< [Scratch] scratch array
|
||||
integer lpe !< [Input] length of Pe buffer
|
||||
double precision Pe(lpe) !< [Output] kinetic energy integrals
|
||||
c:tex-\end{verbatim}
|
||||
c
|
||||
call int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lpe,scr,scr,Pe,
|
||||
& .false.,.false.,.true.,'int_1epe_beta')
|
||||
c doS doT doV msg
|
||||
return
|
||||
end
|
||||
cc AJL/End
|
||||
C> @}
|
||||
c----------------------------------------------------------------------
|
||||
*
|
||||
C> \ingroup nwint
|
||||
|
|
@ -573,8 +644,11 @@ c:tex-\begin{verbatim}
|
|||
c:tex-\end{verbatim}
|
||||
c
|
||||
call int_1estv(i_basis,ish,j_basis,jsh,lscr-lh1,scr(lh1+1),lh1,
|
||||
& scr,scr,H1,.false.,.true.,.true.)
|
||||
c doS doT doV
|
||||
cc AJL/Begin
|
||||
c & scr,scr,H1,.false.,.true.,.true.)
|
||||
& scr,scr,H1,.false.,.true.,.true.,'int_1eh1')
|
||||
c doS doT doV msg
|
||||
cc AJL/End
|
||||
do i = 1,lh1
|
||||
H1(i) = H1(i)+scr(i)
|
||||
end do
|
||||
|
|
@ -631,8 +705,11 @@ c:tex-\begin{verbatim}
|
|||
c:tex-\end{verbatim}
|
||||
c
|
||||
call int_1estv(i_basis,ish,j_basis,jsh,lscr,scr,lstv,S,T,V,
|
||||
& .true.,.true.,.true.)
|
||||
c doS doT doV
|
||||
cc AJL/Begin
|
||||
c & .true.,.true.,.true.)
|
||||
& .true.,.true.,.true.,'int_1eall')
|
||||
c doS doT doV msg
|
||||
cc AJL/End
|
||||
return
|
||||
end
|
||||
c----------------------------------------------------------------------
|
||||
|
|
@ -1397,3 +1474,406 @@ c
|
|||
end
|
||||
C>
|
||||
C> @}
|
||||
cc AJL/BEGIN/FDE
|
||||
c----------------------------------------------------------------------
|
||||
C> \brief Compute the 1-electron Bq embedding potential
|
||||
C>
|
||||
C> This routine computes the 1 electron integrals \f$V\f$:
|
||||
C> \f{eqnarray*}{
|
||||
C> V & = &
|
||||
C({\mu}|\sum_{\alpha}\frac{-q_{\alpha}}{|r_{1}-r_{\alpha}|}|{\nu}) \\\\
|
||||
C> & = &
|
||||
C\int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac{-q_{\alpha}}
|
||||
C> {|r_{1}-r_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1}
|
||||
C> \f}
|
||||
C> where \f$ q_\alpha \f$ is a point charge situated at position \f$
|
||||
Cr_\alpha \f$.
|
||||
C> The Bq embedding point charges are accessed in int_1eefc1.
|
||||
C>
|
||||
subroutine int_1efde1(i_basis,ish,j_basis,jsh,fde_basis,lscr,scr,
|
||||
& lpe,Pe)
|
||||
#include "util.fh"
|
||||
c
|
||||
integer i_basis !< [Input] basis set handle for ish
|
||||
integer ish !< [Input] i shell/contraction
|
||||
integer j_basis !< [Input] basis set handle for jsh
|
||||
integer jsh !< [Input] j shell/contraction
|
||||
integer fde_basis !< [Input] FDE basis set handle
|
||||
integer lscr !< [Input] length of scratch array
|
||||
double precision scr(lscr) !< [Scratch] scratch array
|
||||
integer lpe !< [Input] length of Pe buffer
|
||||
double precision Pe(lpe) !< [Output] potential energy integrals
|
||||
c
|
||||
c call int_1epe(i_basis,ish,j_basis,jsh,lscr,scr,lpe,scr,
|
||||
c & scr,Pe,.false.,.false.,.true.,'pot')
|
||||
call int_1efde(i_basis,ish,j_basis,jsh,fde_basis,lscr,scr,lpe,scr,
|
||||
& scr,Pe,.false.,.false.,.true.)
|
||||
c doS doT doV
|
||||
return
|
||||
end
|
||||
c----------------------------------------------------------------------
|
||||
C> \brief Compute the 1-electron integrals but the potential is due to
|
||||
C> Bq embedding charges
|
||||
C>
|
||||
C> This routine computes the 1 electron integrals \f$V\f$:
|
||||
C> \f{eqnarray*}{
|
||||
C> V & = &
|
||||
C({\mu}|\sum_{\alpha}\frac{-q_{\alpha}}{|r_{1}-r_{\alpha}|}|{\nu}) \\\\
|
||||
C> & = &
|
||||
C\int_{-\infty}^{\infty}g_{\mu}(X_{\mu},r_{1})\sum_{\alpha}\frac{-q_{\alpha}}
|
||||
C> {|r_{1}-r_{\alpha}|}g_{\nu}(X_{\nu},r_{1})dr_{1}
|
||||
C> \f}
|
||||
C> where \f$ q_\alpha \f$ is a point charge situated at position \f$
|
||||
Cr_\alpha \f$.
|
||||
C> The location of the Bq embedding point charges is specified in a
|
||||
C> Bq instance that has been activated elsewhere.
|
||||
C>
|
||||
subroutine int_1efde(i_basis,ish,j_basis,jsh,fde_basis,lscr,scr,
|
||||
& lstv,S,T,V,doS,doT,doV)
|
||||
implicit none
|
||||
#include "nwc_const.fh"
|
||||
#include "errquit.fh"
|
||||
#include "basP.fh"
|
||||
#include "basdeclsP.fh"
|
||||
#include "geomP.fh"
|
||||
#include "geom.fh"
|
||||
#include "geobasmapP.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "bas_exndcf_dec.fh"
|
||||
#include "bas_ibs_dec.fh"
|
||||
#include "int_nbf.fh"
|
||||
#include "stdio.fh"
|
||||
#include "apiP.fh"
|
||||
#include "rel_nwc.fh"
|
||||
#include "util.fh"
|
||||
C
|
||||
C The following include block contains the geometry handle for the
|
||||
C geometry and charges of the solvent for COSMO. This is a clear
|
||||
C violation
|
||||
C of the integral API, but for now the cleanest solution
|
||||
C
|
||||
#include "prop.fh"
|
||||
c::external subroutines used
|
||||
c... errquit
|
||||
c::functions
|
||||
logical cando_hnd_1e
|
||||
logical cando_nw_1e
|
||||
logical cando_nw
|
||||
logical int_chk_init
|
||||
logical int_chk_sh
|
||||
external int_chk_init
|
||||
external int_chk_sh
|
||||
external cando_hnd_1e
|
||||
external cando_nw_1e
|
||||
external cando_nw
|
||||
c::passed
|
||||
c:tex-\begin{verbatim}
|
||||
integer i_basis !< [Input] basis set handle for ish
|
||||
integer ish !< [Input] i shell/contraction
|
||||
integer j_basis !< [Input] basis set handle for jsh
|
||||
integer jsh !< [Input] j shell/contraction
|
||||
integer fde_basis !< [Input] FDE basis set handle
|
||||
integer lscr !< [Input] length of scratch array
|
||||
integer lstv !< [Input] length of integral buffer
|
||||
double precision scr(lscr) !< [Scratch] scratch array
|
||||
double precision S(lstv) !< [Output] overlap integrals
|
||||
double precision T(lstv) !< [Output] kinetic energy integrals
|
||||
double precision V(lstv) !< [Output] potential energy integrals
|
||||
logical doS !< [Input] flag for overlap integrals
|
||||
logical doT !< [Input] flag for kinetic energy integrals
|
||||
logical doV !< [Input] flag for potential energy integrals
|
||||
c:tex-\end{verbatim}
|
||||
c::local
|
||||
logical ohnd_ok, onw_ok, onw1e_ok
|
||||
logical shells_ok, orel, oirel, ojrel, oNR, canAB
|
||||
integer i_geom, j_geom, ibas, jbas, ucont, uconts
|
||||
integer lbas, sbas, abas, isbas, jsbas
|
||||
integer Li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_icfpS
|
||||
integer Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_icfpS
|
||||
c
|
||||
integer fde_bas, fde_geom
|
||||
c
|
||||
logical any_spherical
|
||||
integer i_nbf_x, j_nbf_x
|
||||
integer i_nbf_s, j_nbf_s
|
||||
c
|
||||
integer WarnP
|
||||
save WarnP
|
||||
data WarnP /0/
|
||||
c
|
||||
#include "bas_exndcf_sfn.fh"
|
||||
#include "bas_ibs_sfn.fh"
|
||||
c
|
||||
c check initialization and shells
|
||||
c
|
||||
if (.not.int_chk_init('int_1efde'))
|
||||
& call errquit('int_1efde: int_init was not called' ,0,
|
||||
& INT_ERR)
|
||||
c
|
||||
shells_ok = int_chk_sh(i_basis,ish)
|
||||
shells_ok = shells_ok .and. int_chk_sh(j_basis,jsh)
|
||||
if (.not.shells_ok)
|
||||
& call errquit('int_1efde: invalid contraction/shell',0,
|
||||
& BASIS_ERR)
|
||||
c
|
||||
ibas = i_basis + BASIS_HANDLE_OFFSET
|
||||
jbas = j_basis + BASIS_HANDLE_OFFSET
|
||||
fde_bas = fde_basis + BASIS_HANDLE_OFFSET
|
||||
oNR = .true.
|
||||
oirel = .false.
|
||||
ojrel = .false.
|
||||
orel = .false.
|
||||
canAB = .false.
|
||||
c
|
||||
if (dyall_mod_dir) then
|
||||
c
|
||||
c get basis set handles; relativistic integral option only valid
|
||||
c if both ibas and jbas are the ao basis.
|
||||
c
|
||||
lbas = lc_bsh + BASIS_HANDLE_OFFSET
|
||||
sbas = sc_bsh + BASIS_HANDLE_OFFSET
|
||||
abas = ao_bsh + BASIS_HANDLE_OFFSET
|
||||
orel = ibas .eq. abas .and. jbas .eq. abas
|
||||
end if
|
||||
c
|
||||
c i shell
|
||||
c
|
||||
ucont = (sf_ibs_cn2ucn(ish,ibas))
|
||||
c
|
||||
c check for relativistic shell
|
||||
c
|
||||
if (orel .and. (infbs_cont(CONT_RELLS ,ucont,ibas) .ne. 0)) then
|
||||
oirel = .true.
|
||||
isbas = sbas
|
||||
uconts = ao_to_ls(ucont)
|
||||
if (uconts .eq. 0) call errquit (
|
||||
& 'int_1efde: no relativistic pointer',911, INT_ERR)
|
||||
if (nesc_1e_approx) then
|
||||
ibas = lbas
|
||||
ucont = uconts
|
||||
end if
|
||||
else
|
||||
uconts = ucont
|
||||
isbas = ibas
|
||||
end if
|
||||
c
|
||||
Li = infbs_cont(CONT_TYPE ,ucont,ibas)
|
||||
i_prim = infbs_cont(CONT_NPRIM,ucont,ibas)
|
||||
i_gen = infbs_cont(CONT_NGEN ,ucont,ibas)
|
||||
i_iexp = infbs_cont(CONT_IEXP ,ucont,ibas)
|
||||
i_icfp = infbs_cont(CONT_ICFP ,ucont,ibas)
|
||||
i_cent = (sf_ibs_cn2ce(ish,ibas))
|
||||
i_geom = ibs_geom(ibas)
|
||||
i_icfpS = infbs_cont(CONT_ICFP ,uconts,isbas)
|
||||
c
|
||||
c j shell
|
||||
c
|
||||
ucont = (sf_ibs_cn2ucn(jsh,jbas))
|
||||
c
|
||||
c check for relativistic shell
|
||||
c
|
||||
if (orel .and. (infbs_cont(CONT_RELLS ,ucont,jbas) .ne. 0)) then
|
||||
ojrel = .true.
|
||||
jsbas = sbas
|
||||
uconts = ao_to_ls(ucont)
|
||||
if (uconts .eq. 0) call errquit (
|
||||
& 'int_1efde: no relativistic pointer',911, INT_ERR)
|
||||
if (nesc_1e_approx) then
|
||||
jbas = lbas
|
||||
ucont = uconts
|
||||
end if
|
||||
else
|
||||
uconts = ucont
|
||||
jsbas = jbas
|
||||
end if
|
||||
c
|
||||
Lj = infbs_cont(CONT_TYPE ,ucont,jbas)
|
||||
j_prim = infbs_cont(CONT_NPRIM,ucont,jbas)
|
||||
j_gen = infbs_cont(CONT_NGEN ,ucont,jbas)
|
||||
j_iexp = infbs_cont(CONT_IEXP ,ucont,jbas)
|
||||
j_icfp = infbs_cont(CONT_ICFP ,ucont,jbas)
|
||||
j_cent = (sf_ibs_cn2ce(jsh,jbas))
|
||||
j_geom = ibs_geom(jbas)
|
||||
j_icfpS = infbs_cont(CONT_ICFP ,uconts,jsbas)
|
||||
c
|
||||
oNR = .not.(oirel.and.ojrel)
|
||||
orel = oirel.or.ojrel
|
||||
c
|
||||
if (i_geom.ne.j_geom.and.WarnP.eq.0) then
|
||||
write(luout,*)
|
||||
& 'int_1efde: WARNING: possible geometry inconsistency'
|
||||
write(luout,*)'i_basis geometry handle:',i_geom
|
||||
write(luout,*)'j_basis geometry handle:',j_geom
|
||||
WarnP = 1
|
||||
endif
|
||||
c
|
||||
c ----- current cosmo restriction ... -----
|
||||
c
|
||||
if(i_geom.ne.j_geom) then
|
||||
write(luout,*) 'i_geom and j_geom is not eq'
|
||||
call errquit('int_1efde: i_geom and j_geom must be the same',
|
||||
& 0,INT_ERR)
|
||||
endif
|
||||
|
||||
ohnd_ok = cando_hnd_1e(i_basis,ish,0)
|
||||
& .and. cando_hnd_1e(j_basis,jsh,0)
|
||||
& .and. (.not.geom_any_finuc (i_geom))
|
||||
& .and. (.not.geom_any_finuc (j_geom))
|
||||
onw_ok = cando_nw(i_basis,ish,0) .and. cando_nw(j_basis,jsh,0)
|
||||
onw1e_ok = cando_nw_1e(i_basis,ish,0)
|
||||
& .and. cando_nw_1e(j_basis,jsh,0)
|
||||
c
|
||||
c get external charges here (MV)
|
||||
c ------------------------------
|
||||
c if(.not.geom_extbq_on())
|
||||
c > call errquit('int_1eefc1:no active bqs',0,0)
|
||||
c bq_ncent = geom_extbq_ncenter()
|
||||
c i_cbq = geom_extbq_coord()
|
||||
c i_qbq = geom_extbq_charge()
|
||||
c if(.not.ma_push_get(mt_dbl,bq_ncent,'xbq',h_xbq0,i_xbq0))
|
||||
c + call errquit( 'int_1eefc1',0,0)
|
||||
c call dfill(bq_ncent,0.d0,dbl_mb(i_xbq0),1)
|
||||
c
|
||||
fde_geom = ibs_geom(fde_bas)
|
||||
|
||||
c write(*,*) 'FDE_GEOM ', fde_geom
|
||||
c write(*,*) 'COORDS ', coords(1,1,fde_geom)
|
||||
c write(*,*) 'CHARGE ', charge(1,fde_geom)
|
||||
c write(*,*) 'INVNUCLEARCHARGE', geom_invnucexp(1,fde_geom)
|
||||
c write(*,*) 'NCENTER ', ncenter(fde_geom)
|
||||
|
||||
c call util_flush(luout)
|
||||
|
||||
if (orel) then
|
||||
call rel_onel (
|
||||
& coords(1,i_cent,i_geom),
|
||||
& dbl_mb(mb_exndcf(i_iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(i_icfp,ibas)),
|
||||
& dbl_mb(mb_exndcf(i_icfpS,isbas)),i_prim,i_gen,Li,
|
||||
& coords(1,j_cent,j_geom),
|
||||
& dbl_mb(mb_exndcf(j_iexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(j_icfp,jbas)),
|
||||
& dbl_mb(mb_exndcf(j_icfpS,jsbas)),j_prim,j_gen,Lj,
|
||||
& coords(1,1,fde_geom),charge(1,fde_geom),
|
||||
& geom_invnucexp(1,fde_geom),ncenter(fde_geom),
|
||||
& S,T,V,lstv,doS,doT,doV,canAB,onw_ok,ohnd_ok,oNR,.false.,
|
||||
& scr,lscr,rel_dbg,rel_typ)
|
||||
else if (ohnd_ok) then
|
||||
call hnd_stvint(
|
||||
& coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(i_icfp,ibas)),
|
||||
& i_prim, i_gen, Li,
|
||||
& coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(j_icfp,jbas)),
|
||||
& j_prim, j_gen, Lj,
|
||||
& coords(1,1,fde_geom),charge(1,fde_geom),
|
||||
& ncenter(fde_geom),
|
||||
& S,T,V,lstv,doS,doT,doV,scr,lscr)
|
||||
c
|
||||
elseif (onw1e_ok) then
|
||||
call int_hf1sp(
|
||||
& coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(i_icfp,ibas)),
|
||||
& i_prim, i_gen, Li, i_cent,
|
||||
& coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(j_icfp,jbas)),
|
||||
& j_prim, j_gen, Lj, j_cent,
|
||||
& coords(1,1,fde_geom),charge(1,fde_geom),
|
||||
& geom_invnucexp(1,fde_geom),ncenter(fde_geom),
|
||||
& S,T,V,lstv,doS,doT,doV,canAB,.false.,
|
||||
& scr,lscr,'int_1eefc')
|
||||
elseif (onw_ok) then
|
||||
call hf1(
|
||||
& coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(i_icfp,ibas)), i_prim, i_gen, Li,
|
||||
& coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(j_icfp,jbas)), j_prim, j_gen, Lj,
|
||||
& coords(1,1,fde_geom),charge(1,fde_geom),
|
||||
& geom_invnucexp(1,fde_geom),ncenter(fde_geom),
|
||||
& S,T,V,lstv,doS,doT,doV,canAB,.false.,
|
||||
& scr,lscr)
|
||||
else
|
||||
call errquit('int_1efde: could not do hnd, sp or nw integrals',
|
||||
& 0, INT_ERR)
|
||||
endif
|
||||
c
|
||||
c if(.not.ma_pop_stack(h_xbq0))
|
||||
c + call errquit( 'int_1eefc1',0,0)
|
||||
c
|
||||
* We now have the cartesian integral block(s) (jlo:jhi,ilo:ihi)
|
||||
*
|
||||
any_spherical = bas_spherical(ibas).or.bas_spherical(jbas)
|
||||
if (.not.any_spherical) return
|
||||
c
|
||||
c ... reset general contractions for sp shells to 1 since they are handled
|
||||
c as a block of 4. Since int_nbf_* arrays are set to the appropriate size.
|
||||
c
|
||||
if (li.eq.-1) i_gen = 1
|
||||
if (lj.eq.-1) j_gen = 1
|
||||
c
|
||||
if (bas_spherical(ibas).and.bas_spherical(jbas)) then
|
||||
*... transform both i and j integrals
|
||||
i_nbf_x = int_nbf_x(Li)
|
||||
i_nbf_s = int_nbf_s(Li)
|
||||
j_nbf_x = int_nbf_x(Lj)
|
||||
j_nbf_s = int_nbf_s(Lj)
|
||||
c
|
||||
if (doS) call spcart_tran1e(S,scr,
|
||||
& j_nbf_x,i_nbf_x,Lj,j_gen,
|
||||
& j_nbf_s,i_nbf_s,Li,i_gen,
|
||||
& .false.)
|
||||
if (doT) call spcart_tran1e(T,scr,
|
||||
& j_nbf_x,i_nbf_x,Lj,j_gen,
|
||||
& j_nbf_s,i_nbf_s,Li,i_gen,
|
||||
& .false.)
|
||||
if (doV) call spcart_tran1e(V,scr,
|
||||
& j_nbf_x,i_nbf_x,Lj,j_gen,
|
||||
& j_nbf_s,i_nbf_s,Li,i_gen,
|
||||
& .false.)
|
||||
else if (bas_spherical(ibas)) then
|
||||
*.. transform on i component
|
||||
i_nbf_x = int_nbf_x(Li)
|
||||
i_nbf_s = int_nbf_s(Li)
|
||||
j_nbf_x = int_nbf_x(Lj)
|
||||
j_nbf_s = j_nbf_x
|
||||
if (doS) call spcart_tran1e(S,scr,
|
||||
& j_nbf_x,i_nbf_x,0,j_gen,
|
||||
& j_nbf_s,i_nbf_s,Li,i_gen,
|
||||
& .false.)
|
||||
if (doT) call spcart_tran1e(T,scr,
|
||||
& j_nbf_x,i_nbf_x,0,j_gen,
|
||||
& j_nbf_s,i_nbf_s,Li,i_gen,
|
||||
& .false.)
|
||||
if (doV) call spcart_tran1e(V,scr,
|
||||
& j_nbf_x,i_nbf_x,0,j_gen,
|
||||
& j_nbf_s,i_nbf_s,Li,i_gen,
|
||||
& .false.)
|
||||
else if (bas_spherical(jbas)) then
|
||||
*.. transform on j component
|
||||
i_nbf_x = int_nbf_x(Li)
|
||||
i_nbf_s = i_nbf_x
|
||||
j_nbf_x = int_nbf_x(Lj)
|
||||
j_nbf_s = int_nbf_s(Lj)
|
||||
if (doS) call spcart_tran1e(S,scr,
|
||||
& j_nbf_x,i_nbf_x,Lj,j_gen,
|
||||
& j_nbf_s,i_nbf_s,0,i_gen,
|
||||
& .false.)
|
||||
if (doT) call spcart_tran1e(T,scr,
|
||||
& j_nbf_x,i_nbf_x,Lj,j_gen,
|
||||
& j_nbf_s,i_nbf_s,0,i_gen,
|
||||
& .false.)
|
||||
if (doV) call spcart_tran1e(V,scr,
|
||||
& j_nbf_x,i_nbf_x,Lj,j_gen,
|
||||
& j_nbf_s,i_nbf_s,0,i_gen,
|
||||
& .false.)
|
||||
else
|
||||
call errquit(
|
||||
& 'int_1efde: should never reach transform blocked else',911,
|
||||
& INT_ERR)
|
||||
endif
|
||||
return
|
||||
c
|
||||
end
|
||||
cc AJL/End
|
||||
C>
|
||||
C> @}
|
||||
|
|
|
|||
|
|
@ -294,11 +294,25 @@ c::local
|
|||
& 'int_hf1sp: unknown center for ECPs',
|
||||
& 911, INT_ERR)
|
||||
endif
|
||||
call int_ecp_hf1(
|
||||
& xyzi,expi,coefi,i_nprim,i_ngen,Li,ictri,
|
||||
& xyzj,expj,coefj,j_nprim,j_ngen,Lj,ictrj,
|
||||
& scr,nintV,scr(offset),(lscr-nintV),
|
||||
& dryrun)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (msg.eq.'int_1epe_beta') then
|
||||
cc Call to spin polarised ecp calculation for Beta channel
|
||||
call int_ecp_hf1_beta(
|
||||
& xyzi,expi,coefi,i_nprim,i_ngen,Li,ictri,
|
||||
& xyzj,expj,coefj,j_nprim,j_ngen,Lj,ictrj,
|
||||
& scr,nintV,scr(offset),(lscr-nintV),
|
||||
& dryrun)
|
||||
else
|
||||
cc Original call to spin independent ECP
|
||||
call int_ecp_hf1(
|
||||
& xyzi,expi,coefi,i_nprim,i_ngen,Li,ictri,
|
||||
& xyzj,expj,coefj,j_nprim,j_ngen,Lj,ictrj,
|
||||
& scr,nintV,scr(offset),(lscr-nintV),
|
||||
& dryrun)
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
*... sum ecp into V block
|
||||
call daxpy(nintV, 1.0d00, scr,1,V,1)
|
||||
endif
|
||||
|
|
|
|||
|
|
@ -108,11 +108,14 @@ c
|
|||
offset = offset + 1 ! increment for passing to intd_1eh1P
|
||||
c
|
||||
call intd_1eh1P(i_basis,ish,j_basis,jsh,
|
||||
& scrsize,scr(offset),nint,H1a,scr)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
& scrsize,scr(offset),nint,H1a,scr,1)
|
||||
c ^ ECP Integral (ALPHA)
|
||||
c
|
||||
end
|
||||
subroutine intd_1eh1P(i_basis,ish,j_basis,jsh,lscr,scr,
|
||||
& nint,H1a,Ta)
|
||||
& nint,H1a,Ta,ecp_channel)
|
||||
cc AJL/End
|
||||
implicit none
|
||||
#include "stdio.fh"
|
||||
#include "errquit.fh"
|
||||
|
|
@ -146,6 +149,9 @@ c ! NOTE: nint*3 integral derivatives returned per unique
|
|||
double precision scr(lscr) ! [input] scratch array
|
||||
double precision H1a(nint,3,*) ! [output] derivative integrals (nint,3,n_atoms)
|
||||
double precision Ta(nint,3,2) ! [scratch] space for kinetic integrals
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_channel !< [Input] ECP Channel (1 = Alpha, 2 = Beta)
|
||||
cc AJL/End
|
||||
c::local
|
||||
logical doT
|
||||
integer ucont,uconts
|
||||
|
|
@ -347,23 +353,45 @@ c
|
|||
offset = nintV*3*nat + 1
|
||||
* write(luout,*)' lscr to ecp_hf1:',(lscr-nintV)
|
||||
* call dcopy(nintV*3*nat,0.0d00,0,scr,1) ! buffer zeroed in ecp_gradient
|
||||
call intd_ecp_hf1(
|
||||
& coords(1,iatom,igeom),
|
||||
& dbl_mb(mb_exndcf(iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(icf,ibas)),
|
||||
& inp,igen,itype,iatom,
|
||||
c
|
||||
& coords(1,jatom,jgeom),
|
||||
& dbl_mb(mb_exndcf(jexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(jcf,jbas)),
|
||||
& jnp,jgen,jtype,jatom,
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (ecp_channel.eq.1) then
|
||||
call intd_ecp_hf1(
|
||||
& coords(1,iatom,igeom),
|
||||
& dbl_mb(mb_exndcf(iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(icf,ibas)),
|
||||
& inp,igen,itype,iatom,
|
||||
c
|
||||
& scr,nintV,nat,
|
||||
& scr(offset),(lscr-offset-1),
|
||||
& .false.)
|
||||
c
|
||||
*... sum ecp into derivative H1 block
|
||||
& coords(1,jatom,jgeom),
|
||||
& dbl_mb(mb_exndcf(jexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(jcf,jbas)),
|
||||
& jnp,jgen,jtype,jatom,
|
||||
c
|
||||
& scr,nintV,nat,
|
||||
& scr(offset),(lscr-offset-1),
|
||||
& .false.)
|
||||
c
|
||||
else
|
||||
call intd_ecp_hf1_beta(
|
||||
& coords(1,iatom,igeom),
|
||||
& dbl_mb(mb_exndcf(iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(icf,ibas)),
|
||||
& inp,igen,itype,iatom,
|
||||
c
|
||||
& coords(1,jatom,jgeom),
|
||||
& dbl_mb(mb_exndcf(jexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(jcf,jbas)),
|
||||
& jnp,jgen,jtype,jatom,
|
||||
c
|
||||
& scr,nintV,nat,
|
||||
& scr(offset),(lscr-offset-1),
|
||||
& .false.)
|
||||
c
|
||||
endif
|
||||
c*... sum ecp into derivative H1 block
|
||||
call daxpy(nintV*3*nat,1.0d00,scr,1,H1a,1)
|
||||
cc AJL/End
|
||||
c
|
||||
endif
|
||||
c
|
||||
* H1a now has the cartesian integral block (jlo:jhi,ilo:ihi)
|
||||
|
|
@ -443,4 +471,52 @@ c
|
|||
c
|
||||
|
||||
end
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc This is a wrapper, identical to intd_1epot except we have
|
||||
cc separated out for alpha and beta channels
|
||||
subroutine intd_1eh1_beta(i_basis,ish,j_basis,jsh,lscr,scr,
|
||||
& lH1a,H1a)
|
||||
implicit none
|
||||
#include "stdio.fh"
|
||||
#include "errquit.fh"
|
||||
#include "nwc_const.fh"
|
||||
#include "basP.fh"
|
||||
#include "basdeclsP.fh"
|
||||
#include "geomP.fh"
|
||||
#include "geobasmapP.fh"
|
||||
c::functions
|
||||
integer int_nint_cart
|
||||
external int_nint_cart
|
||||
c::passed
|
||||
integer i_basis !< [Input] ish basis set handle
|
||||
integer ish !< [Input] ``i'' contraction index
|
||||
integer j_basis !< [Input] jsh basis set handle
|
||||
integer jsh !< [Input] ``j'' contraction index
|
||||
integer lscr !< [Input] length of scratch space
|
||||
integer lH1a !< [Input] number of h1 integral derivatives in shells ish and jsh
|
||||
c ! NOTE: nint*3 integral derivatives returned per unique center
|
||||
double precision scr(lscr) !< [Input] scratch array
|
||||
double precision H1a(*) !< [Output] derivative integrals
|
||||
c
|
||||
c::local
|
||||
integer nint, offset, scrsize, nat
|
||||
c
|
||||
nat = ncenter(ibs_geom((i_basis + Basis_Handle_Offset)))
|
||||
c
|
||||
nint = int_nint_cart(i_basis,ish,j_basis,jsh,0,0,0,0)
|
||||
if (nint*3*nat.gt.lH1a) then
|
||||
write(luout,*) 'nint*3*nat = ',nint*3*nat
|
||||
write(luout,*) 'lH1a = ',lH1a
|
||||
call errquit('intd_1eh1: nint>lH1a error',911, INT_ERR)
|
||||
endif
|
||||
c
|
||||
offset = nint*3*2 ! scratch for Ta array in intd_1eh1P routine
|
||||
scrsize = lscr - offset ! new scratch array size
|
||||
offset = offset + 1 ! increment for passing to intd_1eh1P
|
||||
c
|
||||
call intd_1eh1P(i_basis,ish,j_basis,jsh,
|
||||
& scrsize,scr(offset),nint,H1a,scr,2)
|
||||
c ^ ECP Integral (BETA)
|
||||
end
|
||||
cc AJL/End
|
||||
C> @}
|
||||
|
|
|
|||
|
|
@ -107,11 +107,14 @@ c
|
|||
offset = offset + 1 ! increment for passing to intd_1eh1P
|
||||
c
|
||||
call intd_1epotP(i_basis,ish,j_basis,jsh,
|
||||
& scrsize,scr(offset),nint,H1a,scr)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
& scrsize,scr(offset),nint,H1a,scr,1)
|
||||
c ^ ECP Channel (ALPHA)
|
||||
c
|
||||
end
|
||||
subroutine intd_1epotP(i_basis,ish,j_basis,jsh,lscr,scr,
|
||||
& nint,H1a,Ta)
|
||||
& nint,H1a,Ta,ecp_channel)
|
||||
cc AJL/End
|
||||
implicit none
|
||||
#include "stdio.fh"
|
||||
#include "errquit.fh"
|
||||
|
|
@ -145,6 +148,9 @@ c ! NOTE: nint*3 integral derivatives returned per unique
|
|||
double precision scr(lscr) !< [Input] scratch array
|
||||
double precision H1a(nint,3,*) !< [Output] derivative integrals (nint,3,n_atoms)
|
||||
double precision Ta(nint,3,2) !< [Scratch] space for kinetic integrals
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_channel !< [Input] ECP Channel (1 = Alpha, 2 = Beta)
|
||||
cc AJL/End
|
||||
c::local
|
||||
logical doT
|
||||
integer ucont,uconts
|
||||
|
|
@ -395,23 +401,44 @@ c
|
|||
offset = nintV*3*nat + 1
|
||||
* write(luout,*)' lscr to ecp_hf1:',(lscr-nintV)
|
||||
* call dcopy(nintV*3*nat,0.0d00,0,scr,1) ! buffer zeroed in ecp_gradient
|
||||
call intd_ecp_hf1(
|
||||
& coords(1,iatom,igeom),
|
||||
& dbl_mb(mb_exndcf(iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(icf,ibas)),
|
||||
& inp,igen,itype,iatom,
|
||||
c
|
||||
& coords(1,jatom,jgeom),
|
||||
& dbl_mb(mb_exndcf(jexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(jcf,jbas)),
|
||||
& jnp,jgen,jtype,jatom,
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (ecp_channel.eq.1) then
|
||||
call intd_ecp_hf1(
|
||||
& coords(1,iatom,igeom),
|
||||
& dbl_mb(mb_exndcf(iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(icf,ibas)),
|
||||
& inp,igen,itype,iatom,
|
||||
c
|
||||
& scr,nintV,nat,
|
||||
& scr(offset),(lscr-offset-1),
|
||||
& .false.)
|
||||
& coords(1,jatom,jgeom),
|
||||
& dbl_mb(mb_exndcf(jexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(jcf,jbas)),
|
||||
& jnp,jgen,jtype,jatom,
|
||||
c
|
||||
& scr,nintV,nat,
|
||||
& scr(offset),(lscr-offset-1),
|
||||
& .false.)
|
||||
c
|
||||
else
|
||||
call intd_ecp_hf1_beta(
|
||||
& coords(1,iatom,igeom),
|
||||
& dbl_mb(mb_exndcf(iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(icf,ibas)),
|
||||
& inp,igen,itype,iatom,
|
||||
c
|
||||
& coords(1,jatom,jgeom),
|
||||
& dbl_mb(mb_exndcf(jexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(jcf,jbas)),
|
||||
& jnp,jgen,jtype,jatom,
|
||||
c
|
||||
& scr,nintV,nat,
|
||||
& scr(offset),(lscr-offset-1),
|
||||
& .false.)
|
||||
c
|
||||
endif
|
||||
*... sum ecp into derivative H1 block
|
||||
call daxpy(nintV*3*nat,1.0d00,scr,1,H1a,1)
|
||||
cc AJL/End
|
||||
endif
|
||||
c
|
||||
* H1a now has the cartesian integral block (jlo:jhi,ilo:ihi)
|
||||
|
|
@ -490,5 +517,50 @@ c
|
|||
endif
|
||||
c
|
||||
|
||||
end
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc This is a wrapper, identical to intd_1epot except we have separated out
|
||||
cc For alpha and beta channels
|
||||
subroutine intd_1epot_beta(i_basis,ish,j_basis,jsh,lscr,scr,
|
||||
& lH1a,H1a)
|
||||
implicit none
|
||||
#include "stdio.fh"
|
||||
#include "errquit.fh"
|
||||
#include "nwc_const.fh"
|
||||
#include "basP.fh"
|
||||
#include "basdeclsP.fh"
|
||||
#include "geomP.fh"
|
||||
#include "geobasmapP.fh"
|
||||
integer int_nint_cart
|
||||
external int_nint_cart
|
||||
c::passed
|
||||
integer i_basis !< [Input] ish basis set handle
|
||||
integer ish !< [Input] `i` contraction index
|
||||
integer j_basis !< [Input] jsh basis set handle
|
||||
integer jsh !< [Input] `j` contraction index
|
||||
integer lscr !< [Input] length of scratch space
|
||||
integer lH1a !< [Input] number of h1 integral derivatives in shells ish and jsh
|
||||
c ! NOTE: nint*3 integral derivatives returned per unique center
|
||||
double precision scr(lscr) !< [Input] scratch array
|
||||
double precision H1a(*) !< [Output] derivative integrals
|
||||
c::local
|
||||
integer nint, offset, scrsize, nat
|
||||
c
|
||||
nat = ncenter(ibs_geom((i_basis + Basis_Handle_Offset)))
|
||||
c
|
||||
nint = int_nint_cart(i_basis,ish,j_basis,jsh,0,0,0,0)
|
||||
if (nint*3*nat.gt.lH1a) then
|
||||
write(luout,*) 'nint*3*nat = ',nint*3*nat
|
||||
write(luout,*) 'lH1a = ',lH1a
|
||||
call errquit('intd_1eh1: nint>lH1a error',911, INT_ERR)
|
||||
endif
|
||||
c
|
||||
offset = nint*3*2 ! scratch for Ta array in intd_1eh1P routine
|
||||
scrsize = lscr - offset ! new scratch array size
|
||||
offset = offset + 1 ! increment for passing to intd_1eh1P
|
||||
c
|
||||
call intd_1epotP(i_basis,ish,j_basis,jsh,
|
||||
& scrsize,scr(offset),nint,H1a,scr,2)
|
||||
c ^ ECP Channel (BETA)
|
||||
end
|
||||
C> @}
|
||||
|
|
|
|||
|
|
@ -102,11 +102,14 @@ c
|
|||
endif
|
||||
c
|
||||
call intdd_1eecp1P(i_basis,ish,j_basis,jsh,
|
||||
& lscr,scr,nint,H1a,nat)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
& lscr,scr,nint,H1a,nat,1)
|
||||
c ^ Alpha Channel
|
||||
c
|
||||
end
|
||||
subroutine intdd_1eecp1P(i_basis,ish,j_basis,jsh,lscr,scr,
|
||||
& nint,H1a,nat)
|
||||
& nint,H1a,nat,ecp_channel)
|
||||
cc AJL/End
|
||||
implicit none
|
||||
#include "stdio.fh"
|
||||
#include "errquit.fh"
|
||||
|
|
@ -144,6 +147,9 @@ c::local
|
|||
integer zatom, zyx1, zyx2, nintV, offset
|
||||
c
|
||||
logical any_spherical
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_channel
|
||||
cc AJL/End
|
||||
c
|
||||
c Temporary variable that needs to be taken out after testing!
|
||||
c
|
||||
|
|
@ -194,20 +200,40 @@ c
|
|||
c
|
||||
nat = ncenter(igeom) ! needed for both ecp and spherical
|
||||
nintV = int_nbf_x(itype)*int_nbf_x(jtype)
|
||||
call intdd_ecp_hf1(
|
||||
& coords(1,iatom,igeom),
|
||||
& dbl_mb(mb_exndcf(iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(icf,ibas)),
|
||||
& inp,igen,itype,iatom,
|
||||
c
|
||||
& coords(1,jatom,jgeom),
|
||||
& dbl_mb(mb_exndcf(jexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(jcf,jbas)),
|
||||
& jnp,jgen,jtype,jatom,
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (ecp_channel.eq.1) then
|
||||
call intdd_ecp_hf1(
|
||||
& coords(1,iatom,igeom),
|
||||
& dbl_mb(mb_exndcf(iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(icf,ibas)),
|
||||
& inp,igen,itype,iatom,
|
||||
c
|
||||
& H1a,nintV,nat,
|
||||
& scr,lscr,
|
||||
& .false.)
|
||||
& coords(1,jatom,jgeom),
|
||||
& dbl_mb(mb_exndcf(jexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(jcf,jbas)),
|
||||
& jnp,jgen,jtype,jatom,
|
||||
c
|
||||
& H1a,nintV,nat,
|
||||
& scr,lscr,
|
||||
& .false.)
|
||||
else
|
||||
call intdd_ecp_hf1_beta(
|
||||
& coords(1,iatom,igeom),
|
||||
& dbl_mb(mb_exndcf(iexp,ibas)),
|
||||
& dbl_mb(mb_exndcf(icf,ibas)),
|
||||
& inp,igen,itype,iatom,
|
||||
c
|
||||
& coords(1,jatom,jgeom),
|
||||
& dbl_mb(mb_exndcf(jexp,jbas)),
|
||||
& dbl_mb(mb_exndcf(jcf,jbas)),
|
||||
& jnp,jgen,jtype,jatom,
|
||||
c
|
||||
& H1a,nintV,nat,
|
||||
& scr,lscr,
|
||||
& .false.)
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
any_spherical = bas_spherical(ibas).or.bas_spherical(jbas)
|
||||
if (.not.any_spherical) return
|
||||
|
|
@ -289,4 +315,50 @@ c
|
|||
endif
|
||||
c
|
||||
end
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc This is a wrapper, identical to intd_1epot except we have
|
||||
cc separated out for alpha and beta channels
|
||||
subroutine intdd_1eecp1_beta(i_basis,ish,j_basis,jsh,lscr,scr,
|
||||
& lH1a,H1a)
|
||||
implicit none
|
||||
#include "stdio.fh"
|
||||
#include "errquit.fh"
|
||||
#include "nwc_const.fh"
|
||||
#include "basP.fh"
|
||||
#include "basdeclsP.fh"
|
||||
#include "geomP.fh"
|
||||
#include "geobasmapP.fh"
|
||||
c::functions
|
||||
integer int_nint_cart
|
||||
external int_nint_cart
|
||||
c::passed
|
||||
integer i_basis !< [Input] ish basis set handle
|
||||
integer ish !< [Input] `i` contraction index
|
||||
integer j_basis !< [Input] jsh basis set handle
|
||||
integer jsh !< [Input] `j` contraction index
|
||||
integer lscr !< [Input] length of scratch space
|
||||
integer lH1a !< [Input] number of h1 integral derivatives in shells ish and jsh
|
||||
c ! NOTE: nint*3 integral derivatives returned per unique center
|
||||
double precision scr(lscr) !< [Input] scratch array
|
||||
double precision H1a(*) !< [Output] derivative integrals
|
||||
c
|
||||
c::local
|
||||
integer nint, offset, nat
|
||||
c
|
||||
nat = ncenter(ibs_geom((i_basis + Basis_Handle_Offset)))
|
||||
c
|
||||
nint = int_nint_cart(i_basis,ish,j_basis,jsh,0,0,0,0)
|
||||
if (nint*3*3*(nat*(nat-1)/2+nat).gt.lH1a) then
|
||||
write(luout,*) 'nint*3*3*(nat*(nat-1)/2+nat) = ',
|
||||
& nint*3*3*(nat*(nat-1)/2+nat)
|
||||
write(luout,*) 'lH1a = ',lH1a
|
||||
call errquit('intd_1eecp: nint>lH1a error',911, INT_ERR)
|
||||
endif
|
||||
c
|
||||
call intdd_1eecp1P(i_basis,ish,j_basis,jsh,
|
||||
& lscr,scr,nint,H1a,nat,2)
|
||||
c ^ Beta Channel
|
||||
c
|
||||
end
|
||||
cc AJL/End
|
||||
C> @}
|
||||
|
|
|
|||
|
|
@ -20,10 +20,29 @@ c
|
|||
integer k_ecp_c2s ! MA index for c2s routines
|
||||
integer h_ecp_lip ! MA handle for ecp center index
|
||||
integer k_ecp_lip ! MA index for ecp center index
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer h_ecp_e_beta ! MA handle for ecp exponents
|
||||
integer k_ecp_e_beta ! MA index for ecp exponents
|
||||
integer h_ecp_c_beta ! MA handle for ecp coefficients
|
||||
integer k_ecp_c_beta ! MA index for ecp coefficients
|
||||
integer h_ecp_nprim_c_beta ! MA handle for n_prim_C (see ecp_integral)
|
||||
integer k_ecp_nprim_c_beta ! MA index for n_prim_C (see ecp_integral)
|
||||
integer h_ecp_ncoef_c_beta ! MA handle for n_coef_C (see ecp_integral)
|
||||
integer k_ecp_ncoef_c_beta ! MA index for n_coef_C (see ecp_integral)
|
||||
integer h_ecp_ind_c_beta ! MA handle for ind_C (see ecp_integral)
|
||||
integer k_ecp_ind_c_beta ! MA index for ind_C (see ecp_integral)
|
||||
integer h_ecp_ind_z_beta ! MA handle for ind_Z (see ecp_integral)
|
||||
integer k_ecp_ind_z_beta ! MA index for ind_Z (see ecp_integral)
|
||||
integer h_ecp_l_c_beta ! MA handle for l_C (see ecp_integral)
|
||||
integer k_ecp_l_c_beta ! MA index for l_C (see ecp_integral)
|
||||
integer h_ecp_lip_beta ! MA handle for ecp center index
|
||||
integer k_ecp_lip_beta ! MA index for ecp center index
|
||||
cc AJL/End
|
||||
integer n_zeta_c ! length of ecp exp/coef array
|
||||
integer l_ecp ! high ang for ecp basis
|
||||
integer n_ecp ! number of ecp centers (from API)
|
||||
integer mem_c2s ! size of c2s array
|
||||
integer ecp_channels ! Number of ecp channels
|
||||
logical init_ecp_init ! logical saying if ecp is init-ed
|
||||
c
|
||||
common /c_ecp_nwc/ h_xyzecp, k_xyzecp, h_ecp_e, k_ecp_e,
|
||||
|
|
@ -32,5 +51,16 @@ c
|
|||
& h_ecp_ind_c, k_ecp_ind_c, h_ecp_ind_z, k_ecp_ind_z,
|
||||
& h_ecp_l_c, k_ecp_l_c, h_ecp_c2s, k_ecp_c2s,
|
||||
& h_ecp_lip, k_ecp_lip,
|
||||
& n_zeta_c, l_ecp, n_ecp, mem_c2s, init_ecp_init
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c & n_zeta_c, l_ecp, n_ecp, mem_c2s, init_ecp_init
|
||||
& h_ecp_e_beta, k_ecp_e_beta,
|
||||
& h_ecp_c_beta, k_ecp_c_beta,
|
||||
& h_ecp_nprim_c_beta, k_ecp_nprim_c_beta,
|
||||
& h_ecp_ncoef_c_beta, k_ecp_ncoef_c_beta,
|
||||
& h_ecp_ind_c_beta, k_ecp_ind_c_beta,
|
||||
& h_ecp_ind_z_beta, k_ecp_ind_z_beta,
|
||||
& h_ecp_l_c_beta, k_ecp_l_c_beta,
|
||||
& h_ecp_lip_beta, k_ecp_lip_beta,
|
||||
& n_zeta_c, l_ecp, n_ecp, mem_c2s, ecp_channels, init_ecp_init
|
||||
cc AJL/End
|
||||
|
||||
|
|
|
|||
|
|
@ -220,6 +220,19 @@ c
|
|||
endif
|
||||
n_zeta_c = n_zeta_c_e ! use e one
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c
|
||||
c... determine the number of channels needed (Spin-polarised ECPs)
|
||||
ecp_channels = 1
|
||||
if (o_ecp) then
|
||||
if (.not.ecp_get_high_chan(ecpidin,ecp_channels))
|
||||
& call errquit('int_ecp_init error',911, INT_ERR)
|
||||
end if
|
||||
cc AJL/Debug
|
||||
c write(6,*) 'Number of ECP Channels: ', ecp_channels
|
||||
c write(6,*) '(1 = Spin Independent, Default; 2 = Spin Polarised)'
|
||||
cc AJL/End
|
||||
c
|
||||
c.. allocate space for ecp/so exponents and ecp/so coefficients
|
||||
if (.not.ma_alloc_get(mt_dbl,
|
||||
& n_zeta_c,
|
||||
|
|
@ -234,6 +247,25 @@ c.. allocate space for ecp/so exponents and ecp/so coefficients
|
|||
& (' int_ecp_init: ecp coefficients ma failed ',911, BASIS_ERR)
|
||||
call dfill(n_zeta_c,0.0d00,dbl_mb(k_ecp_c),1)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (ecp_channels.gt.1) then
|
||||
if (.not.ma_alloc_get(mt_dbl,
|
||||
& n_zeta_c,
|
||||
& 'beta ecp exponents',
|
||||
& h_ecp_e_beta, k_ecp_e_beta)) call errquit
|
||||
& (' int_ecp_init: beta ecp exponent ma failed ',
|
||||
& 911, BASIS_ERR)
|
||||
call dfill(n_zeta_c,0.0d00,dbl_mb(k_ecp_e_beta),1)
|
||||
if (.not.ma_alloc_get(mt_dbl,
|
||||
& n_zeta_c,
|
||||
& 'beta ecp coefficients',
|
||||
& h_ecp_c_beta, k_ecp_c_beta)) call errquit
|
||||
& (' int_ecp_init: beta ecp coefficients ma failed ',
|
||||
& 911, BASIS_ERR)
|
||||
call dfill(n_zeta_c,0.0d00,dbl_mb(k_ecp_c_beta),1)
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
c... determine maximum angular momentum of ecp basis
|
||||
l_ecp = 0
|
||||
if (o_ecp) then
|
||||
|
|
@ -290,6 +322,55 @@ c... allocate space for ecp center pointer list
|
|||
& h_ecp_lip, k_ecp_lip)) call errquit
|
||||
& (' int_ecp: ecp_lip ma failed',911, BASIS_ERR)
|
||||
call ifill(n_ecp,0,int_mb(k_ecp_lip),1)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c... allocate space for n_prim_C(0:4,-1:l_ecp_max,n_C*2),
|
||||
if (ecp_channels.gt.1) then
|
||||
if (.not.ma_alloc_get(mt_int,
|
||||
& (5*l_ecp_sz*n_ecp*2),
|
||||
& 'beta ecp n_prim_C',
|
||||
& h_ecp_nprim_c_beta, k_ecp_nprim_c_beta)) call errquit
|
||||
& (' int_ecp: beta ecp nprim_c ma failed ',911, BASIS_ERR)
|
||||
call ifill((5*l_ecp_sz*n_ecp*2),
|
||||
& 0,int_mb(k_ecp_nprim_c_beta),1)
|
||||
c... allocate space for n_coef_C(-1:l_ecp_max,n_C*2)
|
||||
if (.not.ma_alloc_get(mt_int,
|
||||
& (l_ecp_sz*n_ecp*2),
|
||||
& 'beta ecp n_coef_C',
|
||||
& h_ecp_ncoef_c_beta, k_ecp_ncoef_c_beta)) call errquit
|
||||
& (' int_ecp: beta ecp ncoef_c ma failed ',911, BASIS_ERR)
|
||||
call ifill((l_ecp_sz*n_ecp*2),0,int_mb(k_ecp_ncoef_c_beta),1)
|
||||
c... allocate space for ind_C
|
||||
if (.not.ma_alloc_get(mt_int,
|
||||
& (l_ecp_sz*n_ecp*2),
|
||||
& 'beta ecp ind_C',
|
||||
& h_ecp_ind_c_beta, k_ecp_ind_c_beta)) call errquit
|
||||
& (' int_ecp: beta ecp ind_c ma failed ',911, BASIS_ERR)
|
||||
call ifill((l_ecp_sz*n_ecp*2),0,int_mb(k_ecp_ind_c_beta),1)
|
||||
c... allocate space for ind_z
|
||||
if (.not.ma_alloc_get(mt_int,
|
||||
& (l_ecp_sz*n_ecp*2),
|
||||
& 'beta ecp ind_z',
|
||||
& h_ecp_ind_z_beta, k_ecp_ind_z_beta)) call errquit
|
||||
& (' int_ecp: beta ecp ind_z ma failed ',911, BASIS_ERR)
|
||||
call ifill((l_ecp_sz*n_ecp*2),0,int_mb(k_ecp_ind_z_beta),1)
|
||||
c... allocate space for l_C
|
||||
if (.not.ma_alloc_get(mt_int,
|
||||
& n_ecp,
|
||||
& 'beta ecp l_C',
|
||||
& h_ecp_l_c_beta, k_ecp_l_c_beta)) call errquit
|
||||
& (' int_ecp: beta ecp l_c ma failed ',911, BASIS_ERR)
|
||||
call ifill(n_ecp,0,int_mb(k_ecp_l_c_beta),1)
|
||||
c... allocate space for ecp center pointer list
|
||||
if (.not.ma_alloc_get(mt_int,
|
||||
& n_ecp,
|
||||
& 'beta ecp lexical indices for ecp centers',
|
||||
& h_ecp_lip_beta, k_ecp_lip_beta)) call errquit
|
||||
& (' int_ecp: beta ecp_lip ma failed',911, BASIS_ERR)
|
||||
call ifill(n_ecp,0,int_mb(k_ecp_lip_beta),1)
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
c
|
||||
call int_ecp_build_ecp_ptrs(ecpidin,soidin,geom,
|
||||
& o_both,o_ecp,o_so,
|
||||
|
|
@ -297,6 +378,10 @@ c
|
|||
& n_ecp,
|
||||
& l_ecp,
|
||||
& n_zeta_c,
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
& 1, ! Channel to be filled. Would be 2 for beta channel
|
||||
& ecp_channels, ! Total no. of channels, for print_ptrs
|
||||
cc AJL/End
|
||||
& int_mb(k_ecp_nprim_c),
|
||||
& int_mb(k_ecp_ncoef_c),
|
||||
& int_mb(k_ecp_ind_c),
|
||||
|
|
@ -306,6 +391,26 @@ c
|
|||
& dbl_mb(k_ecp_e),
|
||||
& dbl_mb(k_ecp_c) )
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (ecp_channels.gt.1) ! Setup beta options
|
||||
& call int_ecp_build_ecp_ptrs(ecpidin,soidin,geom,
|
||||
& o_both,o_ecp,o_so,
|
||||
& ncenters,
|
||||
& n_ecp,
|
||||
& l_ecp,
|
||||
& n_zeta_c,
|
||||
& 2, ! Beta channel
|
||||
& ecp_channels, ! Total no. of channels, for print_ptrs
|
||||
& int_mb(k_ecp_nprim_c_beta),
|
||||
& int_mb(k_ecp_ncoef_c_beta),
|
||||
& int_mb(k_ecp_ind_c_beta),
|
||||
& int_mb(k_ecp_ind_z_beta),
|
||||
& int_mb(k_ecp_l_c_beta),
|
||||
& int_mb(k_ecp_lip_beta),
|
||||
& dbl_mb(k_ecp_e_beta),
|
||||
& dbl_mb(k_ecp_c_beta))
|
||||
cc AJL/End
|
||||
c
|
||||
*... allocate space for c2s and s2c internal ecp transformation routines
|
||||
c
|
||||
c determine lmax among ao basis and ecp basis
|
||||
|
|
@ -464,6 +569,50 @@ c
|
|||
endif
|
||||
c ! nint block used in api routine int_hf1sp
|
||||
lecp_mem = lecp_mem + nV
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc Need to account for scenario where we have two channels so
|
||||
cc double the number of ECPs could be generated
|
||||
if (o_ecp.and.ecp_channels.gt.1) then
|
||||
if (indx_grad.eq.0) then
|
||||
call int_ecp_hf1_beta(
|
||||
& coords(1,ic,ig),
|
||||
& dbl_mb(mb_exndcf(ie,bas_id)),
|
||||
& dbl_mb(mb_exndcf(icf,bas_id)),
|
||||
& inp,ing,il,ic,
|
||||
& coords(1,jc,jg),
|
||||
& dbl_mb(mb_exndcf(je,bas_id)),
|
||||
& dbl_mb(mb_exndcf(jcf,bas_id)),
|
||||
& jnp,jng,jl,jc,
|
||||
& dpdum,nV,dpdum,lecp_mem,.true.)
|
||||
elseif (indx_grad.eq.1) then
|
||||
call intd_ecp_hf1_beta(
|
||||
& coords(1,ic,ig),
|
||||
& dbl_mb(mb_exndcf(ie,bas_id)),
|
||||
& dbl_mb(mb_exndcf(icf,bas_id)),
|
||||
& inp,ing,il,ic,
|
||||
& coords(1,jc,jg),
|
||||
& dbl_mb(mb_exndcf(je,bas_id)),
|
||||
& dbl_mb(mb_exndcf(jcf,bas_id)),
|
||||
& jnp,jng,jl,jc,
|
||||
& dpdum,nV,ncenters,dpdum,lecp_mem,.true.)
|
||||
elseif (indx_grad.eq.2) then
|
||||
call intdd_ecp_hf1_beta(
|
||||
& coords(1,ic,ig),
|
||||
& dbl_mb(mb_exndcf(ie,bas_id)),
|
||||
& dbl_mb(mb_exndcf(icf,bas_id)),
|
||||
& inp,ing,il,ic,
|
||||
& coords(1,jc,jg),
|
||||
& dbl_mb(mb_exndcf(je,bas_id)),
|
||||
& dbl_mb(mb_exndcf(jcf,bas_id)),
|
||||
& jnp,jng,jl,jc,
|
||||
& dpdum,nV,ncenters,dpdum,lecp_mem,.true.)
|
||||
endif
|
||||
c Add data storage space to lecp_mem
|
||||
lecp_mem = lecp_mem + nV
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
* write(6,*)' i_cont ',i_cont
|
||||
* write(6,*)' j_cont ',j_cont
|
||||
mem_ecp = max(mem_ecp,lecp_mem)
|
||||
|
|
@ -537,6 +686,38 @@ c
|
|||
l_ecp = 0
|
||||
n_ecp = 0
|
||||
init_ecp_init = .false.
|
||||
c
|
||||
cc AJL/Begin
|
||||
if (ecp_channels.gt.1) then
|
||||
status = status .and. MA_free_heap(h_ecp_e_beta)
|
||||
status = status .and. MA_free_heap(h_ecp_c_beta)
|
||||
status = status .and. MA_free_heap(h_ecp_nprim_c_beta)
|
||||
status = status .and. MA_free_heap(h_ecp_ncoef_c_beta)
|
||||
status = status .and. MA_free_heap(h_ecp_ind_c_beta)
|
||||
status = status .and. MA_free_heap(h_ecp_ind_z_beta)
|
||||
status = status .and. MA_free_heap(h_ecp_l_c_beta)
|
||||
status = status .and. MA_free_heap(h_ecp_lip_beta)
|
||||
|
||||
h_ecp_e_beta = 0
|
||||
k_ecp_e_beta = 0
|
||||
h_ecp_c_beta = 0
|
||||
k_ecp_c_beta = 0
|
||||
h_ecp_nprim_c_beta = 0
|
||||
k_ecp_nprim_c_beta = 0
|
||||
h_ecp_ncoef_c_beta = 0
|
||||
k_ecp_ncoef_c_beta = 0
|
||||
h_ecp_ind_c_beta = 0
|
||||
k_ecp_ind_c_beta = 0
|
||||
h_ecp_ind_z_beta = 0
|
||||
k_ecp_ind_z_beta = 0
|
||||
h_ecp_l_c_beta = 0
|
||||
k_ecp_l_c_beta = 0
|
||||
h_ecp_lip_beta = 0
|
||||
k_ecp_lip_beta = 0
|
||||
ecp_channels = 0 ! Clear out ECP channels and we are done
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
if (status) return
|
||||
call errquit
|
||||
& ('int_ecp_terminate: error freeing heap',911, MEM_ERR)
|
||||
|
|
@ -547,6 +728,10 @@ c
|
|||
& n_ecp,
|
||||
& l_ecp,
|
||||
& nz_ecp,
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
& channel, ! spin channel
|
||||
& channels, ! total number of channels
|
||||
cc AJL/End
|
||||
& n_prim_C,
|
||||
& n_coef_C,
|
||||
& ind_C,
|
||||
|
|
@ -574,6 +759,10 @@ c
|
|||
integer n_ecp ! [input] number of ecp centers
|
||||
integer l_ecp ! [input] maximum angular momentum in ecp basis
|
||||
integer nz_ecp ! [input] number of prims/coeffs in stored data structure.
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer channel ! [input] Current spin channel used for ECPs
|
||||
integer channels ! [input] Total channels, used for print
|
||||
cc AJL/End
|
||||
integer n_prim_C(0:4,-1:l_ecp,n_ecp,2) ! [output]
|
||||
integer n_coef_C(-1:l_ecp,n_ecp,2) ! [output]
|
||||
integer ind_C(-1:l_ecp,n_ecp,2) ! [output]
|
||||
|
|
@ -603,6 +792,9 @@ c
|
|||
integer n0,n1,n2 ! r-exponent count
|
||||
integer n3,n4 ! r-exponent count
|
||||
integer ie ! ecp/so counter
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer mychan ! channel of contracted ECP
|
||||
cc AJL/End
|
||||
c
|
||||
#include "bas_exndcf_sfn.fh"
|
||||
#include "bas_ibs_sfn.fh"
|
||||
|
|
@ -666,26 +858,32 @@ c
|
|||
iexp = infbs_cont(CONT_IEXP,iucont,ecpid)
|
||||
icfp = infbs_cont(CONT_ICFP,iucont,ecpid)
|
||||
irexp = infbs_cont(CONT_IREXP,iucont,ecpid)
|
||||
if ((p_nprim+nprim).gt.nz_ecp) call errquit
|
||||
& ('int_ecp_build_ecp_ptrs:too many coefficients',911,
|
||||
& BASIS_ERR)
|
||||
call ecp_get_n3(
|
||||
& c_exp((p_nprim+1)),
|
||||
& dbl_mb(mb_exndcf(iexp,ecpid)),
|
||||
& c_coef((p_nprim+1)),
|
||||
& dbl_mb(mb_exndcf(icfp,ecpid)),
|
||||
& dbl_mb(mb_exndcf(irexp,ecpid)),
|
||||
& nprim,n0,n1,n2,n3,n4)
|
||||
n_prim_C(0,type,ie,1) = n0
|
||||
n_prim_C(1,type,ie,1) = n1
|
||||
n_prim_C(2,type,ie,1) = n2
|
||||
n_prim_C(3,type,ie,1) = n3
|
||||
n_prim_C(4,type,ie,1) = n4
|
||||
ind_C(type,ie,1) = p_nprim+1
|
||||
ind_Z(type,ie,1) = p_nprim+1
|
||||
n_coef_c(type,ie,1) = nprim
|
||||
l_C(ie) = max(type,l_C(ie))
|
||||
p_nprim = p_nprim+nprim
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
mychan = infbs_cont(CONT_CHANNEL,iucont,ecpid)
|
||||
if (mychan.eq.channel.or.mychan.eq.0) then
|
||||
cc AJL: If mychan = 0, applies the ECP to all channels
|
||||
if ((p_nprim+nprim).gt.nz_ecp) call errquit
|
||||
& ('int_ecp_build_ecp_ptrs:too many coefficients',911,
|
||||
& BASIS_ERR)
|
||||
call ecp_get_n3(
|
||||
& c_exp((p_nprim+1)),
|
||||
& dbl_mb(mb_exndcf(iexp,ecpid)),
|
||||
& c_coef((p_nprim+1)),
|
||||
& dbl_mb(mb_exndcf(icfp,ecpid)),
|
||||
& dbl_mb(mb_exndcf(irexp,ecpid)),
|
||||
& nprim,n0,n1,n2,n3,n4)
|
||||
n_prim_C(0,type,ie,1) = n0
|
||||
n_prim_C(1,type,ie,1) = n1
|
||||
n_prim_C(2,type,ie,1) = n2
|
||||
n_prim_C(3,type,ie,1) = n3
|
||||
n_prim_C(4,type,ie,1) = n4
|
||||
ind_C(type,ie,1) = p_nprim+1
|
||||
ind_Z(type,ie,1) = p_nprim+1
|
||||
n_coef_c(type,ie,1) = nprim
|
||||
l_C(ie) = max(type,l_C(ie))
|
||||
p_nprim = p_nprim+nprim
|
||||
endif
|
||||
cc AJL/End
|
||||
enddo
|
||||
endif
|
||||
if (on_so) then
|
||||
|
|
@ -706,7 +904,7 @@ c
|
|||
irexp = infbs_cont(CONT_IREXP,iucont,soid)
|
||||
if ((p_nprim+nprim).gt.nz_ecp) call errquit
|
||||
& ('int_ecp_build_ecp_ptrs:too many coefficients',911,
|
||||
& BASIS_ERR)
|
||||
& BASIS_ERR)
|
||||
call ecp_get_n3(
|
||||
& c_exp((p_nprim+1)),
|
||||
& dbl_mb(mb_exndcf(iexp,soid)),
|
||||
|
|
@ -728,19 +926,27 @@ c
|
|||
endif
|
||||
enddo
|
||||
c
|
||||
*debug: call print_ecp_ptrs(n_ecp,
|
||||
*debug: & l_ecp,nz_ecp,
|
||||
*debug: & n_prim_C,
|
||||
*debug: & n_coef_C,
|
||||
*debug: & ind_C,
|
||||
*debug: & ind_Z,
|
||||
*debug: & l_C,
|
||||
*debug: & c_exp,
|
||||
*debug: & c_coef)
|
||||
cc AJL/Begin/SPIN ECPS (Debug)
|
||||
c call print_ecp_ptrs(n_ecp,
|
||||
c & l_ecp,nz_ecp,
|
||||
c & channel, ! current channel [alpha/beta]
|
||||
c & channels, ! total number of channels
|
||||
c & n_prim_C,
|
||||
c & n_coef_C,
|
||||
c & ind_C,
|
||||
c & ind_Z,
|
||||
c & l_C,
|
||||
c & c_exp,
|
||||
c & c_coef)
|
||||
cc AJL/End
|
||||
c
|
||||
end
|
||||
subroutine print_ecp_ptrs(n_ecp,
|
||||
& l_ecp,nz_ecp,
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
& channel, ! current channel [alpha/beta]
|
||||
& channels, ! total number of channels
|
||||
cc AJL/End
|
||||
& n_prim_C,
|
||||
& n_coef_C,
|
||||
& ind_C,
|
||||
|
|
@ -753,13 +959,17 @@ c
|
|||
integer n_ecp ! [input] number of ecp centers
|
||||
integer l_ecp ! [input] maximum angular momentum in ecp basis
|
||||
integer nz_ecp ! [input] number of prims/coeffs in stored data structure.
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer channel ! [input] Current spin channel used for ECPs
|
||||
integer channels ! [input] Total number of spin channels used
|
||||
cc AJL/End
|
||||
integer n_prim_C(0:4,-1:l_ecp,n_ecp,2) ! [output]
|
||||
integer n_coef_C(-1:l_ecp,n_ecp,2) ! [output]
|
||||
integer ind_C(-1:l_ecp,n_ecp,2) ! [output]
|
||||
integer ind_Z(-1:l_ecp,n_ecp,2) ! [output]
|
||||
integer l_C(n_ecp) ! [output]
|
||||
double precision c_exp(nz_ecp) ! [output]
|
||||
double precision c_coef(nz_ecp) ! [output]
|
||||
integer l_C(n_ecp) ! [output]
|
||||
double precision c_exp(nz_ecp) ! [output]
|
||||
double precision c_coef(nz_ecp) ! [output]
|
||||
*
|
||||
integer i, j, k, l
|
||||
* integer low0, high0, low1, high1, low2, high2
|
||||
|
|
@ -767,6 +977,17 @@ c
|
|||
integer pe, pek
|
||||
*
|
||||
write(6,*)' print_ecp_ptrs: start'
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (channels.eq.1) then ! ECP variable for number of channels
|
||||
write(6,*) 'Both Channels'
|
||||
else if (channel.eq.1) then ! Spin-polarised check?
|
||||
write(6,*) 'Alpha Channel'
|
||||
else if (channel.eq.2) then
|
||||
write(6,*) 'Beta Channel'
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
write(6,*)' exponents and coefficients array'
|
||||
do i = 1,nz_ecp
|
||||
write(6,10000)i,c_exp(i),i,c_coef(i)
|
||||
|
|
@ -930,9 +1151,30 @@ c
|
|||
data k_ecp_c2s /0/ ! MA index for c2s routines
|
||||
data h_ecp_lip /0/ ! MA handle for ecp center index
|
||||
data k_ecp_lip /0/ ! MA index for ecp center index
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
data h_ecp_e_beta /0/ ! MA handle for ecp exponents
|
||||
data k_ecp_e_beta /0/ ! MA index for ecp exponents
|
||||
data h_ecp_c_beta /0/ ! MA handle for ecp coefficients
|
||||
data k_ecp_c_beta /0/ ! MA index for ecp coefficients
|
||||
data h_ecp_nprim_c_beta /0/ ! MA handle for n_prim_C (see ecp_integral)
|
||||
data k_ecp_nprim_c_beta /0/ ! MA index for n_prim_C (see ecp_integral)
|
||||
data h_ecp_ncoef_c_beta /0/ ! MA handle for n_coef_C (see ecp_integral)
|
||||
data k_ecp_ncoef_c_beta /0/ ! MA index for n_coef_C (see ecp_integral)
|
||||
data h_ecp_ind_c_beta /0/ ! MA handle for ind_C (see ecp_integral)
|
||||
data k_ecp_ind_c_beta /0/ ! MA index for ind_C (see ecp_integral)
|
||||
data h_ecp_ind_z_beta /0/ ! MA handle for ind_Z (see ecp_integral)
|
||||
data k_ecp_ind_z_beta /0/ ! MA index for ind_Z (see ecp_integral)
|
||||
data h_ecp_l_c_beta /0/ ! MA handle for l_C (see ecp_integral)
|
||||
data k_ecp_l_c_beta /0/ ! MA index for l_C (see ecp_integral)
|
||||
data h_ecp_lip_beta /0/ ! MA handle for ecp center index
|
||||
data k_ecp_lip_beta /0/ ! MA index for ecp center index
|
||||
cc AJL/End
|
||||
data n_zeta_c /0/ ! length of ecp exp/coef array
|
||||
data l_ecp /0/ ! high ang for ecp basis
|
||||
data n_ecp /0/ ! number of ecp centers (from API)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
data ecp_channels /0/ ! number of ecp channels [Alpha/Beta]
|
||||
cc AJL/End
|
||||
data init_ecp_init /.false./ ! logical saying if ecp is init-ed
|
||||
end
|
||||
subroutine int_ecp_hf1(
|
||||
|
|
@ -1247,3 +1489,170 @@ c
|
|||
& 0) ! ibug
|
||||
*
|
||||
end
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc All the subroutines designed for use with the beta ECPs
|
||||
|
||||
subroutine int_ecp_hf1_beta(
|
||||
& xyza,expa,coefa,a_nprim,a_ngen,La,ictra,
|
||||
& xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb,
|
||||
& ecp_ints,sz_ints,scr,lscr,dryrun)
|
||||
implicit none
|
||||
#include "mafdecls.fh"
|
||||
#include "ecp_nwc.fh"
|
||||
*
|
||||
integer a_nprim, a_ngen, La, ictra
|
||||
integer b_nprim, b_ngen, Lb, ictrb
|
||||
double precision expa(a_nprim), expb(b_nprim)
|
||||
double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen)
|
||||
integer sz_ints ! [input] buffer size for ecp_ints
|
||||
integer lscr ! [input] length of scratch array
|
||||
double precision xyza(3), xyzb(3) ! [input] a and b center coords.
|
||||
double precision ecp_ints(sz_ints) ! [output] ecp integrals
|
||||
double precision scr(lscr) ! [scratch] array
|
||||
logical dryrun ! [input] compute vs calculate memory requirements.
|
||||
*
|
||||
* write(6,*)' lscr IN ecp_hf1:',lscr
|
||||
* if (.not.dryrun) then
|
||||
* write(6,*)' int_ecp_hf1: coords a ',xyza
|
||||
* write(6,*)' int_ecp_hf1: coords b ',xyzb
|
||||
* endif
|
||||
c
|
||||
call ecp_integral(
|
||||
& xyza,
|
||||
& expa,
|
||||
& coefa,
|
||||
& a_nprim,a_ngen,La,ictra,
|
||||
& xyzb,
|
||||
& expb,
|
||||
& coefb,
|
||||
& b_nprim,b_ngen,Lb,ictrb,
|
||||
& dbl_mb(k_xyzecp),
|
||||
& dbl_mb(k_ecp_e_beta),
|
||||
& dbl_mb(k_ecp_c_beta),
|
||||
& int_mb(k_ecp_nprim_c_beta),
|
||||
& int_mb(k_ecp_ncoef_c_beta), ! new name is n_colc_C
|
||||
& int_mb(k_ecp_ind_z_beta),
|
||||
& int_mb(k_ecp_ind_c_beta),
|
||||
& n_zeta_c,
|
||||
& n_zeta_c,
|
||||
& int_mb(k_ecp_l_c_beta),
|
||||
& int_mb(k_ecp_lip_beta),
|
||||
& n_ecp,l_ecp,
|
||||
& 0,
|
||||
& dbl_mb(k_ecp_c2s),mem_c2s,
|
||||
& ecp_ints,sz_ints,1, ! nblk 1 for ecp integrals only
|
||||
& dryrun,scr,lscr,
|
||||
& 0) ! ibug
|
||||
*
|
||||
end
|
||||
subroutine intdd_ecp_hf1_beta(
|
||||
& xyza,expa,coefa,a_nprim,a_ngen,La,ictra,
|
||||
& xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb,
|
||||
& ecp_grad,sz_grad,nat,scr,lscr,dryrun)
|
||||
implicit none
|
||||
#include "mafdecls.fh"
|
||||
#include "ecp_nwc.fh"
|
||||
*
|
||||
integer a_nprim, a_ngen, La, ictra
|
||||
integer b_nprim, b_ngen, Lb, ictrb
|
||||
integer nat
|
||||
double precision expa(a_nprim), expb(b_nprim)
|
||||
double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen)
|
||||
integer sz_grad ! [input] buffer size for ecp_grad
|
||||
integer lscr ! [input] length of scratch array
|
||||
double precision xyza(3), xyzb(3) ! [input] a and b center coords.
|
||||
double precision ecp_grad(sz_grad,3,3,(nat*(nat-1)/2+nat))
|
||||
! [output] ecp integrals
|
||||
double precision scr(lscr) ! [scratch] array
|
||||
logical dryrun ! [input] compute vs calculate memory requirements.
|
||||
*
|
||||
* write(6,*)' lscr IN d_ecp_hf1:',lscr
|
||||
* if (.not.dryrun) then
|
||||
* write(6,*)' intd_ecp_hf1: coords a ',xyza
|
||||
* write(6,*)' intd_ecp_hf1: coords b ',xyzb
|
||||
* endif
|
||||
c
|
||||
call ecp_hessian(
|
||||
& xyza,
|
||||
& expa,
|
||||
& coefa,
|
||||
& a_nprim,a_ngen,La,ictra,
|
||||
& xyzb,
|
||||
& expb,
|
||||
& coefb,
|
||||
& b_nprim,b_ngen,Lb,ictrb,
|
||||
& dbl_mb(k_xyzecp),
|
||||
& dbl_mb(k_ecp_e_beta),
|
||||
& dbl_mb(k_ecp_c_beta),
|
||||
& int_mb(k_ecp_nprim_c_beta),
|
||||
& int_mb(k_ecp_ncoef_c_beta),
|
||||
& int_mb(k_ecp_ind_z_beta),
|
||||
& int_mb(k_ecp_ind_c_beta),
|
||||
& n_zeta_c,
|
||||
& n_zeta_c,
|
||||
& int_mb(k_ecp_l_c_beta),
|
||||
& int_mb(k_ecp_lip_beta),
|
||||
& n_ecp,l_ecp,
|
||||
& 0,
|
||||
& dbl_mb(k_ecp_c2s),mem_c2s,
|
||||
& ecp_grad,sz_grad,1,nat, ! nblock = 1 for ECP only
|
||||
& dryrun,scr,lscr,
|
||||
& 0) ! ibug
|
||||
*
|
||||
end
|
||||
subroutine intd_ecp_hf1_beta(
|
||||
& xyza,expa,coefa,a_nprim,a_ngen,La,ictra,
|
||||
& xyzb,expb,coefb,b_nprim,b_ngen,Lb,ictrb,
|
||||
& ecp_grad,sz_grad,nat,scr,lscr,dryrun)
|
||||
implicit none
|
||||
#include "mafdecls.fh"
|
||||
#include "ecp_nwc.fh"
|
||||
*
|
||||
integer a_nprim, a_ngen, La, ictra
|
||||
integer b_nprim, b_ngen, Lb, ictrb
|
||||
integer nat
|
||||
double precision expa(a_nprim), expb(b_nprim)
|
||||
double precision coefa(a_nprim,a_ngen), coefb(b_nprim,b_ngen)
|
||||
integer sz_grad ! [input] buffer size for ecp_grad
|
||||
integer lscr ! [input] length of scratch array
|
||||
double precision xyza(3), xyzb(3) ! [input] a and b center coords.
|
||||
double precision ecp_grad(sz_grad,3,nat) ! [output] ecp integrals
|
||||
double precision scr(lscr) ! [scratch] array
|
||||
logical dryrun ! [input] compute vs calculate memory requirements.
|
||||
*
|
||||
* write(6,*)' lscr IN d_ecp_hf1:',lscr
|
||||
* if (.not.dryrun) then
|
||||
* write(6,*)' intd_ecp_hf1: coords a ',xyza
|
||||
* write(6,*)' intd_ecp_hf1: coords b ',xyzb
|
||||
* endif
|
||||
c
|
||||
call ecp_gradient(
|
||||
& xyza,
|
||||
& expa,
|
||||
& coefa,
|
||||
& a_nprim,a_ngen,La,ictra,
|
||||
& xyzb,
|
||||
& expb,
|
||||
& coefb,
|
||||
& b_nprim,b_ngen,Lb,ictrb,
|
||||
& dbl_mb(k_xyzecp),
|
||||
& dbl_mb(k_ecp_e_beta),
|
||||
& dbl_mb(k_ecp_c_beta),
|
||||
& int_mb(k_ecp_nprim_c_beta),
|
||||
& int_mb(k_ecp_ncoef_c_beta),
|
||||
& int_mb(k_ecp_ind_z_beta),
|
||||
& int_mb(k_ecp_ind_c_beta),
|
||||
& n_zeta_c,
|
||||
& n_zeta_c,
|
||||
& int_mb(k_ecp_l_c_beta),
|
||||
& int_mb(k_ecp_lip_beta),
|
||||
& n_ecp,l_ecp,
|
||||
& 0,
|
||||
& dbl_mb(k_ecp_c2s),mem_c2s,
|
||||
& ecp_grad,sz_grad,1,nat, ! nblock = 1 for ECP only
|
||||
& dryrun,scr,lscr,
|
||||
& 0) ! ibug
|
||||
*
|
||||
end
|
||||
cc AJL/End
|
||||
|
|
|
|||
|
|
@ -61,6 +61,10 @@ c
|
|||
logical ecp_get_num_elec
|
||||
logical ecp_set_parent_handle
|
||||
logical ecp_get_parent_handle
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
logical ecp_set_high_chan
|
||||
logical ecp_get_high_chan
|
||||
cc AJL/End
|
||||
c
|
||||
logical bas_set_so_name
|
||||
logical bas_get_so_name
|
||||
|
|
@ -151,6 +155,9 @@ c
|
|||
external ecp_get_num_elec
|
||||
external ecp_set_parent_handle
|
||||
external ecp_get_parent_handle
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
external ecp_get_high_chan
|
||||
cc AJL/End
|
||||
c
|
||||
external bas_set_so_name
|
||||
external bas_get_so_name
|
||||
|
|
|
|||
|
|
@ -143,7 +143,7 @@ c ! 5 = icoeffpt ! pointer into linearized real*8
|
|||
c ! array for first coefficient of
|
||||
c ! each cont. on each tag
|
||||
c ! 6 = itag ! cunique tag number
|
||||
c ! 7 = irexp ! pointer to R exponents for ECP functions
|
||||
c ! 7 = irexp ! pointer to R exponents for ECP functions
|
||||
c Real*8 exndcf(3) ! linearized real*8 array of
|
||||
c ! exponents and coefficients
|
||||
c ! allocated on heap; exndcf now
|
||||
|
|
@ -174,7 +174,7 @@ c ! 4 = first cont ! first contraction on tag
|
|||
c ! 5 = last cont ! last contraction on tag
|
||||
c ! 6 = nelec ! number of elecs removed by ECP
|
||||
c ! 7 = high ang ! highest angular momentum on tag
|
||||
c Integer infbs_cont(7,num_cont_total,nbasis)
|
||||
c Integer infbs_cont(9,num_cont_total,nbasis)
|
||||
c ! 1 = itype ! type of contraction 0=s,1=p,2=d...
|
||||
c ! later -1=sp,-2=spd etc.
|
||||
c ! 2 = nprim ! num of prims in each cont. on tag
|
||||
|
|
@ -191,6 +191,10 @@ c ! 6 = tag number ! unique center lexical index
|
|||
c ! 7 = irexp ! pointer to R exponents for ECP functions
|
||||
c ! 8 = rel flag ! flag for relativistic shell.
|
||||
c ! 0 = nonrelativistic, 1 = relativistic.
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c ! 9 = channel ! Channel for ECP implementation
|
||||
c ! 0 = default; 1 = alpha; 2 = beta
|
||||
cc AJL/End
|
||||
c Real*8 exndcf(3,nbasis) ! linearized real*8 array of
|
||||
c ! exponents and coefficients
|
||||
c ! allocated on heap; exndcf now
|
||||
|
|
@ -214,7 +218,7 @@ c.. static dimension information for common blocks
|
|||
c
|
||||
c leading dimensions of compressed arrays
|
||||
integer ndbs_tags, ndbs_ucont, ndbs_head
|
||||
parameter(ndbs_tags=7, ndbs_ucont=8, ndbs_head=7)
|
||||
parameter(ndbs_tags=7, ndbs_ucont=9, ndbs_head=7)
|
||||
c.. stored structures
|
||||
character*256 bs_name
|
||||
character*16 bs_tags
|
||||
|
|
|
|||
|
|
@ -6,6 +6,42 @@
|
|||
*.....................................................................
|
||||
logical function bas_add_ucnt(basis, tag, l_value, ngen, nprim,
|
||||
$ rex, expnt, coeffs, ldc, stdtag, oshell_is_rel)
|
||||
cc AJL/Begin: This is a wrapper call to the updated version of this function
|
||||
cc SPIN-POLARISED ECPS
|
||||
implicit none
|
||||
c#include "errquit.fh"
|
||||
c#include "mafdecls.fh"
|
||||
c#include "basdeclsP.fh"
|
||||
c#include "nwc_const.fh"
|
||||
c#include "basP.fh"
|
||||
c#include "bas_exndcf_dec.fh"
|
||||
c#include "ecpso_decP.fh"
|
||||
c#include "stdio.fh"
|
||||
integer basis ! [input] basis handle
|
||||
character*(*) tag ! [input] tag on which to add contraction
|
||||
character*(*) stdtag ! [input] standard basis set name associated
|
||||
* with the tag
|
||||
integer l_value ! [input] type of contraction (s, p, ..., sp)
|
||||
integer ngen ! [input] no. of contractions
|
||||
integer nprim ! [input] no. of primitives
|
||||
integer ldc ! [input] leading dimension of coeffs
|
||||
double precision expnt(nprim) ! [input] exponents
|
||||
double precision coeffs(ldc, 1:*) ! [input] coefficients
|
||||
double precision rex(nprim) ! [input] gaussian R exponents
|
||||
*. . . . . . . . . . . . . . . . . (offset by 2)
|
||||
logical oshell_is_rel ! [input] flag for relativistic shell
|
||||
|
||||
logical bas_add_ucnt0
|
||||
external bas_add_ucnt0
|
||||
|
||||
bas_add_ucnt = bas_add_ucnt0(basis, tag, l_value,ngen, nprim,
|
||||
& rex, expnt, coeffs, ldc, stdtag, oshell_is_rel, 0)
|
||||
|
||||
end
|
||||
c
|
||||
logical function bas_add_ucnt0(basis, tag, l_value, ngen, nprim,
|
||||
$ rex, expnt, coeffs, ldc, stdtag, oshell_is_rel, channel)
|
||||
cc AJL/End
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
#include "mafdecls.fh"
|
||||
|
|
@ -43,6 +79,9 @@ c
|
|||
c
|
||||
logical bas_check_handle
|
||||
external bas_check_handle
|
||||
cc AJL/Begin/SPIN-POLARISED ECPS
|
||||
integer channel ! [input] ecp channel [both/alpha/beta]
|
||||
cc AJL/End
|
||||
c
|
||||
#include "bas_exndcf_sfn.fh"
|
||||
#include "ecpso_sfnP.fh"
|
||||
|
|
@ -67,14 +106,20 @@ c
|
|||
c adds a new general contraction on the specified tag. If the
|
||||
c tag is not present it will also add that by calling bas_add_utag
|
||||
c
|
||||
bas_add_ucnt = bas_check_handle(basis,'bas_add_ucnt')
|
||||
if (.not. bas_add_ucnt) return
|
||||
CC AJL/Begin/SPIN POLARISED ECPS
|
||||
c bas_add_ucnt = bas_check_handle(basis,'bas_add_ucnt')
|
||||
c if (.not. bas_add_ucnt) return
|
||||
bas_add_ucnt0 = bas_check_handle(basis,'bas_add_ucnt0')
|
||||
if (.not. bas_add_ucnt0) return
|
||||
ind = basis + BASIS_HANDLE_OFFSET
|
||||
c
|
||||
c Make sure that the tag is in the list
|
||||
c
|
||||
bas_add_ucnt = bas_add_utag(basis, tag, stdtag, itag)
|
||||
if (.not. bas_add_ucnt) return
|
||||
c bas_add_ucnt = bas_add_utag(basis, tag, stdtag, itag)
|
||||
c if (.not. bas_add_ucnt) return
|
||||
bas_add_ucnt0 = bas_add_utag(basis, tag, stdtag, itag)
|
||||
if (.not. bas_add_ucnt0) return
|
||||
CC AJL/End
|
||||
c
|
||||
c Update header information about all unique contractions on all
|
||||
c tags. Free points to next free word in the exndcf
|
||||
|
|
@ -91,12 +136,12 @@ c
|
|||
if (.not.ma_alloc_get(
|
||||
& mt_dbl,s_new,' input for basis heap ',
|
||||
& h_new, k_new)) then
|
||||
write(LuOut,*)'bas_add_ucnt: too many prims/coeffs'
|
||||
write(LuOut,*)'bas_add_ucnt0: too many prims/coeffs'
|
||||
write(LuOut,*)' allocated size for input is :',
|
||||
& exndcf(SZ_exndcf,ind)
|
||||
write(LuOut,*)' size requested here :',
|
||||
& (free+size_add-1)
|
||||
bas_add_ucnt = .false.
|
||||
bas_add_ucnt0 = .false.
|
||||
return
|
||||
endif
|
||||
call dfill(s_new,0.0d00,dbl_mb(k_new),1)
|
||||
|
|
@ -110,8 +155,8 @@ c
|
|||
endif
|
||||
* if (infbs_head(HEAD_NCONT,ind)+1 .gt. nucont_bsmx) then
|
||||
if (infbs_head(HEAD_NCONT,ind) .gt. nucont_bsmx) then
|
||||
write(LuOut,*) 'bas_add_ucnt: too many contractions '
|
||||
bas_add_ucnt = .false.
|
||||
write(LuOut,*) 'bas_add_ucnt0: too many contractions '
|
||||
bas_add_ucnt0 = .false.
|
||||
return
|
||||
endif
|
||||
c
|
||||
|
|
@ -157,7 +202,7 @@ c
|
|||
$ + nprim*ngen
|
||||
if (infbs_tags(TAG_FCONT,itag,ind).eq.0) then
|
||||
if (itag .ne. ntags) call errquit
|
||||
$ ('bas_add_ucnt: tag error', itag, BASIS_ERR)
|
||||
$ ('bas_add_ucnt0: tag error', itag, BASIS_ERR)
|
||||
infbs_tags(TAG_FCONT,itag,ind) = infbs_head(HEAD_NCONT,ind)
|
||||
infbs_tags(TAG_LCONT,itag,ind) = infbs_head(HEAD_NCONT,ind)
|
||||
else
|
||||
|
|
@ -173,11 +218,14 @@ c
|
|||
infbs_cont(CONT_NGEN, iu_cont,ind) = ngen
|
||||
infbs_cont(CONT_IEXP, iu_cont,ind) = free
|
||||
infbs_cont(CONT_ICFP, iu_cont,ind) = free + nprim
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
infbs_cont(CONT_CHANNEL, iu_cont, ind) = channel
|
||||
cc AJL/End
|
||||
if (oIs_ecp.or.oIs_so) then
|
||||
infbs_cont(Cont_Irexp, iu_cont, ind) =
|
||||
& free + nprim + nprim*ngen
|
||||
else
|
||||
*. . . . . . . . . . . . . . . . . . ! point to exponents for saftey?
|
||||
*. . . . . . . . . . . . . . . . . . ! point to exponents for safety?
|
||||
infbs_cont(Cont_Irexp, iu_cont,ind) = free
|
||||
endif
|
||||
if (oshell_is_rel) then
|
||||
|
|
|
|||
|
|
@ -27,7 +27,7 @@ c [spherical|cartesian] [segment||nosegment] [print] [noprint]
|
|||
c
|
||||
c tag nelec <integer>
|
||||
c tag library [tag in library] <standard set> [file <filename>]
|
||||
c tag <shell type>
|
||||
c tag <shell type> <channel>
|
||||
c <r-exponent> <exponent> <contraction coefficients>
|
||||
c ...
|
||||
c end ecp
|
||||
|
|
@ -343,7 +343,7 @@ c
|
|||
c for ecp or so
|
||||
c tag nelec <integer> ! only for ecp
|
||||
c tag library [<tag in library>] <standard set> [file <filename>]
|
||||
c tag <contraction type>
|
||||
c tag <contraction type> <channel>
|
||||
c <r-exponent> <exponent> <contraction coefficients>
|
||||
c ...
|
||||
c end ecp
|
||||
|
|
@ -369,6 +369,9 @@ c
|
|||
integer spvalues(nsptypes)
|
||||
integer num_elec
|
||||
integer l_value, ind, icount
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
integer channel ! Both (Default)=0; Alpha=1; Beta=2
|
||||
cc AJL/End
|
||||
data ltypes /'s','p','d','f','g','h','i'/
|
||||
data sptypes / 'sp', 'l ', 'ul'/
|
||||
data spvalues/ -1 , -1 , -1 /
|
||||
|
|
@ -430,6 +433,26 @@ c
|
|||
if (oIs_star) goto 10000
|
||||
c
|
||||
l_value = ind - 1
|
||||
|
||||
CC AJL/Begin/SPIN-POLARISED ECPs
|
||||
c
|
||||
c Check if it has a channel specified for spin-polarised ECPs
|
||||
c This is only implemented for user-defined values, not Library
|
||||
c
|
||||
if (oIs_ecp) then
|
||||
if (.not.inp_i(channel)) then
|
||||
channel = 0
|
||||
else if (channel.lt.0.or.channel.gt.2) then
|
||||
call errquit('bas_input_body: channel choice for ECP
|
||||
& must be either 0: Both; 1: Alpha (default); 2: Beta',
|
||||
& 99, INPUT_ERR)
|
||||
end if
|
||||
else
|
||||
channel = 0
|
||||
end if
|
||||
c
|
||||
CC AJL/End
|
||||
|
||||
c
|
||||
c check if it is a relativistic shell
|
||||
c
|
||||
|
|
@ -445,6 +468,25 @@ c
|
|||
c The contraction is an sp-type shell or an ecp Ulocal shell
|
||||
c
|
||||
l_value = spvalues(ind)
|
||||
|
||||
CC AJL/Begin
|
||||
c
|
||||
c check if it has a channel specified for spin-polarised ECPs
|
||||
c
|
||||
if (oIs_ecp) then
|
||||
if (.not.inp_i(channel)) then
|
||||
channel = 0
|
||||
else if (channel.lt.0.or.channel.gt.2) then
|
||||
call errquit('bas_input_body: channel choice for ECP
|
||||
& must be either 0: Both; 1: Alpha (default); 2: Beta',
|
||||
& 99, INPUT_ERR)
|
||||
end if
|
||||
else
|
||||
channel = 0
|
||||
end if
|
||||
c
|
||||
CC AJL/End
|
||||
|
||||
c
|
||||
c check for relativistic shell - this is illegal.
|
||||
c
|
||||
|
|
@ -559,7 +601,7 @@ c
|
|||
c Fall thru to here to read in a set of contraction coefficients
|
||||
c
|
||||
call bas_input_cont(basis, osegment, tag, l_value,
|
||||
& 'user specified', oshell_is_rel)
|
||||
& 'user specified', oshell_is_rel, channel)
|
||||
c
|
||||
c Have already read in the next line ... parse it
|
||||
c
|
||||
|
|
@ -582,7 +624,7 @@ c
|
|||
$ ' [except <exceptions list>] [rel]',/,
|
||||
$ ' tag nelec <number of electrons ',
|
||||
& 'replaced on tag>',/,
|
||||
$ ' tag <contraction type>',/,
|
||||
$ ' tag <contraction type> <channel>',/,
|
||||
$ ' <exponent> <contraction coefficients>',/,
|
||||
$ ' ... ',/,
|
||||
$ ' end ecp',/,/,/,'or',/,
|
||||
|
|
@ -599,7 +641,7 @@ c
|
|||
end
|
||||
*.....................................................................
|
||||
subroutine bas_input_cont(basis, osegment, tag, l_value, stdtag,
|
||||
& oshell_is_rel)
|
||||
& oshell_is_rel, channel)
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
#include "stdio.fh"
|
||||
|
|
@ -612,8 +654,8 @@ c
|
|||
#include "ecpso_decP.fh"
|
||||
#include "util.fh"
|
||||
c::functions
|
||||
logical bas_add_ucnt
|
||||
external bas_add_ucnt
|
||||
logical bas_add_ucnt0
|
||||
external bas_add_ucnt0
|
||||
c
|
||||
integer basis ! [input] basis set handle
|
||||
logical osegment ! [input] segment the basis set?
|
||||
|
|
@ -621,6 +663,9 @@ c
|
|||
integer l_value ! [input] l value of contraction
|
||||
character*(*) stdtag ! [input] standard name
|
||||
logical oshell_is_rel ! [input] flag for relativistic shell
|
||||
cc AJL/Begin/SPIN-POLARISED-ECPs
|
||||
integer channel ! [input] ecp channel [both/alpha/beta]
|
||||
cc AJL/End
|
||||
c
|
||||
integer ngen, iprim, nprim, i, j, nptmp
|
||||
integer cont_max ! Max. no. of contractions
|
||||
|
|
@ -733,6 +778,8 @@ c
|
|||
c
|
||||
c bas_add_ucnt -> adds a new general contraction on the specified
|
||||
c tag.
|
||||
c AJL/SPIN-POLARISED ECPs
|
||||
c bas_add_ucnt0 -> adds the contraction to a specified channel
|
||||
c
|
||||
c If the tag is not present it will also add that.
|
||||
c
|
||||
|
|
@ -772,39 +819,59 @@ c simple shell or all ecp/so shells
|
|||
endif
|
||||
endif
|
||||
enddo
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
c status = status .and.
|
||||
c $ bas_add_ucnt(basis, tag, l_value, 1, nptmp,
|
||||
c $ rex_tmp, etmp, ctmp, prim_max,
|
||||
c $ stdtag, oshell_is_rel)
|
||||
status = status .and.
|
||||
$ bas_add_ucnt(basis, tag, l_value, 1, nptmp,
|
||||
$ bas_add_ucnt0(basis, tag, l_value, 1, nptmp,
|
||||
$ rex_tmp, etmp, ctmp, prim_max,
|
||||
$ stdtag, oshell_is_rel)
|
||||
$ stdtag, oshell_is_rel, channel)
|
||||
cc AJL/End
|
||||
enddo
|
||||
else
|
||||
c sp shell
|
||||
cc AJL/Begin
|
||||
c status = status .and.
|
||||
c $ bas_add_ucnt(basis, tag, 0, 1, nprim,
|
||||
c $ rex, expnt, coeff(1,1), prim_max,
|
||||
c & stdtag, oshell_is_rel)
|
||||
c status = status .and.
|
||||
c $ bas_add_ucnt(basis, tag, 1, 1, nprim,
|
||||
c $ rex, expnt, coeff(1,2), prim_max,
|
||||
c & stdtag, oshell_is_rel)
|
||||
status = status .and.
|
||||
$ bas_add_ucnt(basis, tag, 0, 1, nprim,
|
||||
$ bas_add_ucnt0(basis, tag, 0, 1, nprim,
|
||||
$ rex, expnt, coeff(1,1), prim_max,
|
||||
& stdtag, oshell_is_rel)
|
||||
& stdtag, oshell_is_rel,channel) ! Added Channel
|
||||
status = status .and.
|
||||
$ bas_add_ucnt(basis, tag, 1, 1, nprim,
|
||||
$ bas_add_ucnt0(basis, tag, 1, 1, nprim,
|
||||
$ rex, expnt, coeff(1,2), prim_max,
|
||||
& stdtag, oshell_is_rel)
|
||||
& stdtag, oshell_is_rel,channel) ! Added Channel
|
||||
endif
|
||||
else
|
||||
c
|
||||
c Add as a single general contraction or sp shell
|
||||
c
|
||||
c status = status .and.
|
||||
c $ bas_add_ucnt(basis, tag, l_value, ngen, nprim,
|
||||
c $ rex, expnt, coeff, prim_max,
|
||||
c & stdtag, oshell_is_rel)
|
||||
status = status .and.
|
||||
$ bas_add_ucnt(basis, tag, l_value, ngen, nprim,
|
||||
$ bas_add_ucnt0(basis, tag, l_value, ngen, nprim,
|
||||
$ rex, expnt, coeff, prim_max,
|
||||
& stdtag, oshell_is_rel)
|
||||
& stdtag, oshell_is_rel, channel)
|
||||
endif
|
||||
cc AJL/End
|
||||
if (.not. status) call errquit
|
||||
$ ('bas_input_cont: bas_add_ucnt failed!!', 0, BASIS_ERR)
|
||||
$ ('bas_input_cont: bas_add_ucnt0 failed!!', 0, BASIS_ERR)
|
||||
c
|
||||
return
|
||||
c
|
||||
10000 write(LuOut,1)
|
||||
1 format(' basis/ecp contraction format is:',/,
|
||||
$ ' tag <contraction type> [rel]',/,
|
||||
$ ' tag <contraction type> [rel] [channel]',/,
|
||||
$ ' <exponent> <contraction coefficients>',/,
|
||||
$ ' ... ',/,
|
||||
$ ' end basis/ecp')
|
||||
|
|
@ -1085,7 +1152,13 @@ c
|
|||
c Have tag and l_value ... read in the contraction coeffcients
|
||||
c
|
||||
call bas_input_cont(basis, osegment, tag, l_value,
|
||||
& standard, oshell_is_rel)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c & standard, oshell_is_rel)
|
||||
& standard, oshell_is_rel, 0)
|
||||
c ^ Indicates that this applies to both
|
||||
c channels i.e. no spin-polarised ecps
|
||||
c when reading from libraries
|
||||
cc AJL/End
|
||||
goto 20
|
||||
c
|
||||
9000 close(33)
|
||||
|
|
|
|||
|
|
@ -30,12 +30,15 @@ c define CONT_ICFP 5
|
|||
c define CONT_TAG 6
|
||||
c define CONT_IREXP 7
|
||||
c define CONT_RELLS 8
|
||||
cc AJL/Begin/SPIN-POLAIRSED ECPs
|
||||
c define CONT_CHANNEL 9
|
||||
cc AJL/End
|
||||
integer HEAD_NTAGS, HEAD_NCONT, HEAD_NPRIM, HEAD_NCOEF,
|
||||
& HEAD_EXCFPTR, HEAD_SPH, HEAD_ECP
|
||||
integer TAG_NCONT, TAG_NPRIM, TAG_NCOEF, TAG_FCONT, TAG_LCONT,
|
||||
& TAG_NELEC, TAG_HIGH_ANG
|
||||
integer CONT_TYPE, CONT_NPRIM, CONT_NGEN, CONT_IEXP, CONT_ICFP,
|
||||
& CONT_TAG, CONT_IREXP, CONT_RELLS
|
||||
& CONT_TAG, CONT_IREXP, CONT_RELLS, CONT_CHANNEL
|
||||
parameter (HEAD_NTAGS = 1)
|
||||
parameter (HEAD_NCONT = 2)
|
||||
parameter (HEAD_NPRIM = 3)
|
||||
|
|
@ -60,6 +63,9 @@ c
|
|||
parameter (CONT_TAG = 6)
|
||||
parameter (CONT_IREXP = 7)
|
||||
parameter (CONT_RELLS = 8)
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
parameter (CONT_CHANNEL = 9)
|
||||
cc AJL/End
|
||||
c
|
||||
c bas_normalization specification values
|
||||
c
|
||||
|
|
|
|||
|
|
@ -1906,9 +1906,21 @@ c
|
|||
c::local
|
||||
integer basis ! Actual index into basis set arrays
|
||||
integer size_ex
|
||||
c
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
integer channels
|
||||
logical ecp_get_high_chan
|
||||
external ecp_get_high_chan
|
||||
cc AJL/End
|
||||
c
|
||||
c:: statement functions
|
||||
#include "errquit.fh"
|
||||
#include "bas_exndcf.fh"
|
||||
#include "ecpso_sfnP.fh"
|
||||
c
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
channels = 1
|
||||
cc AJL/End
|
||||
c
|
||||
bas_rtdb_store = bas_check_handle(basisin,'bas_rtdb_store')
|
||||
if (.not. bas_rtdb_store) return
|
||||
|
|
@ -1917,6 +1929,21 @@ c
|
|||
if (Is_ECP(basis).or.Is_SO(basis)) then
|
||||
size_ex = (2*infbs_head(HEAD_NPRIM,basis))+
|
||||
& infbs_head(HEAD_NCOEF,basis)
|
||||
c
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
cc Use this to put a note of spin_polarised_ecps in rtdb
|
||||
cc For use in nwdft/dft_fockbld.F
|
||||
cc
|
||||
cc By using this function we deal with anomalies for all
|
||||
cc applications of the spin-polarised ECPs being for both channels,
|
||||
cc as ecp_get_high_chan will treat that as spin-independent
|
||||
cc calculation of the fock matrices
|
||||
if (Is_ECP(basis)) then
|
||||
if (.not.ecp_get_high_chan(basisin,channels))
|
||||
& call errquit('bas_rtdb_store error',911, BASIS_ERR)
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
else
|
||||
size_ex = infbs_head(HEAD_NPRIM,basis)+
|
||||
& infbs_head(HEAD_NCOEF,basis)
|
||||
|
|
@ -1931,6 +1958,7 @@ c
|
|||
$ size_ex,
|
||||
& bs_stdname(1,basis),
|
||||
& name_assoc(1,basis),
|
||||
& channels, ! AJL/SPIN-POLARISED ECPs
|
||||
& star_nr_tags, star_tag, star_in_lib, star_bas_typ,
|
||||
& star_file, star_nr_excpt, star_excpt, star_tot_excpt,
|
||||
& star_rel, star_segment)
|
||||
|
|
@ -1940,7 +1968,7 @@ C> @}
|
|||
*.....................................................................
|
||||
logical function bas_rtdb_do_store(rtdb, name, tagsin,
|
||||
& head_array, tags_array, ucont_array, excfin, ntagsin,
|
||||
& nucontin, nexcf, stdnamein, ecp_name,
|
||||
& nucontin, nexcf, stdnamein, ecp_name, ecp_channels,
|
||||
& star_nr_tags, star_tag, star_in_lib, star_bas_typ,
|
||||
& star_file, star_nr_excpt, star_excpt, star_tot_excpt,
|
||||
& star_rel, star_segment)
|
||||
|
|
@ -1995,6 +2023,12 @@ c. . . . . . . . . . . . . . . . . . . contractions coeffs.
|
|||
integer star_tot_excpt ! [input] star total number of except tags
|
||||
integer star_nr_excpt(star_nr_tags) ! [input] number of except per tag
|
||||
character*16 star_excpt(star_nr_tags) ! [input] except tag list
|
||||
c
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
integer ecp_channels ! [input] number of channels for ECP (1 or 2).
|
||||
c............................ *Should* be 1 in all other scenarios
|
||||
cc AJL/End
|
||||
c
|
||||
c::: local
|
||||
character*256 tmp
|
||||
integer len_name, lentmp
|
||||
|
|
@ -2108,6 +2142,15 @@ c reset the stag tag data array for reuse
|
|||
c
|
||||
star_nr_tags = 0
|
||||
c
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
cc Add if we are using spin_polarised_ecps
|
||||
cc There may be a better place for this information
|
||||
if (ecp_channels.gt.1) then
|
||||
status = status .and. rtdb_put(
|
||||
& rtdb,'dft:spin_polarised_ecps',mt_int,1,ecp_channels)
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
c read the basis now get check status of read operations
|
||||
c
|
||||
if (.not.status) then
|
||||
|
|
@ -2173,6 +2216,67 @@ c
|
|||
return
|
||||
end
|
||||
C> @}
|
||||
*.....................................................................
|
||||
cc AJL/Begin
|
||||
C> \ingroup ecp
|
||||
C> @{
|
||||
c
|
||||
C> \brief Retrieve the number of channels in the ecp basis
|
||||
c
|
||||
C> \return Return .true. if successfull, and .false. otherwise
|
||||
c
|
||||
logical function ecp_get_high_chan(ecpidin,channels)
|
||||
implicit none
|
||||
c
|
||||
c calculate and return high channel for given ecp.
|
||||
c
|
||||
#include "nwc_const.fh"
|
||||
#include "basP.fh"
|
||||
#include "basdeclsP.fh"
|
||||
#include "stdio.fh"
|
||||
c::functions
|
||||
logical ecp_check_handle
|
||||
external ecp_check_handle
|
||||
c:: passed
|
||||
integer ecpidin !< [Input] the ECP set handle
|
||||
integer channels !< [Output] the number of channels
|
||||
c:local
|
||||
integer ecpid, i, j, mychannel
|
||||
integer myucont, mytags
|
||||
c
|
||||
ecpid = ecpidin + Basis_Handle_Offset
|
||||
c
|
||||
ecp_get_high_chan = ecp_check_handle(ecpidin,'ecp_get_high_chan')
|
||||
if (.not. ecp_get_high_chan) then
|
||||
write (6,*) 'This is not an ECP Basis'
|
||||
return
|
||||
endif
|
||||
c
|
||||
mytags = infbs_head(HEAD_NTAGS,ecpid)
|
||||
if (mytags.le.0) return
|
||||
c
|
||||
channels = 0
|
||||
do i=1,mytags
|
||||
myucont = infbs_tags(TAG_NCONT,i,ecpid)
|
||||
do j = 1,myucont
|
||||
mychannel = infbs_cont(CONT_CHANNEL,j,ecpid)
|
||||
channels = max(channels,mychannel)
|
||||
enddo
|
||||
enddo
|
||||
c If all values are 0 we have a spin independent calculation
|
||||
if (channels.eq.0) then
|
||||
channels = 1
|
||||
c Somewhere we have initialised spin in the alpha of beta channels
|
||||
else if (channels.gt.0) then
|
||||
channels = 2
|
||||
end if
|
||||
|
||||
c Save value process to be added
|
||||
c
|
||||
return
|
||||
end
|
||||
C> @}
|
||||
cc AJL/End
|
||||
*.....................................................................
|
||||
logical function gbs_map_clear(basisin)
|
||||
implicit none
|
||||
|
|
@ -3663,17 +3767,31 @@ c
|
|||
#include "geom.fh"
|
||||
#include "stdio.fh"
|
||||
#include "bas_starP.fh"
|
||||
#include "errquit.fh"
|
||||
c
|
||||
c function declarations
|
||||
c
|
||||
logical ecp_check_handle
|
||||
external ecp_check_handle
|
||||
c
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
logical ecp_get_high_chan
|
||||
external ecp_get_high_chan
|
||||
cc AJL/End
|
||||
c
|
||||
c:: passed
|
||||
integer ecpidin !< [Input] the basis set handle
|
||||
c:: local
|
||||
integer mytags, myucont, myprim, mycoef, ecpid
|
||||
integer i,j,k,l, ifcont, mygen, mytype, irexptr, iexptr, icfptr
|
||||
integer atn, len_tag, len_ele
|
||||
c
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
integer mychannel
|
||||
integer channels
|
||||
character*5 channel_type(0:2)
|
||||
cc AJL/End
|
||||
c
|
||||
character*2 symbol
|
||||
character*16 element
|
||||
character*9 cartorsph
|
||||
|
|
@ -3691,6 +3809,25 @@ c
|
|||
ctype(4) ='U-g'
|
||||
ctype(5) ='U-h'
|
||||
ctype(6) ='U-i'
|
||||
c
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
cc Aesthetic clean up for print. Is this necessary?
|
||||
C Can I just set the array without the if statement? To test.
|
||||
channels = 1
|
||||
if (.not.ecp_get_high_chan(ecpidin,channels))
|
||||
& call errquit('ecp_print error',911, BASIS_ERR)
|
||||
|
||||
if (channels.eq.1) then
|
||||
channel_type(0) = 'Both'
|
||||
channel_type(1) = ''
|
||||
channel_type(2) = ''
|
||||
else
|
||||
channel_type(0) = 'Both'
|
||||
channel_type(1) = 'Alpha'
|
||||
channel_type(2) = 'Beta'
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
ecp_print = .true.
|
||||
ecpid = ecpidin + Basis_Handle_Offset
|
||||
c
|
||||
|
|
@ -3743,8 +3880,12 @@ c
|
|||
ifcont = infbs_tags(TAG_FCONT,i,ecpid)
|
||||
c
|
||||
write(LuOut,6)
|
||||
6 format(
|
||||
$ ' R-exponent Exponent Coefficients '/
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
c 6 format(
|
||||
c $ ' R-exponent Exponent Coefficients '/
|
||||
6 format(13(' '),
|
||||
$ 'Channel R-exponent Exponent Coefficients'/
|
||||
cc AJL/End
|
||||
$ ' ------------ ',57('-'))
|
||||
do 00200 j=1,myucont
|
||||
myprim = infbs_cont(CONT_NPRIM,ifcont,ecpid)
|
||||
|
|
@ -3755,8 +3896,11 @@ c
|
|||
iexptr = infbs_cont(CONT_IEXP, ifcont,ecpid) - 1
|
||||
icfptr = infbs_cont(CONT_ICFP, ifcont,ecpid) - 1
|
||||
irexptr = infbs_cont(CONT_IREXP,ifcont,ecpid) - 1
|
||||
mychannel = infbs_cont(CONT_CHANNEL,ifcont,ecpid)
|
||||
c
|
||||
do 00300 k=1,myprim
|
||||
write(LuOut,7) j, shell_type,
|
||||
& channel_type(mychannel),
|
||||
& sf_exndcf((irexptr+k),ecpid),
|
||||
& sf_exndcf((iexptr+k),ecpid),
|
||||
& (sf_exndcf((icfptr+k+(l-1)*myprim),ecpid),l=1,mygen)
|
||||
|
|
@ -3794,7 +3938,10 @@ c
|
|||
c ... not done yet
|
||||
c
|
||||
return
|
||||
7 format(1x,i2,1x,a3,1x,f9.2,2x,f14.6,20f15.6)
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
c 7 format(1x,i2,1x,a3,1x,f9.2,2x,f14.6,20f15.6)
|
||||
7 format(1x,i2,1x,a3,7x,a5,3x,f9.2,2x,f14.6,20f15.6)
|
||||
cc AJL/End
|
||||
end
|
||||
*.....................................................................
|
||||
C>
|
||||
|
|
|
|||
|
|
@ -3,5 +3,8 @@
|
|||
*
|
||||
logical ocache
|
||||
integer bas1
|
||||
integer g_cache(3)
|
||||
common / cint1ecache / ocache, bas1, g_cache
|
||||
integer g_cache(4)
|
||||
c
|
||||
integer channels ! AJL: For the alpha/beta ECP channels
|
||||
c
|
||||
common / cint1ecache / ocache, bas1, g_cache, channels
|
||||
|
|
|
|||
|
|
@ -25,6 +25,11 @@ C> * "cos_chg_pot": the integrals with the COSMO charges
|
|||
C>
|
||||
C> * "bq_pot": the integrals with BQ-charges
|
||||
C>
|
||||
cc AJL/Begin
|
||||
C> * "potential_beta": the nuclear attraction with ECPs for the
|
||||
C> beta channel
|
||||
cc AJL/End
|
||||
C>
|
||||
C> The routine can exploit symmetry. If oskel is true then the
|
||||
C> "petite-list" symmetry will be used, although this requires both
|
||||
C> basis sets to be the same.
|
||||
|
|
@ -74,6 +79,11 @@ c
|
|||
if (inp_compare(.false., integ_type1, 'potential0')) then
|
||||
integ_type1='potential'
|
||||
dobq=.false.
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
elseif (inp_compare(.false., integ_type1, 'potential_beta0')) then
|
||||
integ_type1='potential_beta'
|
||||
dobq=.false.
|
||||
cc AJL/End
|
||||
end if
|
||||
c
|
||||
if (oskel) then
|
||||
|
|
@ -96,7 +106,11 @@ c
|
|||
type = 6
|
||||
else if (inp_compare(.false., integ_type1, 'bq_pot'))then
|
||||
if(.not.dobq) return
|
||||
else
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
else if (inp_compare(.false., integ_type1, 'potential_beta')) then
|
||||
type = 8
|
||||
cc AJL/End
|
||||
else
|
||||
write(6,*) ' integ_type1 = ', integ_type1
|
||||
call errquit('int_1e_ga: unknown integ_type', 0, INT_ERR)
|
||||
end if
|
||||
|
|
@ -105,7 +119,10 @@ c Check if Douglas-Kroll is required, if so:
|
|||
c getting the Douglas-Kroll kinetic and potential energy
|
||||
c integrals
|
||||
c
|
||||
if (doug_kroll .and. (type.le.2)) then
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
c if (doug_kroll .and. (type.le.2)) then
|
||||
if (doug_kroll .and. ((type.le.2) .or. (type.eq.8))) then
|
||||
cc AJL/End
|
||||
if (ibas.ne.jbas) call errquit
|
||||
$ ('int_1e_ga: Douglas-Kroll requires ibas=jbas', ibas,
|
||||
& BASIS_ERR)
|
||||
|
|
@ -127,8 +144,12 @@ c
|
|||
c
|
||||
if (dobq) then
|
||||
call int_1e_ooldga(ibas, jbas, g, integ_type1, oskel)
|
||||
if(inp_compare(.false., integ_type1, 'potential'))
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
cc if(inp_compare(.false., integ_type1, 'potential'))
|
||||
cc & call int_1e_ooldga(ibas, jbas, g, 'bq_pot', oskel)
|
||||
if (type.eq.2.or.type.eq.8) ! potential or potential_beta
|
||||
& call int_1e_ooldga(ibas, jbas, g, 'bq_pot', oskel)
|
||||
cc AJL/End
|
||||
else
|
||||
call int_1e_oldga(ibas, jbas, g, integ_type1, oskel)
|
||||
end if
|
||||
|
|
@ -219,7 +240,11 @@ c
|
|||
else if (inp_compare(.false., integ_type, 'cos_chg_pot'))then
|
||||
type = 6
|
||||
else if (inp_compare(.false., integ_type, 'bq_pot'))then
|
||||
type = 7
|
||||
type = 7
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
else if (inp_compare(.false., integ_type, 'potential_beta')) then
|
||||
type = 8
|
||||
cc AJL/End
|
||||
else
|
||||
write(6,*) ' integ_type = ', integ_type,ga_nodeid()
|
||||
call errquit('int_1e_ga: unknown integ_type', 0, INT_ERR)
|
||||
|
|
@ -245,7 +270,15 @@ c ocache = .false.
|
|||
* call util_flush(6)
|
||||
* endif
|
||||
* call ga_sync()
|
||||
call ga_dadd(1.0d0, g_cache(type), 1.0d0, g, g)
|
||||
c
|
||||
cc AJL/Begin/SPIN-POLARISED ECPs
|
||||
if (type.ne.8) then
|
||||
call ga_dadd(1.0d0, g_cache(type), 1.0d0, g, g)
|
||||
else ! accomodate for potential_beta, stored in g_cache(4)
|
||||
call ga_dadd(1.0d0, g_cache(4), 1.0d0, g, g)
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
if (oscfps) call pstat_off(ps_int_1e)
|
||||
return
|
||||
endif
|
||||
|
|
@ -335,6 +368,11 @@ c
|
|||
else if (type .eq. 7) then
|
||||
call int_1epot1 (jbas, jshell, ibas, ishell,
|
||||
$ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf))
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
else if (type .eq. 8) then
|
||||
call int_1epe_beta (jbas, jshell, ibas, ishell,
|
||||
$ mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf))
|
||||
cc AJL/End
|
||||
else
|
||||
call errquit('int_1e_ga: invalid type?', type,
|
||||
& GA_ERR)
|
||||
|
|
@ -403,6 +441,11 @@ c
|
|||
integer geom
|
||||
integer ga_create_atom_blocked
|
||||
external ga_create_atom_blocked
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
#include "global.fh"
|
||||
integer ecp_handle
|
||||
cc AJL/End
|
||||
c
|
||||
if (ocache) call int_1e_uncache_ga()
|
||||
c
|
||||
|
|
@ -423,6 +466,24 @@ c
|
|||
call int_1e_ga(basis, basis, g_cache(1), 'kinetic', oskel)
|
||||
call int_1e_ga(basis, basis, g_cache(2), 'potential0', oskel)
|
||||
call int_1e_ga(basis, basis, g_cache(3), 'overlap', .false.)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (bas_get_ecp_handle(basis,ecp_handle)) then
|
||||
if (.not.ecp_get_high_chan(ecp_handle,channels)) channels = 1
|
||||
if (channels.gt.1) then
|
||||
if (ga_nodeid().eq.0) then
|
||||
write(6,*) 'ECP Channels : ', channels
|
||||
call util_flush(6)
|
||||
endif
|
||||
g_cache(4) =
|
||||
& ga_create_atom_blocked(geom, basis,'int1e: v_beta')
|
||||
call ga_zero(g_cache(4))
|
||||
call int_1e_ga(basis, basis, g_cache(4),
|
||||
& 'potential_beta0', oskel)
|
||||
end if
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
bas1 = basis
|
||||
ocache = .true.
|
||||
c
|
||||
|
|
@ -442,6 +503,14 @@ c
|
|||
& BASIS_ERR)
|
||||
if (.not. ga_destroy(g_cache(3))) call errquit('i1ega: gad?',0,
|
||||
& BASIS_ERR)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (channels.gt.1) then
|
||||
if (.not. ga_destroy(g_cache(4))) call errquit('i1ega: gad?',0,
|
||||
& BASIS_ERR)
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
bas1 = -1
|
||||
ocache = .false.
|
||||
end
|
||||
|
|
@ -458,7 +527,10 @@ c
|
|||
logical oskel ! [input] If true generate symmetry unique list
|
||||
c
|
||||
external int_1eke,int_1epe,int_1eov,int_1epvpe,
|
||||
I int_1epot,int_1epot1
|
||||
I int_1epot,int_1epot1,
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
I int_1epe_beta
|
||||
cc AJL/End
|
||||
integer type
|
||||
c
|
||||
|
||||
|
|
@ -489,6 +561,11 @@ cso
|
|||
else if (inp_compare(.false., integ_type, 'bq_pot'))then
|
||||
type = 7
|
||||
call int_1e_oldga0(ibas, g, type,oskel,int_1epot1)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
else if (inp_compare(.false., integ_type, 'potential_beta')) then
|
||||
type = 8
|
||||
call int_1e_oldga0(ibas, g, type, oskel, int_1epe_beta)
|
||||
cc AJL/End/SPIN ECPs
|
||||
else
|
||||
write(6,*) ' integ_type = ', integ_type,ga_nodeid()
|
||||
call errquit('int_1e_oldga: unknown integ_type', 0, INT_ERR)
|
||||
|
|
@ -578,7 +655,13 @@ c ocache = .false.
|
|||
call util_flush(6)
|
||||
endif
|
||||
if (ocache .and. (ibas.eq.bas1)) then
|
||||
call ga_dadd(1.0d0, g_cache(type), 1.0d0, g, g)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (type.ne.8) then ! everything except potential_beta
|
||||
call ga_dadd(1.0d0, g_cache(type), 1.0d0, g, g)
|
||||
else
|
||||
call ga_dadd(1.0d0, g_cache(4), 1.0d0, g, g)
|
||||
end if
|
||||
cc AJL/End
|
||||
if (oscfps) call pstat_off(ps_int_1e)
|
||||
return
|
||||
endif
|
||||
|
|
|
|||
|
|
@ -665,7 +665,10 @@ c
|
|||
integer d(4), f(6), nfock
|
||||
integer g_a_dens, g_a_coul, g_a_exch, g_a_xc
|
||||
integer g_b_dens, g_b_coul, g_b_exch, g_b_xc
|
||||
integer g_hcore
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c integer g_hcore
|
||||
integer g_a_hcore, g_b_hcore
|
||||
c AJL/End
|
||||
integer g_tmp(6),ifock
|
||||
double precision jfac(4), kfac(4), one, zero, mone
|
||||
parameter (one=1.0d0, zero=0.0d0, mone=-1.0d0)
|
||||
|
|
@ -683,6 +686,15 @@ c DIM/QM JEM
|
|||
logical ldimqm
|
||||
integer g_vdim
|
||||
double precision edimqm
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_handle
|
||||
integer ecp_channels
|
||||
logical bas_get_ecp_handle
|
||||
external bas_get_ecp_handle
|
||||
logical ecp_get_high_chan
|
||||
external ecp_get_high_chan
|
||||
cc AJL/End
|
||||
c
|
||||
c Check
|
||||
c
|
||||
|
|
@ -896,13 +908,39 @@ c
|
|||
c
|
||||
c reuse g_a_exch to hold the 1-e integrals
|
||||
c
|
||||
g_hcore = g_a_exch
|
||||
call ga_zero(g_hcore)
|
||||
call int_1e_ga(basis, basis, g_hcore, 'kinetic', oskel) ! kinetic
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c g_hcore = g_a_exch
|
||||
c call ga_zero(g_hcore)
|
||||
|
||||
ecp_channels = 1
|
||||
if (bas_get_ecp_handle(basis,ecp_handle)) then
|
||||
if (ecp_get_high_chan(ecp_handle,ecp_channels))
|
||||
& continue
|
||||
end if
|
||||
|
||||
g_a_hcore = g_a_exch
|
||||
g_b_hcore = g_b_exch
|
||||
call ga_zero(g_a_hcore)
|
||||
call ga_zero(g_b_hcore)
|
||||
c
|
||||
call int_1e_ga(basis, basis, g_a_hcore, 'kinetic', oskel) ! kinetic
|
||||
if (do_zora .and. .not.(do_NonRel)) then
|
||||
call ga_dadd(1.d0,g_hcore,1.d0,g_zora_Kinetic(1),g_hcore) ! zora kinetic
|
||||
call ga_dadd(1.d0,g_a_hcore,1.d0,g_zora_Kinetic(1),g_a_hcore) ! zora kinetic
|
||||
endif
|
||||
call int_1e_ga(basis, basis, g_hcore, 'potential', oskel) ! potential
|
||||
|
||||
cc If spin polarised ECP, split g_hcore
|
||||
if (ecp_channels.gt.1)
|
||||
& call ga_copy(g_a_hcore,g_b_hcore)
|
||||
|
||||
cc Put normal potential in to alpha hcore
|
||||
c call int_1e_ga(basis, basis, g_hcore, 'potential', oskel) !potential
|
||||
call int_1e_ga(basis, basis, g_a_hcore, 'potential', oskel) !potential
|
||||
if (ecp_channels.gt.1) then ! Put beta potential in to beta hcore
|
||||
call int_1e_ga(basis, basis, g_b_hcore, 'potential_beta', oskel)!potential_beta
|
||||
else ! it is spin independent, so copy alpha to beta
|
||||
call ga_copy(g_a_hcore,g_b_hcore)
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
c cosmo charges, potential and energy contribution
|
||||
c
|
||||
|
|
@ -918,7 +956,10 @@ c
|
|||
& 2, ! 2 = open shell
|
||||
& g_dens, ! input density
|
||||
& cosmo_file) ! cosmo charges file name
|
||||
call int_1e_ga(basis,basis,g_hcore,'cos_chg_pot',.false.)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
call int_1e_ga(basis,basis,g_a_hcore,'cos_chg_pot',.false.)
|
||||
call int_1e_ga(basis,basis,g_b_hcore,'cos_chg_pot',.false.)
|
||||
cc AJL/End
|
||||
do i = 1, 2
|
||||
if (.not.ga_destroy(g_dens(i)))
|
||||
& call errquit('uhf_energy: ga_destroy failed g_dens',0,GA_ERR)
|
||||
|
|
@ -942,15 +983,18 @@ c We need to combine the alpha and beta densities to calculate the total D
|
|||
c
|
||||
c add in frozen embedding
|
||||
if (frozemb) then
|
||||
call ga_dadd(1.d0,g_hcore,1.d0,g_frozemb,g_hcore)
|
||||
cc AJL/Begin/Spin ECPs
|
||||
call ga_dadd(1.d0,g_a_hcore,1.d0,g_frozemb,g_a_hcore)
|
||||
call ga_dadd(1.d0,g_b_hcore,1.d0,g_frozemb,g_b_hcore)
|
||||
cc AJL/End
|
||||
end if
|
||||
c
|
||||
call ga_sync()
|
||||
c
|
||||
eone =
|
||||
$ (ga_ddot(g_a_dens,g_hcore) + ga_ddot(g_b_dens,g_hcore))
|
||||
call ga_dadd(one, g_hcore, one, g_a_coul, g_a_coul)
|
||||
call ga_dadd(one, g_hcore, one, g_b_coul, g_b_coul)
|
||||
$ (ga_ddot(g_a_dens,g_a_hcore) + ga_ddot(g_b_dens,g_b_hcore))
|
||||
call ga_dadd(one, g_a_hcore, one, g_a_coul, g_a_coul)
|
||||
call ga_dadd(one, g_b_hcore, one, g_b_coul, g_b_coul)
|
||||
c
|
||||
c DIM/QM JEM
|
||||
c Add DIM potential
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
subroutine grad1 ( H, lbuf, scr, lscr, dens, wdens, frc_nuc,
|
||||
$ frc_kin, frc_wgh, g_force,
|
||||
$ g_dens, g_wdens, basis, geom, nproc, nat,
|
||||
$ max_at_bf, rtdb, oskel )
|
||||
$ max_at_bf, rtdb, oskel, ndens )
|
||||
c$Id$
|
||||
|
||||
C one electron contribution to RHF, ROHF and UHF gradients
|
||||
|
|
@ -22,7 +22,10 @@ C-------------------------parameters--------------------------------
|
|||
$ g_dens, ! density matrix (summed if ROHF, UHF)
|
||||
$ g_wdens, ! weighted density (Lagrangian)
|
||||
$ g_force, ! global force array
|
||||
$ basis, geom, nproc, nat, max_at_bf, rtdb
|
||||
$ basis, geom, nproc, nat, max_at_bf, rtdb,
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
$ ndens ! Number of density matrices in g_dens
|
||||
cc AJL/End
|
||||
|
||||
double precision H, ! integral derivatives
|
||||
$ scr,
|
||||
|
|
@ -32,7 +35,11 @@ C-------------------------parameters--------------------------------
|
|||
|
||||
dimension H ( lbuf ), frc_nuc(3, nat), frc_kin(3, nat),
|
||||
$ frc_wgh(3, nat), scr(lscr),
|
||||
$ dens(max_at_bf,max_at_bf), wdens(max_at_bf,max_at_bf)
|
||||
$ dens(max_at_bf,max_at_bf),
|
||||
$ wdens(max_at_bf,max_at_bf),
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
$ g_dens(ndens)
|
||||
cc AJL/End
|
||||
|
||||
logical oskel ! symmetry?
|
||||
|
||||
|
|
@ -62,6 +69,35 @@ c
|
|||
integer nxtask, task_size
|
||||
external nxtask
|
||||
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_channels
|
||||
integer iecp
|
||||
double precision H_beta (lbuf)
|
||||
double precision dens_beta (max_at_bf,max_at_bf)
|
||||
|
||||
cc Read this value from rtdb vvvv
|
||||
if (.not.rtdb_get(rtdb, 'dft:spin_polarised_ecps',
|
||||
& mt_int, 1, ecp_channels)) then
|
||||
ecp_channels = 1
|
||||
endif
|
||||
|
||||
cc AJL: With spin-polarised ECPs Hcore will be spin dependent
|
||||
cc See Szabo and Ostlund pg. 215
|
||||
cc So we need to separate out the densities
|
||||
cc
|
||||
c if (ecp_channels.gt.1) then
|
||||
c
|
||||
c Restore alpha and beta densities to calculate spin-polarised
|
||||
c derivatives
|
||||
c
|
||||
c call ga_print(g_dens(1))
|
||||
c call ga_print(g_dens(2))
|
||||
c call ga_dadd(1d0, g_dens(1), -1d0, g_dens(2), g_dens(1))
|
||||
c call ga_print(g_dens(1))
|
||||
c call ga_print(g_dens(2))
|
||||
c end if
|
||||
cc AJL/End
|
||||
|
||||
task_size = 1
|
||||
status = rtdb_parallel(.true.) ! Broadcast reads to all processes
|
||||
|
||||
|
|
@ -96,10 +132,26 @@ c
|
|||
|
||||
status = bas_ce2cnr(basis,iat1,iac1f,iac1l)
|
||||
status = bas_ce2cnr(basis,iat2,iac2f,iac2l)
|
||||
|
||||
call ga_get (g_dens, iab1f,iab1l,iab2f,iab2l,dens,max_at_bf)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c call ga_get(g_dens,iab1f,iab1l,iab2f,iab2l,dens,max_at_bf)
|
||||
do 75, iecp = 1, ecp_channels
|
||||
if (iecp.eq.1) then
|
||||
call ga_get(g_dens(iecp),iab1f,iab1l,
|
||||
$ iab2f,iab2l,dens,max_at_bf)
|
||||
else
|
||||
call ga_get(g_dens(iecp),iab1f,iab1l,
|
||||
$ iab2f,iab2l,dens_beta,max_at_bf)
|
||||
endif
|
||||
75 continue
|
||||
cc Recombine g_dens, as it is not used again
|
||||
c if (ecp_channels.gt.1) then
|
||||
c call ga_dadd(1d0, g_dens(1), 1d0, g_dens(2), g_dens(1))
|
||||
c end if
|
||||
cc g_wdens is not dependent on spin, so can leave this
|
||||
call ga_get(g_wdens,iab1f,iab1l,iab2f,iab2l,wdens,max_at_bf)
|
||||
|
||||
cc AJL/End
|
||||
c
|
||||
do 70, ish1 = iac1f, iac1l
|
||||
if ( iat1.eq.iat2 ) iac2l = ish1
|
||||
do 60, ish2 = iac2f, iac2l
|
||||
|
|
@ -147,7 +199,23 @@ C 1el. derivatives
|
|||
call intd_1epot(basis,ish1,basis,ish2,lscr,scr,
|
||||
& lbuf,H)
|
||||
end if
|
||||
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc With spin-polarised ECPs Hcore will be spin dependent
|
||||
cc See Szabo and Ostlund pg. 215
|
||||
if (ecp_channels.gt.1) then
|
||||
C 1el. derivatives
|
||||
if(.not.dobq) then
|
||||
cc For now this will do, but this could be more efficiently done
|
||||
call intd_1eh1_beta(basis,ish1,basis,ish2,lscr,scr,
|
||||
& lbuf,H_beta)
|
||||
else
|
||||
call intd_1epot_beta(basis,ish1,basis,ish2,lscr,scr,
|
||||
& lbuf,H_beta)
|
||||
end if
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
C D x H
|
||||
|
||||
ic=1
|
||||
|
|
@ -157,6 +225,11 @@ C D x H
|
|||
do 31, ip1 = if1, il1
|
||||
do 30, ip2 = if2, il2
|
||||
dE = dE + dens(ip1,ip2) * H(ic)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (ecp_channels.gt.1) then
|
||||
dE = dE + dens_beta(ip1,ip2) * H_beta(ic)
|
||||
end if
|
||||
cc AJL/End
|
||||
ic = ic + 1
|
||||
30 continue
|
||||
31 continue
|
||||
|
|
|
|||
|
|
@ -175,6 +175,11 @@ c
|
|||
logical lcdft ! Qin
|
||||
integer ipop
|
||||
logical out1
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_channels
|
||||
logical osep
|
||||
cc AJL/End
|
||||
c
|
||||
nproc = ga_nnodes()
|
||||
c
|
||||
|
|
@ -524,11 +529,38 @@ c
|
|||
frac_occ = .false.
|
||||
if (has_frac_occ(rtdb)) frac_occ = .true.
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc These need separating in the case of the spin-polarised calculation
|
||||
if (.not.rtdb_get(rtdb, 'dft:spin_polarised_ecps',
|
||||
& mt_int, 1, ecp_channels)) then
|
||||
ecp_channels = 1
|
||||
endif
|
||||
|
||||
if (ecp_channels.eq.1) then
|
||||
osep = .false.
|
||||
else
|
||||
c Split density between alpha and beta.
|
||||
osep = .true.
|
||||
endif
|
||||
|
||||
c call grad_dens(geom, basis, g_dens, g_wdens, g_eigen_diag,
|
||||
c $ dbl_mb(k_evals), dbl_mb(k_occ), ndens, nbf, nbf2,
|
||||
c $ nopen, nclosed, scftype, movecs_in, omp2, odft,
|
||||
c $ dbl_mb(k_pdm2), dbl_mb(k_coeff),frac_occ,.false.)
|
||||
|
||||
call grad_dens(geom, basis, g_dens, g_wdens, g_eigen_diag,
|
||||
$ dbl_mb(k_evals), dbl_mb(k_occ), ndens, nbf, nbf2,
|
||||
$ nopen, nclosed, scftype, movecs_in, omp2, odft,
|
||||
$ dbl_mb(k_pdm2), dbl_mb(k_coeff),frac_occ,.false.)
|
||||
$ dbl_mb(k_pdm2), dbl_mb(k_coeff),frac_occ,osep)
|
||||
|
||||
c We will combine these, and then separate them when we need them
|
||||
c This is what is expected in all subroutines except the spin-polarised
|
||||
c ECPs, within grad1, and must be adhered to.
|
||||
|
||||
if (osep) call ga_dadd(1d0, g_dens(1), 1d0, g_dens(2), g_dens(1))
|
||||
|
||||
cc AJL/End
|
||||
c
|
||||
C free temporary arrays
|
||||
if (.not. ga_destroy(g_eigen_diag))
|
||||
$ call errquit('grad: could not destroy g_eigen_diag',1,
|
||||
|
|
@ -602,11 +634,27 @@ c Qin
|
|||
c Qin / end
|
||||
cpu_tim(1) = util_cpusec()
|
||||
wall_tim(1) = util_wallsec()
|
||||
c
|
||||
cc AJL/Begin
|
||||
cc Added ndens to grad1
|
||||
c
|
||||
if (ecp_channels.gt.1) then
|
||||
call ga_dadd(1d0, g_dens(1), -1d0, g_dens(2), g_dens(1))
|
||||
end if
|
||||
|
||||
call grad1( dbl_mb(k_buf), lbuf, dbl_mb(k_scr), lscratch,
|
||||
$ dbl_mb(k_dens), dbl_mb(k_wdens), dbl_mb(k_frc_nuc),
|
||||
$ dbl_mb(k_frc_kin), dbl_mb(k_frc_wgh), g_force,
|
||||
$ g_dens, g_wdens, basis, geom, nproc, nat, max_at_bf,
|
||||
$ rtdb, oskel)
|
||||
$ rtdb, oskel, ndens)
|
||||
|
||||
cc Recombine g_dens, as it is not used again
|
||||
if (ecp_channels.gt.1) then
|
||||
call ga_dadd(1d0, g_dens(1), 1d0, g_dens(2), g_dens(1))
|
||||
end if
|
||||
c
|
||||
cc AJL/End
|
||||
c
|
||||
c
|
||||
c ----- cosmo contribution -----
|
||||
c
|
||||
|
|
|
|||
|
|
@ -2,7 +2,7 @@ c
|
|||
c $Id$
|
||||
c
|
||||
subroutine oned_calc(nxyz,g_rhs,Ibuf,lbuf,Iscr,lscr,
|
||||
* skelsym,doS,doT,doV,geom,basis,nat)
|
||||
* skelsym,doS,doT,doV,geom,basis,nat,ecp_channel)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -35,6 +35,9 @@ c
|
|||
c
|
||||
integer nxtask
|
||||
external nxtask
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_channel
|
||||
cc AJL/End
|
||||
c
|
||||
nproc = ga_nnodes()
|
||||
task_size = 1
|
||||
|
|
@ -114,8 +117,19 @@ c
|
|||
C dH
|
||||
c
|
||||
if (doV.and.doT) then
|
||||
call intd_1eh1(basis,ish1,basis,ish2,lscr,Iscr,
|
||||
* lbuf,Ibuf)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c
|
||||
if (ecp_channel.eq.1) then
|
||||
call intd_1eh1(basis,ish1,basis,ish2,lscr,Iscr,
|
||||
* lbuf,Ibuf)
|
||||
else
|
||||
call intd_1eh1_beta(basis,ish1,basis,ish2,lscr,Iscr,
|
||||
* lbuf,Ibuf)
|
||||
end if
|
||||
c
|
||||
cc AJL/End
|
||||
c
|
||||
ic = 1
|
||||
do atx = 1, nxyz
|
||||
lo(1) = atx
|
||||
|
|
|
|||
|
|
@ -3,7 +3,8 @@ c $Id$
|
|||
c
|
||||
subroutine onedd_calc(Hxy,nxyz,Ibuf,lbuf,Iscr,lscr,
|
||||
* skelsym,doS,doT,doV,doECP,
|
||||
* g_wdens,g_dens,geom,basis,nat,ndens,dens,wdens,matb)
|
||||
* g_wdens,g_dens,geom,basis,nat,ndens,dens,wdens,matb,
|
||||
* ecp_channels)
|
||||
c
|
||||
implicit none
|
||||
c
|
||||
|
|
@ -43,6 +44,12 @@ c
|
|||
c
|
||||
integer nxtask
|
||||
external nxtask
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_channels
|
||||
cc Buffer for Beta Channel ECP integrals
|
||||
double precision Ibuf_beta(lbuf)
|
||||
cc AJL/End
|
||||
c
|
||||
nproc = ga_nnodes()
|
||||
task_size = 1
|
||||
|
|
@ -80,7 +87,15 @@ c
|
|||
if (ndens.eq.2) then
|
||||
call ga_get (g_dens(2), iab1f,iab1l,iab2f,iab2l,
|
||||
* dens(1,1,2),matb)
|
||||
call daxpy(matb*matb,1.0d+00,dens(1,1,2),1,dens,1)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc We only want to sum this if ecp_channels are one
|
||||
cc Otherwise we are going to apply different ddH to each density channel
|
||||
if (ecp_channels.eq.1) then
|
||||
call daxpy(matb*matb,1.0d+00,dens(1,1,2),1,dens,1)
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
endif
|
||||
call ga_get(g_wdens,iab1f,iab1l,iab2f,iab2l,wdens,matb)
|
||||
|
||||
|
|
@ -154,6 +169,14 @@ c
|
|||
do ip1 = if1, il1
|
||||
do ip2 = if2, il2
|
||||
dE = dE + dens(ip1,ip2,1)*Ibuf(ic)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc We haven't summed the local density so need to add this
|
||||
if (ecp_channels.gt.1) then
|
||||
dE = dE + dens(ip1,ip2,2)*Ibuf(ic)
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
ic = ic + 1
|
||||
enddo
|
||||
enddo
|
||||
|
|
@ -180,6 +203,14 @@ c
|
|||
do 31 ip1 = if1, il1
|
||||
do 30 ip2 = if2, il2
|
||||
dE = dE + dens(ip1,ip2,1) * Ibuf(ic)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc We haven't summed the local density so need to add this
|
||||
if (ecp_channels.gt.1) then
|
||||
dE = dE + dens(ip1,ip2,2)*Ibuf(ic)
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
c write(6,*) 'dens and int ',dens(ip1,ip2,1),
|
||||
c * Ibuf(ic)
|
||||
ic = ic + 1
|
||||
|
|
@ -216,6 +247,14 @@ c
|
|||
do ip1 = if1, il1
|
||||
do ip2 = if2, il2
|
||||
dE = dE + dens(ip1,ip2,1) * Ibuf(ic)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc We haven't summed the local density so need to add this
|
||||
if (ecp_channels.gt.1) then
|
||||
dE = dE + dens(ip1,ip2,2)*Ibuf(ic)
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
ic = ic + 1
|
||||
enddo
|
||||
enddo
|
||||
|
|
@ -242,6 +281,15 @@ c
|
|||
call intdd_1eecp1(basis,ish1,basis,ish2,
|
||||
& lscr,Iscr,
|
||||
& lbuf,Ibuf)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (ecp_channels.gt.1) then
|
||||
call intdd_1eecp1_beta(basis,ish1,basis,ish2,
|
||||
& lscr,Iscr,
|
||||
& lbuf,Ibuf_beta)
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
ic = 0
|
||||
do at2 = 1, nat
|
||||
do at1 = 1, at2
|
||||
|
|
@ -254,6 +302,14 @@ c do ip1 = if1, il1
|
|||
do ip2 = if2, il2
|
||||
ic = ic + 1
|
||||
dE = dE + dens(ip1,ip2,1) * Ibuf(ic)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc Now we compute the Hecp_beta x Dens_beta
|
||||
if (ecp_channels.gt.1) then
|
||||
dE = dE + dens(ip1,ip2,2) * Ibuf_beta(ic)
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
enddo
|
||||
enddo
|
||||
if (iat1.ne.iat2 .or. ish1.ne.ish2) dE = dE + dE
|
||||
|
|
|
|||
|
|
@ -34,6 +34,16 @@ c
|
|||
integer ga_create_atom_blocked
|
||||
external ga_create_atom_blocked
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_channels
|
||||
cc Read this value from rtdb vvvv
|
||||
if (.not.rtdb_get(rtdb, 'dft:spin_polarised_ecps',
|
||||
& mt_int, 1, ecp_channels)) then
|
||||
ecp_channels = 1
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
c
|
||||
c Set up print info
|
||||
c
|
||||
oprint = util_print('hess_follow',print_high)
|
||||
|
|
@ -107,7 +117,10 @@ c
|
|||
call onedd_calc(dbl_mb(k_exy),n3xyz,dbl_mb(k_buf),
|
||||
* max1e,dbl_mb(k_scr),mscratch_1e,skelsym,
|
||||
* doS,doT,doV,doECP,g_wdens,g_dens,geom,basis,ncent,ndens,
|
||||
* dbl_mb(k_dens),dbl_mb(k_wdens),matb)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c * dbl_mb(k_dens),dbl_mb(k_wdens),matb)
|
||||
* dbl_mb(k_dens),dbl_mb(k_wdens),matb,ecp_channels)
|
||||
cc AJL/End
|
||||
c
|
||||
if(skelsym)call sym_hess_symmetrize(geom,dbl_mb(k_exy))
|
||||
if (oneprint.and.(ga_nodeid().eq.0)) then
|
||||
|
|
@ -170,11 +183,31 @@ c
|
|||
doS = .false.
|
||||
doV = .true.
|
||||
doT = .true.
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if (ndens.gt.1.and.ecp_channels.gt.1) then
|
||||
call oned_calc(n3xyz,g_rhs,dbl_mb(k_buf),
|
||||
* max1e,dbl_mb(k_scr),mscratch_1e,.false.,
|
||||
* doS,doT,doV,geom,basis,ncent,2)
|
||||
C ^ ECP Channel, Beta
|
||||
cc We should be OK just copying, as the commands in oned_calc are puts
|
||||
cc Is this necessary at all? To be tested.
|
||||
call nga_copy_patch('N', g_rhs, alo, ahi, g_rhs, blo, bhi)
|
||||
end if
|
||||
c
|
||||
call oned_calc(n3xyz,g_rhs,dbl_mb(k_buf),
|
||||
* max1e,dbl_mb(k_scr),mscratch_1e,.false.,
|
||||
* doS,doT,doV,geom,basis,ncent)
|
||||
if (ndens.gt.1) call nga_copy_patch('N', g_rhs,
|
||||
* alo, ahi, g_rhs, blo, bhi)
|
||||
* max1e,dbl_mb(k_scr),mscratch_1e,.false.,
|
||||
* doS,doT,doV,geom,basis,ncent,1)
|
||||
C ^ ECP Channel, Alpha
|
||||
c
|
||||
c if (ndens.gt.1) call nga_copy_patch('N', g_rhs,
|
||||
c * alo, ahi, g_rhs, blo, bhi)
|
||||
c
|
||||
if (ndens.gt.1.and.ecp_channels.eq.1) then
|
||||
call nga_copy_patch('N', g_rhs, alo, ahi, g_rhs, blo, bhi)
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
c
|
||||
c Print out the rhs values if requested
|
||||
c
|
||||
|
|
|
|||
|
|
@ -84,6 +84,10 @@ c DIM/QM variables
|
|||
c
|
||||
double precision g_vdim
|
||||
double precision Edim
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_channels
|
||||
cc AJL/End
|
||||
c
|
||||
do_cosmo_smd = .false.
|
||||
do_cosmo_vem = 0
|
||||
|
|
@ -165,25 +169,63 @@ c
|
|||
& ' Kinetic energy =', f20.12/)
|
||||
endif
|
||||
endif
|
||||
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc Read this value from rtdb vvvv
|
||||
if (.not.rtdb_get(rtdb, 'dft:spin_polarised_ecps',
|
||||
& mt_int, 1, ecp_channels)) then
|
||||
ecp_channels = 1
|
||||
endif
|
||||
|
||||
if (ecp_channels.eq.2.and.ipol.eq.1)
|
||||
& call errquit('dft_fockbld: spin_polarised_ecps with mult = 1',
|
||||
& 0, BASIS_ERR)
|
||||
|
||||
c == Copy Fock matrix to individual components for spin ECPs ==
|
||||
call ga_copy(g_fock, g_focks(1))
|
||||
c Is this redundant? Need to test
|
||||
if (ecp_channels.gt.1) call ga_copy(g_fock, g_focks(2))
|
||||
c
|
||||
c if requested separate calculation of bq-electron energy (MV)
|
||||
c -----------------------------------------------------------
|
||||
if (rtdb_get(rtdb, 'dft:bq_energy', mt_dbl, 1,Ebq)) then
|
||||
Ecore = ga_ddot(g_dens(1), g_fock)
|
||||
call int_1e_ga(ao_bas_han, ao_bas_han, g_fock,
|
||||
c Ecore = ga_ddot(g_dens(1), g_fock)
|
||||
c call int_1e_ga(ao_bas_han, ao_bas_han, g_fock,
|
||||
c & 'bq_pot', oskel)
|
||||
|
||||
Ecore = ga_ddot(g_dens(1), g_focks(1))
|
||||
call int_1e_ga(ao_bas_han, ao_bas_han, g_focks(1),
|
||||
& 'bq_pot', oskel)
|
||||
Ebq = ga_ddot(g_dens(1), g_fock)-Ecore
|
||||
call int_1e_ga(ao_bas_han, ao_bas_han, g_fock,
|
||||
c Copy to g_focks(2), as this is identical for both channels
|
||||
if (ecp_channels.gt.1) call ga_copy(g_focks(1), g_focks(2))
|
||||
c
|
||||
c Ebq = ga_ddot(g_dens(1), g_fock)-Ecore
|
||||
c call int_1e_ga(ao_bas_han, ao_bas_han, g_fock,
|
||||
c & 'potential0', oskel)
|
||||
|
||||
Ebq = ga_ddot(g_dens(1), g_focks(1))-Ecore
|
||||
call int_1e_ga(ao_bas_han, ao_bas_han, g_focks(1),
|
||||
& 'potential0', oskel)
|
||||
c Calculate Beta channel 1e potential
|
||||
if (ecp_channels.gt.1)
|
||||
& call int_1e_ga(ao_bas_han, ao_bas_han, g_focks(2),
|
||||
& 'potential_beta0', oskel)
|
||||
cc AJL/End
|
||||
c
|
||||
if (.not. rtdb_put(rtdb,'dft:bq_energy',mt_dbl,1,Ebq))
|
||||
$ call errquit('dft: writing bq_energy failed',
|
||||
$ 0, RTDB_ERR)
|
||||
else
|
||||
if(util_print('dft timings', print_high))
|
||||
& time1_pot=util_cpusec() ! start time
|
||||
call int_1e_ga(ao_bas_han, ao_bas_han, g_fock,
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c call int_1e_ga(ao_bas_han, ao_bas_han, g_fock,
|
||||
call int_1e_ga(ao_bas_han, ao_bas_han, g_focks(1),
|
||||
& 'potential', oskel)
|
||||
if (ecp_channels.gt.1)
|
||||
& call int_1e_ga(ao_bas_han, ao_bas_han, g_focks(2),
|
||||
& 'potential_beta', oskel)
|
||||
cc AJL/End
|
||||
if(util_print('dft timings', print_high))
|
||||
& time2_pot=util_cpusec() ! end time
|
||||
c
|
||||
|
|
@ -198,13 +240,50 @@ c print potential matrix time
|
|||
c
|
||||
c ----- Add cosmo potential when needed -----
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
if(cosmo_on.and.cosmo_phase.eq.2) then
|
||||
call int_1e_ga(ao_bas_han,ao_bas_han,g_fock,
|
||||
c call int_1e_ga(ao_bas_han,ao_bas_han,g_fock,
|
||||
call int_1e_ga(ao_bas_han,ao_bas_han,g_focks(1),
|
||||
& 'cos_chg_pot',.false.)
|
||||
if (ecp_channels.gt.1)
|
||||
& call int_1e_ga(ao_bas_han,ao_bas_han,g_focks(2),
|
||||
& 'cos_chg_pot',.false.)
|
||||
endif
|
||||
|
||||
c == Copy Fock matrix to individual components for normal calcs ==
|
||||
if (ecp_channels.eq.1.and.ipol.gt.1)
|
||||
& call ga_copy(g_focks(1), g_focks(2))
|
||||
|
||||
call ga_sync()
|
||||
c
|
||||
Ecore = ga_ddot(g_dens(1), g_fock)
|
||||
cc Ecore = ga_ddot(g_dens(1), g_fock)
|
||||
cc AJL: Ecore = ga_ddot(g_dens(1), g_focks(1))
|
||||
cc write(6,*)'Ecore(1:T):',ga_ddot(g_dens(1), g_focks(1))
|
||||
|
||||
cc With spin-polarised ECPs Ecore will be spin dependent
|
||||
cc See Szabo and Ostlund pg. 215
|
||||
cc
|
||||
if (ecp_channels.gt.1) then
|
||||
c
|
||||
c Restore alpha and beta densities to calculate energy
|
||||
call ga_dadd(1d0, g_dens(1), -1d0, g_dens(2), g_dens(1))
|
||||
|
||||
Ecore = ga_ddot(g_dens(1), g_focks(1))
|
||||
c write(6,*)'Ecore(1:A):',Ecore
|
||||
Ecore = Ecore + ga_ddot(g_dens(2), g_focks(2))
|
||||
c write(6,*)'Ecore(2:B):',ga_ddot(g_dens(2), g_focks(2))
|
||||
c
|
||||
c And then recombine the densities for now.
|
||||
call ga_dadd(1d0, g_dens(1), 1d0, g_dens(2), g_dens(1))
|
||||
c
|
||||
c call ga_sync()
|
||||
else
|
||||
|
||||
Ecore = ga_ddot(g_dens(1), g_focks(1))
|
||||
c write(6,*)'Ecore(1:T):',ga_ddot(g_dens(1), g_focks(1))
|
||||
|
||||
endif
|
||||
cc AJL/End
|
||||
c
|
||||
c add Bonacic-Fantucci repulsive term
|
||||
c
|
||||
|
|
@ -258,7 +337,13 @@ c
|
|||
c
|
||||
c Add V coul to Fock Matrix
|
||||
c
|
||||
call ga_dadd(1d0, g_tmp, 1d0, g_fock, g_fock)
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
c call ga_dadd(1d0, g_tmp, 1d0, g_fock, g_fock)
|
||||
call ga_dadd(1d0, g_tmp, 1d0, g_focks(1), g_focks(1))
|
||||
if (ipol.gt.1)
|
||||
& call ga_dadd(1d0, g_tmp, 1d0, g_focks(2), g_focks(2))
|
||||
cc AJL/End
|
||||
c
|
||||
if (CDFIT)then
|
||||
if (.not.ma_pop_stack(lcd_coef))
|
||||
& call errquit('dft_scf: cannot pop stack',0, MA_ERR)
|
||||
|
|
@ -310,7 +395,11 @@ c == Compute the XC potential and energy ==
|
|||
g_vxc(1) = g_tmp
|
||||
call ga_zero(g_vxc(1))
|
||||
rho_n = 0.0d0
|
||||
if (ipol.eq.2)call ga_zero(g_vxc(2))
|
||||
cc AJL/Begin/FDE
|
||||
cc Can't remember why I implemented this? To be tested.
|
||||
cc if (ipol.eq.2)call ga_zero(g_vxc(2))
|
||||
if (ipol.gt.1) call ga_copy(g_vxc(1), g_vxc(2))
|
||||
cc AJL/End
|
||||
if (oprint_time)call dft_tstamp(' Before call to GETVXC. ')
|
||||
c
|
||||
if (odftps) call pstat_on(ps_xc)
|
||||
|
|
@ -378,10 +467,12 @@ c Jorge Garza/Begin
|
|||
end if
|
||||
c Jorge Garza/End
|
||||
c
|
||||
c == Copy Fock matrix to individual components ==
|
||||
call ga_copy(g_fock, g_focks(1))
|
||||
if (ipol.gt.1)
|
||||
& call ga_copy(g_fock, g_focks(2))
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc == Copy Fock matrix to individual components ==
|
||||
c call ga_copy(g_fock, g_focks(1))
|
||||
c if (ipol.gt.1)
|
||||
c & call ga_copy(g_fock, g_focks(2))
|
||||
cc AJL/End
|
||||
c
|
||||
c == Add exchange-correlation to the Fock matrix ==
|
||||
call ga_dadd(1.d0, g_vxc(1), 1.d0, g_focks(1), g_focks(1))
|
||||
|
|
@ -403,6 +494,18 @@ c DIM potential is from total density, so same for both alpha and beta
|
|||
c
|
||||
if (odftps) call pstat_off(ps_fockb)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs (Debug)
|
||||
c if(util_print('Intermediate Fock Matrix', print_debug)) then
|
||||
c if (ga_nodeid() .eq. 0) then
|
||||
c write(6,*)
|
||||
c call util_print_centered(6,
|
||||
c $ 'Intermediate Fock Matrix',40,.true.)
|
||||
c write(6,*)
|
||||
c call util_flush(6)
|
||||
c end if
|
||||
c call ga_print(g_focks(1))
|
||||
c if (ipol.eq.2) call ga_print(g_focks(2))
|
||||
cc AJL/End
|
||||
c
|
||||
c MN solvation models -->
|
||||
c
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
subroutine tce_grad1 ( H, lbuf, scr, lscr, dens, wdens, frc_nuc,
|
||||
$ frc_kin, frc_wgh, g_force,
|
||||
$ g_dens, g_wdens, basis, geom, nproc, nat,
|
||||
$ max_at_bf, rtdb, oskel )
|
||||
$ max_at_bf, rtdb, oskel, ndens )
|
||||
c$Id$
|
||||
|
||||
C one electron contribution to RHF, ROHF and UHF gradients
|
||||
|
|
@ -15,13 +15,17 @@ C now also UMP2
|
|||
#include "bas.fh"
|
||||
#include "rtdb.fh"
|
||||
#include "sym.fh"
|
||||
#include "stdio.fh"
|
||||
|
||||
C-------------------------parameters--------------------------------
|
||||
integer lbuf, lscr,
|
||||
$ g_dens, ! density matrix (summed if ROHF, UHF)
|
||||
$ g_wdens, ! weighted density (Lagrangian)
|
||||
$ g_force, ! global force array
|
||||
$ basis, geom, nproc, nat, max_at_bf, rtdb
|
||||
$ basis, geom, nproc, nat, max_at_bf, rtdb,
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
$ ndens ! Number of density matrices in g_dens
|
||||
cc AJL/End
|
||||
|
||||
double precision H, ! integral derivatives
|
||||
$ scr,
|
||||
|
|
@ -31,7 +35,10 @@ C-------------------------parameters--------------------------------
|
|||
|
||||
dimension H ( lbuf ), frc_nuc(3, nat), frc_kin(3, nat),
|
||||
$ frc_wgh(3, nat), scr(lscr),
|
||||
$ dens(max_at_bf,max_at_bf), wdens(max_at_bf,max_at_bf)
|
||||
$ dens(max_at_bf,max_at_bf), wdens(max_at_bf,max_at_bf),
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
$ g_dens(ndens)
|
||||
cc AJL/End
|
||||
|
||||
logical oskel ! symmetry?
|
||||
|
||||
|
|
@ -60,6 +67,50 @@ C-------------------------local variables--------------------------
|
|||
integer nxtask, task_size
|
||||
external nxtask
|
||||
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
integer ecp_channels
|
||||
integer iecp
|
||||
double precision H_beta (lbuf)
|
||||
double precision dens_beta (max_at_bf,max_at_bf)
|
||||
|
||||
cc Read this value from rtdb vvvv
|
||||
if (.not.rtdb_get(rtdb, 'dft:spin_polarised_ecps',
|
||||
& mt_int, 1, ecp_channels)) then
|
||||
ecp_channels = 1
|
||||
H_beta = 0
|
||||
endif
|
||||
cc AJL: With spin-polarised ECPs Hcore will be spin dependent
|
||||
cc See Szabo and Ostlund pg. 215
|
||||
cc So we need to separate out the densities
|
||||
cc
|
||||
if (ecp_channels.gt.1) then
|
||||
c
|
||||
c Restore alpha and beta densities to calculate spin-polarised
|
||||
c derivatives
|
||||
c
|
||||
cc Currently, this will not work, as g_dens is formulated only in the spin-paired
|
||||
cc manner, and not separable. Therefore, for now I will just put an
|
||||
cc an error. In the long term we need to be able to supply the separable
|
||||
cc alpha and beta densities
|
||||
|
||||
write(LuOut,*)'-------------- WARNING ---------------'
|
||||
write(LuOut,*)'Spin-polarised ECPs have not been implemented'
|
||||
write(LuOut,*)'for TCE gradients. Any results received will be'
|
||||
write(LuOut,*)'calculated assuming all density is in the alpha'
|
||||
write(LuOut,*)'channel, and therefore must be assumed INCORRECT'
|
||||
write(LuOut,*)''
|
||||
|
||||
stop
|
||||
|
||||
c Note: all of the framework is already in place, for future reference
|
||||
c we just need to pass in the alpha and beta densities combined in
|
||||
c g_dens(1), and the beta densities in g_dens(2), and comment out
|
||||
c the above error.
|
||||
c
|
||||
call ga_dadd(1d0, g_dens(1), -1d0, g_dens(2), g_dens(1))
|
||||
end if
|
||||
cc AJL/End
|
||||
|
||||
task_size = 1
|
||||
status = rtdb_parallel(.true.) ! Broadcast reads to all processes
|
||||
c
|
||||
|
|
@ -95,9 +146,22 @@ c
|
|||
|
||||
status = bas_ce2cnr(basis,iat1,iac1f,iac1l)
|
||||
status = bas_ce2cnr(basis,iat2,iac2f,iac2l)
|
||||
|
||||
call ga_get(g_dens, iab1f,iab1l,iab2f,iab2l,dens,max_at_bf)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc call ga_get(g_dens,iab1f,iab1l,iab2f,iab2l,dens,max_at_bf)
|
||||
do 75, iecp = 1, ecp_channels
|
||||
if (iecp.eq.1) then
|
||||
call ga_get(g_dens(iecp), iab1f,iab1l,
|
||||
$ iab2f,iab2l,dens,max_at_bf)
|
||||
else
|
||||
call ga_get(g_dens(iecp), iab1f,iab1l,
|
||||
$ iab2f,iab2l,dens_beta,max_at_bf)
|
||||
endif
|
||||
75 continue
|
||||
cc g_wdens is not dependent on spin, so can leave this
|
||||
call ga_get(g_wdens,iab1f,iab1l,iab2f,iab2l,wdens,max_at_bf)
|
||||
cc AJL/End
|
||||
c
|
||||
c->debug
|
||||
c call ga_print(g_dens)
|
||||
c call ga_print(g_wdens)
|
||||
|
|
@ -154,8 +218,23 @@ C 1el. derivatives
|
|||
call intd_1epot(basis,ish1,basis,ish2,lscr,scr,
|
||||
& lbuf,H)
|
||||
end if
|
||||
|
||||
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc AJL: With spin-polarised ECPs Hcore will be spin dependent
|
||||
cc See Szabo and Ostlund pg. 215
|
||||
if (ecp_channels.gt.1) then
|
||||
C 1el. derivatives
|
||||
if(.not.dobq) then
|
||||
cc For now this will do, but this could be more efficiently done
|
||||
call intd_1eh1_beta(basis,ish1,basis,ish2,lscr,scr,
|
||||
& lbuf,H_beta)
|
||||
else
|
||||
call intd_1epot_beta(basis,ish1,basis,ish2,lscr,scr,
|
||||
& lbuf,H_beta)
|
||||
end if
|
||||
end if
|
||||
cc AJL/End
|
||||
c
|
||||
C D x H
|
||||
c write(6,*) 'D x H'
|
||||
ic=1
|
||||
|
|
@ -165,6 +244,11 @@ c write(6,*) 'D x H'
|
|||
do 31, ip1 = if1, il1
|
||||
do 30, ip2 = if2, il2
|
||||
dE = dE + dens(ip1,ip2) * H(ic)
|
||||
cc AJL/Begin
|
||||
if (ecp_channels.gt.1) then
|
||||
dE = dE + dens_beta(ip1,ip2) * H_beta(ic)
|
||||
end if
|
||||
cc AJL/End
|
||||
c write(6,'(a,i4,i4,f20.16)') 'dens',ip1,ip2,
|
||||
c 1 dens(ip1,ip2)
|
||||
c write(6,'(a,i4,f20.16)') 'H',ic,H(ic)
|
||||
|
|
|
|||
|
|
@ -569,11 +569,20 @@ c $ dbl_mb(k_dens), dbl_mb(k_wdens), dbl_mb(k_frc_nuc),
|
|||
c $ dbl_mb(k_frc_kin), dbl_mb(k_frc_wgh), g_force,
|
||||
c $ g_dens, g_wdens, basis, geom, nproc, nat, max_at_bf,
|
||||
c $ rtdb, oskel)
|
||||
c
|
||||
cc AJL/Begin/SPIN ECPs
|
||||
cc Added ndens to grad1
|
||||
cc Work is still needed to separate g_dens for alpha and beta
|
||||
cc so spin-polarised will not work at this time. It throws
|
||||
cc an error in tce_grad1. To be looked at later...
|
||||
c
|
||||
call tce_grad1(dbl_mb(k_buf),lbuf,dbl_mb(k_scr),lscratch,
|
||||
$ dbl_mb(k_dens), dbl_mb(k_wdens), dbl_mb(k_frc_nuc),
|
||||
$ dbl_mb(k_frc_kin), dbl_mb(k_frc_wgh), g_force,
|
||||
$ g_1pdm, g_1pwdm, basis, geom, nproc, nat, max_at_bf,
|
||||
$ rtdb, oskel)
|
||||
$ rtdb, oskel, ndens)
|
||||
c
|
||||
cc AJL/End
|
||||
c
|
||||
c
|
||||
c ----- cosmo contribution -----
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue