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 upf_projector_integrals, ONLY: upf_projector_init,&
47 : upf_projector_overlap,&
48 : upf_projector_type
49 : USE virial_methods, ONLY: virial_pair_force
50 : USE virial_types, ONLY: virial_type
51 :
52 : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
53 : !$ omp_init_lock, omp_set_lock, &
54 : !$ omp_unset_lock, omp_destroy_lock
55 :
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 :
60 : PRIVATE
61 :
62 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'core_ppnl'
63 :
64 : PUBLIC :: build_core_ppnl
65 :
66 : CONTAINS
67 :
68 : ! **************************************************************************************************
69 : !> \brief ...
70 : !> \param matrix_h ...
71 : !> \param matrix_p ...
72 : !> \param force ...
73 : !> \param virial ...
74 : !> \param calculate_forces ...
75 : !> \param use_virial ...
76 : !> \param nder ...
77 : !> \param qs_kind_set ...
78 : !> \param atomic_kind_set ...
79 : !> \param particle_set ...
80 : !> \param sab_orb ...
81 : !> \param sap_ppnl ...
82 : !> \param eps_ppnl ...
83 : !> \param nimages ...
84 : !> \param cell_to_index ...
85 : !> \param basis_type ...
86 : !> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
87 : !> \param matrix_l ...
88 : !> \param atcore ...
89 : ! **************************************************************************************************
90 17370 : SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
91 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
92 17370 : nimages, cell_to_index, basis_type, deltaR, matrix_l, atcore)
93 :
94 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
95 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
96 : TYPE(virial_type), POINTER :: virial
97 : LOGICAL, INTENT(IN) :: calculate_forces
98 : LOGICAL :: use_virial
99 : INTEGER :: nder
100 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
101 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
102 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
103 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
104 : POINTER :: sab_orb, sap_ppnl
105 : REAL(KIND=dp), INTENT(IN) :: eps_ppnl
106 : INTEGER, INTENT(IN) :: nimages
107 : INTEGER, DIMENSION(:, :, :), OPTIONAL, POINTER :: cell_to_index
108 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
109 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
110 : OPTIONAL :: deltaR
111 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
112 : POINTER :: matrix_l
113 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
114 : OPTIONAL :: atcore
115 :
116 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_ppnl'
117 :
118 : INTEGER :: atom_a, first_col, handle, i, i_dim, iab, iac, iatom, ib, ibc, icol, ikind, &
119 : ilist, img, irow, iset, j, jatom, jb, jkind, jneighbor, kac, katom, kbc, kkind, l, &
120 : lc_max, lc_min, ldai, ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, maxppnl, &
121 : maxsgf, na, natom, nb, ncoa, ncoc, nkind, nlist, nneighbor, nnl, np, nppnl, nprjc, nseta, &
122 : nsgfa, prjc, sgfa, slot
123 17370 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
124 : INTEGER, DIMENSION(3) :: cell_b, cell_c
125 17370 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
126 17370 : nsgf_seta
127 17370 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
128 : LOGICAL :: do_dR, do_gth, do_kp, do_soc, do_upf, &
129 : doat, found, ppnl_present
130 : REAL(KIND=dp) :: atk, dac, f0, ppnl_radius
131 17370 : REAL(KIND=dp), ALLOCATABLE :: upf_int(:, :, :)
132 17370 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radp
133 17370 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
134 17370 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work, lab, work_l
135 : REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc
136 : REAL(KIND=dp), DIMENSION(3) :: fa, fb, rab, rac, rbc
137 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_thread
138 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
139 17370 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
140 : TYPE(gth_potential_type), POINTER :: gth_potential
141 17370 : TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
142 : TYPE(clist_type), POINTER :: clist
143 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
144 34740 : REAL(KIND=dp), DIMENSION(SIZE(particle_set)) :: at_thread
145 17370 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, alkint, bchint, bcint, &
146 17370 : blkint
147 17370 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_block, l_block_x, l_block_y, &
148 17370 : l_block_z, p_block, r_2block, &
149 17370 : r_3block, rpgfa, sphi_a, vprj_ppnl, &
150 17370 : wprj_ppnl, zeta
151 17370 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
152 34740 : REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
153 17370 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
154 17370 : TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
155 : TYPE(sgp_potential_type), POINTER :: sgp_potential
156 : TYPE(upf_projector_type), POINTER :: upf
157 :
158 : !$ INTEGER(kind=omp_lock_kind), &
159 17370 : !$ ALLOCATABLE, DIMENSION(:) :: locks
160 : !$ INTEGER(KIND=int_8) :: iatom8
161 : !$ INTEGER :: lock_num, hash
162 : !$ INTEGER, PARAMETER :: nlock = 501
163 :
164 : MARK_USED(int_8)
165 :
166 17370 : do_dR = .FALSE.
167 72 : IF (PRESENT(deltaR)) do_dR = .TRUE.
168 17370 : doat = .FALSE.
169 17370 : IF (PRESENT(atcore)) doat = .TRUE.
170 17370 : IF ((calculate_forces .OR. doat) .AND. do_dR) THEN
171 0 : CPABORT("core_ppl: incompatible options")
172 : END IF
173 :
174 17370 : IF (calculate_forces) THEN
175 6863 : CALL timeset(routineN//"_forces", handle)
176 : ELSE
177 10507 : CALL timeset(routineN, handle)
178 : END IF
179 :
180 17370 : do_soc = PRESENT(matrix_l)
181 :
182 17370 : ppnl_present = ASSOCIATED(sap_ppnl)
183 :
184 17370 : IF (ppnl_present) THEN
185 :
186 17370 : nkind = SIZE(atomic_kind_set)
187 17370 : natom = SIZE(particle_set)
188 :
189 17370 : do_kp = (nimages > 1)
190 :
191 17370 : IF (do_kp) THEN
192 1008 : IF (PRESENT(cell_to_index)) THEN
193 1008 : CPASSERT(ASSOCIATED(cell_to_index))
194 : ELSE
195 0 : CPABORT("Missing cell_to_index for k-point calculation")
196 : END IF
197 : END IF
198 :
199 17370 : IF (calculate_forces .OR. doat) THEN
200 6925 : IF (SIZE(matrix_p, 1) == 2) THEN
201 3042 : DO img = 1, nimages
202 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
203 2332 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
204 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
205 3042 : alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
206 : END DO
207 : END IF
208 : END IF
209 :
210 17370 : maxder = ncoset(nder)
211 :
212 : CALL get_qs_kind_set(qs_kind_set, &
213 : maxco=maxco, &
214 : maxlgto=maxlgto, &
215 : maxsgf=maxsgf, &
216 : maxlppnl=maxlppnl, &
217 : maxppnl=maxppnl, &
218 17370 : basis_type=basis_type)
219 :
220 17370 : maxl = MAX(maxlgto, maxlppnl)
221 17370 : CALL init_orbital_pointers(maxl + nder + 1)
222 :
223 17370 : ldsab = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
224 17370 : ldai = ncoset(maxl + nder + 1)
225 :
226 : ! sap_int needs to be shared as multiple threads need to access this
227 118026 : ALLOCATE (sap_int(nkind*nkind))
228 83286 : DO i = 1, nkind*nkind
229 65916 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
230 83286 : sap_int(i)%nalist = 0
231 : END DO
232 :
233 : ! Set up direct access to basis and potential
234 184266 : ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind))
235 49842 : DO ikind = 1, nkind
236 32472 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
237 32472 : IF (ASSOCIATED(orb_basis_set)) THEN
238 32472 : basis_set(ikind)%gto_basis_set => orb_basis_set
239 : ELSE
240 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
241 : END IF
242 32472 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
243 32472 : NULLIFY (gpotential(ikind)%gth_potential)
244 32472 : NULLIFY (spotential(ikind)%sgp_potential)
245 49842 : IF (ASSOCIATED(gth_potential)) THEN
246 32204 : gpotential(ikind)%gth_potential => gth_potential
247 32204 : IF (do_soc .AND. (.NOT. gth_potential%soc)) THEN
248 0 : CPABORT("Spin-orbit coupling selected, but GTH potential without SOC parameters provided")
249 : END IF
250 268 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
251 34 : spotential(ikind)%sgp_potential => sgp_potential
252 34 : IF (ALLOCATED(sgp_potential%upf_data)) THEN
253 34 : IF (sgp_potential%upf_data%number_of_proj > 0) THEN
254 18 : IF (sgp_potential%upf_projectors%degree < maxlgto + nder) THEN
255 18 : CALL upf_projector_init(sgp_potential%upf_data, maxlgto + nder, sgp_potential%upf_projectors)
256 : END IF
257 : END IF
258 : END IF
259 34 : IF (do_soc .AND. ALLOCATED(sgp_potential%upf_data)) THEN
260 0 : CPABORT("Spin-orbit UPF projectors are not supported")
261 : END IF
262 : END IF
263 : END DO
264 :
265 : ! Allocate sap int
266 1155614 : DO slot = 1, sap_ppnl(1)%nl_size
267 :
268 1138244 : ikind = sap_ppnl(1)%nlist_task(slot)%ikind
269 1138244 : kkind = sap_ppnl(1)%nlist_task(slot)%jkind
270 1138244 : iatom = sap_ppnl(1)%nlist_task(slot)%iatom
271 1138244 : katom = sap_ppnl(1)%nlist_task(slot)%jatom
272 1138244 : nlist = sap_ppnl(1)%nlist_task(slot)%nlist
273 1138244 : ilist = sap_ppnl(1)%nlist_task(slot)%ilist
274 1138244 : nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
275 :
276 1138244 : iac = ikind + nkind*(kkind - 1)
277 1138244 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
278 1138244 : IF (.NOT. ASSOCIATED(gpotential(kkind)%gth_potential) .AND. &
279 : .NOT. ASSOCIATED(spotential(kkind)%sgp_potential)) CYCLE
280 1138244 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
281 36080 : sap_int(iac)%a_kind = ikind
282 36080 : sap_int(iac)%p_kind = kkind
283 36080 : sap_int(iac)%nalist = nlist
284 181000 : ALLOCATE (sap_int(iac)%alist(nlist))
285 108840 : DO i = 1, nlist
286 72760 : NULLIFY (sap_int(iac)%alist(i)%clist)
287 72760 : sap_int(iac)%alist(i)%aatom = 0
288 108840 : sap_int(iac)%alist(i)%nclist = 0
289 : END DO
290 : END IF
291 1155614 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
292 72702 : sap_int(iac)%alist(ilist)%aatom = iatom
293 72702 : sap_int(iac)%alist(ilist)%nclist = nneighbor
294 1792562 : ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
295 1210946 : DO i = 1, nneighbor
296 1210946 : sap_int(iac)%alist(ilist)%clist(i)%catom = 0
297 : END DO
298 : END IF
299 : END DO
300 :
301 : ! Calculate the overlap integrals <a|p>
302 : !$OMP PARALLEL &
303 : !$OMP DEFAULT (NONE) &
304 : !$OMP SHARED (basis_set, gpotential, spotential, maxder, ncoset, &
305 : !$OMP sap_ppnl, sap_int, nkind, ldsab, ldai, nder, nco, do_soc ) &
306 : !$OMP PRIVATE (ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
307 : !$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, &
308 : !$OMP slot, sphi_a, zeta, cprj, hprj, lppnl, nppnl, nprj_ppnl, &
309 : !$OMP do_upf, upf, upf_int, &
310 : !$OMP clist, iset, ncoa, sgfa, prjc, work, work_l, sab, lab, ai_work, nprjc, &
311 : !$OMP ppnl_radius, ncoc, rpgfa, first_col, vprj_ppnl, wprj_ppnl, i, j, l, do_gth, &
312 : !$OMP set_radius_a, rprjc, dac, lc_max, lc_min, zetc, alpha_ppnl, &
313 17370 : !$OMP na, nb, np, nnl, a_nl, radp, i_dim, ib, jb)
314 :
315 : ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder))
316 : sab = 0.0_dp
317 : ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1)))
318 : ai_work = 0.0_dp
319 : IF (do_soc) THEN
320 : ALLOCATE (lab(ldsab, ldsab, 3), work_l(ldsab, ldsab, 3))
321 : lab = 0.0_dp
322 : END IF
323 :
324 : !$OMP DO SCHEDULE(GUIDED)
325 : DO slot = 1, sap_ppnl(1)%nl_size
326 :
327 : ikind = sap_ppnl(1)%nlist_task(slot)%ikind
328 : kkind = sap_ppnl(1)%nlist_task(slot)%jkind
329 : iatom = sap_ppnl(1)%nlist_task(slot)%iatom
330 : katom = sap_ppnl(1)%nlist_task(slot)%jatom
331 : nlist = sap_ppnl(1)%nlist_task(slot)%nlist
332 : ilist = sap_ppnl(1)%nlist_task(slot)%ilist
333 : nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
334 : jneighbor = sap_ppnl(1)%nlist_task(slot)%inode
335 : cell_c(:) = sap_ppnl(1)%nlist_task(slot)%cell(:)
336 : rac(1:3) = sap_ppnl(1)%nlist_task(slot)%r(1:3)
337 :
338 : iac = ikind + nkind*(kkind - 1)
339 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
340 : ! Get definition of basis set
341 : first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf
342 : la_max => basis_set(ikind)%gto_basis_set%lmax
343 : la_min => basis_set(ikind)%gto_basis_set%lmin
344 : npgfa => basis_set(ikind)%gto_basis_set%npgf
345 : nseta = basis_set(ikind)%gto_basis_set%nset
346 : nsgfa = basis_set(ikind)%gto_basis_set%nsgf
347 : nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set
348 : rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius
349 : set_radius_a => basis_set(ikind)%gto_basis_set%set_radius
350 : sphi_a => basis_set(ikind)%gto_basis_set%sphi
351 : zeta => basis_set(ikind)%gto_basis_set%zet
352 : ! Get definition of PP projectors
353 : do_upf = .FALSE.
354 : IF (ASSOCIATED(gpotential(kkind)%gth_potential)) THEN
355 : ! GTH potential
356 : do_gth = .TRUE.
357 : alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl
358 : cprj => gpotential(kkind)%gth_potential%cprj
359 : lppnl = gpotential(kkind)%gth_potential%lppnl
360 : nppnl = gpotential(kkind)%gth_potential%nppnl
361 : nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
362 : ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
363 : vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
364 : wprj_ppnl => gpotential(kkind)%gth_potential%wprj_ppnl
365 : ELSE IF (ASSOCIATED(spotential(kkind)%sgp_potential)) THEN
366 : ! SGP potential
367 : do_gth = .FALSE.
368 : nprjc = spotential(kkind)%sgp_potential%nppnl
369 : IF (nprjc == 0) CYCLE
370 : do_upf = ALLOCATED(spotential(kkind)%sgp_potential%upf_data)
371 : IF (do_upf) THEN
372 : upf => spotential(kkind)%sgp_potential%upf_projectors
373 : nppnl = upf%nproj
374 : ELSE
375 : nnl = spotential(kkind)%sgp_potential%n_nonlocal
376 : lppnl = spotential(kkind)%sgp_potential%lmax
377 : a_nl => spotential(kkind)%sgp_potential%a_nonlocal
378 : ppnl_radius = spotential(kkind)%sgp_potential%ppnl_radius
379 : ALLOCATE (radp(nnl))
380 : radp(:) = ppnl_radius
381 : cprj => spotential(kkind)%sgp_potential%cprj_ppnl
382 : hprj => spotential(kkind)%sgp_potential%vprj_ppnl
383 : nppnl = SIZE(cprj, 2)
384 : END IF
385 : ELSE
386 : CYCLE
387 : END IF
388 :
389 : dac = SQRT(SUM(rac*rac))
390 : clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
391 : clist%catom = katom
392 : clist%cell = cell_c
393 : clist%rac = rac
394 : ALLOCATE (clist%acint(nsgfa, nppnl, maxder), &
395 : clist%achint(nsgfa, nppnl, maxder), &
396 : clist%alint(nsgfa, nppnl, 3), &
397 : clist%alkint(nsgfa, nppnl, 3))
398 : clist%acint = 0.0_dp
399 : clist%achint = 0.0_dp
400 : clist%alint = 0.0_dp
401 : clist%alkint = 0.0_dp
402 :
403 : clist%nsgf_cnt = 0
404 : NULLIFY (clist%sgf_list)
405 : DO iset = 1, nseta
406 : ncoa = npgfa(iset)*ncoset(la_max(iset))
407 : sgfa = first_sgfa(1, iset)
408 : IF (do_gth) THEN
409 : ! GTH potential
410 : prjc = 1
411 : work = 0.0_dp
412 : ! Screened projector channels must not retain SOC integrals
413 : ! from the preceding shell set or atom/image pair.
414 : IF (do_soc) work_l = 0.0_dp
415 : DO l = 0, lppnl
416 : nprjc = nprj_ppnl(l)*nco(l)
417 : IF (nprjc == 0) CYCLE
418 : rprjc(1) = ppnl_radius
419 : IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE
420 : lc_max = l + 2*(nprj_ppnl(l) - 1)
421 : lc_min = l
422 : zetc(1) = alpha_ppnl(l)
423 : ncoc = ncoset(lc_max)
424 :
425 : ! Calculate the primitive overlap integrals
426 : CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
427 : lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, nder, .TRUE., ai_work, ldai)
428 : ! Transformation step projector functions (Cartesian -> spherical)
429 : na = ncoa
430 : nb = nprjc
431 : np = ncoc
432 : DO i = 1, maxder
433 : first_col = (i - 1)*ldsab
434 : ! CALL dgemm("N", "N", ncoa, nprjc, ncoc, 1.0_dp, sab(1, first_col + 1), SIZE(sab, 1), &
435 : ! cprj(1, prjc), SIZE(cprj, 1), 0.0_dp, work(1, first_col + prjc), ldsab)
436 : work(1:na, first_col + prjc:first_col + prjc + nb - 1) = &
437 : MATMUL(sab(1:na, first_col + 1:first_col + np), cprj(1:np, prjc:prjc + nb - 1))
438 : END DO
439 :
440 : IF (do_soc) THEN
441 : ! Calculate the primitive angular momentum integrals needed for spin-orbit coupling
442 : lab = 0.0_dp
443 : CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
444 : lc_max, 1, zetc, rprjc, -rac, [0._dp, 0._dp, 0._dp], lab)
445 : DO i_dim = 1, 3
446 : work_l(1:na, prjc:prjc + nb - 1, i_dim) = &
447 : MATMUL(lab(1:na, 1:np, i_dim), cprj(1:np, prjc:prjc + nb - 1))
448 : END DO
449 : END IF
450 :
451 : prjc = prjc + nprjc
452 :
453 : END DO
454 : na = nsgf_seta(iset)
455 : nb = nppnl
456 : np = ncoa
457 : DO i = 1, maxder
458 : first_col = (i - 1)*ldsab + 1
459 : ! Contraction step (basis functions)
460 : ! CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
461 : ! work(1, first_col), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
462 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
463 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, first_col:first_col + nb - 1))
464 : ! Multiply with interaction matrix(h)
465 : ! CALL dgemm("N", "N", nsgf_seta(iset), nppnl, nppnl, 1.0_dp, clist%acint(sgfa, 1, i), nsgfa, &
466 : ! vprj_ppnl(1, 1), SIZE(vprj_ppnl, 1), 0.0_dp, clist%achint(sgfa, 1, i), nsgfa)
467 : clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
468 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), vprj_ppnl(1:nb, 1:nb))
469 : END DO
470 : IF (do_soc) THEN
471 : DO i_dim = 1, 3
472 : clist%alint(sgfa:sgfa + na - 1, 1:nb, i_dim) = &
473 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_l(1:np, 1:nb, i_dim))
474 : clist%alkint(sgfa:sgfa + na - 1, 1:nb, i_dim) = &
475 : MATMUL(clist%alint(sgfa:sgfa + na - 1, 1:nb, i_dim), wprj_ppnl(1:nb, 1:nb))
476 : END DO
477 : END IF
478 : ELSE IF (do_upf) THEN
479 : na = nsgf_seta(iset)
480 : nb = nppnl
481 : np = ncoa
482 : ALLOCATE (upf_int(np, nb, maxder))
483 : CALL upf_projector_overlap(upf, la_max(iset), zeta(:npgfa(iset), iset), &
484 : rac, nder, upf_int)
485 : DO i = 1, maxder
486 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
487 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), upf_int(:, :, i))
488 : clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
489 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), upf%d)
490 : END DO
491 : DEALLOCATE (upf_int)
492 : ELSE
493 : ! SGP potential
494 : ! Calculate the primitive overlap integrals
495 : CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
496 : lppnl, 0, nnl, radp, a_nl, rac, dac, sab, nder, .TRUE., ai_work, ldai)
497 : na = nsgf_seta(iset)
498 : nb = nppnl
499 : np = ncoa
500 : DO i = 1, maxder
501 : first_col = (i - 1)*ldsab + 1
502 : ! Transformation step projector functions (cartesian->spherical)
503 : ! CALL dgemm("N", "N", ncoa, nppnl, nprjc, 1.0_dp, sab(1, first_col), ldsab, &
504 : ! cprj(1, 1), SIZE(cprj, 1), 0.0_dp, work(1, 1), ldsab)
505 : work(1:np, 1:nb) = MATMUL(sab(1:np, first_col:first_col + nprjc - 1), cprj(1:nprjc, 1:nb))
506 : ! Contraction step (basis functions)
507 : ! CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
508 : ! work(1, 1), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
509 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
510 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, 1:nb))
511 : ! *** Multiply with interaction matrix(h) ***
512 : ncoc = sgfa + nsgf_seta(iset) - 1
513 : DO j = 1, nppnl
514 : clist%achint(sgfa:ncoc, j, i) = clist%acint(sgfa:ncoc, j, i)*hprj(j)
515 : END DO
516 : END DO
517 : END IF
518 : END DO
519 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
520 : clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
521 : IF (.NOT. do_gth .AND. .NOT. do_upf) DEALLOCATE (radp)
522 : END DO
523 :
524 : DEALLOCATE (sab, ai_work, work)
525 : IF (do_soc) DEALLOCATE (lab, work_l)
526 : !$OMP END PARALLEL
527 :
528 : ! Set up a sorting index
529 17370 : CALL sap_sort(sap_int)
530 : ! All integrals needed have been calculated and stored in sap_int
531 : ! We now calculate the Hamiltonian matrix elements
532 :
533 277450 : force_thread = 0.0_dp
534 82390 : at_thread = 0.0_dp
535 17370 : pv_thread = 0.0_dp
536 :
537 : !$OMP PARALLEL &
538 : !$OMP DEFAULT (NONE) &
539 : !$OMP SHARED (do_kp, basis_set, matrix_h, matrix_l, cell_to_index,&
540 : !$OMP sab_orb, matrix_p, sap_int, nkind, eps_ppnl, force, &
541 : !$OMP doat, do_dR, deltaR, maxder, nder, &
542 : !$OMP locks, virial, use_virial, calculate_forces, do_soc, natom) &
543 : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
544 : !$OMP slot, iab, atom_a, f0, irow, icol, h_block, &
545 : !$OMP l_block_x, l_block_y, l_block_z, lock_num, &
546 : !$OMP r_2block, r_3block, atk, &
547 : !$OMP found,p_block, iac, ibc, alist_ac, alist_bc, acint, bcint, &
548 : !$OMP achint, bchint, alkint, blkint, &
549 : !$OMP na, np, nb, katom, j, fa, fb, rbc, rac, &
550 : !$OMP kkind, kac, kbc, i, img, hash, iatom8) &
551 17370 : !$OMP REDUCTION (+ : at_thread, pv_thread, force_thread )
552 :
553 : !$OMP SINGLE
554 : !$ ALLOCATE (locks(nlock))
555 : !$OMP END SINGLE
556 :
557 : !$OMP DO
558 : !$ DO lock_num = 1, nlock
559 : !$ call omp_init_lock(locks(lock_num))
560 : !$ END DO
561 : !$OMP END DO
562 :
563 : !$OMP DO SCHEDULE(GUIDED)
564 : DO slot = 1, sab_orb(1)%nl_size
565 :
566 : ikind = sab_orb(1)%nlist_task(slot)%ikind
567 : jkind = sab_orb(1)%nlist_task(slot)%jkind
568 : iatom = sab_orb(1)%nlist_task(slot)%iatom
569 : jatom = sab_orb(1)%nlist_task(slot)%jatom
570 : cell_b(:) = sab_orb(1)%nlist_task(slot)%cell(:)
571 : rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
572 :
573 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
574 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
575 :
576 : iab = ikind + nkind*(jkind - 1)
577 :
578 : ! Use the symmetry of the first derivatives
579 : IF (iatom == jatom) THEN
580 : f0 = 1.0_dp
581 : ELSE
582 : f0 = 2.0_dp
583 : END IF
584 :
585 : IF (do_kp) THEN
586 : img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
587 : ELSE
588 : img = 1
589 : END IF
590 :
591 : ! Create matrix blocks for a new matrix block column
592 : IF (iatom <= jatom) THEN
593 : irow = iatom
594 : icol = jatom
595 : ELSE
596 : irow = jatom
597 : icol = iatom
598 : END IF
599 : NULLIFY (h_block)
600 : CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found)
601 : IF (do_soc) THEN
602 : NULLIFY (l_block_x, l_block_y, l_block_z)
603 : CALL dbcsr_get_block_p(matrix_l(1, img)%matrix, irow, icol, l_block_x, found)
604 : CALL dbcsr_get_block_p(matrix_l(2, img)%matrix, irow, icol, l_block_y, found)
605 : CALL dbcsr_get_block_p(matrix_l(3, img)%matrix, irow, icol, l_block_z, found)
606 : END IF
607 :
608 : IF (do_dR) THEN
609 : NULLIFY (r_2block, r_3block)
610 : CALL dbcsr_get_block_p(matrix_h(2, img)%matrix, irow, icol, r_2block, found)
611 : CALL dbcsr_get_block_p(matrix_h(3, img)%matrix, irow, icol, r_3block, found)
612 : END IF
613 :
614 : IF (calculate_forces .OR. doat) THEN
615 : NULLIFY (p_block)
616 : CALL dbcsr_get_block_p(matrix_p(1, img)%matrix, irow, icol, p_block, found)
617 : END IF
618 :
619 : ! loop over all kinds for projector atom
620 : IF (ASSOCIATED(h_block)) THEN
621 : !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
622 : !$ hash = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
623 :
624 : DO kkind = 1, nkind
625 : iac = ikind + nkind*(kkind - 1)
626 : ibc = jkind + nkind*(kkind - 1)
627 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
628 : IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
629 : CALL get_alist(sap_int(iac), alist_ac, iatom)
630 : CALL get_alist(sap_int(ibc), alist_bc, jatom)
631 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
632 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
633 : DO kac = 1, alist_ac%nclist
634 : DO kbc = 1, alist_bc%nclist
635 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
636 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
637 : IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
638 : acint => alist_ac%clist(kac)%acint
639 : bcint => alist_bc%clist(kbc)%acint
640 : achint => alist_ac%clist(kac)%achint
641 : bchint => alist_bc%clist(kbc)%achint
642 : IF (do_soc) THEN
643 : alkint => alist_ac%clist(kac)%alkint
644 : blkint => alist_bc%clist(kbc)%alkint
645 : END IF
646 : na = SIZE(acint, 1)
647 : np = SIZE(acint, 2)
648 : nb = SIZE(bcint, 1)
649 : !$ CALL omp_set_lock(locks(hash))
650 : IF (.NOT. do_dR) THEN
651 : IF (iatom <= jatom) THEN
652 : h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
653 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1)))
654 : ELSE
655 : h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
656 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1)))
657 : END IF
658 : END IF
659 : IF (do_soc) THEN
660 : IF (iatom <= jatom) THEN
661 : l_block_x(1:na, 1:nb) = l_block_x(1:na, 1:nb) + &
662 : MATMUL(alkint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1)))
663 : l_block_y(1:na, 1:nb) = l_block_y(1:na, 1:nb) + &
664 : MATMUL(alkint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1)))
665 : l_block_z(1:na, 1:nb) = l_block_z(1:na, 1:nb) + &
666 : MATMUL(alkint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 1)))
667 :
668 : ELSE
669 : l_block_x(1:nb, 1:na) = l_block_x(1:nb, 1:na) + &
670 : MATMUL(blkint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1)))
671 : l_block_y(1:nb, 1:na) = l_block_y(1:nb, 1:na) + &
672 : MATMUL(blkint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 1)))
673 : l_block_z(1:nb, 1:na) = l_block_z(1:nb, 1:na) + &
674 : MATMUL(blkint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 1)))
675 : END IF
676 : END IF
677 : !$ CALL omp_unset_lock(locks(hash))
678 : IF (calculate_forces) THEN
679 : IF (ASSOCIATED(p_block)) THEN
680 : katom = alist_ac%clist(kac)%catom
681 : DO i = 1, 3
682 : j = i + 1
683 : IF (iatom <= jatom) THEN
684 : fa(i) = SUM(p_block(1:na, 1:nb)* &
685 : MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1))))
686 : fb(i) = SUM(p_block(1:na, 1:nb)* &
687 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j))))
688 : ELSE
689 : fa(i) = SUM(p_block(1:nb, 1:na)* &
690 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j))))
691 : fb(i) = SUM(p_block(1:nb, 1:na)* &
692 : MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1))))
693 : END IF
694 : force_thread(i, iatom) = force_thread(i, iatom) + f0*fa(i)
695 : force_thread(i, katom) = force_thread(i, katom) - f0*fa(i)
696 : force_thread(i, jatom) = force_thread(i, jatom) + f0*fb(i)
697 : force_thread(i, katom) = force_thread(i, katom) - f0*fb(i)
698 : END DO
699 :
700 : IF (use_virial) THEN
701 : rac = alist_ac%clist(kac)%rac
702 : rbc = alist_bc%clist(kbc)%rac
703 : CALL virial_pair_force(pv_thread, f0, fa, rac)
704 : CALL virial_pair_force(pv_thread, f0, fb, rbc)
705 : END IF
706 : END IF
707 : END IF
708 :
709 : IF (do_dR) THEN
710 : i = 1; j = 2
711 : katom = alist_ac%clist(kac)%catom
712 : IF (iatom <= jatom) THEN
713 : h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
714 : (deltaR(i, iatom) - deltaR(i, katom))* &
715 : MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1)))
716 :
717 : h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
718 : (deltaR(i, jatom) - deltaR(i, katom))* &
719 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j)))
720 : ELSE
721 : h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
722 : (deltaR(i, iatom) - deltaR(i, katom))* &
723 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j)))
724 : h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
725 : (deltaR(i, jatom) - deltaR(i, katom))* &
726 : MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
727 : END IF
728 :
729 : i = 2; j = 3
730 : katom = alist_ac%clist(kac)%catom
731 : IF (iatom <= jatom) THEN
732 : r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
733 : (deltaR(i, iatom) - deltaR(i, katom))* &
734 : MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1)))
735 :
736 : r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
737 : (deltaR(i, jatom) - deltaR(i, katom))* &
738 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j)))
739 : ELSE
740 : r_2block(1:nb, 1:na) = r_2block(1:nb, 1:na) + &
741 : (deltaR(i, iatom) - deltaR(i, katom))* &
742 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j)))
743 : r_2block(1:nb, 1:na) = r_2block(1:nb, 1:na) + &
744 : (deltaR(i, jatom) - deltaR(i, katom))* &
745 : MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
746 : END IF
747 :
748 : i = 3; j = 4
749 : katom = alist_ac%clist(kac)%catom
750 : IF (iatom <= jatom) THEN
751 : r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &
752 : (deltaR(i, iatom) - deltaR(i, katom))* &
753 : MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1)))
754 :
755 : r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &
756 : (deltaR(i, jatom) - deltaR(i, katom))* &
757 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j)))
758 : ELSE
759 : r_3block(1:nb, 1:na) = r_3block(1:nb, 1:na) + &
760 : (deltaR(i, iatom) - deltaR(i, katom))* &
761 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j)))
762 : r_3block(1:nb, 1:na) = r_3block(1:nb, 1:na) + &
763 : (deltaR(i, jatom) - deltaR(i, katom))* &
764 : MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
765 : END IF
766 :
767 : END IF
768 : IF (doat) THEN
769 : IF (ASSOCIATED(p_block)) THEN
770 : katom = alist_ac%clist(kac)%catom
771 : IF (iatom <= jatom) THEN
772 : atk = SUM(p_block(1:na, 1:nb)* &
773 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1))))
774 : ELSE
775 : atk = SUM(p_block(1:nb, 1:na)* &
776 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1))))
777 : END IF
778 : at_thread(katom) = at_thread(katom) + f0*atk
779 : END IF
780 : END IF
781 : EXIT ! We have found a match and there can be only one single match
782 : END IF
783 : END DO
784 : END DO
785 : END DO
786 : END IF
787 : END DO
788 :
789 : !$OMP DO
790 : !$ DO lock_num = 1, nlock
791 : !$ call omp_destroy_lock(locks(lock_num))
792 : !$ END DO
793 : !$OMP END DO
794 :
795 : !$OMP SINGLE
796 : !$ DEALLOCATE (locks)
797 : !$OMP END SINGLE NOWAIT
798 :
799 : !$OMP END PARALLEL
800 :
801 17370 : CALL release_sap_int(sap_int)
802 :
803 17370 : DEALLOCATE (basis_set, gpotential, spotential)
804 17370 : IF (calculate_forces) THEN
805 6863 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
806 : !$OMP DO
807 : DO iatom = 1, natom
808 24865 : atom_a = atom_of_kind(iatom)
809 24865 : ikind = kind_of(iatom)
810 99460 : force(ikind)%gth_ppnl(:, atom_a) = force(ikind)%gth_ppnl(:, atom_a) + force_thread(:, iatom)
811 : END DO
812 : !$OMP END DO
813 6863 : DEALLOCATE (atom_of_kind, kind_of)
814 : END IF
815 :
816 17370 : IF (calculate_forces .AND. use_virial) THEN
817 11284 : virial%pv_ppnl = virial%pv_ppnl + pv_thread
818 11284 : virial%pv_virial = virial%pv_virial + pv_thread
819 : END IF
820 :
821 17370 : IF (doat) THEN
822 280 : atcore(1:natom) = atcore(1:natom) + at_thread
823 : END IF
824 :
825 34740 : IF (calculate_forces .OR. doat) THEN
826 : ! If LSD, then recover alpha density and beta density
827 : ! from the total density (1) and the spin density (2)
828 6925 : IF (SIZE(matrix_p, 1) == 2) THEN
829 3042 : DO img = 1, nimages
830 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
831 2332 : alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
832 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
833 3042 : alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
834 : END DO
835 : END IF
836 : END IF
837 :
838 : END IF !ppnl_present
839 :
840 17370 : CALL timestop(handle)
841 :
842 34740 : END SUBROUTINE build_core_ppnl
843 :
844 : END MODULE core_ppnl
|
