1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Split and build its own idependent core_core SE interaction module
8!> \author Teodoro Laino [tlaino] - 05.2009
9!> \par History
10!>      Teodoro Laino (05.2009) [tlaino] - create
11! **************************************************************************************************
12MODULE se_core_core
13   USE atomic_kind_types,               ONLY: atomic_kind_type,&
14                                              get_atomic_kind_set
15   USE atprop_types,                    ONLY: atprop_array_init,&
16                                              atprop_type
17   USE cell_types,                      ONLY: cell_type
18   USE cp_control_types,                ONLY: dft_control_type,&
19                                              semi_empirical_control_type
20   USE cp_para_types,                   ONLY: cp_para_env_type
21   USE ewald_environment_types,         ONLY: ewald_env_get,&
22                                              ewald_environment_type
23   USE ewald_pw_types,                  ONLY: ewald_pw_get,&
24                                              ewald_pw_type
25   USE input_constants,                 ONLY: &
26        do_method_am1, do_method_mndo, do_method_mndod, do_method_pdg, do_method_pm3, &
27        do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1
28   USE kinds,                           ONLY: dp
29   USE message_passing,                 ONLY: mp_sum
30   USE particle_types,                  ONLY: particle_type
31   USE qs_energy_types,                 ONLY: qs_energy_type
32   USE qs_environment_types,            ONLY: get_qs_env,&
33                                              qs_environment_type
34   USE qs_force_types,                  ONLY: qs_force_type
35   USE qs_kind_types,                   ONLY: get_qs_kind,&
36                                              qs_kind_type
37   USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
38                                              neighbor_list_iterate,&
39                                              neighbor_list_iterator_create,&
40                                              neighbor_list_iterator_p_type,&
41                                              neighbor_list_iterator_release,&
42                                              neighbor_list_set_p_type
43   USE semi_empirical_int_arrays,       ONLY: rij_threshold
44   USE semi_empirical_integrals,        ONLY: corecore,&
45                                              dcorecore
46   USE semi_empirical_types,            ONLY: get_se_param,&
47                                              se_int_control_type,&
48                                              se_taper_type,&
49                                              semi_empirical_p_type,&
50                                              semi_empirical_type,&
51                                              setup_se_int_control_type
52   USE semi_empirical_utils,            ONLY: finalize_se_taper,&
53                                              get_se_type,&
54                                              initialize_se_taper
55   USE virial_methods,                  ONLY: virial_pair_force
56   USE virial_types,                    ONLY: virial_type
57#include "./base/base_uses.f90"
58
59   IMPLICIT NONE
60   PRIVATE
61
62   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_core_core'
63   LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
64
65   PUBLIC :: se_core_core_interaction
66
67CONTAINS
68
69! **************************************************************************************************
70!> \brief Evaluates the core-core interactions for NDDO methods
71!> \param qs_env ...
72!> \param para_env ...
73!> \param calculate_forces ...
74!> \date 04.2008 [tlaino]
75!> \author Teodoro Laino [tlaino] - University of Zurich
76! **************************************************************************************************
77   SUBROUTINE se_core_core_interaction(qs_env, para_env, calculate_forces)
78      TYPE(qs_environment_type), POINTER                 :: qs_env
79      TYPE(cp_para_env_type), POINTER                    :: para_env
80      LOGICAL, INTENT(in)                                :: calculate_forces
81
82      CHARACTER(len=*), PARAMETER :: routineN = 'se_core_core_interaction', &
83         routineP = moduleN//':'//routineN
84
85      INTEGER                                            :: atom_a, atom_b, handle, iab, iatom, &
86                                                            ikind, itype, jatom, jkind, natom, &
87                                                            nkind
88      INTEGER, DIMENSION(:), POINTER                     :: atom_of_kind
89      LOGICAL                                            :: anag, atener, defined, use_virial
90      LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: se_defined
91      REAL(KIND=dp)                                      :: delta, dr1, dr3inv(3), enuc, enucij, &
92                                                            enuclear, r2inv, r3inv, rinv
93      REAL(KIND=dp), DIMENSION(3)                        :: force_ab, rij
94      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
95      TYPE(atprop_type), POINTER                         :: atprop
96      TYPE(cell_type), POINTER                           :: cell
97      TYPE(dft_control_type), POINTER                    :: dft_control
98      TYPE(ewald_environment_type), POINTER              :: ewald_env
99      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
100      TYPE(neighbor_list_iterator_p_type), &
101         DIMENSION(:), POINTER                           :: nl_iterator
102      TYPE(neighbor_list_set_p_type), DIMENSION(:), &
103         POINTER                                         :: sab_se
104      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
105      TYPE(qs_energy_type), POINTER                      :: energy
106      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
107      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
108      TYPE(se_int_control_type)                          :: se_int_control
109      TYPE(se_taper_type), POINTER                       :: se_taper
110      TYPE(semi_empirical_control_type), POINTER         :: se_control
111      TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_param
112      TYPE(semi_empirical_type), POINTER                 :: se_kind_a, se_kind_b
113      TYPE(virial_type), POINTER                         :: virial
114
115      enuclear = 0.0_dp
116      NULLIFY (dft_control, cell, force, particle_set, se_control, se_taper, atomic_kind_set, &
117               virial, atprop)
118
119      CALL timeset(routineN, handle)
120      CPASSERT(ASSOCIATED(qs_env))
121      CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper, &
122                      virial=virial, atprop=atprop, energy=energy)
123
124      CALL initialize_se_taper(se_taper, coulomb=.TRUE.)
125      ! Parameters
126      se_control => dft_control%qs_control%se_control
127      anag = se_control%analytical_gradients
128      use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
129      CALL setup_se_int_control_type(se_int_control, do_ewald_r3=se_control%do_ewald_r3, &
130                                     do_ewald_gks=se_control%do_ewald_gks, integral_screening=se_control%integral_screening, &
131                                     shortrange=(se_control%do_ewald .OR. se_control%do_ewald_gks), &
132                                     max_multipole=se_control%max_multipole, pc_coulomb_int=.FALSE.)
133
134      ! atomic energy decomposition
135      atener = atprop%energy
136      IF (atener) THEN
137         CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
138         natom = SIZE(particle_set)
139         CALL atprop_array_init(atprop%atecc, natom)
140      END IF
141
142      ! Retrieve some information if GKS ewald scheme is used
143      IF (se_control%do_ewald_gks) THEN
144         CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw)
145         CALL ewald_env_get(ewald_env, alpha=se_int_control%ewald_gks%alpha)
146         CALL ewald_pw_get(ewald_pw, pw_big_pool=se_int_control%ewald_gks%pw_pool, &
147                           dg=se_int_control%ewald_gks%dg)
148         ! Virial not implemented
149         CPASSERT(.NOT. use_virial)
150      END IF
151
152      CALL get_qs_env(qs_env=qs_env, sab_se=sab_se, atomic_kind_set=atomic_kind_set, &
153                      qs_kind_set=qs_kind_set)
154
155      nkind = SIZE(atomic_kind_set)
156      ! Possibly compute forces
157      IF (calculate_forces) THEN
158         CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, force=force)
159         natom = SIZE(particle_set)
160         ALLOCATE (atom_of_kind(natom))
161         delta = se_control%delta
162         CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
163      END IF
164
165      itype = get_se_type(dft_control%qs_control%method_id)
166
167      ALLOCATE (se_kind_param(nkind), se_defined(nkind))
168      DO ikind = 1, nkind
169         CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
170         se_kind_param(ikind)%se_param => se_kind_a
171         CALL get_se_param(se_kind_a, defined=defined)
172         se_defined(ikind) = defined
173      END DO
174      CALL neighbor_list_iterator_create(nl_iterator, sab_se)
175      DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
176         CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
177         IF (.NOT. se_defined(ikind)) CYCLE
178         IF (.NOT. se_defined(jkind)) CYCLE
179         se_kind_a => se_kind_param(ikind)%se_param
180         se_kind_b => se_kind_param(jkind)%se_param
181         iab = ikind + nkind*(jkind - 1)
182         dr1 = DOT_PRODUCT(rij, rij)
183         enucij = 0._dp
184         IF (dr1 > rij_threshold) THEN
185            SELECT CASE (dft_control%qs_control%method_id)
186            CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, &
187                  do_method_rm1, do_method_mndod, do_method_pnnl)
188
189               ! Core-Core energy term
190               CALL corecore(se_kind_a, se_kind_b, rij, enuc=enuc, itype=itype, anag=anag, &
191                             se_int_control=se_int_control, se_taper=se_taper)
192               enucij = enucij + enuc
193               ! Residual integral (1/R^3) correction
194               IF (se_int_control%do_ewald_r3) THEN
195                  r2inv = 1.0_dp/dr1
196                  rinv = SQRT(r2inv)
197                  r3inv = rinv**3
198                  ! Core-Core term
199                  enucij = enucij + se_kind_a%expns3_int(jkind)%expns3%core_core*r3inv
200               END IF
201
202               ! Core-Core Derivatives
203               IF (calculate_forces) THEN
204                  atom_a = atom_of_kind(iatom)
205                  atom_b = atom_of_kind(jatom)
206
207                  CALL dcorecore(se_kind_a, se_kind_b, rij, denuc=force_ab, itype=itype, delta=delta, &
208                                 anag=anag, se_int_control=se_int_control, se_taper=se_taper)
209
210                  ! Residual integral (1/R^3) correction
211                  IF (se_int_control%do_ewald_r3) THEN
212                     dr3inv = -3.0_dp*rij*r3inv*r2inv
213                     ! Derivatives of core-core terms
214                     force_ab = force_ab + se_kind_a%expns3_int(jkind)%expns3%core_core*dr3inv
215                  END IF
216                  IF (use_virial) THEN
217                     CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
218                  END IF
219
220                  ! Sum up force components
221                  force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1)
222                  force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1)
223
224                  force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2)
225                  force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2)
226
227                  force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3)
228                  force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3)
229               END IF
230            CASE DEFAULT
231               CPABORT("")
232            END SELECT
233         ELSE
234            IF (se_int_control%do_ewald_gks) THEN
235               ! Core-Core energy term (self term in periodic systems)
236               CALL corecore(se_kind_a, se_kind_b, rij, enuc=enuc, itype=itype, anag=anag, &
237                             se_int_control=se_int_control, se_taper=se_taper)
238               enucij = enucij + 0.5_dp*enuc
239            END IF
240         END IF
241         IF (atener) THEN
242            atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*enucij
243            atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*enucij
244         END IF
245         enuclear = enuclear + enucij
246      END DO
247      CALL neighbor_list_iterator_release(nl_iterator)
248
249      DEALLOCATE (se_kind_param, se_defined)
250
251      IF (calculate_forces) THEN
252         DEALLOCATE (atom_of_kind)
253      END IF
254
255      CALL mp_sum(enuclear, para_env%group)
256      energy%core_overlap = enuclear
257      energy%core_overlap0 = enuclear
258
259      CALL finalize_se_taper(se_taper)
260      CALL timestop(handle)
261   END SUBROUTINE se_core_core_interaction
262
263END MODULE se_core_core
264
265