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