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 Working with the semi empirical parameter types.
10 : !> \author JGH (14.08.2004)
11 : ! **************************************************************************************************
12 : MODULE semi_empirical_utils
13 : USE basis_set_types, ONLY: allocate_sto_basis_set,&
14 : create_gto_from_sto_basis,&
15 : gto_basis_set_type,&
16 : set_sto_basis_set
17 : USE cell_types, ONLY: cell_type,&
18 : plane_distance
19 : USE cp_control_types, ONLY: semi_empirical_control_type
20 : USE input_constants, ONLY: &
21 : do_method_am1, do_method_mndo, do_method_mndod, do_method_pchg, do_method_pdg, &
22 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1
23 : USE input_section_types, ONLY: section_vals_type,&
24 : section_vals_val_get
25 : USE kinds, ONLY: dp
26 : USE semi_empirical_mpole_methods, ONLY: semi_empirical_mpole_p_setup
27 : USE semi_empirical_par_utils, ONLY: get_se_basis,&
28 : setup_1c_2el_int
29 : USE semi_empirical_parameters, ONLY: &
30 : am1_default_parameter, mndo_default_parameter, pcharge_default_parameter, &
31 : pdg_default_parameter, pm3_default_parameter, pm6_default_parameter, &
32 : pm6fm_default_parameter, pnnl_default_parameter, rm1_default_parameter
33 : USE semi_empirical_types, ONLY: se_taper_type,&
34 : semi_empirical_type
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_utils'
42 :
43 : PUBLIC :: init_se_param, se_param_set_default, get_se_type, &
44 : initialize_se_taper, finalize_se_taper, se_cutoff_compatible
45 :
46 : CONTAINS
47 : ! **************************************************************************************************
48 : !> \brief Reset cutoffs trying to be somehow a bit smarter
49 : !> \param se_control ...
50 : !> \param se_section ...
51 : !> \param cell ...
52 : !> \param output_unit ...
53 : !> \author Teodoro Laino [tlaino] - 03.2009
54 : ! **************************************************************************************************
55 3992 : SUBROUTINE se_cutoff_compatible(se_control, se_section, cell, output_unit)
56 : TYPE(semi_empirical_control_type), POINTER :: se_control
57 : TYPE(section_vals_type), POINTER :: se_section
58 : TYPE(cell_type), POINTER :: cell
59 : INTEGER, INTENT(IN) :: output_unit
60 :
61 : LOGICAL :: explicit1, explicit2
62 : REAL(KIND=dp) :: rc
63 :
64 : ! Coulomb Cutoff Taper
65 :
66 998 : CALL section_vals_val_get(se_section, "COULOMB%CUTOFF", explicit=explicit1)
67 998 : CALL section_vals_val_get(se_section, "COULOMB%RC_TAPER", explicit=explicit2)
68 998 : IF ((.NOT. explicit1) .AND. se_control%do_ewald_gks) THEN
69 : rc = MAX(0.5*plane_distance(1, 0, 0, cell), &
70 : 0.5*plane_distance(0, 1, 0, cell), &
71 2 : 0.5*plane_distance(0, 0, 1, cell))
72 2 : IF (rc /= se_control%cutoff_cou) THEN
73 2 : IF (output_unit > 0) THEN
74 1 : WRITE (output_unit, *)
75 1 : WRITE (output_unit, '(A,T37,A)') " SEMIEMPIRICAL|", &
76 2 : " Coulomb Integral cutoff/taper was redefined"
77 1 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Old value [a.u.]", &
78 2 : se_control%cutoff_cou
79 1 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| New value [a.u.]", rc
80 1 : WRITE (output_unit, *)
81 : END IF
82 : END IF
83 2 : se_control%cutoff_cou = rc
84 2 : IF (.NOT. explicit2) se_control%taper_cou = rc
85 2184 : ELSE IF ((.NOT. explicit1) .AND. (ALL(cell%perd == 0))) THEN
86 : rc = MAX(plane_distance(1, 0, 0, cell), &
87 : plane_distance(0, 1, 0, cell), &
88 374 : plane_distance(0, 0, 1, cell))
89 374 : IF (rc /= se_control%cutoff_cou) THEN
90 374 : IF (output_unit > 0) THEN
91 188 : WRITE (output_unit, *)
92 188 : WRITE (output_unit, '(A,T37,A)') " SEMIEMPIRICAL|", &
93 376 : " Coulomb Integral cutoff/taper was redefined"
94 188 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Old value [a.u.]", &
95 376 : se_control%cutoff_cou
96 188 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| New value [a.u.]", rc
97 188 : WRITE (output_unit, *)
98 : END IF
99 : END IF
100 374 : se_control%cutoff_cou = rc
101 374 : IF (.NOT. explicit2) se_control%taper_cou = rc
102 : END IF
103 998 : IF (output_unit > 0) THEN
104 500 : WRITE (output_unit, *)
105 500 : WRITE (output_unit, '(A,T44,A)') " SEMIEMPIRICAL|", &
106 1000 : " Coulomb Integral cutoff/taper values"
107 500 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Cutoff [a.u.]", &
108 1000 : se_control%cutoff_cou
109 500 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Taper [a.u.]", &
110 1000 : se_control%taper_cou
111 500 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Range [a.u.]", &
112 1000 : se_control%range_cou
113 500 : WRITE (output_unit, *)
114 : END IF
115 : ! Exchange Cutoff Taper
116 998 : CALL section_vals_val_get(se_section, "EXCHANGE%CUTOFF", explicit=explicit1)
117 998 : CALL section_vals_val_get(se_section, "EXCHANGE%RC_TAPER", explicit=explicit2)
118 998 : rc = se_control%cutoff_exc
119 998 : IF (.NOT. explicit1) THEN
120 : rc = MIN(rc, MAX(0.25_dp*plane_distance(1, 0, 0, cell), &
121 : 0.25_dp*plane_distance(0, 1, 0, cell), &
122 954 : 0.25_dp*plane_distance(0, 0, 1, cell)))
123 :
124 954 : IF (rc /= se_control%cutoff_exc) THEN
125 952 : IF (output_unit > 0) THEN
126 477 : WRITE (output_unit, *)
127 477 : WRITE (output_unit, '(A,T36,A)') " SEMIEMPIRICAL|", &
128 954 : " Exchange Integral cutoff/taper was redefined"
129 477 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Default value [a.u.]", &
130 954 : se_control%cutoff_exc
131 477 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| New value [a.u.]", rc
132 477 : WRITE (output_unit, *)
133 : END IF
134 : END IF
135 : END IF
136 998 : se_control%cutoff_exc = rc
137 998 : IF (.NOT. explicit2) se_control%taper_exc = rc
138 :
139 998 : IF (output_unit > 0) THEN
140 500 : WRITE (output_unit, *)
141 500 : WRITE (output_unit, '(A,T43,A)') " SEMIEMPIRICAL|", &
142 1000 : " Exchange Integral cutoff/taper values"
143 500 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Cutoff [a.u.]", &
144 1000 : se_control%cutoff_exc
145 500 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Taper [a.u.]", &
146 1000 : se_control%taper_exc
147 500 : WRITE (output_unit, '(A,T71,F10.3)') " SEMIEMPIRICAL| Range [a.u.]", &
148 1000 : se_control%range_exc
149 500 : WRITE (output_unit, *)
150 : END IF
151 :
152 998 : END SUBROUTINE se_cutoff_compatible
153 :
154 : ! **************************************************************************************************
155 : !> \brief Initializes the semi-empirical taper for a chunk calculation
156 : !> \param se_taper ...
157 : !> \param coulomb ...
158 : !> \param exchange ...
159 : !> \param lr_corr ...
160 : !> \author Teodoro Laino [tlaino] - 03.2009
161 : ! **************************************************************************************************
162 91028 : SUBROUTINE initialize_se_taper(se_taper, coulomb, exchange, lr_corr)
163 : TYPE(se_taper_type), POINTER :: se_taper
164 : LOGICAL, INTENT(IN), OPTIONAL :: coulomb, exchange, lr_corr
165 :
166 : LOGICAL :: check, l_coulomb, l_exchange, l_lrc
167 :
168 91028 : check = .NOT. ASSOCIATED(se_taper%taper)
169 91028 : CPASSERT(check)
170 91028 : l_coulomb = .FALSE.
171 91028 : l_exchange = .FALSE.
172 91028 : l_lrc = .FALSE.
173 91028 : IF (PRESENT(coulomb)) l_coulomb = coulomb
174 91028 : IF (PRESENT(exchange)) l_exchange = exchange
175 91028 : IF (PRESENT(lr_corr)) l_lrc = lr_corr
176 91028 : IF (l_coulomb) THEN
177 48432 : check = (.NOT. l_exchange) .AND. (.NOT. l_lrc)
178 48432 : CPASSERT(check)
179 48432 : se_taper%taper => se_taper%taper_cou
180 : END IF
181 91028 : IF (l_exchange) THEN
182 41094 : check = (.NOT. l_coulomb) .AND. (.NOT. l_lrc)
183 41094 : CPASSERT(check)
184 41094 : se_taper%taper => se_taper%taper_exc
185 : END IF
186 91028 : IF (l_lrc) THEN
187 1502 : check = (.NOT. l_coulomb) .AND. (.NOT. l_exchange)
188 1502 : CPASSERT(check)
189 1502 : se_taper%taper => se_taper%taper_lrc
190 : END IF
191 91028 : END SUBROUTINE initialize_se_taper
192 :
193 : ! **************************************************************************************************
194 : !> \brief Finalizes the semi-empirical taper for a chunk calculation
195 : !> \param se_taper ...
196 : !> \author Teodoro Laino [tlaino] - 03.2009
197 : ! **************************************************************************************************
198 91028 : SUBROUTINE finalize_se_taper(se_taper)
199 : TYPE(se_taper_type), POINTER :: se_taper
200 :
201 : LOGICAL :: check
202 :
203 91028 : check = ASSOCIATED(se_taper%taper)
204 91028 : CPASSERT(check)
205 91028 : NULLIFY (se_taper%taper)
206 91028 : END SUBROUTINE finalize_se_taper
207 :
208 : ! **************************************************************************************************
209 : !> \brief Initialize semi_empirical type
210 : !> \param sep ...
211 : !> \param orb_basis_set ...
212 : !> \param ngauss ...
213 : ! **************************************************************************************************
214 2240 : SUBROUTINE init_se_param(sep, orb_basis_set, ngauss)
215 : TYPE(semi_empirical_type), POINTER :: sep
216 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
217 : INTEGER, INTENT(IN) :: ngauss
218 :
219 2240 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: symbol
220 : INTEGER :: l, nshell
221 2240 : INTEGER, DIMENSION(:), POINTER :: lq, nq
222 2240 : REAL(KIND=dp), DIMENSION(:), POINTER :: zet
223 :
224 2240 : IF (ASSOCIATED(sep)) THEN
225 2240 : CALL allocate_sto_basis_set(sep%basis)
226 2240 : nshell = 0
227 2240 : IF (sep%natorb == 1) nshell = 1
228 2240 : IF (sep%natorb == 4) nshell = 2
229 2240 : IF (sep%natorb == 9) nshell = 3
230 2240 : ALLOCATE (nq(0:3), lq(0:3), zet(0:3))
231 :
232 2240 : ALLOCATE (symbol(0:3))
233 :
234 11200 : symbol = ""
235 11200 : nq = 0
236 11200 : lq = 0
237 11200 : zet = 0._dp
238 6990 : DO l = 0, nshell - 1
239 4750 : nq(l) = get_se_basis(sep, l)
240 4750 : lq(l) = l
241 4750 : zet(l) = sep%sto_exponents(l)
242 4750 : IF (l == 0) WRITE (symbol(0), '(I1,A1)') nq(l), "S"
243 4750 : IF (l == 1) WRITE (symbol(1), '(I1,A1)') nq(l), "P"
244 6990 : IF (l == 2) WRITE (symbol(2), '(I1,A1)') nq(l), "D"
245 : END DO
246 :
247 2240 : IF (nshell > 0) THEN
248 2240 : sep%ngauss = ngauss
249 : CALL set_sto_basis_set(sep%basis, name=sep%name, nshell=nshell, symbol=symbol, &
250 2240 : nq=nq, lq=lq, zet=zet)
251 2240 : CALL create_gto_from_sto_basis(sep%basis, orb_basis_set, sep%ngauss)
252 : END IF
253 :
254 2240 : DEALLOCATE (nq)
255 2240 : DEALLOCATE (lq)
256 2240 : DEALLOCATE (zet)
257 2240 : DEALLOCATE (symbol)
258 : ELSE
259 0 : CPABORT("The pointer sep is not associated")
260 : END IF
261 :
262 2240 : END SUBROUTINE init_se_param
263 :
264 : ! **************************************************************************************************
265 : !> \brief Initialize parameter for a semi_empirival type
266 : !> \param sep ...
267 : !> \param z ...
268 : !> \param method ...
269 : ! **************************************************************************************************
270 3964 : SUBROUTINE se_param_set_default(sep, z, method)
271 :
272 : TYPE(semi_empirical_type), POINTER :: sep
273 : INTEGER, INTENT(IN) :: z, method
274 :
275 3964 : IF (ASSOCIATED(sep)) THEN
276 3964 : IF (z < 0) THEN
277 0 : CPABORT("Atomic number < 0")
278 : END IF
279 4232 : SELECT CASE (method)
280 : CASE (do_method_am1)
281 268 : CALL am1_default_parameter(sep, z)
282 : CASE (do_method_rm1)
283 4 : CALL rm1_default_parameter(sep, z)
284 : CASE (do_method_pm3)
285 98 : CALL pm3_default_parameter(sep, z)
286 : CASE (do_method_pm6)
287 1700 : CALL pm6_default_parameter(sep, z)
288 : CASE (do_method_pm6fm)
289 0 : CALL pm6fm_default_parameter(sep, z)
290 : CASE (do_method_pdg)
291 4 : CALL pdg_default_parameter(sep, z)
292 : CASE (do_method_mndo)
293 106 : CALL mndo_default_parameter(sep, z, do_method_mndo)
294 : CASE (do_method_mndod)
295 32 : CALL mndo_default_parameter(sep, z, do_method_mndod)
296 : CASE (do_method_pnnl)
297 28 : CALL pnnl_default_parameter(sep, z)
298 : CASE (do_method_pchg)
299 1724 : CALL pcharge_default_parameter(sep, z)
300 : CASE DEFAULT
301 3964 : CPABORT("Semiempirical method unknown")
302 : END SELECT
303 : ELSE
304 0 : CPABORT("The pointer sep is not associated")
305 : END IF
306 :
307 : ! Check if the element has been defined..
308 3964 : IF (.NOT. sep%defined) &
309 : CALL cp_abort(__LOCATION__, &
310 : "Semiempirical type ("//TRIM(sep%name)//") cannot be defined for "// &
311 0 : "the requested parameterization.")
312 :
313 : ! Fill 1 center - 2 electron integrals
314 3964 : CALL setup_1c_2el_int(sep)
315 :
316 : ! Fill multipolar expansion of atomic orbitals charge distributions
317 3964 : CALL semi_empirical_mpole_p_setup(sep%w_mpole, sep, method)
318 :
319 : ! Get the value of the size of CORE integral array
320 3964 : sep%core_size = 0
321 3964 : IF (sep%natorb > 0) sep%core_size = 1
322 3964 : IF (sep%natorb > 1) sep%core_size = 4
323 3964 : IF (sep%dorb) sep%core_size = 10
324 :
325 : ! Get size of the all possible combinations of atomic orbitals
326 3964 : sep%atm_int_size = (sep%natorb + 1)*sep%natorb/2
327 :
328 3964 : END SUBROUTINE se_param_set_default
329 :
330 : ! **************************************************************************************************
331 : !> \brief Gives back the unique semi_empirical METHOD type
332 : !> \param se_method ...
333 : !> \return ...
334 : ! **************************************************************************************************
335 51658 : FUNCTION get_se_type(se_method) RESULT(se_type)
336 :
337 : INTEGER, INTENT(IN) :: se_method
338 : INTEGER :: se_type
339 :
340 83166 : SELECT CASE (se_method)
341 : CASE DEFAULT
342 31508 : se_type = se_method
343 : CASE (do_method_am1, do_method_rm1)
344 51658 : se_type = do_method_am1
345 : END SELECT
346 :
347 51658 : END FUNCTION get_se_type
348 :
349 : END MODULE semi_empirical_utils
350 :
|
