added RT-TDDFT molecular orbital complex absorbing potential

This commit is contained in:
Ken Lopata 2016-05-21 05:33:52 +00:00
parent 21edee3d13
commit 20ea77a650
12 changed files with 265 additions and 43 deletions

View file

@ -188,10 +188,15 @@ C
call zfock_cs_build (params, g_zdens_ao, energies, g_zfock_ao)
call canorg_trans (params,"F","AO->MO", g_zfock_ao, g_zfock_mo)
if (.not. mat_is_purereal (g_zfock_mo, params%tol_zero))
$ call errquit (pname//
$ "initial Fock matrix in MO basis is not pure real", 0, 0)
C call rt_tddft_print_notice ("Disabled init MO Fock check")
if (params%mocap_active) then
call rt_tddft_print_notice (
$ "MO CAP active: Disabled initial MO Fock check")
else
if (.not. mat_is_purereal (g_zfock_mo, params%tol_zero))
$ call errquit (pname//
$ "initial Fock matrix in MO basis is not pure real", 0, 0)
endif
call rt_tddft_calc_Etot (params, energies)
energies%gs = energies%tot
@ -229,7 +234,7 @@ C if (params%print_moocc .or. do_zora) then
if (.not. mat_is_purereal (g_zfock_mo, params%tol_zero))
$ call errquit (pname//
$ "initial Fock matrix in MO basis is not pure real", 0, 0)
$ "MO occup needs initial MO Fock mat to be pure real",0,0)
call rt_tddft_movecs_gs (params, g_zfock_mo, g_movecs_mo_gs)
endif

View file

