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