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 Utility functions for 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_utils
16 : USE ao_util, ONLY: trace_r_AxB
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE kinds, ONLY: dp
20 : USE memory_utilities, ONLY: reallocate
21 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
22 : harmonics_atom_type
23 : #include "./base/base_uses.f90"
24 :
25 : IMPLICIT NONE
26 :
27 : PRIVATE
28 :
29 : ! *** Global parameters ***
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cneo_utils'
31 :
32 : ! *** Public subroutines ***
33 : PUBLIC :: atom_solve_cneo, cneo_gather, cneo_scatter, create_harmonics_atom_cneo, &
34 : get_maxl_CG_cneo
35 :
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief Mostly copied from qs_harmonics_atom::create_harmonics_atom
40 : !> \param harmonics ...
41 : !> \param my_CG ...
42 : !> \param llmax ...
43 : !> \param maxs ...
44 : !> \param max_s_harm ...
45 : ! **************************************************************************************************
46 8 : SUBROUTINE create_harmonics_atom_cneo(harmonics, my_CG, llmax, maxs, max_s_harm)
47 :
48 : TYPE(harmonics_atom_type), POINTER :: harmonics
49 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_CG
50 : INTEGER, INTENT(IN) :: llmax, maxs, max_s_harm
51 :
52 : INTEGER :: i, is
53 :
54 8 : CPASSERT(ASSOCIATED(harmonics))
55 :
56 8 : harmonics%max_s_harm = max_s_harm
57 8 : harmonics%llmax = llmax
58 :
59 8 : NULLIFY (harmonics%my_CG, harmonics%my_CG_dxyz, harmonics%my_CG_dxyz_asym)
60 8 : CALL reallocate(harmonics%my_CG, 1, maxs, 1, maxs, 1, max_s_harm)
61 :
62 208 : DO i = 1, max_s_harm
63 2008 : DO is = 1, maxs
64 36200 : harmonics%my_CG(1:maxs, is, i) = my_CG(1:maxs, is, i)
65 : END DO
66 : END DO
67 :
68 8 : END SUBROUTINE create_harmonics_atom_cneo
69 :
70 : ! **************************************************************************************************
71 : !> \brief Mostly copied from qs_harmonics_atom::get_maxl_CG
72 : !> \param harmonics ...
73 : !> \param orb_basis ...
74 : !> \param llmax ...
75 : !> \param max_s_harm ...
76 : ! **************************************************************************************************
77 16 : SUBROUTINE get_maxl_CG_cneo(harmonics, orb_basis, llmax, max_s_harm)
78 :
79 : TYPE(harmonics_atom_type), POINTER :: harmonics
80 : TYPE(gto_basis_set_type), POINTER :: orb_basis
81 : INTEGER, INTENT(IN) :: llmax, max_s_harm
82 :
83 : INTEGER :: is1, is2, itmp, max_iso_not0, nset
84 8 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin
85 :
86 0 : CPASSERT(ASSOCIATED(harmonics))
87 :
88 8 : CALL get_gto_basis_set(gto_basis_set=orb_basis, lmax=lmax, lmin=lmin, nset=nset)
89 :
90 : ! *** Assign indices for the non null CG coefficients ***
91 8 : max_iso_not0 = 0
92 80 : DO is1 = 1, nset
93 728 : DO is2 = 1, nset
94 : CALL get_none0_cg_list(harmonics%my_CG, &
95 : lmin(is1), lmax(is1), lmin(is2), lmax(is2), &
96 648 : max_s_harm, llmax, max_iso_not0=itmp)
97 720 : max_iso_not0 = MAX(max_iso_not0, itmp)
98 : END DO ! is2
99 : END DO ! is1
100 8 : harmonics%max_iso_not0 = max_iso_not0
101 :
102 8 : END SUBROUTINE get_maxl_CG_cneo
103 :
104 : ! **************************************************************************************************
105 : !> \brief Mostly copied from atom_utils::atom_solve
106 : !> \param hmat ...
107 : !> \param f ...
108 : !> \param umat ...
109 : !> \param orb ...
110 : !> \param ener ...
111 : !> \param pmat ...
112 : !> \param r ...
113 : !> \param dist ...
114 : !> \param nb ...
115 : !> \param nv ...
116 : ! **************************************************************************************************
117 645 : SUBROUTINE atom_solve_cneo(hmat, f, umat, orb, ener, pmat, r, dist, nb, nv)
118 :
119 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hmat
120 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: f
121 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: umat
122 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: orb
123 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: ener
124 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pmat
125 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: r
126 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dist
127 : INTEGER, INTENT(IN) :: nb, nv
128 :
129 : INTEGER :: info, lwork, m, n
130 645 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
131 645 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, h_fx
132 :
133 645 : CPASSERT(nb >= nv)
134 :
135 356685 : orb = 0._dp
136 645 : n = nb
137 645 : m = nv
138 645 : IF (n > 0 .AND. m > 0) THEN
139 645 : lwork = MAX(n*n, n + 100)
140 7095 : ALLOCATE (a(m, m), b(n, m), w(m), work(lwork))
141 2580 : IF (DOT_PRODUCT(f, f) /= 0.0_dp) THEN
142 2396 : ALLOCATE (h_fx(n, n))
143 : h_fx(1:n, 1:n) = hmat(1:n, 1:n) + f(1)*dist(1:n, 1:n, 1) + &
144 331247 : f(2)*dist(1:n, 1:n, 2) + f(3)*dist(1:n, 1:n, 3)
145 599 : CALL dgemm("N", "N", n, m, n, 1.0_dp, h_fx, n, umat, n, 0.0_dp, b, n)
146 599 : DEALLOCATE (h_fx)
147 : ELSE
148 46 : CALL dgemm("N", "N", n, m, n, 1.0_dp, hmat, n, umat, n, 0.0_dp, b, n)
149 : END IF
150 645 : CALL dgemm("T", "N", m, m, n, 1.0_dp, umat, n, b, n, 0.0_dp, a, m)
151 645 : CALL dsyev("V", "U", m, a, m, w, work, lwork, info)
152 645 : CALL dgemm("N", "N", n, m, m, 1.0_dp, umat, n, a, m, 0.0_dp, b, n)
153 :
154 645 : m = MIN(m, SIZE(orb, 2))
155 356685 : orb(1:n, 1:m) = b(1:n, 1:m)
156 15480 : ener(1:m) = w(1:m)
157 :
158 645 : DEALLOCATE (a, b, w, work)
159 :
160 : ! calculate the density matrix using the orbital with the lowest orbital energy
161 356685 : pmat = 0.0_dp
162 645 : CALL dger(n, n, 1.0_dp, orb(:, 1), 1, orb(:, 1), 1, pmat, n)
163 : ! calculate the expectation position (basis center as the origin)
164 : r = [trace_r_AxB(dist(1:n, 1:n, 1), n, pmat, n, n, n), &
165 : trace_r_AxB(dist(1:n, 1:n, 2), n, pmat, n, n, n), &
166 3225 : trace_r_AxB(dist(1:n, 1:n, 3), n, pmat, n, n, n)]
167 : END IF
168 :
169 645 : END SUBROUTINE atom_solve_cneo
170 :
171 : ! **************************************************************************************************
172 : !> \brief Mostly copied from qs_oce_methods::prj_gather
173 : !> \param ain ...
174 : !> \param aout ...
175 : !> \param nbas ...
176 : !> \param n2oindex ...
177 : ! **************************************************************************************************
178 92 : SUBROUTINE cneo_gather(ain, aout, nbas, n2oindex)
179 :
180 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ain
181 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aout
182 : INTEGER, INTENT(IN) :: nbas
183 : INTEGER, DIMENSION(:), POINTER :: n2oindex
184 :
185 : INTEGER :: i, ip, j, jp
186 :
187 2208 : DO i = 1, nbas
188 2116 : ip = n2oindex(i)
189 50876 : DO j = 1, nbas
190 48668 : jp = n2oindex(j)
191 50784 : aout(j, i) = ain(jp, ip)
192 : END DO
193 : END DO
194 :
195 92 : END SUBROUTINE cneo_gather
196 :
197 : ! **************************************************************************************************
198 : !> \brief Mostly copied from qs_oce_methods::prj_scatter
199 : !> \param ain ...
200 : !> \param aout ...
201 : !> \param nbas ...
202 : !> \param n2oindex ...
203 : ! **************************************************************************************************
204 156 : SUBROUTINE cneo_scatter(ain, aout, nbas, n2oindex)
205 :
206 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ain
207 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aout
208 : INTEGER, INTENT(IN) :: nbas
209 : INTEGER, DIMENSION(:), POINTER :: n2oindex
210 :
211 : INTEGER :: i, ip, j, jp
212 :
213 3744 : DO i = 1, nbas
214 3588 : ip = n2oindex(i)
215 86268 : DO j = 1, nbas
216 82524 : jp = n2oindex(j)
217 86112 : aout(jp, ip) = aout(jp, ip) + ain(j, i)
218 : END DO
219 : END DO
220 :
221 156 : END SUBROUTINE cneo_scatter
222 :
223 : END MODULE qs_cneo_utils
|