@ -54,7 +54,7 @@ C == Variables ==
double precision dftExc(2)
logical status
double precision elapsed, elapsed_tot
integer g_zfock_mo, g_zcapmo ! for MO damping
C == External routines ==
logical xc_gotxc
@ -258,6 +258,34 @@ C not have changed. Perhaps double check this??
if (params%prof) call prof_end (elapsed, "Fock CS d2z")
C
C Apply molecular orbital complex absorbing potential (MO CAP)
C
if (params%mocap_active) then
if (.not. ga_create(mt_dcpl, params%ns_mo, params%ns_mo,
$ "Fock MO tmp" , 0, 0, g_zfock_mo))
$ call errquit(pname//"GA alloc failed",0, GA_ERR)
if (.not. ga_create(mt_dcpl, params%ns_mo, params%ns_mo,
$ "Fock MO tmp" , 0, 0, g_zcapmo))
$ call errquit(pname//"GA alloc failed",0, GA_ERR)
call canorg_trans (params,"F","AO->MO", g_zfock, g_zfock_mo)
call rt_tddft_mocap (params, g_zfock_mo, g_zcapmo)
call ga_zadd (z1, g_zfock_mo, z1, g_zcapmo, g_zfock_mo)
call canorg_trans (params,"F","MO->AO", g_zfock_mo, g_zfock)
call rt_tddft_print_notice ("Applied closed shell MO CAP")
if (.not. ga_destroy(g_zfock_mo))
$ call errquit(pname//"GA destroy failed",0, GA_ERR)
if (.not. ga_destroy(g_zcapmo))
$ call errquit(pname//"GA destroy failed",0, GA_ERR)
endif
C
C Extra checks, disabled for speed.
C

View file

@ -351,6 +351,15 @@ C (zora stuff)
integer :: g_zora_sf(2)
integer :: g_zora_scale_sf(2)
C (MO complex absorbing potential stuff)
logical :: mocap_active
logical :: mocap_print
logical :: mocap_checks
double precision :: mocap_emin
double precision :: mocap_maxval
double precision :: mocap_prefac
double precision :: mocap_expconst
end type rt_params_t
c $Id$

View file

@ -206,6 +206,39 @@ C
$ params%excite(iexcite))
enddo
C
C MO CAP
C
if (.not.rtdb_get(params%rtdb, "rt_tddft:mocap_active",
$ mt_log, 1, params%mocap_active)) call errquit (pname//
$ "failed to read mocap_active from rtdb", 0, RTDB_ERR)
if (params%mocap_active) then
if (.not. rtdb_get (params%rtdb, "rt_tddft:mocap_emin",
$ mt_dbl, 1, params%mocap_emin)) call errquit (pname//
$ "failed to read mocap_emin from rtdb", 0, RTDB_ERR)
if (.not. rtdb_get (params%rtdb, "rt_tddft:mocap_maxval",
$ mt_dbl, 1, params%mocap_maxval)) call errquit (pname//
$ "failed to read mocap_maxval from rtdb", 0, RTDB_ERR)
if (.not. rtdb_get (params%rtdb, "rt_tddft:mocap_prefac",
$ mt_dbl, 1, params%mocap_prefac)) call errquit (pname//
$ "failed to read mocap_prefac from rtdb", 0, RTDB_ERR)
if (.not. rtdb_get (params%rtdb, "rt_tddft:mocap_expconst",
$ mt_dbl, 1, params%mocap_expconst)) call errquit (pname//
$ "failed to read mocap_expconst from rtdb", 0, RTDB_ERR)
if (.not.rtdb_get(params%rtdb, "rt_tddft:mocap_checks",
$ mt_log, 1, params%mocap_checks)) call errquit (pname//
$ "failed to read mocap_checks from rtdb", 0, RTDB_ERR)
if (.not.rtdb_get(params%rtdb, "rt_tddft:mocap_print",
$ mt_log, 1, params%mocap_print)) call errquit (pname//
$ "failed to read mocap_print from rtdb", 0, RTDB_ERR)
endif
end subroutine
c $Id$

View file

@ -5,6 +5,7 @@ rt_tddft_input_print.o \
rt_tddft_input_excite.o \
rt_tddft_input_tolerances.o \
rt_tddft_input_visualization.o \
rt_tddft_input_mocap.o \
rt_tddft_input.o

View file

@ -182,6 +182,13 @@ C visualization sub-block
C
elseif (inp_compare(.false.,test,'visualization')) then
call rt_tddft_input_visualization (rtdb)
C
C MO CAP sub-block
C
elseif (inp_compare(.false.,test,'mocap')) then
call rt_tddft_input_mocap (rtdb)
C

View file

@ -41,6 +41,7 @@ C integer, parameter :: nsnapshots = 0
logical, parameter :: nodisk = .false. !XXX DISABLED NOW--ie, just use DFT settings
logical, parameter :: matrix_checks = .false.
logical, parameter :: use_dmat = .false.
logical, parameter :: mocap = .false.
character(len=*), parameter :: tag = "<rt_tddft>: "
@ -189,6 +190,9 @@ C
if (.not.rtdb_put(rtdb,'rt_tddft:use_dmat',mt_log,1,use_dmat))
$ call errquit(pname//'Write failed to rtdb',0,RTDB_ERR)
if (.not.rtdb_put(rtdb,'rt_tddft:mocap_active',mt_log,1,mocap))
$ call errquit(pname//'Write failed to rtdb',0,RTDB_ERR)
end subroutine

View file

@ -218,14 +218,22 @@ C
call canorg_trans(params,"F","AO->MO",g_zfock_ao(1),g_zfock_mo(1))
call canorg_trans(params,"F","AO->MO",g_zfock_ao(2),g_zfock_mo(2))
if (.not. mat_is_purereal (g_zfock_mo(1), params%tol_zero))
$ call errquit (pname//
$ "initial Fock alpha matrix in MO basis is not pure real",0,0)
if (.not. mat_is_purereal (g_zfock_mo(2), params%tol_zero))
$ call errquit (pname//
$ "initial Fock beta matrix in MO basis is not pure real",0,0)
if (params%mocap_active) then
call rt_tddft_print_notice (
$ "MO CAP active: Disabled initial MO Fock check")
else
if (.not. mat_is_purereal (g_zfock_mo(1), params%tol_zero))
$ call errquit (pname//
$ "initial Fock alpha matrix in MO basis not pure real",0,0)
if (.not. mat_is_purereal (g_zfock_mo(2), params%tol_zero))
$ call errquit (pname//
$ "initial Fock beta matrix in MO basis not pure real",0,0)
endif
call rt_tddft_calc_Etot (params, energies)
energies%gs = energies%tot

View file

@ -62,6 +62,7 @@ C == Variables ==
logical status
double precision elapsed, elapsed_tot
integer is !spin counter
integer g_zcapmo, g_zfock_mo ! tmp GAs for MO CAP
C == External routines ==
@ -252,6 +253,50 @@ C
if (params%prof) call prof_end (elapsed, "Fock OS d2z")
C
C Apply molecular orbital complex absorbing potential (MO CAP).
C Potential is different for alpha and beta spins, as it is
C built from their respective Fock matrices.
C
if (params%mocap_active) then
if (.not. ga_create(mt_dcpl, params%ns_mo, params%ns_mo,
$ "Fock MO tmp" , 0, 0, g_zfock_mo))
$ call errquit(pname//"GA alloc failed",0, GA_ERR)
if (.not. ga_create(mt_dcpl, params%ns_mo, params%ns_mo,
$ "Fock MO tmp" , 0, 0, g_zcapmo))
$ call errquit(pname//"GA alloc failed",0, GA_ERR)
C (alpha spin)
call canorg_trans (params,"F","AO->MO",
$ g_zfock(1), g_zfock_mo)
call rt_tddft_mocap (params, g_zfock_mo, g_zcapmo)
call ga_zadd (z1, g_zfock_mo,
$ z1, g_zcapmo, g_zfock_mo)
call canorg_trans (params,"F","MO->AO",
$ g_zfock_mo, g_zfock(1))
call rt_tddft_print_notice ("Applied MO CAP to alpha spin")
C (beta spin)
call canorg_trans (params,"F","AO->MO",
$ g_zfock(2), g_zfock_mo)
call rt_tddft_mocap (params, g_zfock_mo, g_zcapmo)
call ga_zadd (z1, g_zfock_mo,
$ z1, g_zcapmo, g_zfock_mo)
call canorg_trans (params,"F","MO->AO",
$ g_zfock_mo, g_zfock(2))
call rt_tddft_print_notice ("Applied MO CAP to beta spin")
if (.not. ga_destroy(g_zfock_mo))
$ call errquit(pname//"GA destroy failed",0, GA_ERR)
if (.not. ga_destroy(g_zcapmo))
$ call errquit(pname//"GA destroy failed",0, GA_ERR)
endif
C
C Extra checks, disabled for speed.
C

View file

@ -15,7 +15,8 @@ rt_tddft_excite.o \
rt_tddft_checkpoint.o \
rt_tddft_zora_scale.o \
rt_tddft_project.o \
rt_tddft_zdens_tot_movecs.o
rt_tddft_zdens_tot_movecs.o \
rt_tddft_mocap.o
OBJ_OPTIMIZE =

View file

@ -8,7 +8,8 @@ C where C' is the eigenvector matrix of the matrix B.
C
C XXX CHECK: only valid when B is the MO Fock matrix?
C
subroutine rt_tddft_project (params, dir, g_za, g_zb, g_zaproj)
subroutine rt_tddft_project (params, dir, g_za, g_zb,
$ g_zaproj, lherm)
implicit none
#include "errquit.fh"
@ -26,6 +27,7 @@ C == Inputs ==
character, intent(in) :: dir !="F" for forward, "R" for reverse
integer, intent(in) :: g_za
integer, intent(in) :: g_zb
logical, intent(in) :: lherm !T= g_za is hermitian, F=g_zb is non hermitian
C == Outputs ==
integer, intent(in) :: g_zaproj
@ -37,6 +39,12 @@ C == Variables ==
integer n
integer ivals, lvals
integer g_zevecs, g_zevecs_hconj, g_ztmp
integer g_zevecs_inv, g_zevecs_hconj_inv
integer me
double complex zval
integer i
me = ga_nodeid()
C XXX ASSUME SIZE: nmo
@ -52,39 +60,103 @@ C XXX ASSUME SIZE: nmo
if (.not. ga_duplicate (g_za, g_zevecs_hconj, "zvecs_hconj"))
$ call errquit (pname//"failed to create ztmp1", 0, GA_ERR)
if (.not. ga_duplicate (g_za, g_zevecs_hconj_inv, "zvecs_hconj"))
$ call errquit (pname//"failed to create ztmp1", 0, GA_ERR)
if (.not. ga_duplicate (g_za, g_zevecs_inv, "zvecs_inv"))
$ call errquit (pname//"failed to create inv", 0, GA_ERR)
if (.not. ga_duplicate (g_za, g_ztmp, "ztmp"))
$ call errquit (pname//"failed to create ztmp", 0, GA_ERR)
if (.not.ma_push_get(mt_dcpl, params%ns_ao,'vals',lvals,ivals))
& call errquit(pname//'cannot allocate vals',0, MA_ERR)
& call errquit(pname//'cannot allocate vals',0, MA_ERR) !xxxx size???
C Diagonalize B matrix to get C matrix and eigenvalues
call ga_zero (g_zevecs)
C call zmat_diag (g_zb, g_zevecs, dcpl_mb(ivals), .true.) !true = reorder evals by energy
call zmat_diag (g_zb, g_zevecs, dcpl_mb(ivals), .false., .false.) !true = reorder evals by energy
call zmat_diag (g_zb, g_zevecs, dcpl_mb(ivals), .true., .false.) !true = reorder evals by energy, false = dont check
call ga_zero (g_zevecs_hconj)
call mat_hermconj (g_zevecs, g_zevecs_hconj)
if ( (dir.eq."F").or.(dir.eq."f")) then ! C^+ A C
call ga_zero (g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_za, g_zevecs, z0, g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_zevecs_hconj, g_ztmp, z0, g_zaproj)
elseif ( (dir.eq."R").or.(dir.eq."r")) then ! C A C^+
call ga_zero (g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_za, g_zevecs_hconj, z0, g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_zevecs, g_ztmp, z0, g_zaproj)
else
call errquit (pname//"invalid flag: "//dir, 0,0)
endif
C
C Print all eigenvalues to screen is turned off for now
C
c$$$ if (me.eq.0) then
c$$$ write(luout, *) ""
c$$$ write(luout, *) "Eigenvalues"
c$$$ do i = 1, params%ns_ao !size??
c$$$ zval = dcpl_mb(ivals + i - 1)
c$$$ write(luout, "(1x,a,i8,3x,2e16.6)")
c$$$ $ "MO", i, dble(zval), aimag(zval)
c$$$ enddo
c$$$ endif
call ga_sync()
if (lherm) then
if (.not. mat_is_unitary(g_zevecs, 1d-8)) then
if (me.eq.0) then
write(luout, *) pname//"Eigenvector matrix non-unitary"
call halt()
endif
else
if (me.eq.0) write(luout, *)
$ pname//"Eigenvector matrix is unitary. Good."
endif
endif
if (lherm) then
call mat_hermconj (g_zevecs, g_zevecs_hconj) !this assumes eigenvector matrix is unitary: U^+ = U^-1
if ( (dir.eq."F").or.(dir.eq."f")) then ! C^+ A C
call ga_zero (g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_za, g_zevecs, z0, g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_zevecs_hconj, g_ztmp, z0, g_zaproj)
elseif ( (dir.eq."R").or.(dir.eq."r")) then ! C A C^+
call ga_zero (g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_za, g_zevecs_hconj, z0, g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_zevecs, g_ztmp, z0, g_zaproj)
else
call errquit (pname//"invalid flag: "//dir, 0,0)
endif
else ! non hermitian / non-unitary eigenvector matrix
call mat_hermconj (g_zevecs, g_zevecs_hconj) ! = C^+
call zmat_inv(g_zevecs_hconj, g_zevecs_hconj_inv) ! = (C^+)^-1
call zmat_inv_check(g_zevecs_hconj, g_zevecs_hconj_inv)
call zmat_inv(g_zevecs, g_zevecs_inv) ! = C^-1
call zmat_inv_check(g_zevecs, g_zevecs_inv) ! = C^-1
if ( (dir.eq."F").or.(dir.eq."f")) then ! C^+ A C
call ga_zero (g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_za, g_zevecs, z0, g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_zevecs_hconj, g_ztmp, z0, g_zaproj)
elseif ( (dir.eq."R").or.(dir.eq."r")) then ! (C^+)^-1 A C^-1
call ga_zero (g_ztmp)
call ga_zgemm ("N", "N", n, n, n,
$ z1, g_za, g_zevecs_inv, z0, g_ztmp)
call ga_zgemm("N", "N", n, n, n,
$ z1, g_zevecs_hconj_inv, g_ztmp, z0, g_zaproj)
else !
call errquit (pname//"invalid flag: "//dir, 0,0)
endif
endif
if (.not. ma_chop_stack (lvals))
& call errquit(pname//'cannot chop stack',0, MA_ERR)
@ -95,6 +167,12 @@ C call zmat_diag (g_zb, g_zevecs, dcpl_mb(ivals), .true.) !true = reorder e
if (.not. ga_destroy (g_zevecs_hconj))
$ call errquit (pname//"failed to destrory ztmp1", 0, GA_ERR)
if (.not. ga_destroy (g_zevecs_hconj_inv))
$ call errquit (pname//"failed to destrory ztmp1", 0, GA_ERR)
if (.not. ga_destroy (g_zevecs_inv))
$ call errquit (pname//"failed to destrory inv", 0, GA_ERR)
if (.not. ga_destroy (g_ztmp))
$ call errquit (pname//"failed to destrory ztmp", 0, GA_ERR)
@ -105,7 +183,7 @@ C====================================================================
C
C Check the projecting A->B and reversing yields A again.
C
subroutine rt_tddft_project_check (params, g_za, g_zb)
subroutine rt_tddft_project_check (params, g_za, g_zb, lherm)
implicit none
#include "errquit.fh"
@ -122,6 +200,7 @@ C == Inputs ==
type (rt_params_t), intent(in) :: params
integer, intent(in) :: g_za
integer, intent(in) :: g_zb
logical, intent(in) :: lherm ! hermitian mats?
C == Parameters ==
character(*),parameter :: pname = "rt_tddft_project_check: "
@ -139,8 +218,9 @@ C == Variables ==
call ga_zero (g_zaproj)
call rt_tddft_project (params, "F", g_za, g_zb, g_zaproj)
call rt_tddft_project (params, "R", g_zaproj, g_zb, g_zaback)
call rt_tddft_project (params, "F", g_za, g_zb, g_zaproj, lherm)
call rt_tddft_project (params, "R",
$ g_zaproj, g_zb, g_zaback, lherm)
if (.not. mat_is_same (g_za, g_zaback, params%tol_zero))
$ call errquit (pname//"failed",0,0)
@ -305,4 +385,3 @@ C call ga_print (g_zaproj)
$ call errquit (pname//"failed to destrory zaback", 0, GA_ERR)
end subroutine
c $Id$

View file

@ -184,10 +184,12 @@ C (make diag occup matrix)
enddo
C (project onto TD complex movecs)
call rt_tddft_project_check (params, g_ztmp3, g_zdens_mo)
C Note: only for Hermitian matrices (.true. flag)
call rt_tddft_project_check (params, g_ztmp3, g_zdens_mo, .true.)
if (me.eq.0) write (luout,*) pname//"projection check passed"
call ga_zero (g_ztmp1)
call rt_tddft_project (params, "R", g_ztmp3, g_zdens_mo, g_ztmp1) ! g_ztmp1 now has projected occup mask
call rt_tddft_project (params, "R", g_ztmp3,
$ g_zdens_mo, g_ztmp1, .true.) ! g_ztmp1 now has projected occup mask
call ga_zero (g_ztmp4)
call ga_elem_multiply (g_ztmp1, g_zevecs_ao, g_ztmp4) ! g_ztmp4 now holds masked movecs in AO basis