1!! Copyright (C) 2017 N. Tancogne-Dejean 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! --------------------------------------------------------- 20! TODO: Merge this with the two_body routine in system/output_me_inc.F90 21subroutine compute_complex_coulomb_integrals (this, mesh, der, st, psolver, namespace) 22 type(lda_u_t), intent(inout) :: this 23 type(mesh_t), intent(in) :: mesh 24 type(derivatives_t), intent(in) :: der 25 type(states_elec_t), intent(in) :: st 26 type(poisson_t), intent(in) :: psolver 27 type(namespace_t), intent(in) :: namespace 28 29 integer :: ist, jst, kst, lst, ijst, klst 30 integer :: is1, is2 31 integer :: norbs, np_sphere, ios, ip 32 integer :: idone, ntodo 33 CMPLX, allocatable :: tmp(:), vv(:,:), nn(:,:) 34 type(orbitalset_t), pointer :: os 35 type(profile_t), save :: prof 36 37 call profiling_in(prof, "DFTU_COMPEX_COULOMB_INTEGRALS") 38 39 PUSH_SUB(compute_complex_coulomb_integrals) 40 41 ASSERT(.not. st%parallel_in_states) 42 if(mesh%parallel_in_domains) then 43 call messages_not_implemented("Coulomb integrals parallel in domains", namespace=namespace) 44 end if 45 46 SAFE_ALLOCATE(nn(1:this%max_np,st%d%dim)) 47 SAFE_ALLOCATE(vv(1:this%max_np,st%d%dim)) 48 SAFE_ALLOCATE(tmp(1:this%max_np)) 49 SAFE_ALLOCATE(this%zcoulomb(1:this%maxnorbs,1:this%maxnorbs,1:this%maxnorbs,1:this%maxnorbs,1:st%d%dim,1:st%d%dim,1:this%norbsets)) 50 this%zcoulomb(1:this%maxnorbs, 1:this%maxnorbs, 1:this%maxnorbs, 1:this%maxnorbs, & 51 1:st%d%dim, 1:st%d%dim, 1:this%norbsets) = M_ZERO 52 53 !Lets counts the number of orbital to treat, to display a progress bar 54 ntodo = 0 55 do ios = this%orbs_dist%start, this%orbs_dist%end 56 norbs = this%orbsets(ios)%norbs 57 ntodo= ntodo + norbs**4*4!((norbs+1)*norbs/2)*((norbs+1)*norbs/2+1)/2 58 end do 59 idone = 0 60 if(mpi_world%rank == 0) call loct_progress_bar(-1, ntodo) 61 62 63 do ios = this%orbs_dist%start, this%orbs_dist%end 64 os => this%orbsets(ios) 65 norbs = os%norbs 66 np_sphere = os%sphere%np 67 68 call submesh_build_global(os%sphere) 69 70 call poisson_init_sm(os%poisson, namespace, psolver, der, os%sphere) 71 72 ijst=0 73 do ist = 1, norbs 74 75 do jst = 1, norbs 76 ijst=ijst+1 77 78 do is1 = 1, st%d%dim 79 !$omp parallel do 80 do ip = 1,np_sphere 81 nn(ip,is1) = conjg(os%zorb(ip,is1,ist))*os%zorb(ip,is1,jst) 82 end do 83 !$omp end parallel do 84 85 !Here it is important to use a non-periodic poisson solver, e.g. the direct solver 86 call zpoisson_solve_sm(os%poisson, os%sphere, vv(1:np_sphere,is1), nn(1:np_sphere,is1)) 87 end do !is1 88 89 klst=0 90 do kst = 1, norbs 91 do lst = 1, norbs 92 klst=klst+1 93 94 do is1 = 1, st%d%dim 95 do is2 = 1, st%d%dim 96 97 !$omp parallel do 98 do ip = 1,np_sphere 99 tmp(ip) = vv(ip,is1)*conjg(os%zorb(ip,is2,kst))*os%zorb(ip,is2,lst) 100 end do 101 !$omp end parallel do 102 103 this%zcoulomb(ist,jst,kst,lst,is1,is2,ios) = zsm_integrate(mesh, os%sphere, tmp(1:np_sphere)) 104 end do !is2 105 end do !is1 106 107 do is1 = 1, st%d%dim 108 do is2 = 1, st%d%dim 109 if(abs(this%zcoulomb(ist,jst,kst,lst,is1,is2,ios))<CNST(1.0e-12)) then 110 this%zcoulomb(ist,jst,kst,lst,is1,is2,ios) = M_ZERO 111 end if 112 113 end do !is2 114 end do !is1 115 116 idone = idone + 1 117 if(mpi_world%rank == 0) call loct_progress_bar(idone, ntodo) 118 end do !lst 119 end do !kst 120 end do !jst 121 end do !ist 122 call poisson_end(os%poisson) 123 124 call submesh_end_global(os%sphere) 125 end do !iorb 126 127 if(this%orbs_dist%parallel) then 128 do ios = 1, this%norbsets 129 do is2 = 1, st%d%dim 130 do is1 = 1, st%d%dim 131 call comm_allreduce(this%orbs_dist%mpi_grp%comm, this%zcoulomb(:,:,:,:,is1,is2,ios)) 132 end do 133 end do 134 end do 135 end if 136 137 SAFE_DEALLOCATE_A(nn) 138 SAFE_DEALLOCATE_A(vv) 139 SAFE_DEALLOCATE_A(tmp) 140 141 POP_SUB(compute_complex_coulomb_integrals) 142 call profiling_out(prof) 143end subroutine compute_complex_coulomb_integrals 144 145! --------------------------------------------------------- 146!> This routine computes the effective U in the non-collinear case 147! --------------------------------------------------------- 148subroutine compute_ACBNO_U_noncollinear(this, ios, namespace) 149 type(lda_u_t), intent(inout) :: this 150 integer, intent(in) :: ios 151 type(namespace_t), intent(in) :: namespace 152 153 integer :: im, imp, impp, imppp, ispin1, ispin2, norbs 154 CMPLX :: numU, numJ, tmpU, tmpJ, denomU, denomJ 155 156 PUSH_SUB(compute_ACBNO_U_noncollinear) 157 158 norbs = this%orbsets(ios)%norbs 159 numU = M_z0 160 numJ = M_z0 161 denomU = M_z0 162 denomJ = M_z0 163 164 if(norbs > 1) then 165 166 do im = 1, norbs 167 do imp = 1,norbs 168 do impp = 1, norbs 169 do imppp = 1, norbs 170 ! We first compute the terms 171 ! sum_{alpha,beta} P^alpha_{mmp}P^beta_{mpp,mppp} 172 ! sum_{alpha} P^alpha_{mmp}P^alpha_{mpp,mppp} 173 tmpU = M_z0 174 tmpJ = M_z0 175 176 do ispin1 = 1, this%spin_channels 177 do ispin2 = 1, this%spin_channels 178 tmpU = tmpU + this%zn_alt(im,imp,ispin1,ios)*this%zn_alt(impp,imppp,ispin2,ios)& 179 * this%zcoulomb(im,imp,impp,imppp,ispin1,ispin2,ios) 180 end do 181 tmpJ = tmpJ + this%zn_alt(im,imp,ispin1,ios)*this%zn_alt(impp,imppp,ispin1,ios)& 182 * this%zcoulomb(im,imppp,impp,imp,ispin1,ispin1,ios) 183 end do 184 tmpJ = tmpJ + this%zn_alt(im,imp,3,ios)*this%zn_alt(impp,imppp,4,ios) & 185 * this%zcoulomb(im,imppp,impp,imp,1,2,ios) & 186 +this%zn_alt(im,imp,4,ios)*this%zn_alt(impp,imppp,3,ios) & 187 * this%zcoulomb(im,imppp,impp,imp,2,1,ios) 188 ! These are the numerator of the ACBN0 U and J 189 numU = numU + tmpU 190 numJ = numJ + tmpJ 191 end do 192 end do 193 194 ! We compute the term 195 ! sum_{alpha} sum_{m,mp/=m} N^alpha_{m}N^alpha_{mp} 196 tmpJ = M_z0 197 tmpU = M_z0 198 if(imp/=im) then 199 do ispin1 = 1, this%spin_channels 200 tmpJ = tmpJ + this%zn(im,im,ispin1,ios)*this%zn(imp,imp,ispin1,ios) 201 tmpU = tmpU + this%zn(im,im,ispin1,ios)*this%zn(imp,imp,ispin1,ios) 202 end do 203 tmpJ = tmpJ + this%zn(im,im,3,ios)*this%zn(imp,imp,4,ios) & 204 + this%zn(im,im,4,ios)*this%zn(imp,imp,3,ios) 205 end if 206 denomJ = denomJ + tmpJ 207 208 ! We compute the term 209 ! sum_{alpha,beta} sum_{m,mp} N^alpha_{m}N^beta_{mp} 210 do ispin1 = 1, this%spin_channels 211 do ispin2 = 1, this%spin_channels 212 if(ispin1 /= ispin2) then 213 tmpU = tmpU + this%zn(im,im,ispin1,ios)*this%zn(imp,imp,ispin2,ios) 214 end if 215 end do 216 end do 217 218 if(im == imp) then 219 tmpU = tmpU - (this%zn(im,im,3,ios)*this%zn(im,im,4,ios) & 220 +this%zn(im,im,4,ios)*this%zn(im,im,3,ios)) 221 end if 222 denomU = denomU + tmpU 223 end do 224 end do 225 this%orbsets(ios)%Ueff = TOFLOAT(numU)/TOFLOAT(denomU) - TOFLOAT(numJ)/TOFLOAT(denomJ) 226 this%orbsets(ios)%Ubar = TOFLOAT(numU)/TOFLOAT(denomU) 227 this%orbsets(ios)%Jbar = TOFLOAT(numJ)/TOFLOAT(denomJ) 228 229 else !In the case of s orbitals, the expression is different 230 ! sum_{alpha/=beta} P^alpha_{mmp}P^beta_{mpp,mppp} 231 ! sum_{alpha,beta} sum_{m,mp} N^alpha_{m}N^beta_{mp} 232 numU = M_z0 233 denomU = M_z0 234 do ispin1 = 1, this%spin_channels 235 do ispin2 = 1, this%spin_channels 236 if(ispin1 /= ispin2) then 237 numU = numU + this%zn_alt(1,1,ispin1,ios)*this%zn_alt(1,1,ispin2,ios) & 238 *this%zcoulomb(1,1,1,1,ispin1,ispin2,ios) 239 denomU = denomU + this%zn(1,1,ispin1,ios)*this%zn(1,1,ispin2,ios) 240 end if 241 end do 242 end do 243 denomU = denomU + this%zn(1,1,3,ios)*this%zn(1,1,4,ios) & 244 +this%zn(1,1,4,ios)*this%zn(1,1,3,ios) 245 246 ! We have to be careful in the case of hydrogen atom for instance 247 if(abs(denomU)> CNST(1.0e-3)) then 248 this%orbsets(ios)%Ubar = (TOFLOAT(numU)/TOFLOAT(denomU)) 249 else 250 write(message(1),'(a,a)')' Small denominator value for the s orbital ', this%orbsets(ios)%Ubar 251 write(message(2),'(a)')' U is set to zero ' 252 call messages_warning(2, namespace=namespace) 253 this%orbsets(ios)%Ubar = M_ZERO 254 end if 255 this%orbsets(ios)%Jbar = 0 256 this%orbsets(ios)%Ueff = this%orbsets(ios)%Ubar 257 end if 258 259 POP_SUB(compute_ACBNO_U_noncollinear) 260end subroutine compute_ACBNO_U_noncollinear 261 262 263