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