1!! Copyright (C) 2017 Johannes Flick 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18!> This file handles the evaluation of the photon-OEP potential, 19!! as described in J. Flick et al. ACS Photonics 2018, 5, 3, 992-1005 20! --------------------------------------------------------- 21subroutine X(xc_oep_pt_phi) (namespace, gr, hm, st, is, oep, phi1) 22 type(namespace_t), intent(in) :: namespace 23 type(grid_t), intent(in) :: gr 24 type(hamiltonian_elec_t), intent(in) :: hm 25 type(states_elec_t), intent(in) :: st 26 integer, intent(in) :: is 27 type(xc_oep_t), intent(inout) :: oep 28 R_TYPE, intent(inout) :: phi1(:,:,:) 29 30 integer :: ist, kst, iter_used 31 FLOAT :: rhs_kkbar, residue, kkopii 32 R_TYPE, allocatable :: rhs(:,:), psiii(:, :), psikk(:, :), pol_dip_psiii(:, :) 33 FLOAT, allocatable :: rhs2(:) 34 35 PUSH_SUB(X(xc_oep_pt_phi)) 36 37 SAFE_ALLOCATE(rhs(1:gr%mesh%np, 1:1)) 38 SAFE_ALLOCATE(psiii(1:gr%mesh%np, 1:st%d%dim)) 39 SAFE_ALLOCATE(pol_dip_psiii(1:gr%mesh%np, 1:st%d%dim)) 40 SAFE_ALLOCATE(psikk(1:gr%mesh%np, 1:st%d%dim)) 41 SAFE_ALLOCATE(rhs2(1:gr%mesh%np)) 42 43 if (is == 1) then !only initialize for first spin channel 44 oep%pt%number = M_ZERO 45 oep%pt%correlator = M_ZERO 46 oep%pt%ex = M_ZERO 47 end if 48 49 do ist = st%st_start, oep%noccst 50 call states_elec_get_state(st, gr%mesh, ist, is, psiii) 51 pol_dip_psiii(:, 1) = oep%pt%pol_dipole(:, 1)*psiii(:, 1) 52 rhs(:,1) = -oep%pt%lambda(1)*sqrt(M_HALF*oep%pt%omega(1))*pol_dip_psiii(:, 1) 53 54 do kst = st%st_start, oep%noccst 55 call states_elec_get_state(st, gr%mesh, kst, is, psikk) 56 rhs_kkbar = X(mf_dotp)(gr%mesh, psikk(:, 1), rhs(:, 1), dotu = .true.) 57 call lalg_axpy(gr%mesh%np, -rhs_kkbar, psikk(:, 1), rhs(:, 1)) 58 end do 59 60 call X(states_elec_orthogonalize_single)(st, gr%mesh, st%nst, is, rhs, normalize = .false.) 61 62 call X(linear_solver_solve_HXeY)(oep%solver, namespace, hm, gr, st, ist, is, oep%photon_lr%X(dl_psi)(:, :, ist, is), rhs, & 63 R_TOTYPE(-st%eigenval(ist, is)) + R_REAL(oep%pt%omega(1)), oep%scftol%final_tol, residue, iter_used) 64 65 call X(states_elec_orthogonalize_single)(st, gr%mesh, st%nst, is, & 66 oep%photon_lr%X(dl_psi)(:, :, ist, is), normalize = .false.) 67 68 phi1(1:gr%mesh%np, 1:st%d%dim, ist) = oep%photon_lr%X(dl_psi)(1:gr%mesh%np, 1:st%d%dim, ist, is) 69 70 oep%pt%number(1) = oep%pt%number(1) + st%occ(ist, is)*X(mf_dotp)(gr%mesh, phi1(:, 1, ist), phi1(:, 1, ist), & 71 dotu = .true.) 72 73 oep%pt%ex = oep%pt%ex + st%occ(ist, is)*oep%pt%lambda(1)*sqrt(M_HALF*oep%pt%omega(1))* & 74 X(mf_dotp)(gr%mesh, phi1(:, 1, ist), pol_dip_psiii(:, 1), dotu = .true.) 75 76 oep%pt%ex = oep%pt%ex + st%occ(ist, is)*M_HALF*oep%pt%lambda(1)**2*X(mf_dotp)(gr%mesh, & 77 pol_dip_psiii(:, 1), pol_dip_psiii(:, 1), dotu = .true.) 78 79 do kst = st%st_start, oep%noccst 80 call states_elec_get_state(st, gr%mesh, kst, is, psikk) 81 kkopii = oep%pt%lambda(1)*X(mf_dotp)(gr%mesh, psiii(:, 1), oep%pt%pol_dipole(:, 1)*psikk(:, 1), & 82 dotu = .true.) 83 oep%pt%ex = oep%pt%ex - st%occ(kst, is)*M_HALF*kkopii**2 84 end do 85 86 ! calculate correlator function 87 rhs2(:) = psiii(:,1)*phi1(:, 1, ist) 88 call lalg_axpy(gr%mesh%np, st%occ(ist, is), rhs2(:), oep%pt%correlator(:, 1)) 89 90 end do 91 92 SAFE_DEALLOCATE_A(rhs) 93 SAFE_DEALLOCATE_A(rhs2) 94 SAFE_DEALLOCATE_A(psiii) 95 SAFE_DEALLOCATE_A(psikk) 96 97 POP_SUB(X(xc_oep_pt_phi)) 98end subroutine X(xc_oep_pt_phi) 99 100! --------------------------------------------------------- 101subroutine X(xc_oep_pt_rhs) (gr, st, is, oep, phi1, ist, rhs) 102 type(grid_t), intent(in) :: gr 103 type(states_elec_t), intent(in) :: st 104 integer, intent(in) :: is 105 type(xc_oep_t), intent(inout) :: oep 106 R_TYPE, intent(in) :: phi1(:,:,:) 107 integer, intent(in) :: ist 108 R_TYPE, intent(inout) :: rhs(:,:) 109 110 integer :: kst 111 FLOAT :: abar, kkopii 112 R_TYPE, allocatable :: aa(:,:), psiii(:, :), psikk(:, :) 113 114 PUSH_SUB(X(xc_oep_pt_rhs)) 115 116 SAFE_ALLOCATE(aa(1:gr%mesh%np, 1:1)) 117 SAFE_ALLOCATE(psiii(1:gr%mesh%np, 1:st%d%dim)) 118 SAFE_ALLOCATE(psikk(1:gr%mesh%np, 1:st%d%dim)) 119 120 call states_elec_get_state(st, gr%mesh, ist, is, psiii) 121 122 aa(:,1) = M_HALF*oep%pt%lambda(1)**2*oep%pt%pol_dipole(:, 1)**2 * R_CONJ(psiii(:, 1)) 123 call lalg_axpy(gr%mesh%np, sqrt(M_HALF*oep%pt%omega(1))*oep%pt%lambda(1), & 124 oep%pt%pol_dipole(:, 1)*R_CONJ(phi1(:, 1, ist)), aa(:, 1)) 125 126 do kst = st%st_start, oep%noccst 127 call states_elec_get_state(st, gr%mesh, kst, is, psikk) 128 kkopii = oep%pt%lambda(1)*X(mf_dotp)(gr%mesh, psiii(:, 1), oep%pt%pol_dipole(:, 1)*psikk(:, 1), & 129 dotu = .true.) 130 call lalg_axpy(gr%mesh%np, - oep%pt%lambda(1)*kkopii, oep%pt%pol_dipole(:, 1)*R_CONJ(psikk(:, 1)), aa(:, 1)) 131 call lalg_axpy(gr%mesh%np, - sqrt(M_HALF*oep%pt%omega(1))*kkopii, R_CONJ(phi1(:, 1, kst)), aa(:,1)) 132 end do 133 134 if (ist /= oep%eigen_n + 1 .or. oep%level == XC_OEP_FULL) then 135 abar = X(mf_dotp)(gr%mesh, aa(:, 1), psiii(:, 1)) 136 call lalg_axpy(gr%mesh%np, -abar, R_CONJ(psiii(:, 1)), aa(:, 1)) 137 end if 138 139 call lalg_axpy(gr%mesh%np, M_ONE, aa(:, 1), rhs(:, 1)) 140 141 SAFE_DEALLOCATE_A(aa) 142 SAFE_DEALLOCATE_A(psiii) 143 SAFE_DEALLOCATE_A(psikk) 144 145 POP_SUB(X(xc_oep_pt_rhs)) 146end subroutine X(xc_oep_pt_rhs) 147 148 149! --------------------------------------------------------- 150subroutine X(xc_oep_pt_inhomog) (gr, st, is, phi1, ist, ss) 151 type(grid_t), intent(in) :: gr 152 type(states_elec_t), intent(in) :: st 153 integer, intent(in) :: is 154 R_TYPE, intent(in) :: phi1(:,:,:) 155 integer, intent(in) :: ist 156 FLOAT, intent(inout) :: ss(:) 157 158 FLOAT :: phi_bar 159 R_TYPE, allocatable :: psiii(:, :) 160 FLOAT, allocatable :: rhs(:) 161 162 PUSH_SUB(X(xc_oep_pt_inhomog)) 163 164 SAFE_ALLOCATE(psiii(1:gr%mesh%np, 1:st%d%dim)) 165 SAFE_ALLOCATE(rhs(1:gr%mesh%np)) 166 167 call states_elec_get_state(st, gr%mesh, ist, is, psiii) 168 169 phi_bar = X(mf_dotp)(gr%mesh, phi1(:, 1, ist), phi1(:, 1, ist), dotu = .true.) 170 rhs = phi1(:, 1, ist)*phi1(:, 1, ist) - phi_bar*(R_CONJ(psiii(:, 1))*psiii(:, 1)) 171 call lalg_axpy(gr%mesh%np, -M_ONE, rhs, ss) 172 173 SAFE_DEALLOCATE_A(psiii) 174 SAFE_DEALLOCATE_A(rhs) 175 176 POP_SUB(X(xc_oep_pt_inhomog)) 177end subroutine X(xc_oep_pt_inhomog) 178 179! --------------------------------------------------------- 180subroutine X(xc_oep_pt_uxcbar) (gr, st, is, oep, phi1, ist, vxbar) 181 type(grid_t), intent(in) :: gr 182 type(states_elec_t), intent(in) :: st 183 integer, intent(in) :: is 184 type(xc_oep_t), intent(in) :: oep 185 R_TYPE, intent(in) :: phi1(:,:,:) 186 integer, intent(in) :: ist 187 FLOAT, intent(inout) :: vxbar 188 189 integer :: kst 190 FLOAT :: kkopii, result1 191 R_TYPE, allocatable :: aa(:,:), psiii(:,:), psikk(:,:) 192 193 PUSH_SUB(X(xc_oep_pt_uxcbar)) 194 195 SAFE_ALLOCATE(aa(1:gr%mesh%np, 1:1)) 196 SAFE_ALLOCATE(psiii(1:gr%mesh%np, 1:st%d%dim)) 197 SAFE_ALLOCATE(psikk(1:gr%mesh%np, 1:st%d%dim)) 198 199 call states_elec_get_state(st, gr%mesh, ist, is, psiii) 200 201 aa(:, 1) = oep%pt%lambda(1)*oep%pt%pol_dipole(:, 1)*psiii(:, 1) 202 result1 = sqrt(M_HALF*oep%pt%omega(1))*X(mf_dotp)(gr%mesh, phi1(:, 1, ist), aa(:, 1), dotu = .true.) 203 result1 = result1 + M_HALF*dmf_dotp(gr%mesh, R_REAL(aa(:, 1)), R_REAL(aa(:, 1))) 204 205 do kst = st%st_start, oep%noccst 206 call states_elec_get_state(st, gr%mesh, kst, is, psikk) 207 kkopii = oep%pt%lambda(1)*X(mf_dotp)(gr%mesh, psiii(:, 1), oep%pt%pol_dipole(:, 1)*psikk(:, 1), & 208 dotu = .true.) 209 result1 = result1 - M_HALF*kkopii**2 210 end do 211 212 vxbar = vxbar - result1 213 214 SAFE_DEALLOCATE_A(aa) 215 SAFE_DEALLOCATE_A(psiii) 216 SAFE_DEALLOCATE_A(psikk) 217 218 POP_SUB(X(xc_oep_pt_uxcbar)) 219end subroutine X(xc_oep_pt_uxcbar) 220 221!! Local Variables: 222!! mode: f90 223!! coding: utf-8 224!! End: 225