PROGRESS  master
prg_charges_mod.F90
Go to the documentation of this file.
1 
5 module prg_charges_mod
6 
7  use prg_openfiles_mod
8  use bml
9  use iso_c_binding
10  use prg_parallel_mod
11  use prg_graph_mod
12  use prg_system_mod
13 
14 #ifdef USE_NVTX
15  use prg_nvtx_mod
16 #endif
17 
18  implicit none
19 
20  private
21 
22  integer, parameter :: dp = kind(1.0d0)
23 
24  public :: prg_get_charges, prg_get_hscf, prg_get_hscf_v2, prg_get_hscf_maxorb
25 
26 contains
27 
35  subroutine prg_get_charges(rho_bml,over_bml,hindex,charges,numel,spindex,mdimin,threshold)
36  character(20) :: bml_type, bml_mode
37  integer :: i, j, nats, norb, nsp
38  integer :: mdim, ld
39  integer, intent(in) :: hindex(:,:), mdimin, spindex(:)
40  real(dp) :: znuc, this_charge, this_diag
41  real(dp), allocatable :: rho_diag(:)
42  real(dp), allocatable, intent(inout) :: charges(:)
43  real(dp), intent(in) :: numel(:), threshold
44 #ifdef USE_OFFLOAD
45  real(c_double), pointer :: rho_bml_ptr(:,:), over_bml_ptr(:,:)
46  type(c_ptr) :: rho_bml_c_ptr, over_bml_c_ptr
47 #endif
48  type(bml_matrix_t) :: aux_bml
49  type(bml_matrix_t), intent(inout) :: over_bml, rho_bml
50 
51  norb = bml_get_n(rho_bml)
52  bml_type = bml_get_type(rho_bml)
53  bml_mode = bml_get_distribution_mode(rho_bml)
54  nats = size(hindex,dim=2)
55  nsp = size(spindex)
56 
57  if(mdimin.lt.0)then
58  mdim = norb
59  else
60  mdim = mdimin
61  endif
62 
63  if(.not.allocated(charges)) allocate(charges(nats))
64 
65  if(.not.allocated(rho_diag))allocate(rho_diag(norb))
66 
67 #ifdef USE_OFFLOAD
68  rho_bml_c_ptr = bml_get_data_ptr_dense(rho_bml)
69  over_bml_c_ptr = bml_get_data_ptr_dense(over_bml)
70  ld = bml_get_ld_dense(rho_bml)
71  call c_f_pointer(rho_bml_c_ptr,rho_bml_ptr,shape=[ld,norb])
72  call c_f_pointer(over_bml_c_ptr,over_bml_ptr,shape=[ld,norb])
73  charges = 0.0_dp
74  !$acc enter data copyin(charges(1:nats),numel(1:nsp),hindex(1:2,1:nats)) &
75  !$acc copyin(rho_diag(1:norb))
76 
77  !$acc parallel loop deviceptr(rho_bml_ptr, over_bml_ptr) &
78  !$acc present(rho_diag) &
79  !$acc private(i,j,this_diag)
80  do i = 1,norb
81  this_diag = 0.0_dp
82  !$acc loop reduction(+:this_diag)
83  do j = 1,norb
84  this_diag = this_diag + over_bml_ptr(i,j)*rho_bml_ptr(j,i) ! Transpose for C pointer
85  enddo
86  !$acc end loop
87  rho_diag(i) = this_diag
88  enddo
89  !$acc end parallel loop
90 
91  !$acc parallel loop present(rho_diag) &
92  !$acc present(charges,numel,hindex) &
93  !$acc private(i,j,this_charge)
94  do i = 1,nats
95  this_charge = -numel(spindex(i))
96  !$acc loop reduction(+:this_charge)
97  do j = hindex(1,i),hindex(2,i)
98  this_charge = this_charge + rho_diag(j)
99  enddo
100  !$acc end loop
101  charges(i) = this_charge
102  enddo
103  !$acc end parallel loop
104  !$acc exit data copyout(charges(1:nats)) &
105  !$acc delete(numel(1:nsp),hindex(1:2,1:nats),rho_diag(1:norb))
106 #else
107  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux_bml, &
108  bml_mode)
109 
110  call bml_multiply(rho_bml,over_bml,aux_bml,1.0_dp,0.0_dp,threshold)
111 
112 #ifdef DO_MPI
113  if (getnranks() .gt. 1 .and. &
114  bml_get_distribution_mode(aux_bml) == bml_dmode_distributed) then
115  call prg_allgatherparallel(aux_bml)
116  endif
117 #endif
118 
119  call bml_get_diagonal(aux_bml,rho_diag)
120  !$omp parallel do default(none) private(i) &
121  !$omp private(j,znuc) &
122  !$omp shared(spindex,numel,charges,hindex,rho_diag,nats)
123  do i = 1,nats
124  znuc = numel(spindex(i))
125  charges(i)=0.0_dp
126  do j = hindex(1,i),hindex(2,i)
127  charges(i) = charges(i) + rho_diag(j)
128  enddo
129  charges(i) = charges(i) - znuc
130  enddo
131  !$omp end parallel do
132  call bml_deallocate(aux_bml)
133 #endif
134  deallocate(rho_diag)
135  end subroutine prg_get_charges
136 
151  subroutine prg_get_hscf(ham0_bml,over_bml,ham_bml,spindex,hindex,hubbardu,charges,&
152  coulomb_pot_r,coulomb_pot_k,mdimin,threshold)
153  character(20) :: bml_type
154  integer :: i, j, nats, norb, mdim
155  integer, intent(in) :: hindex(:,:), mdimin, spindex(:)
156  real(dp), allocatable :: coulomb_pot(:), diagonal(:)
157  real(dp), intent(in) :: charges(:), coulomb_pot_k(:), coulomb_pot_r(:), hubbardu(:)
158  real(dp), intent(in) :: threshold
159  type(bml_matrix_t), save :: aux_bml
160  type(bml_matrix_t), intent(in) :: ham0_bml, over_bml
161  type(bml_matrix_t), intent(inout) :: ham_bml
162  integer, save :: maxnorb = 0
163 #ifdef USE_OFFLOAD
164  type(c_ptr) :: aux_bml_c_ptr
165  integer :: ld
166  real(c_double), pointer :: aux_bml_ptr(:,:)
167 #endif
168 
169  nats = size(coulomb_pot_r)
170 
171  norb = bml_get_n(ham0_bml)
172  if(hindex(2,nats) .ne. norb)then
173  write(*,*)"ERROR in prg_get_hscf. hindex incompatible with ham0"
174  stop
175  endif
176 
177  if(mdimin.lt.0.or.mdimin.gt.norb)then
178  mdim = norb
179  else
180  mdim = mdimin
181  endif
182 
183  bml_type = bml_get_type(ham0_bml)
184 
185  ! if(bml_allocated(ham_bml))then
186  ! call bml_deallocate(ham_bml)
187  ! endif
188  ! call bml_copy_new(ham0_bml,ham_bml)
189 
190  if(bml_allocated(ham_bml))then
191  if(norb.ne.bml_get_n(ham_bml))then
192  call bml_deallocate(ham_bml)
193  call bml_copy_new(ham0_bml,ham_bml)
194  else
195  call bml_copy(ham0_bml,ham_bml)
196  endif
197  else
198  call bml_copy_new(ham0_bml,ham_bml)
199  endif
200 
201 
202 
203 #ifdef USE_OFFLOAD
204  if(.not.bml_allocated(aux_bml))then
205  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux_bml)
206  elseif(norb.gt.maxnorb)then
207  call bml_deallocate(aux_bml)
208  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux_bml)
209  maxnorb = norb
210  else
211  call bml_set_n_dense(aux_bml,norb)
212  endif
213 
214  aux_bml_c_ptr = bml_get_data_ptr_dense(aux_bml)
215  ld = bml_get_ld_dense(aux_bml)
216 
217  call c_f_pointer(aux_bml_c_ptr,aux_bml_ptr,shape=[ld,norb])
218 
219  !$acc enter data copyin(hindex(:,:),hubbardu(:),spindex(:),charges(:)) &
220  !$acc copyin(coulomb_pot_r(:),coulomb_pot_k(:))
221 
222  !$acc parallel loop gang vector collapse(2) &
223  !$acc deviceptr(aux_bml_ptr)
224  do i = 1,norb
225  do j = 1,norb
226  aux_bml_ptr(i,j) = 0.0_dp
227  enddo
228  enddo
229 
230  !$acc parallel loop gang vector &
231  !$acc deviceptr(aux_bml_ptr) &
232  !$acc present(hindex,hubbardu,spindex,charges,coulomb_pot_r,coulomb_pot_k)
233  do i = 1,nats
234  do j = hindex(1,i),hindex(2,i)
235  aux_bml_ptr(j,j) = hubbardu(spindex(i))*charges(i) &
236  + coulomb_pot_r(i) + coulomb_pot_k(i)
237  enddo
238  enddo
239  !$acc exit data delete(hindex(:,:),hubbardu(:),spindex(:),charges(:)) &
240  !$acc delete(coulomb_pot_r(:),coulomb_pot_k(:))
241 
242  call bml_multiply(over_bml,aux_bml,ham_bml,0.5_dp,1.0_dp,threshold) ! h = h + 0.5*s*h1
243 
244  call bml_multiply(aux_bml,over_bml,ham_bml,0.5_dp,1.0_dp,threshold) ! h = h + 0.5*h1*s
245 
246 #else
247  call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux_bml)
248 
249  allocate(coulomb_pot(nats))
250 
251  coulomb_pot = coulomb_pot_k + coulomb_pot_r
252 
253  allocate(diagonal(norb))
254 
255  do i = 1,nats
256  do j = hindex(1,i),hindex(2,i)
257  diagonal(j) = hubbardu(spindex(i))*charges(i) + coulomb_pot(i)
258  enddo
259  enddo
260 
261  call bml_set_diagonal(aux_bml,diagonal)
262 
263  deallocate(diagonal)
264  deallocate(coulomb_pot)
265 
266  call bml_multiply(over_bml,aux_bml,ham_bml,0.5_dp,1.0_dp,threshold) ! h = h + 0.5*s*h1
267 
268  call bml_multiply(aux_bml,over_bml,ham_bml,0.5_dp,1.0_dp,threshold) ! h = h + 0.5*h1*s
269 
270  call bml_deallocate(aux_bml)
271 
272 #endif
273 
274  end subroutine prg_get_hscf
275 
276 
277  subroutine prg_get_hscf_v2(ham0_bml,over_bml,ham_bml,spindex,hindex,hubbardu,charges,&
278  coulomb_pot_r,coulomb_pot_k,mdimin,threshold)
279  character(20) :: bml_type
280  integer :: i, j, nats, norb, mdim
281  integer, intent(in) :: hindex(:,:), mdimin, spindex(:)
282  real(dp), allocatable :: coulomb_pot(:), diagonal(:)
283  real(dp), intent(in) :: charges(:), coulomb_pot_k(:),coulomb_pot_r(:), hubbardu(:)
284  real(dp), intent(in) :: threshold
285  type(bml_matrix_t) :: aux_bml
286  type(bml_matrix_t), intent(in) :: ham0_bml, over_bml
287  type(bml_matrix_t), intent(inout) :: ham_bml
288 
289  nats = size(coulomb_pot_r)
290  allocate(coulomb_pot(nats))
291 
292  norb = bml_get_n(ham0_bml)
293  if(hindex(2,nats) .ne. norb)then
294  write(*,*)"ERROR in prg_get_hscf. hindex incompatible with ham0"
295  stop
296  endif
297 
298  if(mdimin.lt.0)then
299  mdim = norb
300  else
301  mdim = mdimin
302  endif
303 
304  allocate(diagonal(norb))
305 
306  if (bml_allocated(ham_bml))then
307  call bml_deallocate(ham_bml)
308  endif
309 
310  bml_type = bml_get_type(ham0_bml)
311 
312  coulomb_pot = coulomb_pot_k + coulomb_pot_r
313 
314  if(bml_allocated(ham_bml))then
315  call bml_deallocate(ham_bml)
316  endif
317 
318  call bml_zero_matrix(bml_type,bml_element_real,dp,mdim,norb,ham_bml)
319  call bml_zero_matrix(bml_type,bml_element_real,dp,mdim,norb,aux_bml)
320 
321  do i = 1,nats
322  do j = hindex(1,i),hindex(2,i)
323  diagonal(j) = hubbardu(spindex(i))*charges(i) + coulomb_pot(i)
324  enddo
325  enddo
326 
327  call bml_set_diagonal(aux_bml,diagonal)
328  call bml_multiply(over_bml,aux_bml,ham_bml,0.5_dp,1.0_dp,threshold) ! h = h+ 0.5*s*h1
329 
330  call bml_multiply(aux_bml,over_bml,ham_bml,0.5_dp,1.0_dp,threshold) ! h = h+ 0.5*h1*s
331 
332  deallocate(diagonal)
333  deallocate(coulomb_pot)
334  call bml_deallocate(aux_bml)
335 
336  end subroutine prg_get_hscf_v2
337 
338  subroutine prg_get_hscf_maxorb(ham0_bml,over_bml,ham_bml,spindex,hindex,hubbardu,charges,&
339  coulomb_pot_r,coulomb_pot_k,mdimin,threshold,maxorb)
340  character(20) :: bml_type
341  integer :: i, j, nats, norb, mdim
342  integer, intent(in) :: hindex(:,:), mdimin, spindex(:), maxorb
343  real(dp), allocatable :: coulomb_pot(:), diagonal(:)
344  real(dp), intent(in) :: charges(:), coulomb_pot_k(:),coulomb_pot_r(:), hubbardu(:)
345  real(dp), intent(in) :: threshold
346  real(dp), allocatable, save :: zero(:,:)
347  type(bml_matrix_t), save :: zero_bml, aux_bml
348  type(bml_matrix_t), intent(in) :: ham0_bml, over_bml
349  type(bml_matrix_t), intent(inout) :: ham_bml
350 
351  nats = size(coulomb_pot_r)
352  allocate(coulomb_pot(nats))
353 
354  bml_type = bml_get_type(ham_bml)
355 
356  if(.not.bml_allocated(zero_bml))then
357  allocate(zero(maxorb,maxorb))
358  zero = 0.0_dp
359  call bml_import_from_dense(bml_type,zero,zero_bml)
360  call bml_import_from_dense(bml_type,zero,aux_bml)
361  endif
362 
363  norb = bml_get_n(ham0_bml)
364  if(hindex(2,nats) .ne. norb)then
365  write(*,*)"ERROR in prg_get_hscf. hindex incompatible with ham0"
366  stop
367  endif
368 
369  if(mdimin.ne.norb)then
370  write(*,*)"ERROR in prg_ge_hscf_maxorb: mdim must equal norb"
371  stop
372  endif
373 
374  call bml_set_n_dense(aux_bml,norb)
375  call bml_set_n_dense(zero_bml,norb)
376  call bml_copy(zero_bml,aux_bml)
377 
378  allocate(diagonal(norb))
379 
380  call bml_copy(ham0_bml,ham_bml)
381 
382  bml_type = bml_get_type(ham_bml)
383 
384  coulomb_pot = coulomb_pot_k + coulomb_pot_r
385 
386 
387  do i = 1,nats
388  do j = hindex(1,i),hindex(2,i)
389  diagonal(j) = hubbardu(spindex(i))*charges(i) + coulomb_pot(i)
390  enddo
391  enddo
392 
393  call bml_set_diagonal(aux_bml,diagonal)
394  call bml_multiply(over_bml,aux_bml,ham_bml,0.5_dp,1.0_dp,threshold) ! h = h+ 0.5*s*h1
395 
396  call bml_multiply(aux_bml,over_bml,ham_bml,0.5_dp,1.0_dp,threshold) ! h = h+ 0.5*h1*s
397 
398  deallocate(diagonal)
399  deallocate(coulomb_pot)
400  !call bml_deallocate(aux_bml)
401 
402  end subroutine prg_get_hscf_maxorb
403 
404 
405 end module prg_charges_mod
prg_charges_mod
A module to compute the Mulliken charges of a chemical system.
Definition: prg_charges_mod.F90:5
prg_charges_mod::prg_get_hscf
subroutine, public prg_get_hscf(ham0_bml, over_bml, ham_bml, spindex, hindex, hubbardu, charges, coulomb_pot_r, coulomb_pot_k, mdimin, threshold)
Constructs the SCF Hamiltonian given H0, HubbardU and charges. This routine does: ,...
Definition: prg_charges_mod.F90:153
prg_charges_mod::prg_get_hscf_v2
subroutine, public prg_get_hscf_v2(ham0_bml, over_bml, ham_bml, spindex, hindex, hubbardu, charges, coulomb_pot_r, coulomb_pot_k, mdimin, threshold)
Definition: prg_charges_mod.F90:279
prg_charges_mod::prg_get_charges
subroutine, public prg_get_charges(rho_bml, over_bml, hindex, charges, numel, spindex, mdimin, threshold)
Constructs the charges from the density matrix.
Definition: prg_charges_mod.F90:36
prg_charges_mod::prg_get_hscf_maxorb
subroutine, public prg_get_hscf_maxorb(ham0_bml, over_bml, ham_bml, spindex, hindex, hubbardu, charges, coulomb_pot_r, coulomb_pot_k, mdimin, threshold, maxorb)
Definition: prg_charges_mod.F90:340
prg_charges_mod::dp
integer, parameter dp
Definition: prg_charges_mod.F90:22
prg_graph_mod
The graph module.
Definition: prg_graph_mod.F90:6
prg_nvtx_mod
Definition: prg_nvtx_mod.F90:1
prg_openfiles_mod
Module to handle input output files for the PROGRESS lib.
Definition: prg_openfiles_mod.F90:4
prg_parallel_mod
The parallel module.
Definition: prg_parallel_mod.F90:6
prg_parallel_mod::prg_allgatherparallel
subroutine, public prg_allgatherparallel(a)
Definition: prg_parallel_mod.F90:774
prg_parallel_mod::getnranks
integer function, public getnranks()
Definition: prg_parallel_mod.F90:90
prg_system_mod
A module to read and handle chemical systems.
Definition: prg_system_mod.F90:5