1!! Copyright (C) 2012-2013 M. Gruning, P. Melo, M. Oliveira 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 19! --------------------------------------------------------- 20subroutine xc_kli_pauli_solve(mesh, namespace, st, oep) 21 type(mesh_t), intent(in) :: mesh 22 type(namespace_t), intent(in) :: namespace 23 type(states_elec_t), intent(in) :: st 24 type(xc_oep_t), intent(inout) :: oep 25 ! 26 integer :: is, ip, ii, jj, ist, eigen_n, it, kssi 27 FLOAT, allocatable :: rho(:,:), lambda(:), n(:), t_rho(:,:) 28 FLOAT, allocatable :: t_v(:,:), vloc(:,:), rhov(:), v_m1(:,:), p_i(:,:,:), t_vi(:,:,:), delta_v(:), vs(:,:) 29 CMPLX, allocatable :: weighted_hf(:,:,:), rho_i(:,:,:,:), psii(:), psij(:) 30 FLOAT :: reached_threshold(4) 31 32 call profiling_in(C_PROFILING_XC_KLI, 'XC_KLI') 33 PUSH_SUB(xc_kli_pauli_solve) 34 35 ! Density related quantities 36 SAFE_ALLOCATE(rho(1:mesh%np, 4)) 37 SAFE_ALLOCATE(n(1:mesh%np)) 38 SAFE_ALLOCATE(lambda(1:mesh%np)) 39 rho(1:mesh%np, 1:4) = st%rho(1:mesh%np, 1:4) 40 do is = 1, 2 41 do ip = 1,mesh%np 42 if (rho(ip, is) < CNST(1e-20)) rho(ip, is) = CNST(1e-20) 43 end do 44 end do 45 n(1:mesh%np) = rho(1:mesh%np,1) + rho(1:mesh%np,2) 46 lambda(1:mesh%np) = rho(1:mesh%np,1)*rho(1:mesh%np,2) - (rho(1:mesh%np,3)**2 + rho(1:mesh%np,4)**2) 47 48 ! Potential related quantities 49 ! (Built from HF potentials weighted with orbital densities) 50 SAFE_ALLOCATE(weighted_hf(1:mesh%np, 1:st%d%dim,st%d%dim)) 51 weighted_hf = M_Z0 52 53 SAFE_ALLOCATE(psij(1:mesh%np)) 54 55 ! w_{up,down} = \sum_i \phi_{i,down}^* u_x^{i,up}^* \phi_{i,up} 56 do ii = 1,st%d%dim 57 do jj = 1,st%d%dim 58 do ist = st%st_start,st%st_end 59 call states_elec_get_state(st, mesh, jj, ist, 1, psij) 60 weighted_hf(1:mesh%np,ii,jj) = weighted_hf(1:mesh%np,ii,jj) + & 61 oep%socc*st%occ(ist,1)*conjg(psij(1:mesh%np)*oep%zlxc(1:mesh%np,ist,ii)) ! oep%zlxc => (\phi_j)^*u_x^j 62 end do 63 end do 64 end do 65 66 SAFE_DEALLOCATE_A(psij) 67 68 ! reduce over states 69 if(st%parallel_in_states) then 70 call comm_allreduce(st%mpi_grp%comm, weighted_hf) 71 end if 72 73 SAFE_ALLOCATE(t_v(1:mesh%np, 1:4)) 74 t_v = M_ZERO 75 t_v(1:mesh%np,1) = TOFLOAT(weighted_hf(1:mesh%np,2,2)) 76 t_v(1:mesh%np,2) = TOFLOAT(weighted_hf(1:mesh%np,1,1)) 77 t_v(1:mesh%np,3) = -TOFLOAT(weighted_hf(1:mesh%np,1,2) + weighted_hf(1:mesh%np,2,1)) 78 t_v(1:mesh%np,4) = -aimag(weighted_hf(1:mesh%np,1,2) - weighted_hf(1:mesh%np,2,1)) 79 SAFE_DEALLOCATE_A(weighted_hf) 80 81 SAFE_ALLOCATE(vloc(1:mesh%np, 1:4)) 82 vloc = M_ZERO 83 vloc(1:mesh%np,1) = t_v(1:mesh%np,2) - t_v(1:mesh%np,1) 84 vloc(1:mesh%np,2) = -vloc(1:mesh%np,1) 85 vloc(1:mesh%np,3) = -t_v(1:mesh%np,3) 86 vloc(1:mesh%np,4) = -t_v(1:mesh%np,4) 87 88 ! Combine them to obtain Slater part 89 SAFE_ALLOCATE(rhov(1:mesh%np)) 90 rhov = M_ZERO 91 do ip = 1, mesh%np 92 rhov(ip) = sum(rho(ip,1:4)*t_v(ip,1:4)) 93 end do 94 rhov(1:mesh%np) = rhov(1:mesh%np)/lambda(1:mesh%np) 95 SAFE_DEALLOCATE_A(t_v) 96 97 SAFE_ALLOCATE(t_rho(1:mesh%np, 1:4)) 98 t_rho(1:mesh%np,1) = rho(1:mesh%np,2) 99 t_rho(1:mesh%np,2) = rho(1:mesh%np,1) 100 t_rho(1:mesh%np,3:4) = -rho(1:mesh%np,3:4) 101 do ip = 1, mesh%np 102 vloc(ip,1:4) = (vloc(ip,1:4) + t_rho(ip,1:4)*rhov(ip))/n(ip) 103 end do 104 105 106 SAFE_ALLOCATE(vs(1:mesh%np, 1:4)) 107 vs = vloc ! Slater part 108 109 ! iteration criteria 110 call scf_tol_init(oep%scftol, namespace, st%qtot, def_maximumiter=50) 111 112 ! get the HOMO state 113 call xc_oep_AnalyzeEigen(oep, st, 1) 114 eigen_n = oep%eigen_n 115 if (eigen_n == 0) then 116 oep%vxc = vs 117 118 else 119 120 ! orbital densities 121 SAFE_ALLOCATE(rho_i(1:mesh%np, 1:st%d%dim,st%d%dim, 1:eigen_n)) 122 rho_i = M_Z0 123 124 SAFE_ALLOCATE(psii(1:mesh%np)) 125 SAFE_ALLOCATE(psij(1:mesh%np)) 126 127 do ii = 1, st%d%dim 128 do jj = ii, st%d%dim 129 do ist = 1, eigen_n 130 kssi = oep%eigen_index(ist) 131 132 call states_elec_get_state(st, mesh, ii, kssi, 1, psii) 133 call states_elec_get_state(st, mesh, jj, kssi, 1, psij) 134 135 rho_i(1:mesh%np,ii,jj,ist) = oep%socc*st%occ(kssi,1)*conjg(psij(1:mesh%np))*psii(1:mesh%np) 136 rho_i(1:mesh%np,jj,ii,ist) = conjg(rho_i(1:mesh%np,ii,jj,ist)) 137 end do 138 end do 139 end do 140 141 SAFE_DEALLOCATE_A(psii) 142 SAFE_DEALLOCATE_A(psij) 143 144 ! arrange them in a 4-vector 145 SAFE_ALLOCATE(p_i(1:mesh%np, 1:4, 1:eigen_n)) 146 p_i = M_ZERO 147 p_i(1:mesh%np,1,:) = TOFLOAT(rho_i(1:mesh%np,1,1,:)) 148 p_i(1:mesh%np,2,:) = TOFLOAT(rho_i(1:mesh%np,2,2,:)) 149 p_i(1:mesh%np,3,:) = M_TWO*TOFLOAT(rho_i(1:mesh%np,1,2,:)) 150 p_i(1:mesh%np,4,:) = M_TWO*aimag(rho_i(1:mesh%np,1,2,:)) 151 152 SAFE_DEALLOCATE_A(rho_i) 153 154 ! Calculate iteratively response part 155 SAFE_ALLOCATE(v_m1(1:mesh%np, 1:4)) 156 SAFE_ALLOCATE(delta_v(1:eigen_n)) 157 SAFE_ALLOCATE(t_vi(1:mesh%np, 1:4, 1:eigen_n)) 158 159 vloc = M_ZERO 160 KLI_iteration: do it = 1,oep%scftol%max_iter 161 v_m1 = vs + vloc 162 163 ! delta_v^KLI 164 delta_v = M_ZERO 165 do ist = 1,eigen_n 166 do is = 1,st%d%nspin 167 delta_v(ist) = delta_v(ist)+ dmf_dotp(mesh,p_i(1:mesh%np,is,ist),v_m1(1:mesh%np,is), reduce = .false.) 168 end do 169 end do 170 if(mesh%parallel_in_domains) call comm_allreduce(mesh%mpi_grp%comm, delta_v, dim = eigen_n) 171 172 do ist = 1,eigen_n 173 kssi = oep%eigen_index(ist) 174 delta_v(ist) = delta_v(ist) - TOFLOAT(sum(oep%uxc_bar(kssi,:))) 175 end do 176 177 ! 178 t_vi(1:mesh%np,1,:) = p_i(1:mesh%np,2,:) 179 t_vi(1:mesh%np,2,:) = p_i(1:mesh%np,1,:) 180 t_vi(1:mesh%np,3,:) =-p_i(1:mesh%np,3,:) 181 t_vi(1:mesh%np,4,:) =-p_i(1:mesh%np,4,:) 182 do is = 1, st%d%nspin 183 do ip = 1, mesh%np 184 t_vi(ip, is, :) = t_vi(ip, is, :)*delta_v(:) 185 end do 186 end do 187 188 vloc = M_ZERO 189 do ip = 1, mesh%np 190 vloc(ip,1) = sum(t_vi(ip,2,:) - t_vi(ip,1,:)) 191 vloc(ip,2) = -vloc(ip,1) 192 vloc(ip,3) = -sum(t_vi(ip,3,:)) 193 vloc(ip,4) = -sum(t_vi(ip,4,:)) 194 end do 195 196 ! 197 rhov = M_ZERO 198 do ist = 1, eigen_n 199 do ip = 1, mesh%np 200 rhov(ip) = rhov(ip) + sum(rho(ip,:)*t_vi(ip,:,ist)) 201 end do 202 end do 203 rhov = rhov/lambda 204 205 do ip = 1, mesh%np 206 vloc(ip,:) = (vloc(ip,:) + t_rho(ip,:)*rhov(ip))/n(ip) 207 end do 208 ! 209 do is = 1, 4 210 reached_threshold(is) = dmf_nrm2(mesh,(vs(1:mesh%np,is) + vloc(1:mesh%np,is) - v_m1(1:mesh%np,is))) 211 end do 212 if (all(reached_threshold(:) <= oep%scftol%conv_abs_dens)) exit 213 214 end do KLI_iteration 215 216 if (.not. all(reached_threshold(:) <= oep%scftol%conv_abs_dens)) & 217 it = it - 1 218 219 write(message(1), '(a,i4,a,es14.6)') & 220 "Info: After ", it, " iterations, KLI converged to ", maxval(reached_threshold(:)) 221 message(2) = '' 222 call messages_info(2) 223 ! 224 oep%vxc = v_m1 225 226 SAFE_DEALLOCATE_A(vs) 227 SAFE_DEALLOCATE_A(v_m1) 228 SAFE_DEALLOCATE_A(delta_v) 229 SAFE_DEALLOCATE_A(t_vi) 230 SAFE_DEALLOCATE_A(p_i) 231 232 end if 233 234 SAFE_DEALLOCATE_A(rho) 235 SAFE_DEALLOCATE_A(lambda) 236 SAFE_DEALLOCATE_A(n) 237 SAFE_DEALLOCATE_A(t_rho) 238 SAFE_DEALLOCATE_A(rhov) 239 240 call profiling_out(C_PROFILING_XC_KLI) 241 242 POP_SUB(xc_kli_pauli_solve) 243end subroutine xc_kli_pauli_solve 244