Line data Source code
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 : ! **************************************************************************************************
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 cp_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 : !> \param atcore ...
86 : ! **************************************************************************************************
87 16998 : 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 16998 : 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 16998 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
121 : INTEGER, DIMENSION(3) :: cell_b, cell_c
122 16998 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
123 16998 : nsgf_seta
124 16998 : 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 16998 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radp
129 16998 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
130 16998 : 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 16998 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
136 : TYPE(gth_potential_type), POINTER :: gth_potential
137 16998 : TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
138 : TYPE(clist_type), POINTER :: clist
139 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
140 33996 : REAL(KIND=dp), DIMENSION(SIZE(particle_set)) :: at_thread
141 16998 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, alkint, bchint, bcint, &
142 16998 : blkint
143 16998 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_block, l_block_x, l_block_y, &
144 16998 : l_block_z, p_block, r_2block, &
145 16998 : r_3block, rpgfa, sphi_a, vprj_ppnl, &
146 16998 : wprj_ppnl, zeta
147 16998 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
148 33996 : REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
149 16998 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
150 16998 : TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
151 : TYPE(sgp_potential_type), POINTER :: sgp_potential
152 :
153 : !$ INTEGER(kind=omp_lock_kind), &
154 16998 : !$ 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 16998 : do_dR = .FALSE.
162 72 : IF (PRESENT(deltaR)) do_dR = .TRUE.
163 16998 : doat = .FALSE.
164 16998 : IF (PRESENT(atcore)) doat = .TRUE.
165 16998 : IF ((calculate_forces .OR. doat) .AND. do_dR) THEN
166 0 : CPABORT("core_ppl: incompatible options")
167 : END IF
168 :
169 16998 : IF (calculate_forces) THEN
170 6911 : CALL timeset(routineN//"_forces", handle)
171 : ELSE
172 10087 : CALL timeset(routineN, handle)
173 : END IF
174 :
175 16998 : do_soc = PRESENT(matrix_l)
176 :
177 16998 : ppnl_present = ASSOCIATED(sap_ppnl)
178 :
179 16998 : IF (ppnl_present) THEN
180 :
181 16998 : nkind = SIZE(atomic_kind_set)
182 16998 : natom = SIZE(particle_set)
183 :
184 16998 : do_kp = (nimages > 1)
185 :
186 16998 : IF (do_kp) THEN
187 578 : IF (PRESENT(cell_to_index)) THEN
188 578 : CPASSERT(ASSOCIATED(cell_to_index))
189 : ELSE
190 0 : CPABORT("Missing cell_to_index for k-point calculation")
191 : END IF
192 : END IF
193 :
194 16998 : IF (calculate_forces .OR. doat) THEN
195 6973 : IF (SIZE(matrix_p, 1) == 2) THEN
196 2114 : DO img = 1, nimages
197 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
198 1358 : 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 2114 : alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
201 : END DO
202 : END IF
203 : END IF
204 :
205 16998 : 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 16998 : basis_type=basis_type)
214 :
215 16998 : maxl = MAX(maxlgto, maxlppnl)
216 16998 : CALL init_orbital_pointers(maxl + nder + 1)
217 :
218 16998 : ldsab = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
219 16998 : ldai = ncoset(maxl + nder + 1)
220 :
221 : ! sap_int needs to be shared as multiple threads need to access this
222 115356 : ALLOCATE (sap_int(nkind*nkind))
223 81360 : DO i = 1, nkind*nkind
224 64362 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
225 81360 : sap_int(i)%nalist = 0
226 : END DO
227 :
228 : ! Set up direct access to basis and potential
229 180120 : ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind))
230 48708 : DO ikind = 1, nkind
231 31710 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
232 31710 : IF (ASSOCIATED(orb_basis_set)) THEN
233 31710 : basis_set(ikind)%gto_basis_set => orb_basis_set
234 : ELSE
235 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
236 : END IF
237 31710 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
238 31710 : NULLIFY (gpotential(ikind)%gth_potential)
239 31710 : NULLIFY (spotential(ikind)%sgp_potential)
240 48708 : IF (ASSOCIATED(gth_potential)) THEN
241 31466 : gpotential(ikind)%gth_potential => gth_potential
242 31466 : IF (do_soc .AND. (.NOT. gth_potential%soc)) THEN
243 0 : CPABORT("Spin-orbit coupling selected, but GTH potential without SOC parameters provided")
244 : END IF
245 244 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
246 10 : spotential(ikind)%sgp_potential => sgp_potential
247 : END IF
248 : END DO
249 :
250 : ! Allocate sap int
251 1132702 : DO slot = 1, sap_ppnl(1)%nl_size
252 :
253 1115704 : ikind = sap_ppnl(1)%nlist_task(slot)%ikind
254 1115704 : kkind = sap_ppnl(1)%nlist_task(slot)%jkind
255 1115704 : iatom = sap_ppnl(1)%nlist_task(slot)%iatom
256 1115704 : katom = sap_ppnl(1)%nlist_task(slot)%jatom
257 1115704 : nlist = sap_ppnl(1)%nlist_task(slot)%nlist
258 1115704 : ilist = sap_ppnl(1)%nlist_task(slot)%ilist
259 1115704 : nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
260 :
261 1115704 : iac = ikind + nkind*(kkind - 1)
262 1115704 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
263 1115704 : IF (.NOT. ASSOCIATED(gpotential(kkind)%gth_potential) .AND. &
264 : .NOT. ASSOCIATED(spotential(kkind)%sgp_potential)) CYCLE
265 1115704 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
266 35292 : sap_int(iac)%a_kind = ikind
267 35292 : sap_int(iac)%p_kind = kkind
268 35292 : sap_int(iac)%nalist = nlist
269 177710 : ALLOCATE (sap_int(iac)%alist(nlist))
270 107126 : DO i = 1, nlist
271 71834 : NULLIFY (sap_int(iac)%alist(i)%clist)
272 71834 : sap_int(iac)%alist(i)%aatom = 0
273 107126 : sap_int(iac)%alist(i)%nclist = 0
274 : END DO
275 : END IF
276 1132702 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
277 71776 : sap_int(iac)%alist(ilist)%aatom = iatom
278 71776 : sap_int(iac)%alist(ilist)%nclist = nneighbor
279 1761688 : ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
280 1187480 : DO i = 1, nneighbor
281 1187480 : 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 16998 : !$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 : ! Screened projector channels must not retain SOC integrals
390 : ! from the preceding shell set or atom/image pair.
391 : IF (do_soc) work_l = 0.0_dp
392 : DO l = 0, lppnl
393 : nprjc = nprj_ppnl(l)*nco(l)
394 : IF (nprjc == 0) CYCLE
395 : rprjc(1) = ppnl_radius
396 : IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE
397 : lc_max = l + 2*(nprj_ppnl(l) - 1)
398 : lc_min = l
399 : zetc(1) = alpha_ppnl(l)
400 : ncoc = ncoset(lc_max)
401 :
402 : ! Calculate the primitive overlap integrals
403 : CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
404 : lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, nder, .TRUE., ai_work, ldai)
405 : ! Transformation step projector functions (Cartesian -> spherical)
406 : na = ncoa
407 : nb = nprjc
408 : np = ncoc
409 : DO i = 1, maxder
410 : first_col = (i - 1)*ldsab
411 : ! CALL dgemm("N", "N", ncoa, nprjc, ncoc, 1.0_dp, sab(1, first_col + 1), SIZE(sab, 1), &
412 : ! cprj(1, prjc), SIZE(cprj, 1), 0.0_dp, work(1, first_col + prjc), ldsab)
413 : work(1:na, first_col + prjc:first_col + prjc + nb - 1) = &
414 : MATMUL(sab(1:na, first_col + 1:first_col + np), cprj(1:np, prjc:prjc + nb - 1))
415 : END DO
416 :
417 : IF (do_soc) THEN
418 : ! Calculate the primitive angular momentum integrals needed for spin-orbit coupling
419 : lab = 0.0_dp
420 : CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
421 : lc_max, 1, zetc, rprjc, -rac, [0._dp, 0._dp, 0._dp], lab)
422 : DO i_dim = 1, 3
423 : work_l(1:na, prjc:prjc + nb - 1, i_dim) = &
424 : MATMUL(lab(1:na, 1:np, i_dim), cprj(1:np, prjc:prjc + nb - 1))
425 : END DO
426 : END IF
427 :
428 : prjc = prjc + nprjc
429 :
430 : END DO
431 : na = nsgf_seta(iset)
432 : nb = nppnl
433 : np = ncoa
434 : DO i = 1, maxder
435 : first_col = (i - 1)*ldsab + 1
436 : ! Contraction step (basis functions)
437 : ! CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
438 : ! work(1, first_col), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
439 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
440 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, first_col:first_col + nb - 1))
441 : ! Multiply with interaction matrix(h)
442 : ! CALL dgemm("N", "N", nsgf_seta(iset), nppnl, nppnl, 1.0_dp, clist%acint(sgfa, 1, i), nsgfa, &
443 : ! vprj_ppnl(1, 1), SIZE(vprj_ppnl, 1), 0.0_dp, clist%achint(sgfa, 1, i), nsgfa)
444 : clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
445 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), vprj_ppnl(1:nb, 1:nb))
446 : END DO
447 : IF (do_soc) THEN
448 : DO i_dim = 1, 3
449 : clist%alint(sgfa:sgfa + na - 1, 1:nb, i_dim) = &
450 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_l(1:np, 1:nb, i_dim))
451 : clist%alkint(sgfa:sgfa + na - 1, 1:nb, i_dim) = &
452 : MATMUL(clist%alint(sgfa:sgfa + na - 1, 1:nb, i_dim), wprj_ppnl(1:nb, 1:nb))
453 : END DO
454 : END IF
455 : ELSE
456 : ! SGP potential
457 : ! Calculate the primitive overlap integrals
458 : CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
459 : lppnl, 0, nnl, radp, a_nl, rac, dac, sab, nder, .TRUE., ai_work, ldai)
460 : na = nsgf_seta(iset)
461 : nb = nppnl
462 : np = ncoa
463 : DO i = 1, maxder
464 : first_col = (i - 1)*ldsab + 1
465 : ! Transformation step projector functions (cartesian->spherical)
466 : ! CALL dgemm("N", "N", ncoa, nppnl, nprjc, 1.0_dp, sab(1, first_col), ldsab, &
467 : ! cprj(1, 1), SIZE(cprj, 1), 0.0_dp, work(1, 1), ldsab)
468 : work(1:np, 1:nb) = MATMUL(sab(1:np, first_col:first_col + nprjc - 1), cprj(1:nprjc, 1:nb))
469 : ! Contraction step (basis functions)
470 : ! CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
471 : ! work(1, 1), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
472 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
473 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, 1:nb))
474 : ! *** Multiply with interaction matrix(h) ***
475 : ncoc = sgfa + nsgf_seta(iset) - 1
476 : DO j = 1, nppnl
477 : clist%achint(sgfa:ncoc, j, i) = clist%acint(sgfa:ncoc, j, i)*hprj(j)
478 : END DO
479 : END DO
480 : END IF
481 : END DO
482 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
483 : clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
484 : IF (.NOT. do_gth) DEALLOCATE (radp)
485 : END DO
486 :
487 : DEALLOCATE (sab, ai_work, work)
488 : IF (do_soc) DEALLOCATE (lab, work_l)
489 : !$OMP END PARALLEL
490 :
491 : ! Set up a sorting index
492 16998 : CALL sap_sort(sap_int)
493 : ! All integrals needed have been calculated and stored in sap_int
494 : ! We now calculate the Hamiltonian matrix elements
495 :
496 273486 : force_thread = 0.0_dp
497 81120 : at_thread = 0.0_dp
498 16998 : pv_thread = 0.0_dp
499 :
500 : !$OMP PARALLEL &
501 : !$OMP DEFAULT (NONE) &
502 : !$OMP SHARED (do_kp, basis_set, matrix_h, matrix_l, cell_to_index,&
503 : !$OMP sab_orb, matrix_p, sap_int, nkind, eps_ppnl, force, &
504 : !$OMP doat, do_dR, deltaR, maxder, nder, &
505 : !$OMP locks, virial, use_virial, calculate_forces, do_soc, natom) &
506 : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
507 : !$OMP slot, iab, atom_a, f0, irow, icol, h_block, &
508 : !$OMP l_block_x, l_block_y, l_block_z, lock_num, &
509 : !$OMP r_2block, r_3block, atk, &
510 : !$OMP found,p_block, iac, ibc, alist_ac, alist_bc, acint, bcint, &
511 : !$OMP achint, bchint, alkint, blkint, &
512 : !$OMP na, np, nb, katom, j, fa, fb, rbc, rac, &
513 : !$OMP kkind, kac, kbc, i, img, hash, iatom8) &
514 16998 : !$OMP REDUCTION (+ : at_thread, pv_thread, force_thread )
515 :
516 : !$OMP SINGLE
517 : !$ ALLOCATE (locks(nlock))
518 : !$OMP END SINGLE
519 :
520 : !$OMP DO
521 : !$ DO lock_num = 1, nlock
522 : !$ call omp_init_lock(locks(lock_num))
523 : !$ END DO
524 : !$OMP END DO
525 :
526 : !$OMP DO SCHEDULE(GUIDED)
527 : DO slot = 1, sab_orb(1)%nl_size
528 :
529 : ikind = sab_orb(1)%nlist_task(slot)%ikind
530 : jkind = sab_orb(1)%nlist_task(slot)%jkind
531 : iatom = sab_orb(1)%nlist_task(slot)%iatom
532 : jatom = sab_orb(1)%nlist_task(slot)%jatom
533 : cell_b(:) = sab_orb(1)%nlist_task(slot)%cell(:)
534 : rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
535 :
536 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
537 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
538 :
539 : iab = ikind + nkind*(jkind - 1)
540 :
541 : ! Use the symmetry of the first derivatives
542 : IF (iatom == jatom) THEN
543 : f0 = 1.0_dp
544 : ELSE
545 : f0 = 2.0_dp
546 : END IF
547 :
548 : IF (do_kp) THEN
549 : img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
550 : ELSE
551 : img = 1
552 : END IF
553 :
554 : ! Create matrix blocks for a new matrix block column
555 : IF (iatom <= jatom) THEN
556 : irow = iatom
557 : icol = jatom
558 : ELSE
559 : irow = jatom
560 : icol = iatom
561 : END IF
562 : NULLIFY (h_block)
563 : CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found)
564 : IF (do_soc) THEN
565 : NULLIFY (l_block_x, l_block_y, l_block_z)
566 : CALL dbcsr_get_block_p(matrix_l(1, img)%matrix, irow, icol, l_block_x, found)
567 : CALL dbcsr_get_block_p(matrix_l(2, img)%matrix, irow, icol, l_block_y, found)
568 : CALL dbcsr_get_block_p(matrix_l(3, img)%matrix, irow, icol, l_block_z, found)
569 : END IF
570 :
571 : IF (do_dR) THEN
572 : NULLIFY (r_2block, r_3block)
573 : CALL dbcsr_get_block_p(matrix_h(2, img)%matrix, irow, icol, r_2block, found)
574 : CALL dbcsr_get_block_p(matrix_h(3, img)%matrix, irow, icol, r_3block, found)
575 : END IF
576 :
577 : IF (calculate_forces .OR. doat) THEN
578 : NULLIFY (p_block)
579 : CALL dbcsr_get_block_p(matrix_p(1, img)%matrix, irow, icol, p_block, found)
580 : END IF
581 :
582 : ! loop over all kinds for projector atom
583 : IF (ASSOCIATED(h_block)) THEN
584 : !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
585 : !$ hash = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
586 :
587 : DO kkind = 1, nkind
588 : iac = ikind + nkind*(kkind - 1)
589 : ibc = jkind + nkind*(kkind - 1)
590 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
591 : IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
592 : CALL get_alist(sap_int(iac), alist_ac, iatom)
593 : CALL get_alist(sap_int(ibc), alist_bc, jatom)
594 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
595 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
596 : DO kac = 1, alist_ac%nclist
597 : DO kbc = 1, alist_bc%nclist
598 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
599 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
600 : IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
601 : acint => alist_ac%clist(kac)%acint
602 : bcint => alist_bc%clist(kbc)%acint
603 : achint => alist_ac%clist(kac)%achint
604 : bchint => alist_bc%clist(kbc)%achint
605 : IF (do_soc) THEN
606 : alkint => alist_ac%clist(kac)%alkint
607 : blkint => alist_bc%clist(kbc)%alkint
608 : END IF
609 : na = SIZE(acint, 1)
610 : np = SIZE(acint, 2)
611 : nb = SIZE(bcint, 1)
612 : !$ CALL omp_set_lock(locks(hash))
613 : IF (.NOT. do_dR) THEN
614 : IF (iatom <= jatom) THEN
615 : h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
616 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1)))
617 : ELSE
618 : h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
619 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1)))
620 : END IF
621 : END IF
622 : IF (do_soc) THEN
623 : IF (iatom <= jatom) THEN
624 : l_block_x(1:na, 1:nb) = l_block_x(1:na, 1:nb) + &
625 : MATMUL(alkint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1)))
626 : l_block_y(1:na, 1:nb) = l_block_y(1:na, 1:nb) + &
627 : MATMUL(alkint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1)))
628 : l_block_z(1:na, 1:nb) = l_block_z(1:na, 1:nb) + &
629 : MATMUL(alkint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 1)))
630 :
631 : ELSE
632 : l_block_x(1:nb, 1:na) = l_block_x(1:nb, 1:na) + &
633 : MATMUL(blkint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1)))
634 : l_block_y(1:nb, 1:na) = l_block_y(1:nb, 1:na) + &
635 : MATMUL(blkint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 1)))
636 : l_block_z(1:nb, 1:na) = l_block_z(1:nb, 1:na) + &
637 : MATMUL(blkint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 1)))
638 : END IF
639 : END IF
640 : !$ CALL omp_unset_lock(locks(hash))
641 : IF (calculate_forces) THEN
642 : IF (ASSOCIATED(p_block)) THEN
643 : katom = alist_ac%clist(kac)%catom
644 : DO i = 1, 3
645 : j = i + 1
646 : IF (iatom <= jatom) THEN
647 : fa(i) = SUM(p_block(1:na, 1:nb)* &
648 : MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1))))
649 : fb(i) = SUM(p_block(1:na, 1:nb)* &
650 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j))))
651 : ELSE
652 : fa(i) = SUM(p_block(1:nb, 1:na)* &
653 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j))))
654 : fb(i) = SUM(p_block(1:nb, 1:na)* &
655 : MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1))))
656 : END IF
657 : force_thread(i, iatom) = force_thread(i, iatom) + f0*fa(i)
658 : force_thread(i, katom) = force_thread(i, katom) - f0*fa(i)
659 : force_thread(i, jatom) = force_thread(i, jatom) + f0*fb(i)
660 : force_thread(i, katom) = force_thread(i, katom) - f0*fb(i)
661 : END DO
662 :
663 : IF (use_virial) THEN
664 : rac = alist_ac%clist(kac)%rac
665 : rbc = alist_bc%clist(kbc)%rac
666 : CALL virial_pair_force(pv_thread, f0, fa, rac)
667 : CALL virial_pair_force(pv_thread, f0, fb, rbc)
668 : END IF
669 : END IF
670 : END IF
671 :
672 : IF (do_dR) THEN
673 : i = 1; j = 2
674 : katom = alist_ac%clist(kac)%catom
675 : IF (iatom <= jatom) THEN
676 : h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
677 : (deltaR(i, iatom) - deltaR(i, katom))* &
678 : MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1)))
679 :
680 : h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
681 : (deltaR(i, jatom) - deltaR(i, katom))* &
682 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j)))
683 : ELSE
684 : h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
685 : (deltaR(i, iatom) - deltaR(i, katom))* &
686 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j)))
687 : h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
688 : (deltaR(i, jatom) - deltaR(i, katom))* &
689 : MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
690 : END IF
691 :
692 : i = 2; j = 3
693 : katom = alist_ac%clist(kac)%catom
694 : IF (iatom <= jatom) THEN
695 : r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
696 : (deltaR(i, iatom) - deltaR(i, katom))* &
697 : MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1)))
698 :
699 : r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
700 : (deltaR(i, jatom) - deltaR(i, katom))* &
701 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j)))
702 : ELSE
703 : r_2block(1:nb, 1:na) = r_2block(1:nb, 1:na) + &
704 : (deltaR(i, iatom) - deltaR(i, katom))* &
705 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j)))
706 : r_2block(1:nb, 1:na) = r_2block(1:nb, 1:na) + &
707 : (deltaR(i, jatom) - deltaR(i, katom))* &
708 : MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
709 : END IF
710 :
711 : i = 3; j = 4
712 : katom = alist_ac%clist(kac)%catom
713 : IF (iatom <= jatom) THEN
714 : r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &
715 : (deltaR(i, iatom) - deltaR(i, katom))* &
716 : MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1)))
717 :
718 : r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &
719 : (deltaR(i, jatom) - deltaR(i, katom))* &
720 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j)))
721 : ELSE
722 : r_3block(1:nb, 1:na) = r_3block(1:nb, 1:na) + &
723 : (deltaR(i, iatom) - deltaR(i, katom))* &
724 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j)))
725 : r_3block(1:nb, 1:na) = r_3block(1:nb, 1:na) + &
726 : (deltaR(i, jatom) - deltaR(i, katom))* &
727 : MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
728 : END IF
729 :
730 : END IF
731 : IF (doat) THEN
732 : IF (ASSOCIATED(p_block)) THEN
733 : katom = alist_ac%clist(kac)%catom
734 : IF (iatom <= jatom) THEN
735 : atk = SUM(p_block(1:na, 1:nb)* &
736 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1))))
737 : ELSE
738 : atk = SUM(p_block(1:nb, 1:na)* &
739 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1))))
740 : END IF
741 : at_thread(katom) = at_thread(katom) + f0*atk
742 : END IF
743 : END IF
744 : EXIT ! We have found a match and there can be only one single match
745 : END IF
746 : END DO
747 : END DO
748 : END DO
749 : END IF
750 : END DO
751 :
752 : !$OMP DO
753 : !$ DO lock_num = 1, nlock
754 : !$ call omp_destroy_lock(locks(lock_num))
755 : !$ END DO
756 : !$OMP END DO
757 :
758 : !$OMP SINGLE
759 : !$ DEALLOCATE (locks)
760 : !$OMP END SINGLE NOWAIT
761 :
762 : !$OMP END PARALLEL
763 :
764 16998 : CALL release_sap_int(sap_int)
765 :
766 16998 : DEALLOCATE (basis_set, gpotential, spotential)
767 16998 : IF (calculate_forces) THEN
768 6911 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
769 : !$OMP DO
770 : DO iatom = 1, natom
771 25021 : atom_a = atom_of_kind(iatom)
772 25021 : ikind = kind_of(iatom)
773 100084 : force(ikind)%gth_ppnl(:, atom_a) = force(ikind)%gth_ppnl(:, atom_a) + force_thread(:, iatom)
774 : END DO
775 : !$OMP END DO
776 6911 : DEALLOCATE (atom_of_kind, kind_of)
777 : END IF
778 :
779 16998 : IF (calculate_forces .AND. use_virial) THEN
780 11336 : virial%pv_ppnl = virial%pv_ppnl + pv_thread
781 11336 : virial%pv_virial = virial%pv_virial + pv_thread
782 : END IF
783 :
784 16998 : IF (doat) THEN
785 280 : atcore(1:natom) = atcore(1:natom) + at_thread
786 : END IF
787 :
788 33996 : IF (calculate_forces .OR. doat) THEN
789 : ! If LSD, then recover alpha density and beta density
790 : ! from the total density (1) and the spin density (2)
791 6973 : IF (SIZE(matrix_p, 1) == 2) THEN
792 2114 : DO img = 1, nimages
793 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
794 1358 : alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
795 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
796 2114 : alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
797 : END DO
798 : END IF
799 : END IF
800 :
801 : END IF !ppnl_present
802 :
803 16998 : CALL timestop(handle)
804 :
805 33996 : END SUBROUTINE build_core_ppnl
806 :
807 : END MODULE core_ppnl
|
