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 : ! **************************************************************************************************
9 : !> \brief A collection of functions used by CNEO-DFT
10 : !> (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
11 : !> \par History
12 : !> 08.2025 created [zc62]
13 : !> \author Zehua Chen
14 : ! **************************************************************************************************
15 : MODULE qs_cneo_methods
16 : USE ai_verfc, ONLY: verfc
17 : USE ao_util, ONLY: trace_r_AxB
18 : USE atom_operators, ONLY: atom_int_release,&
19 : atom_int_setup
20 : USE atom_types, ONLY: CGTO_BASIS,&
21 : atom_basis_type,&
22 : atom_integrals,&
23 : lmat,&
24 : release_atom_basis
25 : USE atomic_kind_types, ONLY: atomic_kind_type,&
26 : get_atomic_kind,&
27 : get_atomic_kind_set
28 : USE basis_set_types, ONLY: get_gto_basis_set,&
29 : gto_basis_set_type
30 : USE bibliography, ONLY: Chen2025,&
31 : cite_reference
32 : USE core_ae, ONLY: verfc_force
33 : USE cp_control_types, ONLY: gapw_control_type
34 : USE cp_log_handling, ONLY: cp_to_string
35 : USE distribution_1d_types, ONLY: distribution_1d_type
36 : USE kinds, ONLY: dp
37 : USE mathconstants, ONLY: dfac,&
38 : fourpi,&
39 : pi
40 : USE mathlib, ONLY: get_pseudo_inverse_svd
41 : USE memory_utilities, ONLY: reallocate
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE orbital_pointers, ONLY: indso,&
44 : indso_inv,&
45 : init_orbital_pointers,&
46 : ncoset,&
47 : nso,&
48 : nsoset
49 : USE particle_types, ONLY: particle_type
50 : USE physcon, ONLY: massunit
51 : USE qs_cneo_types, ONLY: allocate_rhoz_cneo_set,&
52 : cneo_potential_type,&
53 : get_cneo_potential,&
54 : rhoz_cneo_type,&
55 : set_cneo_potential
56 : USE qs_cneo_utils, ONLY: atom_solve_cneo,&
57 : cneo_gather,&
58 : create_harmonics_atom_cneo,&
59 : create_my_CG_cneo,&
60 : get_maxl_CG_cneo
61 : USE qs_environment_types, ONLY: get_qs_env,&
62 : qs_environment_type
63 : USE qs_force_types, ONLY: qs_force_type
64 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
65 : grid_atom_type
66 : USE qs_harmonics_atom, ONLY: allocate_harmonics_atom,&
67 : get_none0_cg_list,&
68 : harmonics_atom_type
69 : USE qs_kind_types, ONLY: get_qs_kind,&
70 : get_qs_kind_set,&
71 : qs_kind_type
72 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
73 : neighbor_list_iterator_create,&
74 : neighbor_list_iterator_p_type,&
75 : neighbor_list_iterator_release,&
76 : neighbor_list_set_p_type,&
77 : nl_set_sub_iterator,&
78 : nl_sub_iterate
79 : USE util, ONLY: get_limit
80 : USE virial_methods, ONLY: virial_pair_force
81 : USE virial_types, ONLY: virial_type
82 : USE whittaker, ONLY: whittaker_c0a,&
83 : whittaker_ci
84 :
85 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
86 :
87 : #include "./base/base_uses.f90"
88 :
89 : IMPLICIT NONE
90 :
91 : PRIVATE
92 :
93 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cneo_methods'
94 :
95 : PUBLIC :: allocate_rhoz_cneo_internals, calculate_rhoz_cneo, cneo_core_matrices, &
96 : init_cneo_potential_internals, Vh_1c_nuc_integrals
97 :
98 : CONTAINS
99 :
100 : ! **************************************************************************************************
101 : !> \brief ...
102 : !> \param potential ...
103 : !> \param nuc_basis ...
104 : !> \param nuc_soft_basis ...
105 : !> \param gapw_control ...
106 : !> \param grid_atom ...
107 : ! **************************************************************************************************
108 8 : SUBROUTINE init_cneo_potential_internals(potential, nuc_basis, nuc_soft_basis, gapw_control, grid_atom)
109 :
110 : TYPE(cneo_potential_type), POINTER :: potential
111 : TYPE(gto_basis_set_type), POINTER :: nuc_basis, nuc_soft_basis
112 : TYPE(gapw_control_type), POINTER :: gapw_control
113 : TYPE(grid_atom_type), POINTER :: grid_atom
114 :
115 : CHARACTER(len=*), PARAMETER :: routineN = 'init_cneo_potential_internals'
116 :
117 : INTEGER :: handle, i, icg, ico, ii, ipgf, ipgf1, ipgf2, ir, is1, is2, iset, iset1, iset2, &
118 : iso, iso1, iso2, iso_pgf, iso_set, j, k, k1, k2, l, l_iso, l_sub, l_sum, ll, llmax, &
119 : lmax12, lmax_expansion, lmax_sphere, lmin12, m, m1, m2, max_iso_not0, max_iso_not0_local, &
120 : max_s, max_s_harm, maxl, maxso, n1, n2, nl, nne, npgf2, npgf_sum, npsgf, nr, ns, nset, &
121 : nsgf, nsotot, nsox
122 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
123 8 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
124 : INTEGER, DIMENSION(0:lmat, 100) :: set_index, shell_index
125 8 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, n2oindex, npgf, npgf_s, &
126 8 : nshell, o2nindex
127 8 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, ls
128 8 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: done_vgg
129 : REAL(KIND=dp) :: c1, c2, gcc_tmp, mass, massinv, &
130 : root_zet12, scal, scal1, zet12
131 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: erf_zet12, g1, g2, gg0, int1, int2
132 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
133 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dist
134 8 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: kin, my_gcc_h, my_gcc_s, oorad2l, ovlp, &
135 8 : rad2l, utrans, zet
136 8 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: distance, gcc_h, gcc_s, gg, my_CG
137 8 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: vgg
138 : TYPE(atom_basis_type), POINTER :: basis
139 : TYPE(atom_integrals), POINTER :: integrals
140 : TYPE(grid_atom_type), POINTER :: grid
141 : TYPE(harmonics_atom_type), POINTER :: harmonics
142 :
143 0 : CPASSERT(ASSOCIATED(potential))
144 8 : CPASSERT(ASSOCIATED(nuc_basis))
145 8 : CPASSERT(ASSOCIATED(nuc_soft_basis))
146 :
147 8 : CALL cite_reference(Chen2025)
148 :
149 8 : CALL timeset(routineN, handle)
150 :
151 : CALL get_cneo_potential(potential, my_gcc_h=my_gcc_h, my_gcc_s=my_gcc_s, &
152 : ovlp=ovlp, kin=kin, utrans=utrans, distance=distance, &
153 : harmonics=harmonics, gg=gg, vgg=vgg, n2oindex=n2oindex, &
154 8 : o2nindex=o2nindex, rad2l=rad2l, oorad2l=oorad2l)
155 8 : CPASSERT(.NOT. ASSOCIATED(my_gcc_h))
156 8 : CPASSERT(.NOT. ASSOCIATED(my_gcc_s))
157 8 : CPASSERT(.NOT. ASSOCIATED(ovlp))
158 8 : CPASSERT(.NOT. ASSOCIATED(kin))
159 8 : CPASSERT(.NOT. ASSOCIATED(utrans))
160 8 : CPASSERT(.NOT. ASSOCIATED(distance))
161 8 : CPASSERT(.NOT. ASSOCIATED(harmonics))
162 8 : CPASSERT(.NOT. ASSOCIATED(gg))
163 8 : CPASSERT(.NOT. ASSOCIATED(vgg))
164 8 : CPASSERT(.NOT. ASSOCIATED(n2oindex))
165 8 : CPASSERT(.NOT. ASSOCIATED(o2nindex))
166 8 : CPASSERT(.NOT. ASSOCIATED(rad2l))
167 8 : CPASSERT(.NOT. ASSOCIATED(oorad2l))
168 :
169 : ! ovlp, kin and utrans parts are mostly copied from atom_kind_orbitals::calculate_atomic_orbitals
170 : ! and atom_set_basis::set_kind_basis_atomic
171 8 : NULLIFY (basis, integrals, grid)
172 1848 : ALLOCATE (basis, integrals)
173 8 : CALL allocate_grid_atom(grid)
174 8 : basis%grid => grid
175 8 : NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
176 : ! fill in the basis data structures
177 8 : basis%basis_type = CGTO_BASIS
178 8 : basis%eps_eig = 1.e-12_dp
179 :
180 8 : NULLIFY (nshell, npgf, lmin, lmax, ls, zet, gcc_h, first_sgf)
181 : CALL get_gto_basis_set(nuc_basis, nset=nset, nshell=nshell, npgf=npgf, lmin=lmin, &
182 : lmax=lmax, l=ls, nsgf=nsgf, zet=zet, gcc=gcc_h, first_sgf=first_sgf, &
183 8 : maxl=maxl, maxso=maxso, npgf_sum=npgf_sum)
184 8 : NULLIFY (npgf_s, gcc_s)
185 8 : CALL get_gto_basis_set(nuc_soft_basis, npgf=npgf_s, gcc=gcc_s)
186 : ! There is such a limitation because we rely on atomic code to build S, T and U.
187 : ! Usually l=5 is more than enough, suppoting PB6H basis.
188 8 : IF (maxl > lmat) THEN
189 : CALL cp_abort(__LOCATION__, "Nuclear basis with angular momentum higher than "// &
190 0 : "atom_types::lmat is not supported yet.")
191 : END IF
192 :
193 8 : set_index = 0
194 8 : shell_index = 0
195 56 : basis%nprim = 0
196 56 : basis%nbas = 0
197 80 : DO i = 1, nset
198 144 : DO j = lmin(i), MIN(lmax(i), lmat)
199 144 : basis%nprim(j) = basis%nprim(j) + npgf(i)
200 : END DO
201 152 : DO j = 1, nshell(i)
202 72 : l = ls(j, i)
203 144 : IF (l <= lmat) THEN
204 72 : basis%nbas(l) = basis%nbas(l) + 1
205 72 : k = basis%nbas(l)
206 72 : CPASSERT(k <= 100)
207 72 : set_index(l, k) = i
208 72 : shell_index(l, k) = j
209 : END IF
210 : END DO
211 : END DO
212 :
213 56 : nl = MAXVAL(basis%nprim)
214 56 : ns = MAXVAL(basis%nbas)
215 24 : ALLOCATE (basis%am(nl, 0:lmat))
216 248 : basis%am = 0._dp
217 40 : ALLOCATE (basis%cm(nl, ns, 0:lmat))
218 1016 : basis%cm = 0._dp
219 56 : DO l = 0, lmat
220 : nl = 0
221 : ns = 0
222 488 : DO i = 1, nset
223 480 : IF (l >= lmin(i) .AND. l <= lmax(i)) THEN
224 144 : DO ipgf = 1, npgf(i)
225 144 : basis%am(nl + ipgf, l) = zet(ipgf, i)
226 : END DO
227 144 : DO ii = 1, nshell(i)
228 144 : IF (ls(ii, i) == l) THEN
229 72 : ns = ns + 1
230 144 : DO ipgf = 1, npgf(i)
231 144 : basis%cm(nl + ipgf, ns, l) = gcc_h(ipgf, ii, i) ! NOTE: not normalized
232 : END DO
233 : END IF
234 : END DO
235 72 : nl = nl + npgf(i)
236 : END IF
237 : END DO
238 : END DO
239 :
240 : ! overlap, kinetic and transformation matrices
241 8 : CALL atom_int_setup(integrals, basis)
242 :
243 : ! make the integrals full matrix form
244 64 : ALLOCATE (ovlp(nsgf, nsgf), kin(nsgf, nsgf), utrans(nsgf, nsgf))
245 4424 : ovlp = 0.0_dp
246 4424 : kin = 0.0_dp
247 4424 : utrans = 0.0_dp
248 8 : CALL get_cneo_potential(potential, mass=mass)
249 8 : mass = mass*massunit
250 8 : massinv = 1._dp/mass
251 8 : nne = 0 ! number of linear-independent spherical basis functions
252 56 : DO l = 0, lmat
253 48 : ll = 2*l
254 120 : DO k2 = 1, integrals%nne(l)
255 304 : DO m = 0, ll
256 184 : nne = nne + 1
257 760 : DO k1 = 1, basis%nbas(l)
258 504 : scal1 = SQRT(integrals%ovlp(k1, k1, l))
259 504 : i = first_sgf(shell_index(l, k1), set_index(l, k1))
260 688 : utrans(i + m, nne) = integrals%utrans(k1, k2, l)*scal1
261 : END DO
262 : END DO
263 : END DO
264 128 : DO k1 = 1, basis%nbas(l)
265 72 : scal1 = 1._dp/SQRT(integrals%ovlp(k1, k1, l))
266 72 : i = first_sgf(shell_index(l, k1), set_index(l, k1))
267 352 : DO k2 = 1, basis%nbas(l)
268 232 : scal = scal1/SQRT(integrals%ovlp(k2, k2, l))
269 232 : j = first_sgf(shell_index(l, k2), set_index(l, k2))
270 808 : DO m = 0, ll
271 : ! normalize the integrals
272 504 : ovlp(i + m, j + m) = integrals%ovlp(k1, k2, l)*scal
273 736 : kin(i + m, j + m) = integrals%kin(k1, k2, l)*scal*massinv
274 : END DO
275 : END DO
276 : END DO
277 : END DO
278 :
279 8 : nsotot = maxso*nset
280 48 : ALLOCATE (my_gcc_h(nsotot, nsgf), my_gcc_s(nsotot, nsgf))
281 15096 : my_gcc_h = 0.0_dp
282 15096 : my_gcc_s = 0.0_dp
283 : ! create gcc that really 3D-normalize the basis functions
284 32 : DO l = 0, MIN(maxl, lmat)
285 24 : ns = 0
286 24 : m = 0
287 24 : ll = 2*l
288 24 : k = nsoset(l - 1) + 1
289 248 : DO i = 1, nset
290 216 : IF (l >= lmin(i) .AND. l <= lmax(i)) THEN
291 72 : nsox = nsoset(lmax(i))
292 144 : DO ii = 1, nshell(i)
293 144 : IF (ls(ii, i) == l) THEN
294 72 : ns = ns + 1
295 72 : k1 = first_sgf(shell_index(l, ns), set_index(l, ns))
296 72 : scal = 1._dp/SQRT(integrals%ovlp(ns, ns, l))
297 144 : DO ipgf = 1, npgf(i)
298 72 : gcc_tmp = gcc_h(ipgf, ii, i)*scal
299 72 : k2 = (ipgf - 1)*nsox + m
300 328 : DO j = 0, ll
301 256 : my_gcc_h(k + k2 + j, k1 + j) = gcc_tmp
302 : END DO
303 : END DO
304 86 : DO ipgf = 1, npgf_s(i)
305 14 : gcc_tmp = gcc_s(ipgf, ii, i)*scal
306 14 : k2 = (ipgf - 1)*nsox + m
307 112 : DO j = 0, ll
308 40 : my_gcc_s(k + k2 + j, k1 + j) = gcc_tmp
309 : END DO
310 : END DO
311 : END IF
312 : END DO
313 : END IF
314 240 : m = m + maxso
315 : END DO
316 : END DO
317 :
318 8 : CALL atom_int_release(integrals)
319 : CALL set_cneo_potential(potential, nsgf=nsgf, nne=nne, nsotot=nsotot, &
320 : my_gcc_h=my_gcc_h, my_gcc_s=my_gcc_s, &
321 8 : ovlp=ovlp, kin=kin, utrans=utrans)
322 8 : CALL release_atom_basis(basis)
323 8 : DEALLOCATE (basis, integrals)
324 :
325 : ! initialize my_CG
326 8 : lmax_sphere = gapw_control%lmax_sphere
327 : ! make sure llmax is at least 1 such that distance matrices can be generated
328 8 : llmax = MAX(1, MIN(lmax_sphere, 2*maxl))
329 8 : max_s_harm = nsoset(llmax)
330 8 : max_s = nsoset(maxl)
331 8 : NULLIFY (my_CG)
332 8 : CALL reallocate(my_CG, 1, max_s, 1, max_s, 1, max_s_harm)
333 8 : CALL create_my_CG_cneo(my_CG, MAX(llmax, 2*maxl, 1), maxl, llmax)
334 :
335 : ! initialize harmonics
336 8 : CALL allocate_harmonics_atom(harmonics)
337 8 : CALL create_harmonics_atom_cneo(harmonics, my_CG, llmax, max_s, max_s_harm)
338 8 : DEALLOCATE (my_CG)
339 8 : CALL get_maxl_CG_cneo(harmonics, nuc_basis, llmax, max_s_harm)
340 :
341 8 : CALL set_cneo_potential(potential, harmonics=harmonics)
342 :
343 : ! initialize my own rad2l and oorad2l
344 : ! copied from qs_grid_atom::create_grid_atom
345 8 : nr = grid_atom%nr
346 8 : NULLIFY (rad2l, oorad2l)
347 8 : CALL reallocate(rad2l, 1, nr, 0, llmax + 1)
348 8 : CALL reallocate(oorad2l, 1, nr, 0, llmax + 1)
349 408 : rad2l(:, 0) = 1._dp
350 408 : oorad2l(:, 0) = 1._dp
351 48 : DO l = 1, llmax + 1
352 4040 : rad2l(:, l) = rad2l(:, l - 1)*grid_atom%rad(:)
353 4048 : oorad2l(:, l) = oorad2l(:, l - 1)/grid_atom%rad(:)
354 : END DO
355 8 : CALL set_cneo_potential(potential, rad2l=rad2l, oorad2l=oorad2l)
356 : ! still need to bump lmax in grid_atom as qs_rho0_types::calculate_g0 uses it
357 8 : IF (SIZE(rad2l, 2) > SIZE(grid_atom%rad2l, 2)) THEN
358 2 : CPASSERT(SIZE(rad2l, 1) == SIZE(grid_atom%rad2l, 1))
359 2 : DEALLOCATE (grid_atom%rad2l)
360 2 : NULLIFY (grid_atom%rad2l)
361 2 : CALL reallocate(grid_atom%rad2l, 1, nr, 0, llmax + 1)
362 1226 : grid_atom%rad2l = rad2l
363 : END IF
364 8 : IF (SIZE(oorad2l, 2) > SIZE(grid_atom%oorad2l, 2)) THEN
365 2 : CPASSERT(SIZE(oorad2l, 1) == SIZE(grid_atom%oorad2l, 1))
366 2 : DEALLOCATE (grid_atom%oorad2l)
367 2 : NULLIFY (grid_atom%oorad2l)
368 2 : CALL reallocate(grid_atom%oorad2l, 1, nr, 0, llmax + 1)
369 1226 : grid_atom%oorad2l = oorad2l
370 : END IF
371 :
372 : ! distance matrices
373 72 : ALLOCATE (distance(nsgf, nsgf, 3), dist(nsotot, nsotot, 3))
374 13280 : distance = 0.0_dp
375 8 : dist = 0.0_dp
376 : ! initialize gg and vgg
377 : ! mostly copied from qs_rho_atom_methods::calculate_rho_atom
378 8 : max_iso_not0 = harmonics%max_iso_not0
379 8 : lmax_expansion = indso(1, max_iso_not0)
380 8 : my_CG => harmonics%my_CG
381 72 : ALLOCATE (g1(nr), g2(nr), gg0(nr), gg(nr, 0:2*maxl, npgf_sum*(npgf_sum + 1)/2))
382 56 : ALLOCATE (erf_zet12(nr), vgg(nr, 0:2*maxl, 0:indso(1, max_iso_not0), npgf_sum*(npgf_sum + 1)/2))
383 32 : ALLOCATE (done_vgg(0:2*maxl, 0:indso(1, max_iso_not0)))
384 24 : ALLOCATE (int1(nr), int2(nr))
385 48 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
386 :
387 8 : j = 0
388 8 : m1 = 0
389 80 : DO iset1 = 1, nset
390 72 : n1 = nsoset(lmax(iset1))
391 72 : m2 = 0
392 432 : DO iset2 = 1, iset1
393 360 : n2 = nsoset(lmax(iset2))
394 :
395 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
396 360 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
397 360 : CPASSERT(max_iso_not0_local <= max_iso_not0)
398 :
399 720 : DO ipgf1 = 1, npgf(iset1)
400 18360 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
401 :
402 360 : IF (iset2 == iset1) THEN
403 : npgf2 = ipgf1
404 : ELSE
405 288 : npgf2 = npgf(iset2)
406 : END IF
407 1080 : DO ipgf2 = 1, npgf2
408 360 : zet12 = zet(ipgf1, iset1) + zet(ipgf2, iset2)
409 :
410 : ! distance part
411 : ! -1 -> y -> 2, 0 -> z -> 3, 1 -> x -> 1
412 1440 : DO m = -1, 1
413 1080 : k = m + 3
414 1080 : IF (m == 1) k = 1
415 1080 : iso = indso_inv(1, m)
416 2256 : DO icg = 1, cg_n_list(iso)
417 816 : is1 = cg_list(1, icg, iso)
418 816 : is2 = cg_list(2, icg, iso)
419 :
420 816 : iso1 = is1 + n1*(ipgf1 - 1) + m1
421 816 : iso2 = is2 + n2*(ipgf2 - 1) + m2
422 :
423 816 : l = indso(1, is1) + indso(1, is2)
424 : dist(iso1, iso2, k) = dist(iso1, iso2, k) + my_CG(is1, is2, iso)* &
425 : pi*dfac(l + 2)/ &
426 816 : ((2.0_dp*zet12)**((l + 3)/2)*SQRT(3.0_dp*zet12))
427 1896 : dist(iso2, iso1, k) = dist(iso1, iso2, k) ! symmetric
428 : END DO !icg
429 : END DO
430 :
431 : ! gg and vgg part
432 360 : j = j + 1
433 18360 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
434 360 : lmin12 = lmin(iset1) + lmin(iset2)
435 360 : lmax12 = lmax(iset1) + lmax(iset2)
436 :
437 360 : root_zet12 = SQRT(zet12)
438 18360 : DO ir = 1, nr
439 18360 : erf_zet12(ir) = erf(root_zet12*grid_atom%rad(ir))
440 : END DO
441 :
442 92160 : gg(:, :, j) = 0.0_dp
443 461160 : vgg(:, :, :, j) = 0.0_dp
444 360 : done_vgg = .FALSE.
445 : ! reduce the number of terms in the expansion local densities
446 720 : IF (lmin12 <= lmax_expansion) THEN
447 360 : IF (lmin12 == 0) THEN
448 4080 : gg(1:nr, lmin12, j) = g1(1:nr)*g2(1:nr)
449 4080 : gg0(1:nr) = gg(1:nr, lmin12, j)
450 : ELSE
451 14280 : gg0(1:nr) = g1(1:nr)*g2(1:nr)
452 28280 : gg(1:nr, lmin12, j) = rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
453 : END IF
454 :
455 : ! reduce the number of terms in the expansion local densities
456 360 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
457 :
458 360 : DO l = lmin12 + 1, lmax12
459 360 : gg(1:nr, l, j) = grid_atom%rad(1:nr)*gg(1:nr, l - 1, j)
460 : END DO
461 :
462 360 : c2 = SQRT(pi*pi*pi/(zet12*zet12*zet12))
463 :
464 3200 : DO iso = 1, max_iso_not0_local
465 2840 : l_iso = indso(1, iso)
466 2840 : c1 = fourpi/(2._dp*REAL(l_iso, dp) + 1._dp)
467 7704 : DO icg = 1, cg_n_list(iso)
468 4504 : iso1 = cg_list(1, icg, iso)
469 4504 : iso2 = cg_list(2, icg, iso)
470 :
471 4504 : l = indso(1, iso1) + indso(1, iso2)
472 4504 : CPASSERT(l <= lmax_expansion)
473 4504 : IF (done_vgg(l, l_iso)) CYCLE
474 504 : L_sum = l + l_iso
475 504 : L_sub = l - l_iso
476 :
477 504 : IF (l_sum == 0) THEN
478 8080 : vgg(1:nr, l, l_iso, j) = erf_zet12(1:nr)*oorad2l(1:nr, 1)*c2
479 : ELSE
480 424 : CALL whittaker_c0a(int1, grid_atom%rad, gg0, erf_zet12, zet12, l, l_iso, nr)
481 424 : CALL whittaker_ci(int2, grid_atom%rad, gg0, zet12, L_sub, nr)
482 :
483 21624 : DO ir = 1, nr
484 21200 : int2(ir) = rad2l(ir, l_iso)*int2(ir)
485 21624 : vgg(ir, l, l_iso, j) = c1*(int1(ir) + int2(ir))
486 : END DO
487 : END IF
488 7344 : done_vgg(l, l_iso) = .TRUE.
489 : END DO
490 : END DO
491 : END IF ! lmax_expansion
492 :
493 : END DO ! ipgf2
494 : END DO ! ipgf1
495 792 : m2 = m2 + maxso
496 : END DO ! iset2
497 80 : m1 = m1 + maxso
498 : END DO ! iset1
499 :
500 8 : DEALLOCATE (g1, g2, gg0, erf_zet12, int1, int2, done_vgg)
501 8 : DEALLOCATE (cg_list, cg_n_list)
502 :
503 24 : ALLOCATE (work(nsotot, nsgf))
504 32 : DO k = 1, 3
505 : CALL dgemm("N", "N", nsotot, nsgf, nsotot, 1.0_dp, dist(:, :, k), nsotot, my_gcc_h, &
506 24 : nsotot, 0.0_dp, work, nsotot)
507 : CALL dgemm("T", "N", nsgf, nsgf, nsotot, 1.0_dp, my_gcc_h, nsotot, work, &
508 32 : nsotot, 0.0_dp, distance(:, :, k), nsgf)
509 : END DO
510 8 : DEALLOCATE (work, dist)
511 8 : CALL set_cneo_potential(potential, distance=distance, gg=gg, vgg=vgg)
512 :
513 : ! Index transformation OLD-NEW
514 : ! copied from paw_proj_set_types::build_projector
515 24 : ALLOCATE (o2nindex(nsotot))
516 16 : ALLOCATE (n2oindex(nsotot))
517 656 : o2nindex = 0
518 656 : n2oindex = 0
519 : ico = 1
520 80 : DO iset = 1, nset
521 72 : iso_set = (iset - 1)*maxso + 1
522 72 : nsox = nsoset(lmax(iset))
523 152 : DO ipgf = 1, npgf(iset)
524 72 : iso_pgf = iso_set + (ipgf - 1)*nsox
525 72 : iso = iso_pgf + nsoset(lmin(iset) - 1)
526 216 : DO l = lmin(iset), lmax(iset)
527 328 : DO k = 1, nso(l)
528 184 : n2oindex(ico) = iso
529 184 : o2nindex(iso) = ico
530 184 : iso = iso + 1
531 256 : ico = ico + 1
532 : END DO
533 : END DO
534 : END DO
535 : END DO
536 8 : npsgf = ico - 1
537 8 : CALL set_cneo_potential(potential, npsgf=npsgf, n2oindex=n2oindex, o2nindex=o2nindex)
538 :
539 8 : CALL timestop(handle)
540 :
541 32 : END SUBROUTINE init_cneo_potential_internals
542 :
543 : ! **************************************************************************************************
544 : !> \brief ...
545 : !> \param rhoz_cneo_set ...
546 : !> \param atomic_kind_set ...
547 : !> \param qs_kind_set ...
548 : !> \param qs_env ...
549 : ! **************************************************************************************************
550 16 : SUBROUTINE allocate_rhoz_cneo_internals(rhoz_cneo_set, atomic_kind_set, &
551 : qs_kind_set, qs_env)
552 :
553 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
554 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
555 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
556 : TYPE(qs_environment_type), POINTER :: qs_env
557 :
558 : CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rhoz_cneo_internals'
559 :
560 : INTEGER :: bo(2), handle, iat, iatom, ikind, &
561 : max_iso_not0, mepos, nat, natom, &
562 : npsgf, nr, nsgf, nsotot, num_pe
563 8 : INTEGER, DIMENSION(:), POINTER :: atom_list
564 : LOGICAL :: paw_atom
565 : TYPE(cneo_potential_type), POINTER :: cneo_potential
566 : TYPE(mp_para_env_type), POINTER :: para_env
567 :
568 8 : CALL timeset(routineN, handle)
569 :
570 8 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
571 :
572 8 : CALL allocate_rhoz_cneo_set(rhoz_cneo_set, natom)
573 :
574 8 : NULLIFY (para_env)
575 8 : CALL get_qs_env(qs_env, para_env=para_env)
576 :
577 22 : DO ikind = 1, SIZE(atomic_kind_set)
578 :
579 14 : NULLIFY (cneo_potential)
580 : CALL get_qs_kind(qs_kind_set(ikind), &
581 : ngrid_rad=nr, &
582 : paw_atom=paw_atom, &
583 14 : cneo_potential=cneo_potential)
584 :
585 22 : IF (ASSOCIATED(cneo_potential)) THEN
586 8 : CPASSERT(paw_atom)
587 :
588 8 : NULLIFY (atom_list)
589 8 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
590 :
591 8 : nsgf = cneo_potential%nsgf
592 8 : npsgf = cneo_potential%npsgf
593 8 : nsotot = cneo_potential%nsotot
594 :
595 22 : DO iat = 1, nat
596 14 : iatom = atom_list(iat)
597 :
598 : ! density matrices, core and soft vmat will be broadcast to all processes
599 : ALLOCATE (rhoz_cneo_set(iatom)%pmat(1:nsgf, 1:nsgf), &
600 : rhoz_cneo_set(iatom)%cpc_h(1:npsgf, 1:npsgf), &
601 : rhoz_cneo_set(iatom)%cpc_s(1:npsgf, 1:npsgf), &
602 : rhoz_cneo_set(iatom)%core(1:nsgf, 1:nsgf), &
603 182 : rhoz_cneo_set(iatom)%vmat(1:nsgf, 1:nsgf))
604 7742 : rhoz_cneo_set(iatom)%pmat = 0.0_dp
605 7742 : rhoz_cneo_set(iatom)%cpc_h = 0.0_dp
606 7742 : rhoz_cneo_set(iatom)%cpc_s = 0.0_dp
607 7742 : rhoz_cneo_set(iatom)%core = 0.0_dp
608 7750 : rhoz_cneo_set(iatom)%vmat = 0.0_dp
609 : END DO
610 :
611 8 : max_iso_not0 = cneo_potential%harmonics%max_iso_not0
612 8 : num_pe = para_env%num_pe
613 8 : mepos = para_env%mepos
614 8 : bo = get_limit(nat, num_pe, mepos)
615 15 : DO iat = bo(1), bo(2)
616 7 : iatom = atom_list(iat)
617 :
618 : ALLOCATE (rhoz_cneo_set(iatom)%fmat(1:nsgf, 1:nsgf), &
619 42 : rhoz_cneo_set(iatom)%wfn(1:nsgf, 1:nsgf))
620 3871 : rhoz_cneo_set(iatom)%fmat = 0.0_dp
621 3871 : rhoz_cneo_set(iatom)%wfn = 0.0_dp
622 :
623 : ALLOCATE (rhoz_cneo_set(iatom)%rho_rad_h(1:nr, 1:max_iso_not0), &
624 : rhoz_cneo_set(iatom)%rho_rad_s(1:nr, 1:max_iso_not0), &
625 : rhoz_cneo_set(iatom)%vrho_rad_h(1:nr, 1:max_iso_not0), &
626 70 : rhoz_cneo_set(iatom)%vrho_rad_s(1:nr, 1:max_iso_not0))
627 8932 : rhoz_cneo_set(iatom)%rho_rad_h = 0.0_dp
628 8932 : rhoz_cneo_set(iatom)%rho_rad_s = 0.0_dp
629 8932 : rhoz_cneo_set(iatom)%vrho_rad_h = 0.0_dp
630 8932 : rhoz_cneo_set(iatom)%vrho_rad_s = 0.0_dp
631 :
632 7 : NULLIFY (rhoz_cneo_set(iatom)%ga_Vlocal_gb_h)
633 7 : CALL reallocate(rhoz_cneo_set(iatom)%ga_Vlocal_gb_h, 1, nsotot, 1, nsotot)
634 46501 : rhoz_cneo_set(iatom)%ga_Vlocal_gb_h = 0.0_dp
635 7 : NULLIFY (rhoz_cneo_set(iatom)%ga_Vlocal_gb_s)
636 7 : CALL reallocate(rhoz_cneo_set(iatom)%ga_Vlocal_gb_s, 1, nsotot, 1, nsotot)
637 46509 : rhoz_cneo_set(iatom)%ga_Vlocal_gb_s = 0.0_dp
638 : END DO ! iat
639 : END IF
640 :
641 : END DO
642 :
643 8 : CALL timestop(handle)
644 :
645 8 : END SUBROUTINE allocate_rhoz_cneo_internals
646 :
647 : ! **************************************************************************************************
648 : !> \brief ...
649 : !> \param qs_env ...
650 : !> \param calculate_forces ...
651 : !> \param nder ...
652 : ! **************************************************************************************************
653 18707 : SUBROUTINE cneo_core_matrices(qs_env, calculate_forces, nder)
654 : TYPE(qs_environment_type), POINTER :: qs_env
655 : LOGICAL, INTENT(IN) :: calculate_forces
656 : INTEGER, INTENT(IN) :: nder
657 :
658 : LOGICAL :: use_virial
659 18707 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
660 : TYPE(distribution_1d_type), POINTER :: distribution_1d
661 : TYPE(mp_para_env_type), POINTER :: para_env
662 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
663 18707 : POINTER :: sab_cneo
664 18707 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
665 18707 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
666 18707 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
667 18707 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
668 : TYPE(virial_type), POINTER :: virial
669 :
670 18707 : NULLIFY (rhoz_cneo_set)
671 18707 : CALL get_qs_env(qs_env=qs_env, rhoz_cneo_set=rhoz_cneo_set)
672 :
673 18707 : IF (ASSOCIATED(rhoz_cneo_set)) THEN
674 14 : NULLIFY (force, virial)
675 : ! force
676 14 : IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force)
677 : ! virial
678 14 : CALL get_qs_env(qs_env=qs_env, virial=virial)
679 14 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
680 :
681 14 : NULLIFY (qs_kind_set, atomic_kind_set, particle_set, distribution_1d, para_env, sab_cneo)
682 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
683 : particle_set=particle_set, local_particles=distribution_1d, &
684 14 : para_env=para_env, sab_cneo=sab_cneo)
685 : CALL build_core_cneo(rhoz_cneo_set, force, virial, calculate_forces, use_virial, nder, &
686 : qs_kind_set, atomic_kind_set, particle_set, distribution_1d, &
687 14 : sab_cneo, para_env)
688 : END IF
689 :
690 18707 : END SUBROUTINE cneo_core_matrices
691 :
692 : ! **************************************************************************************************
693 : !> \brief ...
694 : !> \param rhoz_cneo_set ...
695 : !> \param force ...
696 : !> \param virial ...
697 : !> \param calculate_forces ...
698 : !> \param use_virial ...
699 : !> \param nder ...
700 : !> \param qs_kind_set ...
701 : !> \param atomic_kind_set ...
702 : !> \param particle_set ...
703 : !> \param distribution_1d ...
704 : !> \param sab_cneo ...
705 : !> \param para_env ...
706 : ! **************************************************************************************************
707 14 : SUBROUTINE build_core_cneo(rhoz_cneo_set, force, virial, calculate_forces, use_virial, nder, &
708 : qs_kind_set, atomic_kind_set, particle_set, distribution_1d, &
709 : sab_cneo, para_env)
710 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
711 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
712 : TYPE(virial_type), POINTER :: virial
713 : LOGICAL, INTENT(IN) :: calculate_forces, use_virial
714 : INTEGER, INTENT(IN) :: nder
715 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
716 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
717 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
718 : TYPE(distribution_1d_type), POINTER :: distribution_1d
719 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
720 : POINTER :: sab_cneo
721 : TYPE(mp_para_env_type), POINTER :: para_env
722 :
723 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_cneo'
724 :
725 : INTEGER :: atom_a, handle, iat, iatom, ikind, iset, jatom, jkind, jset, ldai, ldsab, maxco, &
726 : maxl, maxnset, maxsgf, mepos, na_plus, nat, natom, nb_plus, ncoa, ncob, nij, nkind, nset, &
727 : nthread, sgfa, sgfb
728 14 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
729 14 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf, nsgf
730 14 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
731 : REAL(KIND=dp) :: alpha_c, core_charge, core_radius, dab, &
732 : f0, rab2, zeta_i, zeta_j
733 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ff
734 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: habd, work
735 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab, pab, verf, vnuc
736 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, force_i, rab
737 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_thread
738 : TYPE(neighbor_list_iterator_p_type), &
739 14 : DIMENSION(:), POINTER :: ap_iterator
740 : TYPE(gto_basis_set_type), POINTER :: basis_set
741 : TYPE(cneo_potential_type), POINTER :: cneo_potential, cneo_tmp
742 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: core, pmat, rpgf, sphi, zet
743 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius
744 28 : REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
745 :
746 28 : IF (calculate_forces) THEN
747 6 : CALL timeset(routineN//"_forces", handle)
748 : ELSE
749 8 : CALL timeset(routineN, handle)
750 : END IF
751 :
752 14 : nkind = SIZE(atomic_kind_set)
753 14 : natom = SIZE(particle_set)
754 :
755 166 : force_thread = 0.0_dp
756 14 : pv_thread = 0.0_dp
757 :
758 : ! re-initialize core matrices to zero, as later will use para_env%sum to broadcast
759 40 : DO ikind = 1, nkind
760 26 : NULLIFY (cneo_potential)
761 26 : CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
762 :
763 40 : IF (ASSOCIATED(cneo_potential)) THEN
764 14 : NULLIFY (atom_list)
765 14 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
766 40 : DO iat = 1, nat
767 26 : iatom = atom_list(iat)
768 14392 : rhoz_cneo_set(iatom)%core = 0.0_dp
769 : END DO
770 : END IF
771 : END DO
772 :
773 : CALL get_qs_kind_set(qs_kind_set, basis_type="NUC", &
774 14 : maxco=maxco, maxlgto=maxl, maxsgf=maxsgf, maxnset=maxnset)
775 14 : CALL init_orbital_pointers(maxl + nder + 1)
776 14 : ldsab = MAX(maxco, maxsgf)
777 14 : ldai = ncoset(maxl + nder + 1)
778 :
779 : nthread = 1
780 14 : !$ nthread = omp_get_max_threads()
781 :
782 14 : CALL neighbor_list_iterator_create(ap_iterator, sab_cneo, search=.TRUE., nthread=nthread)
783 :
784 : !$OMP PARALLEL &
785 : !$OMP DEFAULT (NONE) &
786 : !$OMP SHARED (rhoz_cneo_set, ap_iterator, distribution_1d, calculate_forces, use_virial, &
787 : !$OMP qs_kind_set, nthread, ncoset, nkind, iat, ldsab, maxnset, ldai, nder, maxl, &
788 : !$OMP maxco, para_env) &
789 : !$OMP PRIVATE (ikind, jkind, iatom, jatom, basis_set, first_sgf, lmax, lmin, npgf, nset, &
790 : !$OMP nsgf, rpgf, sphi, zet, set_radius, zeta_i, zeta_j, alpha_c, core_charge, &
791 : !$OMP core_radius, rab, rab2, dab, core, pmat, iset, ncoa, sgfa, jset, ncob, sgfb, &
792 : !$OMP work, pab, hab, na_plus, nb_plus, verf, vnuc, force_a, force_b, force_i, &
793 : !$OMP mepos, habd, f0, nij, ff, cneo_potential, cneo_tmp) &
794 14 : !$OMP REDUCTION (+ : pv_thread, force_thread )
795 :
796 : mepos = 0
797 : !$ mepos = omp_get_thread_num()
798 :
799 : ALLOCATE (hab(ldsab, ldsab, maxnset*(maxnset + 1)/2), work(ldsab, ldsab))
800 : ALLOCATE (verf(ldai, ldai, 2*maxl + nder + 1), vnuc(ldai, ldai, 2*maxl + nder + 1), ff(0:2*maxl + nder))
801 : IF (calculate_forces) THEN
802 : ALLOCATE (pab(maxco, maxco, maxnset*(maxnset + 1)/2))
803 : END IF
804 :
805 : DO ikind = 1, nkind
806 : NULLIFY (cneo_potential)
807 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type="NUC", &
808 : cneo_potential=cneo_potential, zeff=zeta_i)
809 : IF (ASSOCIATED(cneo_potential)) THEN
810 : CPASSERT(ASSOCIATED(basis_set))
811 : first_sgf => basis_set%first_sgf
812 : lmax => basis_set%lmax
813 : lmin => basis_set%lmin
814 : npgf => basis_set%npgf
815 : nset = basis_set%nset
816 : nsgf => basis_set%nsgf_set
817 : rpgf => basis_set%pgf_radius
818 : set_radius => basis_set%set_radius
819 : sphi => basis_set%sphi
820 : zet => basis_set%zet
821 :
822 : !$OMP DO SCHEDULE(GUIDED)
823 : DO iat = 1, distribution_1d%n_el(ikind)
824 : iatom = distribution_1d%list(ikind)%array(iat)
825 : core => rhoz_cneo_set(iatom)%core
826 : CPASSERT(ASSOCIATED(core))
827 : core = cneo_potential%kin ! copy kinetic matrix to core
828 : IF (calculate_forces) THEN
829 : CPASSERT(rhoz_cneo_set(iatom)%ready)
830 : pmat => rhoz_cneo_set(iatom)%pmat
831 : CPASSERT(ASSOCIATED(pmat))
832 : ! *** Decontract density matrix ***
833 : DO iset = 1, nset
834 : ncoa = npgf(iset)*ncoset(lmax(iset))
835 : sgfa = first_sgf(1, iset)
836 : DO jset = 1, iset
837 : ncob = npgf(jset)*ncoset(lmax(jset))
838 : sgfb = first_sgf(1, jset)
839 : nij = jset + (iset - 1)*iset/2
840 : work(1:ncoa, 1:nsgf(jset)) = MATMUL(sphi(1:ncoa, sgfa:sgfa + nsgf(iset) - 1), &
841 : pmat(sgfa:sgfa + nsgf(iset) - 1, sgfb:sgfb + nsgf(jset) - 1))
842 : pab(1:ncoa, 1:ncob, nij) = MATMUL(work(1:ncoa, 1:nsgf(jset)), &
843 : TRANSPOSE(sphi(1:ncob, sgfb:sgfb + nsgf(jset) - 1)))
844 : END DO
845 : END DO
846 : END IF
847 :
848 : hab = 0._dp
849 : DO jkind = 1, nkind
850 : NULLIFY (cneo_tmp)
851 : CALL get_qs_kind(qs_kind_set(jkind), cneo_potential=cneo_tmp)
852 : IF (.NOT. ASSOCIATED(cneo_tmp)) THEN
853 : CALL get_qs_kind(qs_kind_set(jkind), &
854 : alpha_core_charge=alpha_c, zeff=zeta_j, &
855 : ccore_charge=core_charge, core_charge_radius=core_radius)
856 : CALL nl_set_sub_iterator(ap_iterator, ikind, jkind, iatom, mepos=mepos)
857 :
858 : DO WHILE (nl_sub_iterate(ap_iterator, mepos=mepos) == 0)
859 : CALL get_iterator_info(ap_iterator, jatom=jatom, r=rab, mepos=mepos)
860 : rab2 = SUM(rab*rab)
861 : dab = SQRT(rab2)
862 : IF (MAXVAL(set_radius(:)) + core_radius < dab) CYCLE
863 : DO iset = 1, nset
864 : IF (set_radius(iset) + core_radius < dab) CYCLE
865 : ncoa = npgf(iset)*ncoset(lmax(iset))
866 : sgfa = first_sgf(1, iset)
867 : DO jset = 1, iset ! symmetric
868 : IF (set_radius(jset) + core_radius < dab) CYCLE
869 : ncob = npgf(jset)*ncoset(lmax(jset))
870 : sgfb = first_sgf(1, jset)
871 : nij = jset + (iset - 1)*iset/2
872 : IF (calculate_forces) THEN
873 : IF (jset == iset) THEN
874 : f0 = -zeta_i
875 : ELSE
876 : f0 = -2.0_dp*zeta_i
877 : END IF
878 : na_plus = npgf(iset)*ncoset(lmax(iset) + nder)
879 : nb_plus = npgf(jset)*ncoset(lmax(jset))
880 : ALLOCATE (habd(na_plus, nb_plus))
881 : habd = 0._dp
882 : CALL verfc( &
883 : lmax(iset) + nder, npgf(iset), zet(:, iset), rpgf(:, iset), lmin(iset), &
884 : lmax(jset), npgf(jset), zet(:, jset), rpgf(:, jset), lmin(jset), &
885 : alpha_c, core_radius, zeta_j, core_charge, &
886 : [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp, rab, rab2, rab2, &
887 : hab(:, :, nij), verf, vnuc, ff(0:), nder, habd)
888 :
889 : ! *** The derivatives w.r.t. atomic center b are ***
890 : ! *** calculated using the translational invariance ***
891 : ! *** of the first derivatives ***
892 : CALL verfc_force(habd, pab(:, :, nij), force_a, force_b, nder, &
893 : lmax(iset), lmin(iset), npgf(iset), zet(:, iset), &
894 : lmax(jset), lmin(jset), npgf(jset), zet(:, jset), &
895 : [0.0_dp, 0.0_dp, 0.0_dp])
896 :
897 : DEALLOCATE (habd)
898 : force_i = force_a + force_b
899 :
900 : force_thread(1, iatom) = force_thread(1, iatom) + f0*force_i(1)
901 : force_thread(2, iatom) = force_thread(2, iatom) + f0*force_i(2)
902 : force_thread(3, iatom) = force_thread(3, iatom) + f0*force_i(3)
903 :
904 : force_thread(1, jatom) = force_thread(1, jatom) - f0*force_i(1)
905 : force_thread(2, jatom) = force_thread(2, jatom) - f0*force_i(2)
906 : force_thread(3, jatom) = force_thread(3, jatom) - f0*force_i(3)
907 :
908 : IF (use_virial) THEN
909 : CALL virial_pair_force(pv_thread, f0, force_i, rab)
910 : END IF
911 : ELSE
912 : CALL verfc( &
913 : lmax(iset), npgf(iset), zet(:, iset), rpgf(:, iset), lmin(iset), &
914 : lmax(jset), npgf(jset), zet(:, jset), rpgf(:, jset), lmin(jset), &
915 : alpha_c, core_radius, zeta_j, core_charge, &
916 : [0.0_dp, 0.0_dp, 0.0_dp], 0.0_dp, rab, rab2, rab2, &
917 : hab(:, :, nij), verf, vnuc, ff(0:))
918 : END IF
919 : END DO
920 : END DO
921 : END DO
922 : END IF
923 : END DO
924 : ! *** Contract nuclear repulsion integrals
925 : DO iset = 1, nset
926 : ncoa = npgf(iset)*ncoset(lmax(iset))
927 : sgfa = first_sgf(1, iset)
928 : DO jset = 1, iset
929 : ncob = npgf(jset)*ncoset(lmax(jset))
930 : sgfb = first_sgf(1, jset)
931 : nij = jset + (iset - 1)*iset/2
932 : work(1:ncoa, 1:nsgf(jset)) = MATMUL(hab(1:ncoa, 1:ncob, nij), &
933 : sphi(1:ncob, sgfb:sgfb + nsgf(jset) - 1))
934 : core(sgfa:sgfa + nsgf(iset) - 1, sgfb:sgfb + nsgf(jset) - 1) = &
935 : core(sgfa:sgfa + nsgf(iset) - 1, sgfb:sgfb + nsgf(jset) - 1) - zeta_i* &
936 : MATMUL(TRANSPOSE(sphi(1:ncoa, sgfa:sgfa + nsgf(iset) - 1)), work(1:ncoa, 1:nsgf(jset)))
937 : ! symmetrize core matrix
938 : IF (iset /= jset) THEN
939 : core(sgfb:sgfb + nsgf(jset) - 1, sgfa:sgfa + nsgf(iset) - 1) = &
940 : TRANSPOSE(core(sgfa:sgfa + nsgf(iset) - 1, sgfb:sgfb + nsgf(jset) - 1))
941 : END IF
942 : END DO
943 : END DO
944 : END DO
945 : END IF
946 : END DO
947 :
948 : DEALLOCATE (hab, work, verf, vnuc, ff)
949 : IF (calculate_forces) THEN
950 : DEALLOCATE (pab)
951 : END IF
952 :
953 : !$OMP END PARALLEL
954 :
955 14 : CALL neighbor_list_iterator_release(ap_iterator)
956 :
957 14 : IF (calculate_forces) THEN
958 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, &
959 6 : kind_of=kind_of)
960 : !$OMP DO
961 : DO iatom = 1, natom
962 18 : atom_a = atom_of_kind(iatom)
963 18 : ikind = kind_of(iatom)
964 : force(ikind)%cneo_potential(:, atom_a) = force(ikind)%cneo_potential(:, atom_a) + &
965 72 : force_thread(:, iatom)
966 : END DO
967 : !$OMP END DO
968 : END IF
969 :
970 14 : IF (calculate_forces .AND. use_virial) THEN
971 0 : virial%pv_ppl = virial%pv_ppl + pv_thread
972 0 : virial%pv_virial = virial%pv_virial + pv_thread
973 : END IF
974 :
975 : ! broadcast core matrices
976 40 : DO ikind = 1, nkind
977 26 : NULLIFY (cneo_potential)
978 26 : CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
979 :
980 40 : IF (ASSOCIATED(cneo_potential)) THEN
981 14 : NULLIFY (atom_list)
982 14 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
983 40 : DO iat = 1, nat
984 26 : iatom = atom_list(iat)
985 28744 : CALL para_env%sum(rhoz_cneo_set(iatom)%core)
986 : END DO
987 : END IF
988 : END DO
989 :
990 14 : CALL timestop(handle)
991 :
992 28 : END SUBROUTINE build_core_cneo
993 :
994 : ! **************************************************************************************************
995 : !> \brief ...
996 : !> \param rho ...
997 : !> \param potential ...
998 : !> \param cg_list ...
999 : !> \param cg_n_list ...
1000 : !> \param nset ...
1001 : !> \param npgf ...
1002 : !> \param lmin ...
1003 : !> \param lmax ...
1004 : !> \param maxl ...
1005 : !> \param maxso ...
1006 : ! **************************************************************************************************
1007 46 : SUBROUTINE calculate_rhoz_cneo(rho, potential, cg_list, cg_n_list, nset, npgf, &
1008 : lmin, lmax, maxl, maxso)
1009 :
1010 : TYPE(rhoz_cneo_type), POINTER :: rho
1011 : TYPE(cneo_potential_type), POINTER :: potential
1012 : INTEGER, DIMENSION(:, :, :), INTENT(INOUT) :: cg_list
1013 : INTEGER, DIMENSION(:), INTENT(INOUT) :: cg_n_list
1014 : INTEGER, INTENT(IN) :: nset
1015 : INTEGER, DIMENSION(:), POINTER :: npgf, lmin, lmax
1016 : INTEGER, INTENT(IN) :: maxl, maxso
1017 :
1018 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rhoz_cneo'
1019 :
1020 : INTEGER :: handle, i, i1, i2, icg, ipgf1, ipgf2, iset1, iset2, iso, iso1, iso1_first, &
1021 : iso1_last, iso2, iso2_first, iso2_last, iter, j, l, l1, l2, l_iso, lmax_expansion, m1s, &
1022 : m2s, max_iso_not0, max_iso_not0_local, max_iter, max_s_harm, n1s, n2s, nne, npgf2, npsgf, &
1023 : nsgf, nsotot, size1, size2
1024 46 : INTEGER, DIMENSION(:), POINTER :: n2oindex, o2nindex
1025 : REAL(KIND=dp) :: det, df_norm, factor, g0, g0p, g1, step, &
1026 : zeff
1027 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ener
1028 46 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: CPCH_sphere, CPCS_sphere, work
1029 : REAL(KIND=dp), DIMENSION(3) :: df, f_tmp, r, r_tmp
1030 : REAL(KIND=dp), DIMENSION(3, 3) :: jac, jac_inv
1031 : REAL(KIND=dp), DIMENSION(:), POINTER :: f
1032 46 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: core, cpc_h, cpc_s, fmat, int_local_h, &
1033 46 : int_local_s, my_gcc_h, my_gcc_s, pmat, rho_rad_h, rho_rad_s, utrans, vmat, vrho_rad_h, &
1034 46 : vrho_rad_s, wfn
1035 46 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: distance, gg, my_CG
1036 46 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: vgg
1037 : TYPE(harmonics_atom_type), POINTER :: harmonics
1038 :
1039 0 : CPASSERT(ASSOCIATED(rho))
1040 46 : CPASSERT(ASSOCIATED(potential))
1041 :
1042 46 : CALL timeset(routineN, handle)
1043 :
1044 : ! convert ga_Vlocal_gb to compressed form V_Hartree
1045 : ! use fmat to store V_Hartree
1046 46 : NULLIFY (utrans, my_gcc_h, my_gcc_s, distance, n2oindex, o2nindex)
1047 : CALL get_cneo_potential(potential, zeff=zeff, nsgf=nsgf, nne=nne, npsgf=npsgf, &
1048 : nsotot=nsotot, my_gcc_h=my_gcc_h, my_gcc_s=my_gcc_s, &
1049 : utrans=utrans, distance=distance, n2oindex=n2oindex, &
1050 46 : o2nindex=o2nindex)
1051 46 : fmat => rho%fmat
1052 46 : int_local_h => rho%ga_Vlocal_gb_h
1053 46 : int_local_s => rho%ga_Vlocal_gb_s
1054 184 : ALLOCATE (work(nsotot, nsgf))
1055 : CALL dgemm("N", "N", nsotot, nsgf, nsotot, 1.0_dp, int_local_h, nsotot, my_gcc_h, &
1056 46 : nsotot, 0.0_dp, work, nsotot)
1057 : CALL dgemm("T", "N", nsgf, nsgf, nsotot, 1.0_dp, my_gcc_h, nsotot, work, &
1058 46 : nsotot, 0.0_dp, fmat, nsgf)
1059 : CALL dgemm("N", "N", nsotot, nsgf, nsotot, 1.0_dp, int_local_s, nsotot, my_gcc_s, &
1060 46 : nsotot, 0.0_dp, work, nsotot)
1061 : CALL dgemm("T", "N", nsgf, nsgf, nsotot, -1.0_dp, my_gcc_s, nsotot, work, &
1062 46 : nsotot, 1.0_dp, fmat, nsgf)
1063 : ! add the soft basis FFT grid part
1064 46 : vmat => rho%vmat
1065 50830 : fmat = fmat + vmat
1066 :
1067 46 : core => rho%core
1068 46 : wfn => rho%wfn
1069 46 : pmat => rho%pmat
1070 46 : f => rho%f
1071 138 : ALLOCATE (ener(nne))
1072 : ! build the fock matrix: F = T + V_core + V_Hartree
1073 50830 : fmat = fmat + core
1074 : ! conduct the constrained optimization with F + f*x
1075 : ! initial guess of f is taken from the result of last iteration
1076 46 : CALL atom_solve_cneo(fmat, f, utrans, wfn, ener, pmat, r, distance, nsgf, nne)
1077 : ! test if zero initial guess is better
1078 322 : IF (NORM2(r) > 1.e-12_dp .AND. DOT_PRODUCT(f, f) /= 0.0_dp) THEN
1079 : CALL atom_solve_cneo(fmat, [0.0_dp, 0.0_dp, 0.0_dp], utrans, wfn, &
1080 39 : ener, pmat, r_tmp, distance, nsgf, nne)
1081 273 : IF (DOT_PRODUCT(r_tmp, r_tmp) < DOT_PRODUCT(r, r)) THEN
1082 24 : f = 0.0_dp
1083 6 : r = r_tmp
1084 : END IF
1085 : END IF
1086 46 : max_iter = 20
1087 46 : iter = 0
1088 : ! using Newton's method to solve for f
1089 884 : DO WHILE (NORM2(r) > 1.e-12_dp)
1090 140 : iter = iter + 1
1091 : ! construct numerical Jacobian with one-side finite difference
1092 560 : DO i = 1, 3
1093 1680 : f_tmp = f
1094 420 : f_tmp(i) = f(i) + SIGN(1.e-4_dp, r(i)) ! forward or backward based on the sign of r
1095 420 : CALL atom_solve_cneo(fmat, f_tmp, utrans, wfn, ener, pmat, r_tmp, distance, nsgf, nne)
1096 1820 : DO j = 1, 3
1097 1680 : jac(j, i) = (r_tmp(j) - r(j))*SIGN(1.e4_dp, r(i))
1098 : END DO
1099 : END DO
1100 140 : CALL invert_matrix_3x3(jac, jac_inv, det)
1101 140 : IF (ABS(det) < 1.0E-8_dp) THEN
1102 : CALL cp_warn(__LOCATION__, "Determinant of the CNEO position Jacobian is small! "// &
1103 0 : TRIM(cp_to_string(det))//" Trying central difference.")
1104 : ! construct numerical Jacobian with central finite difference
1105 0 : DO i = 1, 3
1106 0 : f_tmp = f
1107 0 : f_tmp(i) = f(i) - SIGN(1.e-4_dp, r(i))
1108 0 : CALL atom_solve_cneo(fmat, f_tmp, utrans, wfn, ener, pmat, r_tmp, distance, nsgf, nne)
1109 0 : DO j = 1, 3
1110 : jac(j, i) = (jac(j, i)*SIGN(1.e-4_dp, r(i)) + r(j) - r_tmp(j)) &
1111 0 : /SIGN(2.e-4_dp, r(i))
1112 : END DO
1113 : END DO
1114 0 : CALL invert_matrix_3x3(jac, jac_inv, det)
1115 0 : IF (ABS(det) < 1.0E-8_dp) THEN
1116 : CALL cp_warn(__LOCATION__, "Determinant of the CNEO position Jacobian is small! "// &
1117 0 : "(Central difference) "//TRIM(cp_to_string(det))//" Using pseudoinverse.")
1118 : END IF
1119 0 : CALL invert_matrix_3x3(jac, jac_inv, det, try_svd=.TRUE.)
1120 : END IF
1121 2380 : df = -RESHAPE(MATMUL(jac_inv, RESHAPE(r, [3, 1])), [3])
1122 560 : df_norm = NORM2(df)
1123 560 : f_tmp = f
1124 140 : r_tmp = r
1125 560 : g0 = NORM2(r_tmp)
1126 560 : f = f_tmp + df
1127 140 : CALL atom_solve_cneo(fmat, f, utrans, wfn, ener, pmat, r, distance, nsgf, nne)
1128 560 : g1 = NORM2(r)
1129 140 : step = 1.0_dp
1130 140 : DO WHILE (g1 >= g0)
1131 : ! line search
1132 0 : IF (step < 0.0101_dp) THEN
1133 0 : CPWARN("CNEO nuclear position constraint solver line search failure.")
1134 0 : EXIT
1135 : END IF
1136 0 : g0p = -g0/(step*df_norm)
1137 0 : step = step*MAX(-g0p/(2.0_dp*(g1 - g0 - g0p)), 0.1_dp)
1138 0 : f = f_tmp + step*df
1139 0 : CALL atom_solve_cneo(fmat, f, utrans, wfn, ener, pmat, r, distance, nsgf, nne)
1140 140 : g1 = NORM2(r)
1141 : END DO
1142 280 : IF (iter >= max_iter) THEN
1143 : CALL cp_warn(__LOCATION__, "CNEO nuclear position constraint solver failed to "// &
1144 : "converge in "//TRIM(cp_to_string(max_iter))//" steps. "// &
1145 : "Nuclear position error (x,y,z): "//TRIM(cp_to_string(r(1)))// &
1146 : ", "//TRIM(cp_to_string(r(2)))//", "//TRIM(cp_to_string(r(3)))// &
1147 0 : ". This does not hurt as long as it is not the final SCF iteration.")
1148 0 : EXIT
1149 : END IF
1150 : END DO
1151 46 : DEALLOCATE (ener)
1152 46 : rho%e_core = trace_r_AxB(core, nsgf, pmat, nsgf, nsgf, nsgf)
1153 :
1154 : ! decontract the density matrix
1155 : ! first use ga_Vlocal_gb to store the decompressed form
1156 : CALL dgemm("N", "N", nsotot, nsgf, nsgf, 1.0_dp, my_gcc_h, nsotot, pmat, nsgf, &
1157 46 : 0.0_dp, work, nsotot)
1158 : CALL dgemm("N", "T", nsotot, nsotot, nsgf, 1.0_dp, work, nsotot, my_gcc_h, nsotot, &
1159 46 : 0.0_dp, int_local_h, nsotot)
1160 : CALL dgemm("N", "N", nsotot, nsgf, nsgf, 1.0_dp, my_gcc_s, nsotot, pmat, nsgf, &
1161 46 : 0.0_dp, work, nsotot)
1162 : CALL dgemm("N", "T", nsotot, nsotot, nsgf, 1.0_dp, work, nsotot, my_gcc_s, nsotot, &
1163 46 : 0.0_dp, int_local_s, nsotot)
1164 46 : DEALLOCATE (work)
1165 : ! compress the density matrix
1166 46 : cpc_h => rho%cpc_h
1167 46 : cpc_s => rho%cpc_s
1168 46 : CALL cneo_gather(int_local_h, cpc_h, npsgf, n2oindex)
1169 46 : CALL cneo_gather(int_local_s, cpc_s, npsgf, n2oindex)
1170 : ! restore ga_Vlocal_gb to zeros
1171 305578 : int_local_h = 0.0_dp
1172 305578 : int_local_s = 0.0_dp
1173 :
1174 : ! construct the nuclear density and its Hartree potential
1175 : ! rho_rad_h and vrho_rad_h should contain the -Zeff factor
1176 : ! mostly copied from qs_rho_atom_methods::calculate_rho_atom
1177 46 : NULLIFY (harmonics, gg, vgg)
1178 46 : CALL get_cneo_potential(potential, harmonics=harmonics, gg=gg, vgg=vgg)
1179 46 : rho_rad_h => rho%rho_rad_h
1180 46 : rho_rad_s => rho%rho_rad_s
1181 58696 : rho_rad_h = 0.0_dp
1182 58696 : rho_rad_s = 0.0_dp
1183 46 : vrho_rad_h => rho%vrho_rad_h
1184 46 : vrho_rad_s => rho%vrho_rad_s
1185 58696 : vrho_rad_h = 0.0_dp
1186 58696 : vrho_rad_s = 0.0_dp
1187 46 : my_CG => harmonics%my_CG
1188 46 : max_iso_not0 = harmonics%max_iso_not0
1189 46 : max_s_harm = harmonics%max_s_harm
1190 46 : lmax_expansion = indso(1, max_iso_not0)
1191 :
1192 184 : ALLOCATE (CPCH_sphere(nsoset(maxl), nsoset(maxl)))
1193 138 : ALLOCATE (CPCS_sphere(nsoset(maxl), nsoset(maxl)))
1194 46 : j = 0
1195 46 : m1s = 0
1196 460 : DO iset1 = 1, nset
1197 414 : m2s = 0
1198 414 : n1s = nsoset(lmax(iset1))
1199 2484 : DO iset2 = 1, iset1
1200 :
1201 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
1202 2070 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
1203 2070 : CPASSERT(max_iso_not0_local <= max_iso_not0)
1204 :
1205 2070 : n2s = nsoset(lmax(iset2))
1206 4140 : DO ipgf1 = 1, npgf(iset1)
1207 2070 : iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s
1208 2070 : iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s
1209 2070 : size1 = iso1_last - iso1_first + 1
1210 2070 : iso1_first = o2nindex(iso1_first)
1211 2070 : iso1_last = o2nindex(iso1_last)
1212 2070 : i1 = iso1_last - iso1_first + 1
1213 2070 : CPASSERT(size1 == i1)
1214 2070 : i1 = nsoset(lmin(iset1) - 1) + 1
1215 :
1216 2070 : IF (iset2 == iset1) THEN
1217 : npgf2 = ipgf1
1218 : ELSE
1219 1656 : npgf2 = npgf(iset2)
1220 : END IF
1221 6210 : DO ipgf2 = 1, npgf2
1222 2070 : j = j + 1
1223 2070 : iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s
1224 2070 : iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s
1225 2070 : size2 = iso2_last - iso2_first + 1
1226 2070 : iso2_first = o2nindex(iso2_first)
1227 2070 : iso2_last = o2nindex(iso2_last)
1228 2070 : i2 = iso2_last - iso2_first + 1
1229 2070 : CPASSERT(size2 == i2)
1230 2070 : i2 = nsoset(lmin(iset2) - 1) + 1
1231 :
1232 2070 : IF (iset2 == iset1 .AND. ipgf2 == ipgf1) THEN
1233 414 : factor = -zeff
1234 : ELSE
1235 1656 : factor = -2.0_dp*zeff
1236 : END IF
1237 :
1238 2070 : CPCH_sphere = 0.0_dp
1239 2070 : CPCS_sphere = 0.0_dp
1240 19826 : CPCH_sphere(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = cpc_h(iso1_first:iso1_last, iso2_first:iso2_last)
1241 19826 : CPCS_sphere(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = cpc_s(iso1_first:iso1_last, iso2_first:iso2_last)
1242 20470 : DO iso = 1, max_iso_not0_local
1243 16330 : l_iso = indso(1, iso)
1244 44298 : DO icg = 1, cg_n_list(iso)
1245 25898 : iso1 = cg_list(1, icg, iso)
1246 25898 : iso2 = cg_list(2, icg, iso)
1247 :
1248 25898 : l1 = indso(1, iso1)
1249 25898 : l2 = indso(1, iso2)
1250 :
1251 25898 : l = indso(1, iso1) + indso(1, iso2)
1252 25898 : CPASSERT(l <= lmax_expansion)
1253 :
1254 : rho_rad_h(:, iso) = rho_rad_h(:, iso) + gg(:, l, j)* &
1255 2615698 : CPCH_sphere(iso1, iso2)*my_CG(iso1, iso2, iso)*factor
1256 :
1257 : rho_rad_s(:, iso) = rho_rad_s(:, iso) + gg(:, l, j)* &
1258 2615698 : CPCS_sphere(iso1, iso2)*my_CG(iso1, iso2, iso)*factor
1259 :
1260 : vrho_rad_h(:, iso) = vrho_rad_h(:, iso) + vgg(:, l, l_iso, j)* &
1261 2615698 : CPCH_sphere(iso1, iso2)*my_CG(iso1, iso2, iso)*factor
1262 :
1263 : vrho_rad_s(:, iso) = vrho_rad_s(:, iso) + vgg(:, l, l_iso, j)* &
1264 2632028 : CPCS_sphere(iso1, iso2)*my_CG(iso1, iso2, iso)*factor
1265 : END DO ! icg
1266 : END DO ! iso
1267 : END DO ! ipgf2
1268 : END DO ! ipgf1
1269 2484 : m2s = m2s + maxso
1270 : END DO ! iset2
1271 460 : m1s = m1s + maxso
1272 : END DO ! iset1
1273 46 : DEALLOCATE (CPCH_sphere, CPCS_sphere)
1274 :
1275 46 : CALL timestop(handle)
1276 :
1277 92 : END SUBROUTINE calculate_rhoz_cneo
1278 :
1279 : ! **************************************************************************************************
1280 : !> \brief Mostly copied from hartree_local_methods::Vh_1c_atom_integrals
1281 : !> \param rhoz_cneo ...
1282 : !> \param zeff ...
1283 : !> \param aVh1b_hh ...
1284 : !> \param aVh1b_ss ...
1285 : !> \param aVh1b_00 ...
1286 : !> \param Vh1_h ...
1287 : !> \param Vh1_s ...
1288 : !> \param max_iso_not0_elec ...
1289 : !> \param max_iso_not0_nuc ...
1290 : !> \param max_s_harm ...
1291 : !> \param llmax ...
1292 : !> \param cg_list ...
1293 : !> \param cg_n_list ...
1294 : !> \param nset ...
1295 : !> \param npgf ...
1296 : !> \param lmin ...
1297 : !> \param lmax ...
1298 : !> \param nsotot ...
1299 : !> \param maxso ...
1300 : !> \param nchan_0 ...
1301 : !> \param gsph ...
1302 : !> \param g0_h_w ...
1303 : !> \param my_CG ...
1304 : !> \param Qlm_gg ...
1305 : ! **************************************************************************************************
1306 46 : SUBROUTINE Vh_1c_nuc_integrals(rhoz_cneo, zeff, &
1307 46 : aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, &
1308 : max_iso_not0_elec, max_iso_not0_nuc, &
1309 46 : max_s_harm, llmax, cg_list, cg_n_list, &
1310 : nset, npgf, lmin, lmax, nsotot, maxso, nchan_0, gsph, &
1311 46 : g0_h_w, my_CG, Qlm_gg)
1312 :
1313 : TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
1314 : REAL(KIND=dp), INTENT(IN) :: zeff
1315 : REAL(KIND=dp), DIMENSION(:, :) :: aVh1b_hh, aVh1b_ss, aVh1b_00
1316 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s
1317 : INTEGER, INTENT(IN) :: max_iso_not0_elec, max_iso_not0_nuc, &
1318 : max_s_harm, llmax
1319 : INTEGER, DIMENSION(:, :, :) :: cg_list
1320 : INTEGER, DIMENSION(:) :: cg_n_list
1321 : INTEGER, INTENT(IN) :: nset
1322 : INTEGER, DIMENSION(:), POINTER :: npgf, lmin, lmax
1323 : INTEGER, INTENT(IN) :: nsotot, maxso, nchan_0
1324 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gsph
1325 : REAL(KIND=dp), DIMENSION(:, 0:) :: g0_h_w
1326 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg
1327 :
1328 : INTEGER :: icg, ipgf1, ipgf2, ir, is1, is2, iset1, &
1329 : iset2, iso, iso1, iso2, l_ang, m1, m2, &
1330 : max_iso_not0_local, n1, n2, nr
1331 : REAL(KIND=dp) :: gVg_0, gVg_h, gVg_s
1332 :
1333 : ! Calculate the integrals of the potential with 2 primitives
1334 305578 : aVh1b_hh = 0.0_dp
1335 305578 : aVh1b_ss = 0.0_dp
1336 305578 : aVh1b_00 = 0.0_dp
1337 :
1338 46 : nr = SIZE(gsph, 1)
1339 :
1340 46 : m1 = 0
1341 460 : DO iset1 = 1, nset
1342 414 : n1 = nsoset(lmax(iset1))
1343 414 : m2 = 0
1344 4140 : DO iset2 = 1, nset
1345 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
1346 3726 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
1347 :
1348 3726 : n2 = nsoset(lmax(iset2))
1349 7452 : DO ipgf1 = 1, npgf(iset1)
1350 11178 : DO ipgf2 = 1, npgf(iset2)
1351 37260 : DO iso = 1, MIN(max_iso_not0_elec, max_iso_not0_nuc)
1352 62376 : DO icg = 1, cg_n_list(iso)
1353 25116 : is1 = cg_list(1, icg, iso)
1354 25116 : is2 = cg_list(2, icg, iso)
1355 :
1356 25116 : iso1 = is1 + n1*(ipgf1 - 1) + m1
1357 25116 : iso2 = is2 + n2*(ipgf2 - 1) + m2
1358 25116 : gVg_h = 0.0_dp
1359 25116 : gVg_s = 0.0_dp
1360 :
1361 1280916 : DO ir = 1, nr
1362 1255800 : gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso)
1363 1280916 : gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso)
1364 : END DO ! ir
1365 :
1366 25116 : aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso)
1367 58650 : aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso)
1368 :
1369 : END DO !icg
1370 : END DO ! iso
1371 63342 : DO iso = max_iso_not0_elec + 1, max_iso_not0_nuc
1372 81742 : DO icg = 1, cg_n_list(iso)
1373 18400 : is1 = cg_list(1, icg, iso)
1374 18400 : is2 = cg_list(2, icg, iso)
1375 :
1376 18400 : iso1 = is1 + n1*(ipgf1 - 1) + m1
1377 18400 : iso2 = is2 + n2*(ipgf2 - 1) + m2
1378 18400 : gVg_s = 0.0_dp
1379 :
1380 938400 : DO ir = 1, nr
1381 938400 : gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso)
1382 : END DO ! ir
1383 :
1384 78016 : aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso)
1385 :
1386 : END DO !icg
1387 : END DO ! iso
1388 40986 : DO iso = 1, MIN(nchan_0, max_iso_not0_nuc)
1389 33534 : l_ang = indso(1, iso)
1390 1710234 : gVg_0 = SUM(Vh1_s(:, iso)*g0_h_w(:, l_ang))
1391 62376 : DO icg = 1, cg_n_list(iso)
1392 25116 : is1 = cg_list(1, icg, iso)
1393 25116 : is2 = cg_list(2, icg, iso)
1394 :
1395 25116 : iso1 = is1 + n1*(ipgf1 - 1) + m1
1396 25116 : iso2 = is2 + n2*(ipgf2 - 1) + m2
1397 :
1398 58650 : aVh1b_00(iso1, iso2) = aVh1b_00(iso1, iso2) + gVg_0*Qlm_gg(iso1, iso2, iso)
1399 :
1400 : END DO !icg
1401 : END DO ! iso
1402 : END DO ! ipgf2
1403 : END DO ! ipgf1
1404 4140 : m2 = m2 + maxso
1405 : END DO ! iset2
1406 460 : m1 = m1 + maxso
1407 : END DO !iset1
1408 :
1409 46 : CALL daxpy(nsotot*nsotot, -zeff, aVh1b_hh, 1, rhoz_cneo%ga_Vlocal_gb_h, 1)
1410 46 : CALL daxpy(nsotot*nsotot, -zeff, aVh1b_ss, 1, rhoz_cneo%ga_Vlocal_gb_s, 1)
1411 46 : CALL daxpy(nsotot*nsotot, zeff, aVh1b_00, 1, rhoz_cneo%ga_Vlocal_gb_h, 1)
1412 46 : CALL daxpy(nsotot*nsotot, zeff, aVh1b_00, 1, rhoz_cneo%ga_Vlocal_gb_s, 1)
1413 :
1414 46 : END SUBROUTINE Vh_1c_nuc_integrals
1415 :
1416 : ! **************************************************************************************************
1417 : !> \brief Analytical inversion of a 3x3 matrix
1418 : !> \param matrix ...
1419 : !> \param inv_matrix ...
1420 : !> \param det ...
1421 : !> \param try_svd ...
1422 : ! **************************************************************************************************
1423 140 : SUBROUTINE invert_matrix_3x3(matrix, inv_matrix, det, try_svd)
1424 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: matrix
1425 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: inv_matrix
1426 : REAL(KIND=dp), INTENT(OUT) :: det
1427 : LOGICAL, INTENT(IN), OPTIONAL :: try_svd
1428 :
1429 : LOGICAL :: my_try_svd
1430 :
1431 140 : my_try_svd = .FALSE.
1432 140 : IF (PRESENT(try_svd)) my_try_svd = try_svd
1433 :
1434 : det = matrix(1, 1)*(matrix(2, 2)*matrix(3, 3) - matrix(2, 3)*matrix(3, 2)) &
1435 : - matrix(1, 2)*(matrix(2, 1)*matrix(3, 3) - matrix(2, 3)*matrix(3, 1)) &
1436 140 : + matrix(1, 3)*(matrix(2, 1)*matrix(3, 2) - matrix(2, 2)*matrix(3, 1))
1437 140 : IF (ABS(det) < 1.0E-8_dp) THEN
1438 0 : IF (my_try_svd) THEN
1439 : ! pseudo inverse using SVD
1440 0 : CALL get_pseudo_inverse_svd(matrix, inv_matrix, 1.0E-6_dp, det)
1441 : ELSE
1442 0 : inv_matrix = 0.0_dp
1443 : END IF
1444 : ELSE
1445 140 : inv_matrix(1, 1) = matrix(2, 2)*matrix(3, 3) - matrix(2, 3)*matrix(3, 2)
1446 140 : inv_matrix(1, 2) = matrix(1, 3)*matrix(3, 2) - matrix(1, 2)*matrix(3, 3)
1447 140 : inv_matrix(1, 3) = matrix(1, 2)*matrix(2, 3) - matrix(1, 3)*matrix(2, 2)
1448 140 : inv_matrix(2, 1) = matrix(2, 3)*matrix(3, 1) - matrix(2, 1)*matrix(3, 3)
1449 140 : inv_matrix(2, 2) = matrix(1, 1)*matrix(3, 3) - matrix(1, 3)*matrix(3, 1)
1450 140 : inv_matrix(2, 3) = matrix(1, 3)*matrix(2, 1) - matrix(1, 1)*matrix(2, 3)
1451 140 : inv_matrix(3, 1) = matrix(2, 1)*matrix(3, 2) - matrix(2, 2)*matrix(3, 1)
1452 140 : inv_matrix(3, 2) = matrix(1, 2)*matrix(3, 1) - matrix(1, 1)*matrix(3, 2)
1453 140 : inv_matrix(3, 3) = matrix(1, 1)*matrix(2, 2) - matrix(1, 2)*matrix(2, 1)
1454 1820 : inv_matrix = inv_matrix/det
1455 : END IF
1456 140 : END SUBROUTINE invert_matrix_3x3
1457 :
1458 : END MODULE qs_cneo_methods
|
