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