(git:34ef472)
manybody_siepmann.F
Go to the documentation of this file.
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  SUBROUTINE siepmann_energy(pot_loc, siepmann, r_last_update_pbc, atom_a, atom_b, &
69  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  cell)
89  phi_ij = siep_phi_ij(siepmann, r_last_update_pbc, atom_a, atom_b, &
90  cell_v, cell, siepmann%rcutsq, particle_set, nr_oh, nr_h3o, nr_o)
91  f2 = siep_f2(siepmann, drij)
92  d = siepmann%D
93  e = siepmann%E
94 
95  !two-body part --> dipole term
96  pot_loc_v2 = -d*f2*drij**(-3)*phi_ij
97 
98  !three-body part
99  pot_loc_v3 = e*f2*drij**(-siepmann%beta)*a_ij
100 
101  pot_loc = pot_loc_v2 + pot_loc_v3
102 
103  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  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  rcut = sqrt(siepmann%rcutsq)
119  siep_f2 = 0.0_dp
120  IF (r < rcut) THEN
121  siep_f2 = exp(siepmann%B/(r - rcut))
122  END IF
123  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  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  rcut = sqrt(siepmann%rcutsq)
139  b = siepmann%B
140  siep_f2_d = 0.0_dp
141  IF (r < rcut) THEN
142  siep_f2_d = -b*exp(b/(r - rcut))/(r - rcut)**2
143  END IF
144 
145  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  FUNCTION siep_a_ij(siepmann, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
165  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  siep_a_ij = 0.0_dp
182  rab2_max = rcutsq
183  f = siepmann%F
184  CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
185  element_symbol=element_symbol)
186  IF (element_symbol /= "O") RETURN
187  rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
188  drji = sqrt(dot_product(rji, rji))
189  DO ilist = 1, n_loc_size
190  kparticle = full_loc_list(2, ilist)
191  IF (kparticle == jparticle) cycle
192  rjk(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(kparticle)%r(:), cell)
193  drjk = dot_product(rjk, rjk)
194  IF (drjk > rab2_max) cycle
195  drjk = sqrt(drjk)
196  costheta = dot_product(rji, rjk)/(drji*drjk)
197  IF (costheta < -1.0_dp) costheta = -1.0_dp
198  IF (costheta > +1.0_dp) costheta = +1.0_dp
199  theta = acos(costheta)
200  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  SUBROUTINE siep_a_ij_d(siepmann, r_last_update_pbc, iparticle, jparticle, f_nonbond, &
222  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  rab2_max = rcutsq
243  f = siepmann%F
244 
245  rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
246  drji = sqrt(dot_product(rji, rji))
247  rji_hat(:) = rji(:)/drji
248 
249  nparticle = SIZE(r_last_update_pbc)
250  DO ilist = 1, n_loc_size
251  kparticle = full_loc_list(2, ilist)
252  IF (kparticle == jparticle) cycle
253  rjk(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(kparticle)%r(:), cell)
254  drjk = dot_product(rjk, rjk)
255  IF (drjk > rab2_max) cycle
256  drjk = sqrt(drjk)
257  rjk_hat(:) = rjk(:)/drjk
258  costheta = dot_product(rji, rjk)/(drji*drjk)
259  IF (costheta < -1.0_dp) costheta = -1.0_dp
260  IF (costheta > +1.0_dp) costheta = +1.0_dp
261 
262  dcosdri(:) = (1.0_dp/(drji))*(rjk_hat(:) - costheta*rji_hat(:))
263  dcosdrk(:) = (1.0_dp/(drjk))*(rji_hat(:) - costheta*rjk_hat(:))
264  dcosdrj(:) = -(dcosdri(:) + dcosdrk(:))
265 
266  theta = acos(costheta)
267  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  *exp(f*(cos(theta/2.0_dp))**2)
270 
271  dri = d_expterm*dcos_thetahalf*dcosdri
272 
273  drj = d_expterm*dcos_thetahalf*dcosdrj
274 
275  drk = d_expterm*dcos_thetahalf*dcosdrk
276 
277  f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + prefactor*dri(1)
278  f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + prefactor*dri(2)
279  f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + prefactor*dri(3)
280 
281  f_nonbond(1, jparticle) = f_nonbond(1, jparticle) + prefactor*drj(1)
282  f_nonbond(2, jparticle) = f_nonbond(2, jparticle) + prefactor*drj(2)
283  f_nonbond(3, jparticle) = f_nonbond(3, jparticle) + prefactor*drj(3)
284 
285  f_nonbond(1, kparticle) = f_nonbond(1, kparticle) + prefactor*drk(1)
286  f_nonbond(2, kparticle) = f_nonbond(2, kparticle) + prefactor*drk(2)
287  f_nonbond(3, kparticle) = f_nonbond(3, kparticle) + prefactor*drk(3)
288 
289  IF (use_virial) THEN
290  CALL cp_abort(__location__, &
291  "using virial with Siepmann-Sprik"// &
292  " not implemented")
293  END IF
294  END DO
295  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  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  siep_phi_ij = 0.0_dp
332  count_h = 0
333  index_h1 = 0
334  index_h2 = 0
335  rab2_max = rcutsq
336  h_max_dist = 2.27_dp ! 1.2 angstrom
337  natom = SIZE(particle_set)
338  CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
339  element_symbol=element_symbol)
340  IF (element_symbol /= "O") RETURN
341  rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
342  drji = sqrt(dot_product(rji, rji))
343 
344  DO iatom = 1, natom
345  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
346  element_symbol=element_symbol)
347  IF (element_symbol /= "H") cycle
348  rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell)
349  drih = sqrt(dot_product(rih, rih))
350  IF (drih >= h_max_dist) cycle
351  count_h = count_h + 1
352  IF (count_h == 1) THEN
353  index_h1 = iatom
354  ELSEIF (count_h == 2) THEN
355  index_h2 = iatom
356  END IF
357  END DO
358 
359  IF (count_h == 0) THEN
360  IF (siepmann%allow_o_formation) THEN
361  IF (PRESENT(nr_o)) nr_o = nr_o + 1
362  siep_phi_ij = 0.0_dp
363  ELSE
364  cpabort("No H atoms for O found")
365  END IF
366  ELSEIF (count_h == 1) THEN
367  IF (siepmann%allow_oh_formation) THEN
368  IF (PRESENT(nr_oh)) nr_oh = nr_oh + 1
369  siep_phi_ij = 0.0_dp
370  ELSE
371  cpabort("Only one H atom of O atom found")
372  END IF
373  ELSEIF (count_h == 3) THEN
374  IF (siepmann%allow_h3o_formation) THEN
375  IF (PRESENT(nr_h3o)) nr_h3o = nr_h3o + 1
376  siep_phi_ij = 0.0_dp
377  ELSE
378  cpabort("Three H atoms for O atom found")
379  END IF
380  ELSEIF (count_h > 3) THEN
381  cpabort("Error in Siepmann-Sprik part: too many H atoms for O")
382  END IF
383 
384  IF (count_h == 2) THEN
385  !dipole vector rix of the H2O molecule
386  rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
387  rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
388  rix(:) = rih1(:) + rih2(:)
389  drix = sqrt(dot_product(rix, rix))
390  cosphi = dot_product(rji, rix)/(drji*drix)
391  IF (cosphi < -1.0_dp) cosphi = -1.0_dp
392  IF (cosphi > +1.0_dp) cosphi = +1.0_dp
393  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  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  count_h = 0
435  index_h1 = 0
436  index_h2 = 0
437  rab2_max = rcutsq
438  h_max_dist = 2.27_dp ! 1.2 angstrom
439  natom = SIZE(particle_set)
440  phi_ij = siep_phi_ij(siepmann, r_last_update_pbc, iparticle, jparticle, &
441  cell_v, cell, rcutsq, &
442  particle_set)
443  rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
444  drji = sqrt(dot_product(rji, rji))
445  rji_hat(:) = rji(:)/drji
446 
447  DO iatom = 1, natom
448  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
449  element_symbol=element_symbol)
450  IF (element_symbol /= "H") cycle
451  rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell)
452  drih = sqrt(dot_product(rih, rih))
453  IF (drih >= h_max_dist) cycle
454  count_h = count_h + 1
455  IF (count_h == 1) THEN
456  index_h1 = iatom
457  ELSEIF (count_h == 2) THEN
458  index_h2 = iatom
459  END IF
460  END DO
461 
462  IF (count_h == 0 .AND. .NOT. siepmann%allow_o_formation) THEN
463  cpabort("No H atoms for O found")
464  ELSEIF (count_h == 1 .AND. .NOT. siepmann%allow_oh_formation) THEN
465  cpabort("Only one H atom for O atom found")
466  ELSEIF (count_h == 3 .AND. .NOT. siepmann%allow_h3o_formation) THEN
467  cpabort("Three H atoms for O atom found")
468  ELSEIF (count_h > 3) THEN
469  cpabort("Error in Siepmann-Sprik part: too many H atoms for O")
470  END IF
471 
472  IF (count_h == 2) THEN
473  !dipole vector rix of the H2O molecule
474  rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
475  rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
476  rix(:) = rih1(:) + rih2(:)
477  drix = sqrt(dot_product(rix, rix))
478  rix_hat(:) = rix(:)/drix
479  cosphi = dot_product(rji, rix)/(drji*drix)
480  IF (cosphi < -1.0_dp) cosphi = -1.0_dp
481  IF (cosphi > +1.0_dp) cosphi = +1.0_dp
482 
483  dcosdrj(:) = (1.0_dp/(drji))*(-rix_hat(:) + cosphi*rji_hat(:))
484  ! for H atoms:
485  dcosdrh(:) = (1.0_dp/(drix))*(rji_hat(:) - cosphi*rix_hat(:))
486  dcosdri(:) = -dcosdrj - 2.0_dp*dcosdrh
487 
488  dphi = phi_ij*(-8.0_dp)*((cosphi - 1)/4.0_dp)**3
489 
490  dri = dphi*dcosdri
491  drj = dphi*dcosdrj
492  drh = dphi*dcosdrh
493 
494  f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + prefactor*dri(1)
495  f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + prefactor*dri(2)
496  f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + prefactor*dri(3)
497 
498  f_nonbond(1, jparticle) = f_nonbond(1, jparticle) + prefactor*drj(1)
499  f_nonbond(2, jparticle) = f_nonbond(2, jparticle) + prefactor*drj(2)
500  f_nonbond(3, jparticle) = f_nonbond(3, jparticle) + prefactor*drj(3)
501 
502  f_nonbond(1, index_h1) = f_nonbond(1, index_h1) + prefactor*drh(1)
503  f_nonbond(2, index_h1) = f_nonbond(2, index_h1) + prefactor*drh(2)
504  f_nonbond(3, index_h1) = f_nonbond(3, index_h1) + prefactor*drh(3)
505 
506  f_nonbond(1, index_h2) = f_nonbond(1, index_h2) + prefactor*drh(1)
507  f_nonbond(2, index_h2) = f_nonbond(2, index_h2) + prefactor*drh(2)
508  f_nonbond(3, index_h2) = f_nonbond(3, index_h2) + prefactor*drh(3)
509 
510  IF (use_virial) THEN
511  CALL cp_abort(__location__, &
512  "using virial with Siepmann-Sprik"// &
513  " not implemented")
514  END IF
515 
516  END IF
517  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  SUBROUTINE siepmann_forces_v3(siepmann, r_last_update_pbc, cell_v, n_loc_size, &
537  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  beta = siepmann%beta
557  e = siepmann%E
558 
559  CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
560  element_symbol=element_symbol)
561  IF (element_symbol /= "O") RETURN
562 
563  rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
564  drji = sqrt(dot_product(rji, rji))
565  rji_hat(:) = rji(:)/drji
566 
567  fac = -1.0_dp !gradient to force
568  a_ij = siep_a_ij(siepmann, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
569  full_loc_list, cell_v, rcutsq, particle_set, cell)
570  f2 = siep_f2(siepmann, drji)
571  f2_d = siep_f2_d(siepmann, drji)
572 
573  ! Lets do the f_A1 piece derivative of f2
574  f_a1 = e*f2_d*drji**(-beta)*a_ij*fac*(1.0_dp/drji)
575  f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_a1*rji(1)
576  f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_a1*rji(2)
577  f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_a1*rji(3)
578  f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_a1*rji(1)
579  f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_a1*rji(2)
580  f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_a1*rji(3)
581 
582  IF (use_virial) THEN
583  CALL cp_abort(__location__, &
584  "using virial with Siepmann-Sprik"// &
585  " not implemented")
586  END IF
587 
588  ! Lets do the f_A2 piece derivative of rji**(-beta)
589  f_a2 = e*f2*(-beta)*drji**(-beta - 1)*a_ij*fac*(1.0_dp/drji)
590  f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_a2*rji(1)
591  f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_a2*rji(2)
592  f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_a2*rji(3)
593  f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_a2*rji(1)
594  f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_a2*rji(2)
595  f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_a2*rji(3)
596 
597  IF (use_virial) THEN
598  CALL cp_abort(__location__, &
599  "using virial with Siepmann-Sprik"// &
600  " not implemented")
601  END IF
602 
603  ! Lets do the f_A3 piece derivative: of a_ij
604  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  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  SUBROUTINE siepmann_forces_v2(siepmann, r_last_update_pbc, cell_v, cell, &
627  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  d = siepmann%D
643 
644  CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
645  element_symbol=element_symbol)
646  IF (element_symbol /= "O") RETURN
647 
648  rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v)
649  drji = sqrt(dot_product(rji, rji))
650 
651  fac = -1.0_dp
652  phi_ij = siep_phi_ij(siepmann, r_last_update_pbc, iparticle, jparticle, &
653  cell_v, cell, rcutsq, particle_set)
654  f2 = siep_f2(siepmann, drji)
655  f2_d = siep_f2_d(siepmann, drji)
656 
657  ! Lets do the f_A1 piece derivative of f2
658  f_a1 = -d*f2_d*drji**(-3)*phi_ij*fac*(1.0_dp/drji)
659  f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_a1*rji(1)
660  f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_a1*rji(2)
661  f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_a1*rji(3)
662  f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_a1*rji(1)
663  f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_a1*rji(2)
664  f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_a1*rji(3)
665 
666  IF (use_virial) THEN
667  CALL cp_abort(__location__, &
668  "using virial with Siepmann-Sprik"// &
669  " not implemented")
670  END IF
671 
672 ! ! Lets do the f_A2 piece derivative of rji**(-3)
673  f_a2 = -d*f2*(-3.0_dp)*drji**(-4)*phi_ij*fac*(1.0_dp/drji)
674  f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_a2*rji(1)
675  f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_a2*rji(2)
676  f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_a2*rji(3)
677  f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_a2*rji(1)
678  f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_a2*rji(2)
679  f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_a2*rji(3)
680 
681  IF (use_virial) THEN
682  CALL cp_abort(__location__, &
683  "using virial with Siepmann-Sprik"// &
684  " not implemented")
685  END IF
686 
687  ! Lets do the f_A3 piece derivative: of Phi_ij
688  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  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  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  INTEGER, DIMENSION(:), POINTER :: work_list, work_list2
720  INTEGER, DIMENSION(:, :), POINTER :: list
721  REAL(kind=dp), DIMENSION(3) :: cell_v, cvi
722  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  cpassert(.NOT. ASSOCIATED(glob_loc_list))
727  cpassert(.NOT. ASSOCIATED(glob_loc_list_a))
728  cpassert(.NOT. ASSOCIATED(glob_cell_v))
729  CALL timeset(routinen, handle)
730  npairs_tot = 0
731  nkinds = SIZE(potparm%pot, 1)
732  DO ilist = 1, nonbonded%nlists
733  neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
734  npairs = neighbor_kind_pair%npairs
735  IF (npairs == 0) cycle
736  kind_group_loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
737  istart = neighbor_kind_pair%grp_kind_start(igrp)
738  iend = neighbor_kind_pair%grp_kind_end(igrp)
739  ikind = neighbor_kind_pair%ij_kind(1, igrp)
740  jkind = neighbor_kind_pair%ij_kind(2, igrp)
741  pot => potparm%pot(ikind, jkind)%pot
742  npairs = iend - istart + 1
743  IF (pot%no_mb) cycle
744  DO i = 1, SIZE(pot%type)
745  IF (pot%type(i) == siepmann_type) npairs_tot = npairs_tot + npairs
746  END DO
747  END DO kind_group_loop1
748  END DO
749  ALLOCATE (work_list(npairs_tot))
750  ALLOCATE (work_list2(npairs_tot))
751  ALLOCATE (glob_loc_list(2, npairs_tot))
752  ALLOCATE (glob_cell_v(3, npairs_tot))
753  ! Fill arrays with data
754  npairs_tot = 0
755  DO ilist = 1, nonbonded%nlists
756  neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
757  npairs = neighbor_kind_pair%npairs
758  IF (npairs == 0) cycle
759  kind_group_loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
760  istart = neighbor_kind_pair%grp_kind_start(igrp)
761  iend = neighbor_kind_pair%grp_kind_end(igrp)
762  ikind = neighbor_kind_pair%ij_kind(1, igrp)
763  jkind = neighbor_kind_pair%ij_kind(2, igrp)
764  list => neighbor_kind_pair%list
765  cvi = neighbor_kind_pair%cell_vector
766  pot => potparm%pot(ikind, jkind)%pot
767  npairs = iend - istart + 1
768  IF (pot%no_mb) cycle
769  cell_v = matmul(cell%hmat, cvi)
770  DO i = 1, SIZE(pot%type)
771  ! SIEPMANN
772  IF (pot%type(i) == siepmann_type) THEN
773  DO ipair = 1, npairs
774  glob_loc_list(:, npairs_tot + ipair) = list(:, istart - 1 + ipair)
775  glob_cell_v(1:3, npairs_tot + ipair) = cell_v(1:3)
776  END DO
777  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  CALL sort(glob_loc_list(1, :), npairs_tot, work_list)
784  DO ipair = 1, npairs_tot
785  work_list2(ipair) = glob_loc_list(2, work_list(ipair))
786  END DO
787  glob_loc_list(2, :) = work_list2
788  DEALLOCATE (work_list2)
789  ALLOCATE (rwork_list(3, npairs_tot))
790  DO ipair = 1, npairs_tot
791  rwork_list(:, ipair) = glob_cell_v(:, work_list(ipair))
792  END DO
793  glob_cell_v = rwork_list
794  DEALLOCATE (rwork_list)
795  DEALLOCATE (work_list)
796  ALLOCATE (glob_loc_list_a(npairs_tot))
797  glob_loc_list_a = glob_loc_list(1, :)
798  CALL timestop(handle)
799  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  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  IF (ASSOCIATED(glob_loc_list)) THEN
813  DEALLOCATE (glob_loc_list)
814  END IF
815  IF (ASSOCIATED(glob_loc_list_a)) THEN
816  DEALLOCATE (glob_loc_list_a)
817  END IF
818  IF (ASSOCIATED(glob_cell_v)) THEN
819  DEALLOCATE (glob_cell_v)
820  END IF
821 
822  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  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  NULLIFY (logger)
844 
845  CALL para_env%sum(nr_ions)
846  logger => cp_get_default_logger()
847 
848  iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%PROGRAM_RUN_INFO", &
849  extension=".mmLog")
850 
851  IF (iw > 0 .AND. nr_ions > 0 .AND. print_oh) THEN
852  WRITE (iw, '(/,A,T71,I10,/)') " SIEPMANN: number of OH- ions at surface", nr_ions
853  END IF
854  IF (iw > 0 .AND. nr_ions > 0 .AND. print_h3o) THEN
855  WRITE (iw, '(/,A,T71,I10,/)') " SIEPMANN: number of H3O+ ions at surface", nr_ions
856  END IF
857  IF (iw > 0 .AND. nr_ions > 0 .AND. print_o) THEN
858  WRITE (iw, '(/,A,T71,I10,/)') " SIEPMANN: number of O^2- ions at surface", nr_ions
859  END IF
860 
861  CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%PROGRAM_RUN_INFO")
862 
863  END SUBROUTINE print_nr_ions_siepmann
864 
865 END MODULE manybody_siepmann
866 
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
fist_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: fist_neighbor_list_types.F:11
fist_nonbond_env_types
Definition: fist_nonbond_env_types.F:13
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
manybody_siepmann
implementation of dipole and three-body part of Siepmann-Sprik potential dipole term: 3rd term in Eq....
Definition: manybody_siepmann.F:17
manybody_siepmann::siepmann_energy
subroutine, public siepmann_energy(pot_loc, siepmann, r_last_update_pbc, atom_a, atom_b, nloc_size, full_loc_list, cell_v, cell, drij, particle_set, nr_oh, nr_h3o, nr_o)
energy of two-body dipole term and three-body term
Definition: manybody_siepmann.F:71
manybody_siepmann::setup_siepmann_arrays
subroutine, public setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
...
Definition: manybody_siepmann.F:707
manybody_siepmann::siepmann_forces_v2
subroutine, public siepmann_forces_v2(siepmann, r_last_update_pbc, cell_v, cell, iparticle, jparticle, f_nonbond, use_virial, rcutsq, particle_set)
forces generated by the dipole term
Definition: manybody_siepmann.F:628
manybody_siepmann::siepmann_forces_v3
subroutine, public siepmann_forces_v3(siepmann, r_last_update_pbc, cell_v, n_loc_size, full_loc_list, iparticle, jparticle, f_nonbond, use_virial, rcutsq, cell, particle_set)
forces generated by the three-body term
Definition: manybody_siepmann.F:539
manybody_siepmann::print_nr_ions_siepmann
subroutine, public print_nr_ions_siepmann(nr_ions, mm_section, para_env, print_oh, print_h3o, print_o)
prints the number of OH- ions or H3O+ ions near surface
Definition: manybody_siepmann.F:835
manybody_siepmann::destroy_siepmann_arrays
subroutine, public destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
...
Definition: manybody_siepmann.F:808
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
pair_potential_types
Definition: pair_potential_types.F:14
pair_potential_types::siepmann_type
integer, parameter, public siepmann_type
Definition: pair_potential_types.F:35
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
util
All kind of helpful little routines.
Definition: util.F:14