diff --git a/src/NWints/api/int_1estv.F b/src/NWints/api/int_1estv.F index 7effead7fe..b87af6118b 100644 --- a/src/NWints/api/int_1estv.F +++ b/src/NWints/api/int_1estv.F @@ -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> @} diff --git a/src/NWints/api/int_hf1sp.F b/src/NWints/api/int_hf1sp.F index aabb94cb6d..a4f13fd0a1 100644 --- a/src/NWints/api/int_hf1sp.F +++ b/src/NWints/api/int_hf1sp.F @@ -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 diff --git a/src/NWints/api/intd_1eh1.F b/src/NWints/api/intd_1eh1.F index c813f903da..c10964f45a 100644 --- a/src/NWints/api/intd_1eh1.F +++ b/src/NWints/api/intd_1eh1.F @@ -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> @} diff --git a/src/NWints/api/intd_1epot.F b/src/NWints/api/intd_1epot.F index 1826b67f0a..2ac059abe2 100644 --- a/src/NWints/api/intd_1epot.F +++ b/src/NWints/api/intd_1epot.F @@ -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> @} diff --git a/src/NWints/api/intdd_1eecp1.F b/src/NWints/api/intdd_1eecp1.F index 1bdd90598c..04a23bbcdc 100644 --- a/src/NWints/api/intdd_1eecp1.F +++ b/src/NWints/api/intdd_1eecp1.F @@ -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> @} diff --git a/src/NWints/ecp/ecp_nwc.fh b/src/NWints/ecp/ecp_nwc.fh index 197cccda44..a26005636c 100644 --- a/src/NWints/ecp/ecp_nwc.fh +++ b/src/NWints/ecp/ecp_nwc.fh @@ -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 diff --git a/src/NWints/ecp/int_ecp.F b/src/NWints/ecp/int_ecp.F index 05dbf060ba..5ac89b297a 100644 --- a/src/NWints/ecp/int_ecp.F +++ b/src/NWints/ecp/int_ecp.F @@ -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 diff --git a/src/basis/bas.fh b/src/basis/bas.fh index 0d35fd2a80..f4b36d20ae 100644 --- a/src/basis/bas.fh +++ b/src/basis/bas.fh @@ -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 diff --git a/src/basis/basP.fh b/src/basis/basP.fh index dec9a46b1c..c939086150 100644 --- a/src/basis/basP.fh +++ b/src/basis/basP.fh @@ -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 diff --git a/src/basis/bas_blasP.F b/src/basis/bas_blasP.F index 2fc3fd92c0..ec9e555e19 100644 --- a/src/basis/bas_blasP.F +++ b/src/basis/bas_blasP.F @@ -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 diff --git a/src/basis/bas_input.F b/src/basis/bas_input.F index 94ece5c950..a8ff751c14 100644 --- a/src/basis/bas_input.F +++ b/src/basis/bas_input.F @@ -27,7 +27,7 @@ c [spherical|cartesian] [segment||nosegment] [print] [noprint] c c tag nelec c tag library [tag in library] [file ] -c tag +c tag c c ... c end ecp @@ -343,7 +343,7 @@ c c for ecp or so c tag nelec ! only for ecp c tag library [] [file ] -c tag +c tag c 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 ] [rel]',/, $ ' tag nelec ',/, - $ ' tag ',/, + $ ' tag ',/, $ ' ',/, $ ' ... ',/, $ ' 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 [rel]',/, + $ ' tag [rel] [channel]',/, $ ' ',/, $ ' ... ',/, $ ' 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) diff --git a/src/basis/basdeclsP.fh b/src/basis/basdeclsP.fh index 9d419da341..cb83d27f0c 100644 --- a/src/basis/basdeclsP.fh +++ b/src/basis/basdeclsP.fh @@ -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 diff --git a/src/basis/basis.F b/src/basis/basis.F index 86994c89f6..73bbb29e1e 100644 --- a/src/basis/basis.F +++ b/src/basis/basis.F @@ -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> diff --git a/src/ddscf/cint1cache.fh b/src/ddscf/cint1cache.fh index 6a6c00aafa..4abddc45e5 100644 --- a/src/ddscf/cint1cache.fh +++ b/src/ddscf/cint1cache.fh @@ -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 diff --git a/src/ddscf/int_1e_ga.F b/src/ddscf/int_1e_ga.F index 99c5b12b8e..c5c409a0a4 100644 --- a/src/ddscf/int_1e_ga.F +++ b/src/ddscf/int_1e_ga.F @@ -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 diff --git a/src/ddscf/uhf.F b/src/ddscf/uhf.F index 37442661ec..d07e864ecb 100644 --- a/src/ddscf/uhf.F +++ b/src/ddscf/uhf.F @@ -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 diff --git a/src/gradients/grad1.F b/src/gradients/grad1.F index 005147dc61..edad85521b 100644 --- a/src/gradients/grad1.F +++ b/src/gradients/grad1.F @@ -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 diff --git a/src/gradients/grad_force.F b/src/gradients/grad_force.F index 5227683984..1bbe38fa59 100644 --- a/src/gradients/grad_force.F +++ b/src/gradients/grad_force.F @@ -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 diff --git a/src/hessian/analytic/oned_calc.F b/src/hessian/analytic/oned_calc.F index 56369d5ef0..d4498cfe29 100644 --- a/src/hessian/analytic/oned_calc.F +++ b/src/hessian/analytic/oned_calc.F @@ -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 diff --git a/src/hessian/analytic/onedd_calc.F b/src/hessian/analytic/onedd_calc.F index 97169b2cd9..be655c6589 100644 --- a/src/hessian/analytic/onedd_calc.F +++ b/src/hessian/analytic/onedd_calc.F @@ -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 diff --git a/src/hessian/analytic/onedd_cont.F b/src/hessian/analytic/onedd_cont.F index 50d4e19003..36a7e7313b 100644 --- a/src/hessian/analytic/onedd_cont.F +++ b/src/hessian/analytic/onedd_cont.F @@ -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 diff --git a/src/nwdft/scf_dft/dft_fockbld.F b/src/nwdft/scf_dft/dft_fockbld.F index d092ff7de6..83509128e6 100644 --- a/src/nwdft/scf_dft/dft_fockbld.F +++ b/src/nwdft/scf_dft/dft_fockbld.F @@ -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 diff --git a/src/tce/gradients/tce_grad1.F b/src/tce/gradients/tce_grad1.F index d235b20119..8e279a990a 100644 --- a/src/tce/gradients/tce_grad1.F +++ b/src/tce/gradients/tce_grad1.F @@ -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) diff --git a/src/tce/gradients/tce_grad_force.F b/src/tce/gradients/tce_grad_force.F index 5ea57a973a..a08b64cbe5 100644 --- a/src/tce/gradients/tce_grad_force.F +++ b/src/tce/gradients/tce_grad_force.F @@ -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 -----