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 Setup and Methods for semi-empirical multipole types
10 : !> \author Teodoro Laino [tlaino] - 08.2008 Zurich University
11 : ! **************************************************************************************************
12 : MODULE semi_empirical_mpole_methods
13 :
14 : USE input_constants, ONLY: do_method_pnnl
15 : USE kinds, ONLY: dp
16 : USE mathconstants, ONLY: sqrt3
17 : USE semi_empirical_int_arrays, ONLY: alm,&
18 : indexa,&
19 : indexb,&
20 : se_map_alm
21 : USE semi_empirical_mpole_types, ONLY: nddo_mpole_create,&
22 : nddo_mpole_release,&
23 : nddo_mpole_type,&
24 : semi_empirical_mpole_p_create,&
25 : semi_empirical_mpole_p_type,&
26 : semi_empirical_mpole_type
27 : USE semi_empirical_par_utils, ONLY: amn_l
28 : USE semi_empirical_types, ONLY: semi_empirical_type
29 : #include "./base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 :
33 : PRIVATE
34 :
35 : ! *** Global parameters ***
36 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_mpole_methods'
38 :
39 : PUBLIC :: semi_empirical_mpole_p_setup, &
40 : nddo_mpole_setup, &
41 : quadrupole_sph_to_cart
42 :
43 : CONTAINS
44 :
45 : ! **************************************************************************************************
46 : !> \brief Setup semi-empirical mpole type
47 : !> This function setup for each semi-empirical type a structure containing
48 : !> the multipolar expansion for all possible combination on-site of atomic
49 : !> orbitals ( \mu \nu |
50 : !> \param mpoles ...
51 : !> \param se_parameter ...
52 : !> \param method ...
53 : !> \date 09.2008
54 : !> \author Teodoro Laino [tlaino] - University of Zurich
55 : ! **************************************************************************************************
56 3968 : SUBROUTINE semi_empirical_mpole_p_setup(mpoles, se_parameter, method)
57 : TYPE(semi_empirical_mpole_p_type), DIMENSION(:), &
58 : POINTER :: mpoles
59 : TYPE(semi_empirical_type), POINTER :: se_parameter
60 : INTEGER, INTENT(IN) :: method
61 :
62 : CHARACTER(LEN=3), DIMENSION(9), PARAMETER :: &
63 : label_print_orb = [" s", " px", " py", " pz", "dx2", "dzx", "dz2", "dzy", "dxy"]
64 : INTEGER, DIMENSION(9), PARAMETER :: loc_index = [1, 2, 2, 2, 3, 3, 3, 3, 3]
65 :
66 : INTEGER :: a, b, i, ind1, ind2, j, k, k1, k2, mu, &
67 : natorb, ndim, nr
68 : REAL(KIND=dp) :: dlm, tmp, wp, ws, zb, ZP, ZS, zt
69 : REAL(KIND=dp), DIMENSION(3, 3, 45) :: M2
70 : REAL(KIND=dp), DIMENSION(3, 45) :: M1
71 : REAL(KIND=dp), DIMENSION(45) :: M0
72 : REAL(KIND=dp), DIMENSION(6, 0:2) :: amn
73 : TYPE(semi_empirical_mpole_type), POINTER :: mpole
74 :
75 3968 : CPASSERT(.NOT. ASSOCIATED(mpoles))
76 : ! If there are atomic orbitals proceed with the expansion in multipoles
77 3968 : natorb = se_parameter%natorb
78 3968 : IF (natorb /= 0) THEN
79 2244 : ndim = natorb*(natorb + 1)/2
80 2244 : CALL semi_empirical_mpole_p_create(mpoles, ndim)
81 :
82 : ! Select method for multipolar expansion
83 :
84 : ! Fill in information on multipole expansion due to atomic orbitals charge
85 : ! distribution
86 2244 : NULLIFY (mpole)
87 2244 : CALL amn_l(se_parameter, amn)
88 13500 : DO i = 1, natorb
89 57540 : DO j = 1, i
90 44040 : ind1 = indexa(se_map_alm(i), se_map_alm(j))
91 44040 : ind2 = indexb(i, j)
92 : ! the order in the mpoles structure is like the standard one for the
93 : ! integrals: s px py pz dx2-y2 dzx dz2 dzy dxy (lower triangular)
94 : ! which differs from the order of the Hamiltonian in CP2K. But I
95 : ! preferred to keep this order for consistency with the integrals
96 44040 : mpole => mpoles(ind2)%mpole
97 44040 : mpole%indi = i
98 44040 : mpole%indj = j
99 44040 : a = loc_index(i)
100 44040 : b = loc_index(j)
101 44040 : mpole%c = HUGE(0.0_dp)
102 176160 : mpole%d = HUGE(0.0_dp)
103 264240 : mpole%qs = HUGE(0.0_dp)
104 572520 : mpole%qc = HUGE(0.0_dp)
105 :
106 : ! Charge
107 44040 : IF (alm(ind1, 0, 0) /= 0.0_dp) THEN
108 11256 : dlm = 1.0_dp/SQRT(REAL((2*0 + 1), KIND=dp))
109 11256 : tmp = -dlm*amn(indexb(a, b), 0)
110 11256 : mpole%c = tmp*alm(ind1, 0, 0)
111 11256 : mpole%task(1) = .TRUE.
112 : END IF
113 :
114 : ! Dipole
115 149280 : IF (ANY(alm(ind1, 1, -1:1) /= 0.0_dp)) THEN
116 13440 : dlm = 1.0_dp/SQRT(REAL((2*1 + 1), KIND=dp))
117 13440 : tmp = -dlm*amn(indexb(a, b), 1)
118 13440 : mpole%d(1) = tmp*alm(ind1, 1, 1)
119 13440 : mpole%d(2) = tmp*alm(ind1, 1, -1)
120 13440 : mpole%d(3) = tmp*alm(ind1, 1, 0)
121 13440 : mpole%task(2) = .TRUE.
122 : END IF
123 :
124 : ! Quadrupole
125 186694 : IF (ANY(alm(ind1, 2, -2:2) /= 0.0_dp)) THEN
126 23928 : dlm = 1.0_dp/SQRT(REAL((2*2 + 1), KIND=dp))
127 23928 : tmp = -dlm*amn(indexb(a, b), 2)
128 :
129 : ! Spherical components
130 23928 : mpole%qs(1) = tmp*alm(ind1, 2, 0) ! d3z2-r2
131 23928 : mpole%qs(2) = tmp*alm(ind1, 2, 1) ! dzx
132 23928 : mpole%qs(3) = tmp*alm(ind1, 2, -1) ! dzy
133 23928 : mpole%qs(4) = tmp*alm(ind1, 2, 2) ! dx2-y2
134 23928 : mpole%qs(5) = tmp*alm(ind1, 2, -2) ! dxy
135 :
136 : ! Convert into cartesian components
137 23928 : CALL quadrupole_sph_to_cart(mpole%qc, mpole%qs)
138 23928 : mpole%task(3) = .TRUE.
139 : END IF
140 :
141 11256 : IF (debug_this_module) THEN
142 : WRITE (*, '(A,2I6,A)') "Orbitals ", i, j, &
143 : " ("//label_print_orb(i)//","//label_print_orb(j)//")"
144 : IF (mpole%task(1)) WRITE (*, '(9F12.6)') mpole%c
145 : IF (mpole%task(2)) WRITE (*, '(9F12.6)') mpole%d
146 : IF (mpole%task(3)) WRITE (*, '(9F12.6)') mpole%qc
147 : WRITE (*, *)
148 : END IF
149 : END DO
150 : END DO
151 :
152 2244 : IF (method == do_method_pnnl) THEN
153 : ! No d-function for Schenter type integrals
154 28 : CPASSERT(natorb <= 4)
155 :
156 28 : M0 = 0.0_dp
157 28 : M1 = 0.0_dp
158 28 : M2 = 0.0_dp
159 :
160 140 : DO mu = 1, se_parameter%natorb
161 140 : M0(indexb(mu, mu)) = 1.0_dp
162 : END DO
163 :
164 28 : ZS = se_parameter%sto_exponents(0)
165 28 : ZP = se_parameter%sto_exponents(1)
166 28 : nr = se_parameter%nr
167 :
168 28 : ws = REAL((2*nr + 2)*(2*nr + 1), dp)/(24.0_dp*ZS**2)
169 112 : DO k = 1, 3
170 112 : M2(k, k, indexb(1, 1)) = ws
171 : END DO
172 :
173 28 : IF (ZP > 0._dp) THEN
174 28 : zt = SQRT(ZS*ZP)
175 28 : zb = 0.5_dp*(ZS + ZP)
176 112 : DO k = 1, 3
177 112 : M1(k, indexb(1, 1 + k)) = (zt/zb)**(2*nr + 1)*REAL(2*nr + 1, dp)/(2.0_dp*zb*sqrt3)
178 : END DO
179 :
180 28 : wp = REAL((2*nr + 2)*(2*nr + 1), dp)/(40.0_dp*ZP**2)
181 112 : DO k1 = 1, 3
182 364 : DO k2 = 1, 3
183 336 : IF (k1 == k2) THEN
184 84 : M2(k2, k2, indexb(1 + k1, 1 + k1)) = 3.0_dp*wp
185 : ELSE
186 168 : M2(k2, k2, indexb(1 + k1, 1 + k1)) = wp
187 : END IF
188 : END DO
189 : END DO
190 28 : M2(1, 2, indexb(1 + 1, 1 + 2)) = wp
191 28 : M2(2, 1, indexb(1 + 1, 1 + 2)) = wp
192 28 : M2(2, 3, indexb(1 + 2, 1 + 3)) = wp
193 28 : M2(3, 2, indexb(1 + 2, 1 + 3)) = wp
194 28 : M2(3, 1, indexb(1 + 3, 1 + 1)) = wp
195 28 : M2(1, 3, indexb(1 + 3, 1 + 1)) = wp
196 : END IF
197 :
198 140 : DO i = 1, natorb
199 420 : DO j = 1, i
200 280 : ind1 = indexa(se_map_alm(i), se_map_alm(j))
201 280 : ind2 = indexb(i, j)
202 280 : mpole => mpoles(ind2)%mpole
203 280 : mpole%indi = i
204 280 : mpole%indj = j
205 : ! Charge
206 280 : mpole%cs = -M0(indexb(i, j))
207 : ! Dipole
208 1120 : mpole%ds = -M1(1:3, indexb(i, j))
209 : ! Quadrupole
210 3640 : mpole%qq = -3._dp*M2(1:3, 1:3, indexb(i, j))
211 112 : IF (debug_this_module) THEN
212 : WRITE (*, '(A,2I6,A)') "Orbitals ", i, j, &
213 : " ("//label_print_orb(i)//","//label_print_orb(j)//")"
214 : WRITE (*, '(9F12.6)') mpole%cs
215 : WRITE (*, '(9F12.6)') mpole%ds
216 : WRITE (*, '(9F12.6)') mpole%qq
217 : WRITE (*, *)
218 : END IF
219 : END DO
220 : END DO
221 : ELSE
222 2216 : mpole%cs = mpole%c
223 8864 : mpole%ds = mpole%d
224 28808 : mpole%qq = mpole%qc
225 : END IF
226 : END IF
227 :
228 3968 : END SUBROUTINE semi_empirical_mpole_p_setup
229 :
230 : ! **************************************************************************************************
231 : !> \brief Transforms the quadrupole components from sphericals to cartesians
232 : !> \param qcart ...
233 : !> \param qsph ...
234 : !> \date 09.2008
235 : !> \author Teodoro Laino [tlaino] - University of Zurich
236 : ! **************************************************************************************************
237 51930 : SUBROUTINE quadrupole_sph_to_cart(qcart, qsph)
238 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: qcart
239 : REAL(KIND=dp), DIMENSION(5), INTENT(IN) :: qsph
240 :
241 : ! Notation
242 : ! qs(1) - d3z2-r2
243 : ! qs(2) - dzx
244 : ! qs(3) - dzy
245 : ! qs(4) - dx2-y2
246 : ! qs(5) - dxy
247 : ! Cartesian components
248 :
249 51930 : qcart(1, 1) = (qsph(4) - qsph(1)/SQRT(3.0_dp))*SQRT(3.0_dp)/2.0_dp
250 51930 : qcart(2, 1) = qsph(5)*SQRT(3.0_dp)/2.0_dp
251 51930 : qcart(3, 1) = qsph(2)*SQRT(3.0_dp)/2.0_dp
252 51930 : qcart(2, 2) = -(qsph(4) + qsph(1)/SQRT(3.0_dp))*SQRT(3.0_dp)/2.0_dp
253 51930 : qcart(3, 2) = qsph(3)*SQRT(3.0_dp)/2.0_dp
254 51930 : qcart(3, 3) = qsph(1)
255 : ! Symmetrize tensor
256 51930 : qcart(1, 2) = qcart(2, 1)
257 51930 : qcart(1, 3) = qcart(3, 1)
258 51930 : qcart(2, 3) = qcart(3, 2)
259 :
260 51930 : END SUBROUTINE quadrupole_sph_to_cart
261 :
262 : ! **************************************************************************************************
263 : !> \brief Setup NDDO multipole type
264 : !> \param nddo_mpole ...
265 : !> \param natom ...
266 : !> \date 09.2008
267 : !> \author Teodoro Laino [tlaino] - University of Zurich
268 : ! **************************************************************************************************
269 32 : SUBROUTINE nddo_mpole_setup(nddo_mpole, natom)
270 : TYPE(nddo_mpole_type), POINTER :: nddo_mpole
271 : INTEGER, INTENT(IN) :: natom
272 :
273 : CHARACTER(len=*), PARAMETER :: routineN = 'nddo_mpole_setup'
274 :
275 : INTEGER :: handle
276 :
277 32 : CALL timeset(routineN, handle)
278 :
279 32 : IF (ASSOCIATED(nddo_mpole)) THEN
280 0 : CALL nddo_mpole_release(nddo_mpole)
281 : END IF
282 32 : CALL nddo_mpole_create(nddo_mpole)
283 : ! Allocate Global Arrays
284 96 : ALLOCATE (nddo_mpole%charge(natom))
285 96 : ALLOCATE (nddo_mpole%dipole(3, natom))
286 96 : ALLOCATE (nddo_mpole%quadrupole(3, 3, natom))
287 :
288 64 : ALLOCATE (nddo_mpole%efield0(natom))
289 64 : ALLOCATE (nddo_mpole%efield1(3, natom))
290 64 : ALLOCATE (nddo_mpole%efield2(9, natom))
291 :
292 32 : CALL timestop(handle)
293 :
294 32 : END SUBROUTINE nddo_mpole_setup
295 :
296 : END MODULE semi_empirical_mpole_methods
|
