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