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