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