(git:1f285aa)
core_ppnl.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 ! **************************************************************************************************
8 !> \brief Calculation of the non-local pseudopotential contribution to the core Hamiltonian
9 !> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
10 !> \par History
11 !> - refactered from qs_core_hamiltian [Joost VandeVondele, 2008-11-01]
12 !> - full rewrite [jhu, 2009-01-23]
13 !> - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
14 ! **************************************************************************************************
15 MODULE core_ppnl
16  USE ai_angmom, ONLY: angmom
17  USE ai_overlap, ONLY: overlap
18  USE atomic_kind_types, ONLY: atomic_kind_type,&
19  get_atomic_kind_set
20  USE basis_set_types, ONLY: gto_basis_set_p_type,&
21  gto_basis_set_type
22  USE dbcsr_api, ONLY: dbcsr_add,&
23  dbcsr_get_block_p,&
24  dbcsr_p_type
25  USE external_potential_types, ONLY: gth_potential_p_type,&
26  gth_potential_type,&
27  sgp_potential_p_type,&
28  sgp_potential_type
29  USE kinds, ONLY: dp,&
30  int_8
31  USE orbital_pointers, ONLY: init_orbital_pointers,&
32  nco,&
33  ncoset
34  USE particle_types, ONLY: particle_type
35  USE qs_force_types, ONLY: qs_force_type
36  USE qs_kind_types, ONLY: get_qs_kind,&
37  get_qs_kind_set,&
38  qs_kind_type
39  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
40  USE sap_kind_types, ONLY: alist_type,&
41  clist_type,&
42  get_alist,&
43  release_sap_int,&
44  sap_int_type,&
45  sap_sort
46  USE virial_methods, ONLY: virial_pair_force
47  USE virial_types, ONLY: virial_type
48 
49 !$ USE OMP_LIB, ONLY: omp_lock_kind, &
50 !$ omp_init_lock, omp_set_lock, &
51 !$ omp_unset_lock, omp_destroy_lock
52 
53 #include "./base/base_uses.f90"
54 
55  IMPLICIT NONE
56 
57  PRIVATE
58 
59  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'core_ppnl'
60 
61  PUBLIC :: build_core_ppnl
62 
63 CONTAINS
64 
65 ! **************************************************************************************************
66 !> \brief ...
67 !> \param matrix_h ...
68 !> \param matrix_p ...
69 !> \param force ...
70 !> \param virial ...
71 !> \param calculate_forces ...
72 !> \param use_virial ...
73 !> \param nder ...
74 !> \param qs_kind_set ...
75 !> \param atomic_kind_set ...
76 !> \param particle_set ...
77 !> \param sab_orb ...
78 !> \param sap_ppnl ...
79 !> \param eps_ppnl ...
80 !> \param nimages ...
81 !> \param cell_to_index ...
82 !> \param basis_type ...
83 !> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
84 !> \param matrix_l ...
85 ! **************************************************************************************************
86  SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
87  qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
88  nimages, cell_to_index, basis_type, deltaR, matrix_l)
89 
90  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
91  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
92  TYPE(virial_type), POINTER :: virial
93  LOGICAL, INTENT(IN) :: calculate_forces
94  LOGICAL :: use_virial
95  INTEGER :: nder
96  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
97  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
98  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
99  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
100  POINTER :: sab_orb, sap_ppnl
101  REAL(kind=dp), INTENT(IN) :: eps_ppnl
102  INTEGER, INTENT(IN) :: nimages
103  INTEGER, DIMENSION(:, :, :), OPTIONAL, POINTER :: cell_to_index
104  CHARACTER(LEN=*), INTENT(IN) :: basis_type
105  REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
106  OPTIONAL :: deltar
107  TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
108  POINTER :: matrix_l
109 
110  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_core_ppnl'
111 
112  INTEGER :: atom_a, first_col, handle, i, i_dim, iab, iac, iatom, ib, ibc, icol, ikind, &
113  ilist, img, irow, iset, j, jatom, jb, jkind, jneighbor, kac, katom, kbc, kkind, l, &
114  lc_max, lc_min, ldai, ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, maxppnl, &
115  maxsgf, na, natom, nb, ncoa, ncoc, nkind, nlist, nneighbor, nnl, np, nppnl, nprjc, nseta, &
116  nsgfa, prjc, sgfa, slot
117  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
118  INTEGER, DIMENSION(3) :: cell_b, cell_c
119  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
120  nsgf_seta
121  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
122  LOGICAL :: do_dr, do_gth, do_kp, do_soc, found, &
123  ppnl_present
124  REAL(kind=dp) :: dac, f0, ppnl_radius
125  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: radp
126  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
127  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work, lab, work_l
128  REAL(kind=dp), DIMENSION(1) :: rprjc, zetc
129  REAL(kind=dp), DIMENSION(3) :: fa, fb, rab, rac, rbc
130  REAL(kind=dp), DIMENSION(3, 3) :: pv_thread
131  TYPE(gto_basis_set_type), POINTER :: orb_basis_set
132  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
133  TYPE(gth_potential_type), POINTER :: gth_potential
134  TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
135  TYPE(clist_type), POINTER :: clist
136  TYPE(alist_type), POINTER :: alist_ac, alist_bc
137  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, alkint, bchint, bcint, &
138  blkint
139  REAL(kind=dp), DIMENSION(:, :), POINTER :: cprj, h_block, l_block_x, l_block_y, &
140  l_block_z, p_block, r_2block, &
141  r_3block, rpgfa, sphi_a, vprj_ppnl, &
142  wprj_ppnl, zeta
143  REAL(kind=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
144  REAL(kind=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
145  TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
146  TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
147  TYPE(sgp_potential_type), POINTER :: sgp_potential
148 
149 !$ INTEGER(kind=omp_lock_kind), &
150 !$ ALLOCATABLE, DIMENSION(:) :: locks
151 !$ INTEGER(KIND=int_8) :: iatom8
152 !$ INTEGER :: lock_num, hash
153 !$ INTEGER, PARAMETER :: nlock = 501
154 
155  mark_used(int_8)
156 
157  do_dr = .false.
158  IF (PRESENT(deltar)) do_dr = .true.
159 
160  IF (calculate_forces) THEN
161  CALL timeset(routinen//"_forces", handle)
162  ELSE
163  CALL timeset(routinen, handle)
164  END IF
165 
166  do_soc = PRESENT(matrix_l)
167 
168  ppnl_present = ASSOCIATED(sap_ppnl)
169 
170  IF (ppnl_present) THEN
171 
172  nkind = SIZE(atomic_kind_set)
173  natom = SIZE(particle_set)
174 
175  do_kp = (nimages > 1)
176 
177  IF (do_kp) THEN
178  cpassert(PRESENT(cell_to_index) .AND. ASSOCIATED(cell_to_index))
179  END IF
180 
181  IF (calculate_forces) THEN
182  IF (SIZE(matrix_p, 1) == 2) THEN
183  DO img = 1, nimages
184  CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
185  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
186  CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
187  alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
188  END DO
189  END IF
190  END IF
191 
192  maxder = ncoset(nder)
193 
194  CALL get_qs_kind_set(qs_kind_set, &
195  maxco=maxco, &
196  maxlgto=maxlgto, &
197  maxsgf=maxsgf, &
198  maxlppnl=maxlppnl, &
199  maxppnl=maxppnl, &
200  basis_type=basis_type)
201 
202  maxl = max(maxlgto, maxlppnl)
203  CALL init_orbital_pointers(maxl + nder + 1)
204 
205  ldsab = max(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
206  ldai = ncoset(maxl + nder + 1)
207 
208  ! sap_int needs to be shared as multiple threads need to access this
209  ALLOCATE (sap_int(nkind*nkind))
210  DO i = 1, nkind*nkind
211  NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
212  sap_int(i)%nalist = 0
213  END DO
214 
215  ! Set up direct access to basis and potential
216  ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind))
217  DO ikind = 1, nkind
218  CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
219  IF (ASSOCIATED(orb_basis_set)) THEN
220  basis_set(ikind)%gto_basis_set => orb_basis_set
221  ELSE
222  NULLIFY (basis_set(ikind)%gto_basis_set)
223  END IF
224  CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
225  NULLIFY (gpotential(ikind)%gth_potential)
226  NULLIFY (spotential(ikind)%sgp_potential)
227  IF (ASSOCIATED(gth_potential)) THEN
228  gpotential(ikind)%gth_potential => gth_potential
229  IF (do_soc .AND. (.NOT. gth_potential%soc)) THEN
230  cpwarn("Spin-orbit coupling selected, but GTH potential without SOC parameters provided")
231  END IF
232  ELSE IF (ASSOCIATED(sgp_potential)) THEN
233  spotential(ikind)%sgp_potential => sgp_potential
234  END IF
235  END DO
236 
237  ! Allocate sap int
238  DO slot = 1, sap_ppnl(1)%nl_size
239 
240  ikind = sap_ppnl(1)%nlist_task(slot)%ikind
241  kkind = sap_ppnl(1)%nlist_task(slot)%jkind
242  iatom = sap_ppnl(1)%nlist_task(slot)%iatom
243  katom = sap_ppnl(1)%nlist_task(slot)%jatom
244  nlist = sap_ppnl(1)%nlist_task(slot)%nlist
245  ilist = sap_ppnl(1)%nlist_task(slot)%ilist
246  nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
247 
248  iac = ikind + nkind*(kkind - 1)
249  IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
250  IF (.NOT. ASSOCIATED(gpotential(kkind)%gth_potential) .AND. &
251  .NOT. ASSOCIATED(spotential(kkind)%sgp_potential)) cycle
252  IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
253  sap_int(iac)%a_kind = ikind
254  sap_int(iac)%p_kind = kkind
255  sap_int(iac)%nalist = nlist
256  ALLOCATE (sap_int(iac)%alist(nlist))
257  DO i = 1, nlist
258  NULLIFY (sap_int(iac)%alist(i)%clist)
259  sap_int(iac)%alist(i)%aatom = 0
260  sap_int(iac)%alist(i)%nclist = 0
261  END DO
262  END IF
263  IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
264  sap_int(iac)%alist(ilist)%aatom = iatom
265  sap_int(iac)%alist(ilist)%nclist = nneighbor
266  ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
267  DO i = 1, nneighbor
268  sap_int(iac)%alist(ilist)%clist(i)%catom = 0
269  END DO
270  END IF
271  END DO
272 
273  ! Calculate the overlap integrals <a|p>
274 !$OMP PARALLEL &
275 !$OMP DEFAULT (NONE) &
276 !$OMP SHARED (basis_set, gpotential, spotential, maxder, ncoset, &
277 !$OMP sap_ppnl, sap_int, nkind, ldsab, ldai, nder, nco, do_soc ) &
278 !$OMP PRIVATE (ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
279 !$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, &
280 !$OMP slot, sphi_a, zeta, cprj, hprj, lppnl, nppnl, nprj_ppnl, &
281 !$OMP clist, iset, ncoa, sgfa, prjc, work, work_l, sab, lab, ai_work, nprjc, &
282 !$OMP ppnl_radius, ncoc, rpgfa, first_col, vprj_ppnl, wprj_ppnl, i, j, l, do_gth, &
283 !$OMP set_radius_a, rprjc, dac, lc_max, lc_min, zetc, alpha_ppnl, &
284 !$OMP na, nb, np, nnl, a_nl, radp, i_dim, ib, jb)
285 
286  ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder))
287  sab = 0.0_dp
288  ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1)))
289  ai_work = 0.0_dp
290  IF (do_soc) THEN
291  ALLOCATE (lab(ldsab, ldsab, 3), work_l(ldsab, ldsab, 3))
292  lab = 0.0_dp
293  END IF
294 
295 !$OMP DO SCHEDULE(GUIDED)
296  DO slot = 1, sap_ppnl(1)%nl_size
297 
298  ikind = sap_ppnl(1)%nlist_task(slot)%ikind
299  kkind = sap_ppnl(1)%nlist_task(slot)%jkind
300  iatom = sap_ppnl(1)%nlist_task(slot)%iatom
301  katom = sap_ppnl(1)%nlist_task(slot)%jatom
302  nlist = sap_ppnl(1)%nlist_task(slot)%nlist
303  ilist = sap_ppnl(1)%nlist_task(slot)%ilist
304  nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
305  jneighbor = sap_ppnl(1)%nlist_task(slot)%inode
306  cell_c(:) = sap_ppnl(1)%nlist_task(slot)%cell(:)
307  rac(1:3) = sap_ppnl(1)%nlist_task(slot)%r(1:3)
308 
309  iac = ikind + nkind*(kkind - 1)
310  IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
311  ! Get definition of basis set
312  first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf
313  la_max => basis_set(ikind)%gto_basis_set%lmax
314  la_min => basis_set(ikind)%gto_basis_set%lmin
315  npgfa => basis_set(ikind)%gto_basis_set%npgf
316  nseta = basis_set(ikind)%gto_basis_set%nset
317  nsgfa = basis_set(ikind)%gto_basis_set%nsgf
318  nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set
319  rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius
320  set_radius_a => basis_set(ikind)%gto_basis_set%set_radius
321  sphi_a => basis_set(ikind)%gto_basis_set%sphi
322  zeta => basis_set(ikind)%gto_basis_set%zet
323  ! Get definition of PP projectors
324  IF (ASSOCIATED(gpotential(kkind)%gth_potential)) THEN
325  ! GTH potential
326  do_gth = .true.
327  alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl
328  cprj => gpotential(kkind)%gth_potential%cprj
329  lppnl = gpotential(kkind)%gth_potential%lppnl
330  nppnl = gpotential(kkind)%gth_potential%nppnl
331  nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
332  ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
333  vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
334  wprj_ppnl => gpotential(kkind)%gth_potential%wprj_ppnl
335  ELSE IF (ASSOCIATED(spotential(kkind)%sgp_potential)) THEN
336  ! SGP potential
337  do_gth = .false.
338  nprjc = spotential(kkind)%sgp_potential%nppnl
339  IF (nprjc == 0) cycle
340  nnl = spotential(kkind)%sgp_potential%n_nonlocal
341  lppnl = spotential(kkind)%sgp_potential%lmax
342  a_nl => spotential(kkind)%sgp_potential%a_nonlocal
343  ppnl_radius = spotential(kkind)%sgp_potential%ppnl_radius
344  ALLOCATE (radp(nnl))
345  radp(:) = ppnl_radius
346  cprj => spotential(kkind)%sgp_potential%cprj_ppnl
347  hprj => spotential(kkind)%sgp_potential%vprj_ppnl
348  nppnl = SIZE(cprj, 2)
349  ELSE
350  cycle
351  END IF
352 
353  dac = sqrt(sum(rac*rac))
354  clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
355  clist%catom = katom
356  clist%cell = cell_c
357  clist%rac = rac
358  ALLOCATE (clist%acint(nsgfa, nppnl, maxder), &
359  clist%achint(nsgfa, nppnl, maxder), &
360  clist%alint(nsgfa, nppnl, 3), &
361  clist%alkint(nsgfa, nppnl, 3))
362  clist%acint = 0.0_dp
363  clist%achint = 0.0_dp
364  clist%alint = 0.0_dp
365  clist%alkint = 0.0_dp
366 
367  clist%nsgf_cnt = 0
368  NULLIFY (clist%sgf_list)
369  DO iset = 1, nseta
370  ncoa = npgfa(iset)*ncoset(la_max(iset))
371  sgfa = first_sgfa(1, iset)
372  IF (do_gth) THEN
373  ! GTH potential
374  prjc = 1
375  work = 0.0_dp
376  DO l = 0, lppnl
377  nprjc = nprj_ppnl(l)*nco(l)
378  IF (nprjc == 0) cycle
379  rprjc(1) = ppnl_radius
380  IF (set_radius_a(iset) + rprjc(1) < dac) cycle
381  lc_max = l + 2*(nprj_ppnl(l) - 1)
382  lc_min = l
383  zetc(1) = alpha_ppnl(l)
384  ncoc = ncoset(lc_max)
385 
386  ! Calculate the primitive overlap integrals
387  CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
388  lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, nder, .true., ai_work, ldai)
389  ! Transformation step projector functions (Cartesian -> spherical)
390  na = ncoa
391  nb = nprjc
392  np = ncoc
393  DO i = 1, maxder
394  first_col = (i - 1)*ldsab
395  ! CALL dgemm("N", "N", ncoa, nprjc, ncoc, 1.0_dp, sab(1, first_col + 1), SIZE(sab, 1), &
396  ! cprj(1, prjc), SIZE(cprj, 1), 0.0_dp, work(1, first_col + prjc), ldsab)
397  work(1:na, first_col + prjc:first_col + prjc + nb - 1) = &
398  matmul(sab(1:na, first_col + 1:first_col + np), cprj(1:np, prjc:prjc + nb - 1))
399  END DO
400 
401  IF (do_soc) THEN
402  ! Calculate the primitive angular momentum integrals needed for spin-orbit coupling
403  lab = 0.0_dp
404  CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
405  lc_max, 1, zetc, rprjc, -rac, (/0._dp, 0._dp, 0._dp/), lab)
406  DO i_dim = 1, 3
407  work_l(1:na, prjc:prjc + nb - 1, i_dim) = &
408  matmul(lab(1:na, 1:np, i_dim), cprj(1:np, prjc:prjc + nb - 1))
409  END DO
410  END IF
411 
412  prjc = prjc + nprjc
413 
414  END DO
415  na = nsgf_seta(iset)
416  nb = nppnl
417  np = ncoa
418  DO i = 1, maxder
419  first_col = (i - 1)*ldsab + 1
420  ! Contraction step (basis functions)
421  ! CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
422  ! work(1, first_col), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
423  clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
424  matmul(transpose(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, first_col:first_col + nb - 1))
425  ! Multiply with interaction matrix(h)
426  ! CALL dgemm("N", "N", nsgf_seta(iset), nppnl, nppnl, 1.0_dp, clist%acint(sgfa, 1, i), nsgfa, &
427  ! vprj_ppnl(1, 1), SIZE(vprj_ppnl, 1), 0.0_dp, clist%achint(sgfa, 1, i), nsgfa)
428  clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
429  matmul(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), vprj_ppnl(1:nb, 1:nb))
430  END DO
431  IF (do_soc) THEN
432  DO i_dim = 1, 3
433  clist%alint(sgfa:sgfa + na - 1, 1:nb, i_dim) = &
434  matmul(transpose(sphi_a(1:np, sgfa:sgfa + na - 1)), work_l(1:np, 1:nb, i_dim))
435  clist%alkint(sgfa:sgfa + na - 1, 1:nb, i_dim) = &
436  matmul(clist%alint(sgfa:sgfa + na - 1, 1:nb, i_dim), wprj_ppnl(1:nb, 1:nb))
437  END DO
438  END IF
439  ELSE
440  ! SGP potential
441  ! Calculate the primitive overlap integrals
442  CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
443  lppnl, 0, nnl, radp, a_nl, rac, dac, sab, nder, .true., ai_work, ldai)
444  na = nsgf_seta(iset)
445  nb = nppnl
446  np = ncoa
447  DO i = 1, maxder
448  first_col = (i - 1)*ldsab + 1
449  ! Transformation step projector functions (cartesian->spherical)
450  ! CALL dgemm("N", "N", ncoa, nppnl, nprjc, 1.0_dp, sab(1, first_col), ldsab, &
451  ! cprj(1, 1), SIZE(cprj, 1), 0.0_dp, work(1, 1), ldsab)
452  work(1:np, 1:nb) = matmul(sab(1:np, first_col:first_col + nprjc - 1), cprj(1:nprjc, 1:nb))
453  ! Contraction step (basis functions)
454  ! CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
455  ! work(1, 1), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
456  clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
457  matmul(transpose(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, 1:nb))
458  ! *** Multiply with interaction matrix(h) ***
459  ncoc = sgfa + nsgf_seta(iset) - 1
460  DO j = 1, nppnl
461  clist%achint(sgfa:ncoc, j, i) = clist%acint(sgfa:ncoc, j, i)*hprj(j)
462  END DO
463  END DO
464  END IF
465  END DO
466  clist%maxac = maxval(abs(clist%acint(:, :, 1)))
467  clist%maxach = maxval(abs(clist%achint(:, :, 1)))
468  IF (.NOT. do_gth) DEALLOCATE (radp)
469  END DO
470 
471  DEALLOCATE (sab, ai_work, work)
472  IF (do_soc) DEALLOCATE (lab, work_l)
473 !$OMP END PARALLEL
474 
475  ! Set up a sorting index
476  CALL sap_sort(sap_int)
477  ! All integrals needed have been calculated and stored in sap_int
478  ! We now calculate the Hamiltonian matrix elements
479 
480  force_thread = 0.0_dp
481  pv_thread = 0.0_dp
482 
483 !$OMP PARALLEL &
484 !$OMP DEFAULT (NONE) &
485 !$OMP SHARED (do_kp, basis_set, matrix_h, matrix_l, cell_to_index,&
486 !$OMP sab_orb, matrix_p, sap_int, nkind, eps_ppnl, force, &
487 !$OMP do_dR, deltaR, maxder, nder, &
488 !$OMP locks, virial, use_virial, calculate_forces, do_soc, natom) &
489 !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
490 !$OMP slot, iab, atom_a, f0, irow, icol, h_block, &
491 !$OMP l_block_x, l_block_y, l_block_z, &
492 !$OMP r_2block, r_3block, &
493 !$OMP found,p_block, iac, ibc, alist_ac, alist_bc, acint, bcint, &
494 !$OMP achint, bchint, alkint, blkint, &
495 !$OMP na, np, nb, katom, j, fa, fb, rbc, rac, &
496 !$OMP kkind, kac, kbc, i, img, hash, iatom8) &
497 !$OMP REDUCTION (+ : pv_thread, force_thread )
498 
499 !$OMP SINGLE
500 !$ ALLOCATE (locks(nlock))
501 !$OMP END SINGLE
502 
503 !$OMP DO
504 !$ DO lock_num = 1, nlock
505 !$ call omp_init_lock(locks(lock_num))
506 !$ END DO
507 !$OMP END DO
508 
509 !$OMP DO SCHEDULE(GUIDED)
510  DO slot = 1, sab_orb(1)%nl_size
511 
512  ikind = sab_orb(1)%nlist_task(slot)%ikind
513  jkind = sab_orb(1)%nlist_task(slot)%jkind
514  iatom = sab_orb(1)%nlist_task(slot)%iatom
515  jatom = sab_orb(1)%nlist_task(slot)%jatom
516  cell_b(:) = sab_orb(1)%nlist_task(slot)%cell(:)
517  rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
518 
519  IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
520  IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
521 
522  iab = ikind + nkind*(jkind - 1)
523 
524  ! Use the symmetry of the first derivatives
525  IF (iatom == jatom) THEN
526  f0 = 1.0_dp
527  ELSE
528  f0 = 2.0_dp
529  END IF
530 
531  IF (do_kp) THEN
532  img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
533  ELSE
534  img = 1
535  END IF
536 
537  ! Create matrix blocks for a new matrix block column
538  IF (iatom <= jatom) THEN
539  irow = iatom
540  icol = jatom
541  ELSE
542  irow = jatom
543  icol = iatom
544  END IF
545  NULLIFY (h_block)
546  CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found)
547  IF (do_soc) THEN
548  NULLIFY (l_block_x, l_block_y, l_block_z)
549  CALL dbcsr_get_block_p(matrix_l(1, img)%matrix, irow, icol, l_block_x, found)
550  CALL dbcsr_get_block_p(matrix_l(2, img)%matrix, irow, icol, l_block_y, found)
551  CALL dbcsr_get_block_p(matrix_l(3, img)%matrix, irow, icol, l_block_z, found)
552  END IF
553 
554  IF (do_dr) THEN
555  NULLIFY (r_2block, r_3block)
556  CALL dbcsr_get_block_p(matrix_h(2, img)%matrix, irow, icol, r_2block, found)
557  CALL dbcsr_get_block_p(matrix_h(3, img)%matrix, irow, icol, r_3block, found)
558  END IF
559 
560  IF (calculate_forces) THEN
561  NULLIFY (p_block)
562  CALL dbcsr_get_block_p(matrix_p(1, img)%matrix, irow, icol, p_block, found)
563  END IF
564 
565  ! loop over all kinds for projector atom
566  IF (ASSOCIATED(h_block)) THEN
567 !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
568 !$ hash = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
569 
570  DO kkind = 1, nkind
571  iac = ikind + nkind*(kkind - 1)
572  ibc = jkind + nkind*(kkind - 1)
573  IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
574  IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
575  CALL get_alist(sap_int(iac), alist_ac, iatom)
576  CALL get_alist(sap_int(ibc), alist_bc, jatom)
577  IF (.NOT. ASSOCIATED(alist_ac)) cycle
578  IF (.NOT. ASSOCIATED(alist_bc)) cycle
579  DO kac = 1, alist_ac%nclist
580  DO kbc = 1, alist_bc%nclist
581  IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
582  IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
583  IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
584  acint => alist_ac%clist(kac)%acint
585  bcint => alist_bc%clist(kbc)%acint
586  achint => alist_ac%clist(kac)%achint
587  bchint => alist_bc%clist(kbc)%achint
588  IF (do_soc) THEN
589  alkint => alist_ac%clist(kac)%alkint
590  blkint => alist_bc%clist(kbc)%alkint
591  END IF
592  na = SIZE(acint, 1)
593  np = SIZE(acint, 2)
594  nb = SIZE(bcint, 1)
595 !$ CALL omp_set_lock(locks(hash))
596  IF (.NOT. do_dr) THEN
597  IF (iatom <= jatom) THEN
598  h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
599  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 1)))
600  ELSE
601  h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
602  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 1)))
603  END IF
604  END IF
605  IF (do_soc) THEN
606  IF (iatom <= jatom) THEN
607  l_block_x(1:na, 1:nb) = l_block_x(1:na, 1:nb) + &
608  matmul(alkint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 1)))
609  l_block_y(1:na, 1:nb) = l_block_y(1:na, 1:nb) + &
610  matmul(alkint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1)))
611  l_block_z(1:na, 1:nb) = l_block_z(1:na, 1:nb) + &
612  matmul(alkint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1)))
613 
614  ELSE
615  l_block_x(1:nb, 1:na) = l_block_x(1:nb, 1:na) - &
616  matmul(blkint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 1)))
617  l_block_y(1:nb, 1:na) = l_block_y(1:nb, 1:na) - &
618  matmul(blkint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))
619  l_block_z(1:nb, 1:na) = l_block_z(1:nb, 1:na) - &
620  matmul(blkint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))
621  END IF
622  END IF
623 !$ CALL omp_unset_lock(locks(hash))
624  IF (calculate_forces) THEN
625  IF (ASSOCIATED(p_block)) THEN
626  katom = alist_ac%clist(kac)%catom
627  DO i = 1, 3
628  j = i + 1
629  IF (iatom <= jatom) THEN
630  fa(i) = sum(p_block(1:na, 1:nb)* &
631  matmul(acint(1:na, 1:np, j), transpose(bchint(1:nb, 1:np, 1))))
632  fb(i) = sum(p_block(1:na, 1:nb)* &
633  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, j))))
634  ELSE
635  fa(i) = sum(p_block(1:nb, 1:na)* &
636  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, j))))
637  fb(i) = sum(p_block(1:nb, 1:na)* &
638  matmul(bcint(1:nb, 1:np, j), transpose(achint(1:na, 1:np, 1))))
639  END IF
640  force_thread(i, iatom) = force_thread(i, iatom) + f0*fa(i)
641  force_thread(i, katom) = force_thread(i, katom) - f0*fa(i)
642  force_thread(i, jatom) = force_thread(i, jatom) + f0*fb(i)
643  force_thread(i, katom) = force_thread(i, katom) - f0*fb(i)
644  END DO
645 
646  IF (use_virial) THEN
647  rac = alist_ac%clist(kac)%rac
648  rbc = alist_bc%clist(kbc)%rac
649  CALL virial_pair_force(pv_thread, f0, fa, rac)
650  CALL virial_pair_force(pv_thread, f0, fb, rbc)
651  END IF
652  END IF
653  END IF
654 
655  IF (do_dr) THEN
656  i = 1; j = 2;
657  katom = alist_ac%clist(kac)%catom
658  IF (iatom <= jatom) THEN
659  h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
660  (deltar(i, iatom) - deltar(i, katom))* &
661  matmul(acint(1:na, 1:np, j), transpose(bchint(1:nb, 1:np, 1)))
662 
663  h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
664  (deltar(i, jatom) - deltar(i, katom))* &
665  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, j)))
666  ELSE
667  h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
668  (deltar(i, iatom) - deltar(i, katom))* &
669  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, j)))
670  h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
671  (deltar(i, jatom) - deltar(i, katom))* &
672  matmul(bcint(1:nb, 1:np, j), transpose(achint(1:na, 1:np, 1)))
673  END IF
674 
675  i = 2; j = 3;
676  katom = alist_ac%clist(kac)%catom
677  IF (iatom <= jatom) THEN
678  r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
679  (deltar(i, iatom) - deltar(i, katom))* &
680  matmul(acint(1:na, 1:np, j), transpose(bchint(1:nb, 1:np, 1)))
681 
682  r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
683  (deltar(i, jatom) - deltar(i, katom))* &
684  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, j)))
685  ELSE
686  r_2block(1:nb, 1:na) = r_2block(1:nb, 1:na) + &
687  (deltar(i, iatom) - deltar(i, katom))* &
688  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, j)))
689  r_2block(1:nb, 1:na) = r_2block(1:nb, 1:na) + &
690  (deltar(i, jatom) - deltar(i, katom))* &
691  matmul(bcint(1:nb, 1:np, j), transpose(achint(1:na, 1:np, 1)))
692  END IF
693 
694  i = 3; j = 4;
695  katom = alist_ac%clist(kac)%catom
696  IF (iatom <= jatom) THEN
697  r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &
698  (deltar(i, iatom) - deltar(i, katom))* &
699  matmul(acint(1:na, 1:np, j), transpose(bchint(1:nb, 1:np, 1)))
700 
701  r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &
702  (deltar(i, jatom) - deltar(i, katom))* &
703  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, j)))
704  ELSE
705  r_3block(1:nb, 1:na) = r_3block(1:nb, 1:na) + &
706  (deltar(i, iatom) - deltar(i, katom))* &
707  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, j)))
708  r_3block(1:nb, 1:na) = r_3block(1:nb, 1:na) + &
709  (deltar(i, jatom) - deltar(i, katom))* &
710  matmul(bcint(1:nb, 1:np, j), transpose(achint(1:na, 1:np, 1)))
711  END IF
712 
713  END IF
714  EXIT ! We have found a match and there can be only one single match
715  END IF
716  END DO
717  END DO
718  END DO
719  END IF
720  END DO
721 
722 !$OMP DO
723 !$ DO lock_num = 1, nlock
724 !$ call omp_destroy_lock(locks(lock_num))
725 !$ END DO
726 !$OMP END DO
727 
728 !$OMP SINGLE
729 !$ DEALLOCATE (locks)
730 !$OMP END SINGLE NOWAIT
731 
732 !$OMP END PARALLEL
733 
734  CALL release_sap_int(sap_int)
735 
736  DEALLOCATE (basis_set, gpotential, spotential)
737  IF (calculate_forces) THEN
738  CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
739 !$OMP DO
740  DO iatom = 1, natom
741  atom_a = atom_of_kind(iatom)
742  ikind = kind_of(iatom)
743  force(ikind)%gth_ppnl(:, atom_a) = force(ikind)%gth_ppnl(:, atom_a) + force_thread(:, iatom)
744  END DO
745 !$OMP END DO
746  DEALLOCATE (atom_of_kind, kind_of)
747  END IF
748 
749  IF (calculate_forces .AND. use_virial) THEN
750  virial%pv_ppnl = virial%pv_ppnl + pv_thread
751  virial%pv_virial = virial%pv_virial + pv_thread
752  END IF
753 
754  IF (calculate_forces) THEN
755  ! If LSD, then recover alpha density and beta density
756  ! from the total density (1) and the spin density (2)
757  IF (SIZE(matrix_p, 1) == 2) THEN
758  DO img = 1, nimages
759  CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
760  alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
761  CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
762  alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
763  END DO
764  END IF
765  END IF
766 
767  END IF !ppnl_present
768 
769  CALL timestop(handle)
770 
771  END SUBROUTINE build_core_ppnl
772 
773 END MODULE core_ppnl
ai_angmom
Calculation of the angular momentum integrals over Cartesian Gaussian-type functions.
Definition: ai_angmom.F:17
ai_angmom::angmom
subroutine, public angmom(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, rac, rbc, angab)
...
Definition: ai_angmom.F:52
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap.F:18
ai_overlap::overlap
subroutine, public overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, rab, dab, sab, da_max_set, return_derivatives, s, lds, sdab, pab, force_a)
Purpose: Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions...
Definition: ai_overlap.F:73
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
basis_set_types
Definition: basis_set_types.F:15
core_ppnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition: core_ppnl.F:15
core_ppnl::build_core_ppnl
subroutine, public build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, basis_type, deltaR, matrix_l)
...
Definition: core_ppnl.F:89
external_potential_types
Definition of the atomic potential types.
Definition: external_potential_types.F:14
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition: kinds.F:54
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition: orbital_pointers.F:253
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition: orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_force_types
Definition: qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
sap_kind_types
General overlap type integrals containers.
Definition: sap_kind_types.F:12
sap_kind_types::release_sap_int
subroutine, public release_sap_int(sap_int)
...
Definition: sap_kind_types.F:86
sap_kind_types::sap_sort
subroutine, public sap_sort(sap_int)
...
Definition: sap_kind_types.F:247
sap_kind_types::get_alist
subroutine, public get_alist(sap_int, alist, atom)
...
Definition: sap_kind_types.F:140
virial_methods
Definition: virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition: virial_methods.F:139
virial_types
Definition: virial_types.F:13