Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief implementation of dipole and three-body part of Siepmann-Sprik potential
10 : !> dipole term: 3rd term in Eq. (1) in J. Chem. Phys., Vol. 102, p.511
11 : !> three-body term: Eq. (4) in J. Chem. Phys., Vol. 102, p. 511
12 : !> remaining terms of Siepmann-Sprik potential can be given via the GENPOT section
13 : !> \par History
14 : !> 12.2012 created
15 : !> \author Dorothea Golze
16 : ! **************************************************************************************************
17 : MODULE manybody_siepmann
18 :
19 : USE atomic_kind_types, ONLY: get_atomic_kind
20 : USE cell_types, ONLY: cell_type,&
21 : pbc
22 : USE cp_log_handling, ONLY: cp_get_default_logger,&
23 : cp_logger_type
24 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
25 : cp_print_key_unit_nr
26 : USE fist_neighbor_list_types, ONLY: fist_neighbor_type,&
27 : neighbor_kind_pairs_type
28 : USE fist_nonbond_env_types, ONLY: pos_type
29 : USE input_section_types, ONLY: section_vals_type
30 : USE kinds, ONLY: dp
31 : USE message_passing, ONLY: mp_para_env_type
32 : USE pair_potential_types, ONLY: pair_potential_pp_type,&
33 : pair_potential_single_type,&
34 : siepmann_pot_type,&
35 : siepmann_type
36 : USE particle_types, ONLY: particle_type
37 : USE util, ONLY: sort
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 : PUBLIC :: setup_siepmann_arrays, destroy_siepmann_arrays, &
44 : siepmann_energy, siepmann_forces_v2, siepmann_forces_v3, &
45 : print_nr_ions_siepmann
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_siepmann'
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief energy of two-body dipole term and three-body term
52 : !> \param pot_loc ...
53 : !> \param siepmann ...
54 : !> \param r_last_update_pbc ...
55 : !> \param atom_a ...
56 : !> \param atom_b ...
57 : !> \param nloc_size ...
58 : !> \param full_loc_list ...
59 : !> \param cell_v ...
60 : !> \param cell ...
61 : !> \param drij ...
62 : !> \param particle_set ...
63 : !> \param nr_oh number of OH- ions near surface
64 : !> \param nr_h3o number of hydronium ions near surface
65 : !> \param nr_o number of O^(2-) ions near surface
66 : !> \author Dorothea Golze - 11.2012 - University of Zurich
67 : ! **************************************************************************************************
68 318 : SUBROUTINE siepmann_energy(pot_loc, siepmann, r_last_update_pbc, atom_a, atom_b, &
69 318 : nloc_size, full_loc_list, cell_v, cell, drij, particle_set, &
70 : nr_oh, nr_h3o, nr_o)
71 :
72 : REAL(KIND=dp), INTENT(OUT) :: pot_loc
73 : TYPE(siepmann_pot_type), POINTER :: siepmann
74 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
75 : INTEGER, INTENT(IN) :: atom_a, atom_b, nloc_size
76 : INTEGER, DIMENSION(2, 1:nloc_size) :: full_loc_list
77 : REAL(KIND=dp), DIMENSION(3) :: cell_v
78 : TYPE(cell_type), POINTER :: cell
79 : REAL(KIND=dp) :: drij
80 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
81 : INTEGER, INTENT(INOUT) :: nr_oh, nr_h3o, nr_o
82 :
83 : REAL(KIND=dp) :: a_ij, D, E, f2, Phi_ij, pot_loc_v2, &
84 : pot_loc_v3
85 :
86 : a_ij = siep_a_ij(siepmann, r_last_update_pbc, atom_a, atom_b, nloc_size, &
87 : full_loc_list, cell_v, siepmann%rcutsq, particle_set, &
88 318 : cell)
89 : Phi_ij = siep_Phi_ij(siepmann, r_last_update_pbc, atom_a, atom_b, &
90 318 : cell_v, cell, siepmann%rcutsq, particle_set, nr_oh, nr_h3o, nr_o)
91 318 : f2 = siep_f2(siepmann, drij)
92 318 : D = siepmann%D
93 318 : E = siepmann%E
94 :
95 : !two-body part --> dipole term
96 318 : pot_loc_v2 = -D*f2*drij**(-3)*Phi_ij
97 :
98 : !three-body part
99 318 : pot_loc_v3 = E*f2*drij**(-siepmann%beta)*a_ij
100 :
101 318 : pot_loc = pot_loc_v2 + pot_loc_v3
102 :
103 318 : END SUBROUTINE siepmann_energy
104 :
105 : ! **************************************************************************************************
106 : !> \brief f2(r) corresponds to Equation (2) in J. Chem. Phys., Vol. 102, p.511
107 : !> \param siepmann ...
108 : !> \param r distance between oxygen and metal atom
109 : !> \return ...
110 : ! **************************************************************************************************
111 636 : FUNCTION siep_f2(siepmann, r)
112 : TYPE(siepmann_pot_type), POINTER :: siepmann
113 : REAL(KIND=dp), INTENT(IN) :: r
114 : REAL(KIND=dp) :: siep_f2
115 :
116 : REAL(KIND=dp) :: rcut
117 :
118 636 : rcut = SQRT(siepmann%rcutsq)
119 636 : siep_f2 = 0.0_dp
120 636 : IF (r < rcut) THEN
121 636 : siep_f2 = EXP(siepmann%B/(r - rcut))
122 : END IF
123 636 : END FUNCTION siep_f2
124 :
125 : ! **************************************************************************************************
126 : !> \brief f2(r)_d derivative of f2
127 : !> \param siepmann ...
128 : !> \param r distance between oxygen and metal atom
129 : !> \return ...
130 : ! **************************************************************************************************
131 318 : FUNCTION siep_f2_d(siepmann, r)
132 : TYPE(siepmann_pot_type), POINTER :: siepmann
133 : REAL(KIND=dp), INTENT(IN) :: r
134 : REAL(KIND=dp) :: siep_f2_d
135 :
136 : REAL(KIND=dp) :: B, rcut
137 :
138 318 : rcut = SQRT(siepmann%rcutsq)
139 318 : B = siepmann%B
140 318 : siep_f2_d = 0.0_dp
141 318 : IF (r < rcut) THEN
142 318 : siep_f2_d = -B*EXP(B/(r - rcut))/(r - rcut)**2
143 : END IF
144 :
145 318 : END FUNCTION siep_f2_d
146 :
147 : ! **************************************************************************************************
148 : !> \brief exponential part of three-body term, see Equation (4) in J. Chem.
149 : !> Phys., Vol. 102, p.511
150 : !> \param siepmann ...
151 : !> \param r_last_update_pbc ...
152 : !> \param iparticle ...
153 : !> \param jparticle ...
154 : !> \param n_loc_size ...
155 : !> \param full_loc_list ...
156 : !> \param cell_v ...
157 : !> \param rcutsq ...
158 : !> \param particle_set ...
159 : !> \param cell ...
160 : !> \return ...
161 : !> \par History
162 : !> Using a local list of neighbors
163 : ! **************************************************************************************************
164 477 : FUNCTION siep_a_ij(siepmann, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
165 477 : full_loc_list, cell_v, rcutsq, particle_set, cell)
166 : TYPE(siepmann_pot_type), POINTER :: siepmann
167 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
168 : INTEGER, INTENT(IN) :: iparticle, jparticle, n_loc_size
169 : INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list
170 : REAL(KIND=dp), DIMENSION(3) :: cell_v
171 : REAL(KIND=dp), INTENT(IN) :: rcutsq
172 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
173 : TYPE(cell_type), POINTER :: cell
174 : REAL(KIND=dp) :: siep_a_ij
175 :
176 : CHARACTER(LEN=2) :: element_symbol
177 : INTEGER :: ilist, kparticle
178 : REAL(KIND=dp) :: costheta, drji, drjk, F, rab2_max, &
179 : rji(3), rjk(3), theta
180 :
181 477 : siep_a_ij = 0.0_dp
182 477 : rab2_max = rcutsq
183 477 : F = siepmann%F
184 : CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
185 477 : element_symbol=element_symbol)
186 477 : IF (element_symbol /= "O") RETURN
187 1272 : rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
188 1272 : drji = SQRT(DOT_PRODUCT(rji, rji))
189 22266 : DO ilist = 1, n_loc_size
190 21948 : kparticle = full_loc_list(2, ilist)
191 21948 : IF (kparticle == jparticle) CYCLE
192 21630 : rjk(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(kparticle)%r(:), cell)
193 86520 : drjk = DOT_PRODUCT(rjk, rjk)
194 21630 : IF (drjk > rab2_max) CYCLE
195 2862 : drjk = SQRT(drjk)
196 11448 : costheta = DOT_PRODUCT(rji, rjk)/(drji*drjk)
197 2862 : IF (costheta < -1.0_dp) costheta = -1.0_dp
198 2862 : IF (costheta > +1.0_dp) costheta = +1.0_dp
199 2862 : theta = ACOS(costheta)
200 22266 : siep_a_ij = siep_a_ij + EXP(F*(COS(theta/2.0_dp))**2)
201 : END DO
202 : END FUNCTION siep_a_ij
203 :
204 : ! **************************************************************************************************
205 : !> \brief derivative of a_ij
206 : !> \param siepmann ...
207 : !> \param r_last_update_pbc ...
208 : !> \param iparticle ...
209 : !> \param jparticle ...
210 : !> \param f_nonbond ...
211 : !> \param prefactor ...
212 : !> \param n_loc_size ...
213 : !> \param full_loc_list ...
214 : !> \param cell_v ...
215 : !> \param cell ...
216 : !> \param rcutsq ...
217 : !> \param use_virial ...
218 : !> \par History
219 : !> Using a local list of neighbors
220 : ! **************************************************************************************************
221 318 : SUBROUTINE siep_a_ij_d(siepmann, r_last_update_pbc, iparticle, jparticle, f_nonbond, &
222 159 : prefactor, n_loc_size, full_loc_list, &
223 : cell_v, cell, rcutsq, use_virial)
224 : TYPE(siepmann_pot_type), POINTER :: siepmann
225 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
226 : INTEGER, INTENT(IN) :: iparticle, jparticle
227 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond
228 : REAL(KIND=dp), INTENT(IN) :: prefactor
229 : INTEGER, INTENT(IN) :: n_loc_size
230 : INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list
231 : REAL(KIND=dp), DIMENSION(3) :: cell_v
232 : TYPE(cell_type), POINTER :: cell
233 : REAL(KIND=dp), INTENT(IN) :: rcutsq
234 : LOGICAL, INTENT(IN) :: use_virial
235 :
236 : INTEGER :: ilist, kparticle, nparticle
237 : REAL(KIND=dp) :: costheta, d_expterm, dcos_thetahalf, &
238 : drji, drjk, F, rab2_max, theta
239 : REAL(KIND=dp), DIMENSION(3) :: dcosdri, dcosdrj, dcosdrk, dri, drj, &
240 : drk, rji, rji_hat, rjk, rjk_hat
241 :
242 159 : rab2_max = rcutsq
243 159 : F = siepmann%F
244 :
245 636 : rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
246 636 : drji = SQRT(DOT_PRODUCT(rji, rji))
247 636 : rji_hat(:) = rji(:)/drji
248 :
249 11133 : nparticle = SIZE(r_last_update_pbc)
250 11133 : DO ilist = 1, n_loc_size
251 10974 : kparticle = full_loc_list(2, ilist)
252 10974 : IF (kparticle == jparticle) CYCLE
253 10815 : rjk(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(kparticle)%r(:), cell)
254 43260 : drjk = DOT_PRODUCT(rjk, rjk)
255 10815 : IF (drjk > rab2_max) CYCLE
256 1431 : drjk = SQRT(drjk)
257 5724 : rjk_hat(:) = rjk(:)/drjk
258 5724 : costheta = DOT_PRODUCT(rji, rjk)/(drji*drjk)
259 1431 : IF (costheta < -1.0_dp) costheta = -1.0_dp
260 1431 : IF (costheta > +1.0_dp) costheta = +1.0_dp
261 :
262 5724 : dcosdri(:) = (1.0_dp/(drji))*(rjk_hat(:) - costheta*rji_hat(:))
263 5724 : dcosdrk(:) = (1.0_dp/(drjk))*(rji_hat(:) - costheta*rjk_hat(:))
264 5724 : dcosdrj(:) = -(dcosdri(:) + dcosdrk(:))
265 :
266 1431 : theta = ACOS(costheta)
267 1431 : dcos_thetahalf = -1.0_dp/(2.0_dp*SQRT(1 - costheta**2))
268 : d_expterm = -2.0_dp*F*COS(theta/2.0_dp)*SIN(theta/2.0_dp) &
269 1431 : *EXP(F*(COS(theta/2.0_dp))**2)
270 :
271 5724 : dri = d_expterm*dcos_thetahalf*dcosdri
272 :
273 5724 : drj = d_expterm*dcos_thetahalf*dcosdrj
274 :
275 5724 : drk = d_expterm*dcos_thetahalf*dcosdrk
276 :
277 1431 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + prefactor*dri(1)
278 1431 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + prefactor*dri(2)
279 1431 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + prefactor*dri(3)
280 :
281 1431 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) + prefactor*drj(1)
282 1431 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) + prefactor*drj(2)
283 1431 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) + prefactor*drj(3)
284 :
285 1431 : f_nonbond(1, kparticle) = f_nonbond(1, kparticle) + prefactor*drk(1)
286 1431 : f_nonbond(2, kparticle) = f_nonbond(2, kparticle) + prefactor*drk(2)
287 1431 : f_nonbond(3, kparticle) = f_nonbond(3, kparticle) + prefactor*drk(3)
288 :
289 1590 : IF (use_virial) THEN
290 : CALL cp_abort(__LOCATION__, &
291 : "using virial with Siepmann-Sprik"// &
292 0 : " not implemented")
293 : END IF
294 : END DO
295 159 : END SUBROUTINE siep_a_ij_d
296 : ! **************************************************************************************************
297 : !> \brief Phi_ij corresponds to Equation (3) in J. Chem. Phys., Vol. 102, p.511
298 : !> \param siepmann ...
299 : !> \param r_last_update_pbc ...
300 : !> \param iparticle ...
301 : !> \param jparticle ...
302 : !> \param cell_v ...
303 : !> \param cell ...
304 : !> \param rcutsq ...
305 : !> \param particle_set ...
306 : !> \param nr_oh ...
307 : !> \param nr_h3o ...
308 : !> \param nr_o ...
309 : !> \return ...
310 : !> \par History
311 : !> Using a local list of neighbors
312 : ! **************************************************************************************************
313 636 : FUNCTION siep_Phi_ij(siepmann, r_last_update_pbc, iparticle, jparticle, &
314 : cell_v, cell, rcutsq, particle_set, nr_oh, nr_h3o, nr_o)
315 : TYPE(siepmann_pot_type), POINTER :: siepmann
316 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
317 : INTEGER, INTENT(IN) :: iparticle, jparticle
318 : REAL(KIND=dp), DIMENSION(3) :: cell_v
319 : TYPE(cell_type), POINTER :: cell
320 : REAL(KIND=dp), INTENT(IN) :: rcutsq
321 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
322 : INTEGER, INTENT(INOUT), OPTIONAL :: nr_oh, nr_h3o, nr_o
323 : REAL(KIND=dp) :: siep_Phi_ij
324 :
325 : CHARACTER(LEN=2) :: element_symbol
326 : INTEGER :: count_h, iatom, index_h1, index_h2, natom
327 : REAL(KIND=dp) :: cosphi, drih, drix, drji, h_max_dist, &
328 : rab2_max, rih(3), rih1(3), rih2(3), &
329 : rix(3), rji(3)
330 :
331 636 : siep_Phi_ij = 0.0_dp
332 636 : count_h = 0
333 636 : index_h1 = 0
334 636 : index_h2 = 0
335 636 : rab2_max = rcutsq
336 636 : h_max_dist = 2.27_dp ! 1.2 angstrom
337 636 : natom = SIZE(particle_set)
338 : CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
339 636 : element_symbol=element_symbol)
340 636 : IF (element_symbol /= "O") RETURN
341 1908 : rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
342 1908 : drji = SQRT(DOT_PRODUCT(rji, rji))
343 :
344 281817 : DO iatom = 1, natom
345 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
346 281340 : element_symbol=element_symbol)
347 281340 : IF (element_symbol /= "H") CYCLE
348 11214 : rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell)
349 44856 : drih = SQRT(DOT_PRODUCT(rih, rih))
350 11214 : IF (drih >= h_max_dist) CYCLE
351 927 : count_h = count_h + 1
352 1404 : IF (count_h == 1) THEN
353 : index_h1 = iatom
354 468 : ELSEIF (count_h == 2) THEN
355 405 : index_h2 = iatom
356 : END IF
357 : END DO
358 :
359 477 : IF (count_h == 0) THEN
360 18 : IF (siepmann%allow_o_formation) THEN
361 18 : IF (PRESENT(nr_o)) nr_o = nr_o + 1
362 : siep_Phi_ij = 0.0_dp
363 : ELSE
364 0 : CPABORT("No H atoms for O found")
365 : END IF
366 459 : ELSEIF (count_h == 1) THEN
367 54 : IF (siepmann%allow_oh_formation) THEN
368 54 : IF (PRESENT(nr_oh)) nr_oh = nr_oh + 1
369 : siep_Phi_ij = 0.0_dp
370 : ELSE
371 0 : CPABORT("Only one H atom of O atom found")
372 : END IF
373 405 : ELSEIF (count_h == 3) THEN
374 63 : IF (siepmann%allow_h3o_formation) THEN
375 63 : IF (PRESENT(nr_h3o)) nr_h3o = nr_h3o + 1
376 : siep_Phi_ij = 0.0_dp
377 : ELSE
378 0 : CPABORT("Three H atoms for O atom found")
379 : END IF
380 342 : ELSEIF (count_h > 3) THEN
381 0 : CPABORT("Error in Siepmann-Sprik part: too many H atoms for O")
382 : END IF
383 :
384 477 : IF (count_h == 2) THEN
385 : !dipole vector rix of the H2O molecule
386 342 : rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
387 342 : rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
388 1368 : rix(:) = rih1(:) + rih2(:)
389 1368 : drix = SQRT(DOT_PRODUCT(rix, rix))
390 1368 : cosphi = DOT_PRODUCT(rji, rix)/(drji*drix)
391 342 : IF (cosphi < -1.0_dp) cosphi = -1.0_dp
392 342 : IF (cosphi > +1.0_dp) cosphi = +1.0_dp
393 342 : siep_Phi_ij = EXP(-8.0_dp*((cosphi - 1)/4.0_dp)**4)
394 : END IF
395 : END FUNCTION siep_Phi_ij
396 :
397 : ! **************************************************************************************************
398 : !> \brief derivative of Phi_ij
399 : !> \param siepmann ...
400 : !> \param r_last_update_pbc ...
401 : !> \param iparticle ...
402 : !> \param jparticle ...
403 : !> \param f_nonbond ...
404 : !> \param prefactor ...
405 : !> \param cell_v ...
406 : !> \param cell ...
407 : !> \param rcutsq ...
408 : !> \param use_virial ...
409 : !> \param particle_set ...
410 : !> \par History
411 : !> Using a local list of neighbors
412 : ! **************************************************************************************************
413 159 : SUBROUTINE siep_Phi_ij_d(siepmann, r_last_update_pbc, iparticle, jparticle, f_nonbond, &
414 : prefactor, cell_v, cell, rcutsq, use_virial, particle_set)
415 : TYPE(siepmann_pot_type), POINTER :: siepmann
416 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
417 : INTEGER, INTENT(IN) :: iparticle, jparticle
418 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond
419 : REAL(KIND=dp), INTENT(IN) :: prefactor
420 : REAL(KIND=dp), DIMENSION(3) :: cell_v
421 : TYPE(cell_type), POINTER :: cell
422 : REAL(KIND=dp), INTENT(IN) :: rcutsq
423 : LOGICAL, INTENT(IN) :: use_virial
424 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
425 :
426 : CHARACTER(LEN=2) :: element_symbol
427 : INTEGER :: count_h, iatom, index_h1, index_h2, natom
428 : REAL(KIND=dp) :: cosphi, dphi, drih, drix, drji, &
429 : h_max_dist, Phi_ij, rab2_max
430 : REAL(KIND=dp), DIMENSION(3) :: dcosdrh, dcosdri, dcosdrj, drh, dri, &
431 : drj, rih, rih1, rih2, rix, rix_hat, &
432 : rji, rji_hat
433 :
434 159 : count_h = 0
435 159 : index_h1 = 0
436 159 : index_h2 = 0
437 159 : rab2_max = rcutsq
438 159 : h_max_dist = 2.27_dp ! 1.2 angstrom
439 159 : natom = SIZE(particle_set)
440 : Phi_ij = siep_Phi_ij(siepmann, r_last_update_pbc, iparticle, jparticle, &
441 : cell_v, cell, rcutsq, &
442 159 : particle_set)
443 636 : rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
444 636 : drji = SQRT(DOT_PRODUCT(rji, rji))
445 636 : rji_hat(:) = rji(:)/drji
446 :
447 93939 : DO iatom = 1, natom
448 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
449 93780 : element_symbol=element_symbol)
450 93780 : IF (element_symbol /= "H") CYCLE
451 3738 : rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell)
452 14952 : drih = SQRT(DOT_PRODUCT(rih, rih))
453 3738 : IF (drih >= h_max_dist) CYCLE
454 309 : count_h = count_h + 1
455 468 : IF (count_h == 1) THEN
456 : index_h1 = iatom
457 156 : ELSEIF (count_h == 2) THEN
458 135 : index_h2 = iatom
459 : END IF
460 : END DO
461 :
462 159 : IF (count_h == 0 .AND. .NOT. siepmann%allow_o_formation) THEN
463 0 : CPABORT("No H atoms for O found")
464 159 : ELSEIF (count_h == 1 .AND. .NOT. siepmann%allow_oh_formation) THEN
465 0 : CPABORT("Only one H atom for O atom found")
466 159 : ELSEIF (count_h == 3 .AND. .NOT. siepmann%allow_h3o_formation) THEN
467 0 : CPABORT("Three H atoms for O atom found")
468 159 : ELSEIF (count_h > 3) THEN
469 0 : CPABORT("Error in Siepmann-Sprik part: too many H atoms for O")
470 : END IF
471 :
472 159 : IF (count_h == 2) THEN
473 : !dipole vector rix of the H2O molecule
474 114 : rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
475 114 : rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
476 456 : rix(:) = rih1(:) + rih2(:)
477 456 : drix = SQRT(DOT_PRODUCT(rix, rix))
478 456 : rix_hat(:) = rix(:)/drix
479 456 : cosphi = DOT_PRODUCT(rji, rix)/(drji*drix)
480 114 : IF (cosphi < -1.0_dp) cosphi = -1.0_dp
481 114 : IF (cosphi > +1.0_dp) cosphi = +1.0_dp
482 :
483 456 : dcosdrj(:) = (1.0_dp/(drji))*(-rix_hat(:) + cosphi*rji_hat(:))
484 : ! for H atoms:
485 456 : dcosdrh(:) = (1.0_dp/(drix))*(rji_hat(:) - cosphi*rix_hat(:))
486 456 : dcosdri(:) = -dcosdrj - 2.0_dp*dcosdrh
487 :
488 114 : dphi = Phi_ij*(-8.0_dp)*((cosphi - 1)/4.0_dp)**3
489 :
490 456 : dri = dphi*dcosdri
491 456 : drj = dphi*dcosdrj
492 456 : drh = dphi*dcosdrh
493 :
494 114 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + prefactor*dri(1)
495 114 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + prefactor*dri(2)
496 114 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + prefactor*dri(3)
497 :
498 114 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) + prefactor*drj(1)
499 114 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) + prefactor*drj(2)
500 114 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) + prefactor*drj(3)
501 :
502 114 : f_nonbond(1, index_h1) = f_nonbond(1, index_h1) + prefactor*drh(1)
503 114 : f_nonbond(2, index_h1) = f_nonbond(2, index_h1) + prefactor*drh(2)
504 114 : f_nonbond(3, index_h1) = f_nonbond(3, index_h1) + prefactor*drh(3)
505 :
506 114 : f_nonbond(1, index_h2) = f_nonbond(1, index_h2) + prefactor*drh(1)
507 114 : f_nonbond(2, index_h2) = f_nonbond(2, index_h2) + prefactor*drh(2)
508 114 : f_nonbond(3, index_h2) = f_nonbond(3, index_h2) + prefactor*drh(3)
509 :
510 114 : IF (use_virial) THEN
511 : CALL cp_abort(__LOCATION__, &
512 : "using virial with Siepmann-Sprik"// &
513 0 : " not implemented")
514 : END IF
515 :
516 : END IF
517 159 : END SUBROUTINE siep_Phi_ij_d
518 :
519 : ! **************************************************************************************************
520 : !> \brief forces generated by the three-body term
521 : !> \param siepmann ...
522 : !> \param r_last_update_pbc ...
523 : !> \param cell_v ...
524 : !> \param n_loc_size ...
525 : !> \param full_loc_list ...
526 : !> \param iparticle ...
527 : !> \param jparticle ...
528 : !> \param f_nonbond ...
529 : !> \param use_virial ...
530 : !> \param rcutsq ...
531 : !> \param cell ...
532 : !> \param particle_set ...
533 : !> \par History
534 : !> Using a local list of neighbors
535 : ! **************************************************************************************************
536 318 : SUBROUTINE siepmann_forces_v3(siepmann, r_last_update_pbc, cell_v, n_loc_size, &
537 318 : full_loc_list, iparticle, jparticle, f_nonbond, &
538 : use_virial, rcutsq, cell, particle_set)
539 : TYPE(siepmann_pot_type), POINTER :: siepmann
540 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
541 : REAL(KIND=dp), DIMENSION(3) :: cell_v
542 : INTEGER, INTENT(IN) :: n_loc_size
543 : INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list
544 : INTEGER, INTENT(IN) :: iparticle, jparticle
545 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond
546 : LOGICAL, INTENT(IN) :: use_virial
547 : REAL(KIND=dp), INTENT(IN) :: rcutsq
548 : TYPE(cell_type), POINTER :: cell
549 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
550 :
551 : CHARACTER(LEN=2) :: element_symbol
552 : REAL(KIND=dp) :: a_ij, beta, drji, E, f2, f2_d, f_A1, &
553 : f_A2, fac, prefactor, rji(3), &
554 : rji_hat(3)
555 :
556 318 : beta = siepmann%beta
557 318 : E = siepmann%E
558 :
559 : CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
560 318 : element_symbol=element_symbol)
561 318 : IF (element_symbol /= "O") RETURN
562 :
563 636 : rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
564 636 : drji = SQRT(DOT_PRODUCT(rji, rji))
565 636 : rji_hat(:) = rji(:)/drji
566 :
567 159 : fac = -1.0_dp !gradient to force
568 : a_ij = siep_a_ij(siepmann, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
569 159 : full_loc_list, cell_v, rcutsq, particle_set, cell)
570 159 : f2 = siep_f2(siepmann, drji)
571 159 : f2_d = siep_f2_d(siepmann, drji)
572 :
573 : ! Lets do the f_A1 piece derivative of f2
574 159 : f_A1 = E*f2_d*drji**(-beta)*a_ij*fac*(1.0_dp/drji)
575 159 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A1*rji(1)
576 159 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A1*rji(2)
577 159 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A1*rji(3)
578 159 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A1*rji(1)
579 159 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A1*rji(2)
580 159 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A1*rji(3)
581 :
582 159 : IF (use_virial) THEN
583 : CALL cp_abort(__LOCATION__, &
584 : "using virial with Siepmann-Sprik"// &
585 0 : " not implemented")
586 : END IF
587 :
588 : ! Lets do the f_A2 piece derivative of rji**(-beta)
589 159 : f_A2 = E*f2*(-beta)*drji**(-beta - 1)*a_ij*fac*(1.0_dp/drji)
590 159 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A2*rji(1)
591 159 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A2*rji(2)
592 159 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A2*rji(3)
593 159 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A2*rji(1)
594 159 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A2*rji(2)
595 159 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A2*rji(3)
596 :
597 159 : IF (use_virial) THEN
598 : CALL cp_abort(__LOCATION__, &
599 : "using virial with Siepmann-Sprik"// &
600 0 : " not implemented")
601 : END IF
602 :
603 : ! Lets do the f_A3 piece derivative: of a_ij
604 159 : prefactor = E*f2*drji**(-beta)*fac
605 : CALL siep_a_ij_d(siepmann, r_last_update_pbc, iparticle, jparticle, f_nonbond, &
606 : prefactor, n_loc_size, full_loc_list, cell_v, &
607 159 : cell, rcutsq, use_virial)
608 :
609 : END SUBROUTINE siepmann_forces_v3
610 :
611 : ! **************************************************************************************************
612 : !> \brief forces generated by the dipole term
613 : !> \param siepmann ...
614 : !> \param r_last_update_pbc ...
615 : !> \param cell_v ...
616 : !> \param cell ...
617 : !> \param iparticle ...
618 : !> \param jparticle ...
619 : !> \param f_nonbond ...
620 : !> \param use_virial ...
621 : !> \param rcutsq ...
622 : !> \param particle_set ...
623 : !> \par History
624 : !> Using a local list of neighbors
625 : ! **************************************************************************************************
626 318 : SUBROUTINE siepmann_forces_v2(siepmann, r_last_update_pbc, cell_v, cell, &
627 318 : iparticle, jparticle, f_nonbond, use_virial, rcutsq, particle_set)
628 : TYPE(siepmann_pot_type), POINTER :: siepmann
629 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
630 : REAL(KIND=dp), DIMENSION(3) :: cell_v
631 : TYPE(cell_type), POINTER :: cell
632 : INTEGER, INTENT(IN) :: iparticle, jparticle
633 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond
634 : LOGICAL, INTENT(IN) :: use_virial
635 : REAL(KIND=dp), INTENT(IN) :: rcutsq
636 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
637 :
638 : CHARACTER(LEN=2) :: element_symbol
639 : REAL(KIND=dp) :: D, drji, f2, f2_d, f_A1, f_A2, fac, &
640 : Phi_ij, prefactor, rji(3)
641 :
642 318 : D = siepmann%D
643 :
644 : CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
645 318 : element_symbol=element_symbol)
646 318 : IF (element_symbol /= "O") RETURN
647 :
648 636 : rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
649 636 : drji = SQRT(DOT_PRODUCT(rji, rji))
650 :
651 159 : fac = -1.0_dp
652 : Phi_ij = siep_Phi_ij(siepmann, r_last_update_pbc, iparticle, jparticle, &
653 159 : cell_v, cell, rcutsq, particle_set)
654 159 : f2 = siep_f2(siepmann, drji)
655 159 : f2_d = siep_f2_d(siepmann, drji)
656 :
657 : ! Lets do the f_A1 piece derivative of f2
658 159 : f_A1 = -D*f2_d*drji**(-3)*Phi_ij*fac*(1.0_dp/drji)
659 159 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A1*rji(1)
660 159 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A1*rji(2)
661 159 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A1*rji(3)
662 159 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A1*rji(1)
663 159 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A1*rji(2)
664 159 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A1*rji(3)
665 :
666 159 : IF (use_virial) THEN
667 : CALL cp_abort(__LOCATION__, &
668 : "using virial with Siepmann-Sprik"// &
669 0 : " not implemented")
670 : END IF
671 :
672 : ! ! Lets do the f_A2 piece derivative of rji**(-3)
673 159 : f_A2 = -D*f2*(-3.0_dp)*drji**(-4)*Phi_ij*fac*(1.0_dp/drji)
674 159 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A2*rji(1)
675 159 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A2*rji(2)
676 159 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A2*rji(3)
677 159 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A2*rji(1)
678 159 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A2*rji(2)
679 159 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A2*rji(3)
680 :
681 159 : IF (use_virial) THEN
682 : CALL cp_abort(__LOCATION__, &
683 : "using virial with Siepmann-Sprik"// &
684 0 : " not implemented")
685 : END IF
686 :
687 : ! Lets do the f_A3 piece derivative: of Phi_ij
688 159 : prefactor = -D*f2*drji**(-3)*fac
689 : CALL siep_Phi_ij_d(siepmann, r_last_update_pbc, iparticle, jparticle, f_nonbond, &
690 : prefactor, cell_v, cell, &
691 159 : rcutsq, use_virial, particle_set)
692 :
693 : END SUBROUTINE siepmann_forces_v2
694 :
695 : ! **************************************************************************************************
696 : !> \brief ...
697 : !> \param nonbonded ...
698 : !> \param potparm ...
699 : !> \param glob_loc_list ...
700 : !> \param glob_cell_v ...
701 : !> \param glob_loc_list_a ...
702 : !> \param cell ...
703 : !> \par History
704 : ! **************************************************************************************************
705 42 : SUBROUTINE setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, &
706 : glob_loc_list_a, cell)
707 : TYPE(fist_neighbor_type), POINTER :: nonbonded
708 : TYPE(pair_potential_pp_type), POINTER :: potparm
709 : INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list
710 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v
711 : INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a
712 : TYPE(cell_type), POINTER :: cell
713 :
714 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_siepmann_arrays'
715 :
716 : INTEGER :: handle, i, iend, igrp, ikind, ilist, &
717 : ipair, istart, jkind, nkinds, npairs, &
718 : npairs_tot
719 42 : INTEGER, DIMENSION(:), POINTER :: work_list, work_list2
720 42 : INTEGER, DIMENSION(:, :), POINTER :: list
721 : REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi
722 42 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rwork_list
723 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
724 : TYPE(pair_potential_single_type), POINTER :: pot
725 :
726 0 : CPASSERT(.NOT. ASSOCIATED(glob_loc_list))
727 42 : CPASSERT(.NOT. ASSOCIATED(glob_loc_list_a))
728 42 : CPASSERT(.NOT. ASSOCIATED(glob_cell_v))
729 42 : CALL timeset(routineN, handle)
730 42 : npairs_tot = 0
731 42 : nkinds = SIZE(potparm%pot, 1)
732 1176 : DO ilist = 1, nonbonded%nlists
733 1134 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
734 1134 : npairs = neighbor_kind_pair%npairs
735 1134 : IF (npairs == 0) CYCLE
736 1836 : Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
737 1416 : istart = neighbor_kind_pair%grp_kind_start(igrp)
738 1416 : iend = neighbor_kind_pair%grp_kind_end(igrp)
739 1416 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
740 1416 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
741 1416 : pot => potparm%pot(ikind, jkind)%pot
742 1416 : npairs = iend - istart + 1
743 1416 : IF (pot%no_mb) CYCLE
744 1794 : DO i = 1, SIZE(pot%type)
745 1746 : IF (pot%type(i) == siepmann_type) npairs_tot = npairs_tot + npairs
746 : END DO
747 : END DO Kind_Group_Loop1
748 : END DO
749 126 : ALLOCATE (work_list(npairs_tot))
750 84 : ALLOCATE (work_list2(npairs_tot))
751 126 : ALLOCATE (glob_loc_list(2, npairs_tot))
752 126 : ALLOCATE (glob_cell_v(3, npairs_tot))
753 : ! Fill arrays with data
754 42 : npairs_tot = 0
755 1176 : DO ilist = 1, nonbonded%nlists
756 1134 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
757 1134 : npairs = neighbor_kind_pair%npairs
758 1134 : IF (npairs == 0) CYCLE
759 1836 : Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
760 1416 : istart = neighbor_kind_pair%grp_kind_start(igrp)
761 1416 : iend = neighbor_kind_pair%grp_kind_end(igrp)
762 1416 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
763 1416 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
764 1416 : list => neighbor_kind_pair%list
765 5664 : cvi = neighbor_kind_pair%cell_vector
766 1416 : pot => potparm%pot(ikind, jkind)%pot
767 1416 : npairs = iend - istart + 1
768 1416 : IF (pot%no_mb) CYCLE
769 4290 : cell_v = MATMUL(cell%hmat, cvi)
770 1794 : DO i = 1, SIZE(pot%type)
771 : ! SIEPMANN
772 1746 : IF (pot%type(i) == siepmann_type) THEN
773 36786 : DO ipair = 1, npairs
774 218736 : glob_loc_list(:, npairs_tot + ipair) = list(:, istart - 1 + ipair)
775 146154 : glob_cell_v(1:3, npairs_tot + ipair) = cell_v(1:3)
776 : END DO
777 330 : npairs_tot = npairs_tot + npairs
778 : END IF
779 : END DO
780 : END DO Kind_Group_Loop2
781 : END DO
782 : ! Order the arrays w.r.t. the first index of glob_loc_list
783 42 : CALL sort(glob_loc_list(1, :), npairs_tot, work_list)
784 36498 : DO ipair = 1, npairs_tot
785 36498 : work_list2(ipair) = glob_loc_list(2, work_list(ipair))
786 : END DO
787 72996 : glob_loc_list(2, :) = work_list2
788 42 : DEALLOCATE (work_list2)
789 126 : ALLOCATE (rwork_list(3, npairs_tot))
790 36498 : DO ipair = 1, npairs_tot
791 291690 : rwork_list(:, ipair) = glob_cell_v(:, work_list(ipair))
792 : END DO
793 291732 : glob_cell_v = rwork_list
794 42 : DEALLOCATE (rwork_list)
795 42 : DEALLOCATE (work_list)
796 126 : ALLOCATE (glob_loc_list_a(npairs_tot))
797 72996 : glob_loc_list_a = glob_loc_list(1, :)
798 42 : CALL timestop(handle)
799 84 : END SUBROUTINE setup_siepmann_arrays
800 :
801 : ! **************************************************************************************************
802 : !> \brief ...
803 : !> \param glob_loc_list ...
804 : !> \param glob_cell_v ...
805 : !> \param glob_loc_list_a ...
806 : ! **************************************************************************************************
807 42 : SUBROUTINE destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
808 : INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list
809 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v
810 : INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a
811 :
812 42 : IF (ASSOCIATED(glob_loc_list)) THEN
813 42 : DEALLOCATE (glob_loc_list)
814 : END IF
815 42 : IF (ASSOCIATED(glob_loc_list_a)) THEN
816 42 : DEALLOCATE (glob_loc_list_a)
817 : END IF
818 42 : IF (ASSOCIATED(glob_cell_v)) THEN
819 42 : DEALLOCATE (glob_cell_v)
820 : END IF
821 :
822 42 : END SUBROUTINE destroy_siepmann_arrays
823 :
824 : ! **************************************************************************************************
825 : !> \brief prints the number of OH- ions or H3O+ ions near surface
826 : !> \param nr_ions number of ions
827 : !> \param mm_section ...
828 : !> \param para_env ...
829 : !> \param print_oh flag indicating if number OH- is printed
830 : !> \param print_h3o flag indicating if number H3O+ is printed
831 : !> \param print_o flag indicating if number O^(2-) is printed
832 : ! **************************************************************************************************
833 63 : SUBROUTINE print_nr_ions_siepmann(nr_ions, mm_section, para_env, print_oh, &
834 : print_h3o, print_o)
835 : INTEGER, INTENT(INOUT) :: nr_ions
836 : TYPE(section_vals_type), POINTER :: mm_section
837 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
838 : LOGICAL, INTENT(IN) :: print_oh, print_h3o, print_o
839 :
840 : INTEGER :: iw
841 : TYPE(cp_logger_type), POINTER :: logger
842 :
843 63 : NULLIFY (logger)
844 :
845 63 : CALL para_env%sum(nr_ions)
846 63 : logger => cp_get_default_logger()
847 :
848 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%PROGRAM_RUN_INFO", &
849 63 : extension=".mmLog")
850 :
851 63 : IF (iw > 0 .AND. nr_ions > 0 .AND. print_oh) THEN
852 6 : WRITE (iw, '(/,A,T71,I10,/)') " SIEPMANN: number of OH- ions at surface", nr_ions
853 : END IF
854 63 : IF (iw > 0 .AND. nr_ions > 0 .AND. print_h3o) THEN
855 6 : WRITE (iw, '(/,A,T71,I10,/)') " SIEPMANN: number of H3O+ ions at surface", nr_ions
856 : END IF
857 63 : IF (iw > 0 .AND. nr_ions > 0 .AND. print_o) THEN
858 3 : WRITE (iw, '(/,A,T71,I10,/)') " SIEPMANN: number of O^2- ions at surface", nr_ions
859 : END IF
860 :
861 63 : CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%PROGRAM_RUN_INFO")
862 :
863 63 : END SUBROUTINE print_nr_ions_siepmann
864 :
865 : END MODULE manybody_siepmann
866 :
|
