(git:b279b6b)
manybody_gal21.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 the GAL21 potential
10 !>
11 !> \author Clabaut Paul
12 ! **************************************************************************************************
13 MODULE manybody_gal21
14 
15  USE atomic_kind_types, ONLY: get_atomic_kind
16  USE cell_types, ONLY: cell_type,&
17  pbc
18  USE cp_log_handling, ONLY: cp_get_default_logger,&
19  cp_logger_type,&
20  cp_to_string
21  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
22  cp_print_key_unit_nr
23  USE fist_neighbor_list_types, ONLY: fist_neighbor_type,&
24  neighbor_kind_pairs_type
25  USE fist_nonbond_env_types, ONLY: pos_type
26  USE input_section_types, ONLY: section_vals_type
27  USE kinds, ONLY: dp
28  USE message_passing, ONLY: mp_para_env_type
29  USE pair_potential_types, ONLY: gal21_pot_type,&
30  gal21_type,&
31  pair_potential_pp_type,&
32  pair_potential_single_type
33  USE particle_types, ONLY: particle_type
34  USE util, ONLY: sort
35 #include "./base/base_uses.f90"
36 
37  IMPLICIT NONE
38 
39  PRIVATE
40  PUBLIC :: setup_gal21_arrays, destroy_gal21_arrays, &
41  gal21_energy, gal21_forces, &
42  print_nr_ions_gal21
43  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_gal21'
44 
45 CONTAINS
46 
47 ! **************************************************************************************************
48 !> \brief Main part of the energy evaluation of GAL2119
49 !> \param pot_loc value of total potential energy
50 !> \param gal21 all parameters of GAL2119
51 !> \param r_last_update_pbc position of every atoms on previous frame
52 !> \param iparticle first index of the atom of the evaluated pair
53 !> \param jparticle second index of the atom of the evaluated pair
54 !> \param cell dimension of the pbc cell
55 !> \param particle_set full list of atoms of the system
56 !> \param mm_section ...
57 !> \author Clabaut Paul - 2019 - ENS de Lyon
58 ! **************************************************************************************************
59  SUBROUTINE gal21_energy(pot_loc, gal21, r_last_update_pbc, iparticle, jparticle, &
60  cell, particle_set, mm_section)
61 
62  REAL(kind=dp), INTENT(OUT) :: pot_loc
63  TYPE(gal21_pot_type), POINTER :: gal21
64  TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
65  INTEGER, INTENT(IN) :: iparticle, jparticle
66  TYPE(cell_type), POINTER :: cell
67  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
68  TYPE(section_vals_type), POINTER :: mm_section
69 
70  CHARACTER(LEN=2) :: element_symbol
71  INTEGER :: index_outfile
72  REAL(kind=dp) :: anglepart, ao, bo, bxy, bz, cosalpha, &
73  drji2, eps, nvec(3), rji(3), sinalpha, &
74  sum_weight, vang, vgaussian, vh, vtt, &
75  weight
76  TYPE(cp_logger_type), POINTER :: logger
77 
78  pot_loc = 0.0_dp
79  CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
80  element_symbol=element_symbol) !Read the atom type of i
81 
82  IF (element_symbol == "O") THEN !To avoid counting two times each pair
83 
84  rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell) !Vector in pbc from j to i
85 
86  IF (.NOT. ALLOCATED(gal21%n_vectors)) THEN !First calling of the forcefield only
87  ALLOCATE (gal21%n_vectors(3, SIZE(particle_set)))
88  gal21%n_vectors(:, :) = 0.0_dp
89  END IF
90 
91  IF (gal21%express) THEN
92  logger => cp_get_default_logger()
93  index_outfile = cp_print_key_unit_nr(logger, mm_section, &
94  "PRINT%PROGRAM_RUN_INFO", extension=".mmLog")
95  IF (index_outfile > 0) WRITE (index_outfile, *) "GCN", gal21%gcn(jparticle)
96  CALL cp_print_key_finished_output(index_outfile, logger, mm_section, &
97  "PRINT%PROGRAM_RUN_INFO")
98  END IF
99 
100  !Build epsilon attraction and the parameters of the gaussian attraction as a function of gcn
101  eps = gal21%epsilon1*gal21%gcn(jparticle)*gal21%gcn(jparticle) + gal21%epsilon2*gal21%gcn(jparticle) + gal21%epsilon3
102  bxy = gal21%bxy1 + gal21%bxy2*gal21%gcn(jparticle)
103  bz = gal21%bz1 + gal21%bz2*gal21%gcn(jparticle)
104 
105  !Angular dependance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
106  vang = 0.0_dp
107 
108  !Calculation of the normal vector centered on the Me atom of the pair, only the first time that an interaction with the metal atom of the pair is evaluated
109  IF (gal21%n_vectors(1, jparticle) == 0.0_dp .AND. &
110  gal21%n_vectors(2, jparticle) == 0.0_dp .AND. &
111  gal21%n_vectors(3, jparticle) == 0.0_dp) THEN
112  gal21%n_vectors(:, jparticle) = normale(gal21, r_last_update_pbc, jparticle, &
113  particle_set, cell)
114  END IF
115 
116  nvec(:) = gal21%n_vectors(:, jparticle) !Else, retrive it, should not have moved sinc metal is supposed to be frozen
117 
118  !Calculation of the sum of the expontial weights of each Me surrounding the principal one
119  sum_weight = somme(gal21, r_last_update_pbc, iparticle, particle_set, cell)
120 
121  !Exponential damping weight for angular dependance
122  weight = exp(-sqrt(dot_product(rji, rji))/gal21%r1)
123 
124  !Calculation of the truncated fourier series of the water-dipole/surface-normal angle
125  anglepart = 0.0_dp
126  vh = 0.0_dp
127  CALL angular(anglepart, vh, gal21, r_last_update_pbc, iparticle, jparticle, cell, particle_set, nvec, &
128  .true., mm_section)
129 
130  !Build the complete angular potential while avoiding division by 0
131  IF (weight /= 0) THEN
132  vang = weight*weight*anglepart/sum_weight
133  IF (gal21%express) THEN
134  logger => cp_get_default_logger()
135  index_outfile = cp_print_key_unit_nr(logger, mm_section, &
136  "PRINT%PROGRAM_RUN_INFO", extension=".mmLog")
137  IF (index_outfile > 0) WRITE (index_outfile, *) "Fermi", weight*weight/sum_weight
138  CALL cp_print_key_finished_output(index_outfile, logger, mm_section, &
139  "PRINT%PROGRAM_RUN_INFO")
140  END IF
141  END IF
142  !END Angular%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
143 
144  !Attractive Gaussian %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145  vgaussian = 0.0_dp
146  drji2 = dot_product(rji, rji)
147  !Alpha is the angle of the Me-O vector with the normale vector. Used for gaussian attaction
148 
149  cosalpha = dot_product(rji, nvec)/sqrt(drji2)
150  IF (cosalpha < -1.0_dp) cosalpha = -1.0_dp
151  IF (cosalpha > +1.0_dp) cosalpha = +1.0_dp
152  sinalpha = sin(acos(cosalpha))
153 
154  !Gaussian component of the energy
155  vgaussian = -1.0_dp*eps*exp(-bz*drji2*cosalpha*cosalpha &
156  - bxy*drji2*sinalpha*sinalpha)
157  !END Gaussian%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
158 
159  ao = gal21%AO1 + gal21%AO2*gal21%gcn(jparticle)
160  bo = gal21%BO1 + gal21%BO2*gal21%gcn(jparticle)
161 
162  !Tang and toennies potential for physisorption
163  vtt = ao*exp(-bo*sqrt(drji2)) - (1.0 - exp(-bo*sqrt(drji2)) &
164  - bo*sqrt(drji2)*exp(-bo*sqrt(drji2)) &
165  - (((bo*sqrt(drji2))**2)/2)*exp(-bo*sqrt(drji2)) &
166  - (((bo*sqrt(drji2))**3)/6)*exp(-bo*sqrt(drji2)) &
167  - (((bo*sqrt(drji2))**4)/24)*exp(-bo*sqrt(drji2)) &
168  - (((bo*sqrt(drji2))**5)/120)*exp(-bo*sqrt(drji2)) &
169  - (((bo*sqrt(drji2))**6)/720)*exp(-bo*sqrt(drji2))) &
170  *gal21%c/(sqrt(drji2)**6)
171 
172  !For fit purpose only
173  IF (gal21%express) THEN
174  logger => cp_get_default_logger()
175  index_outfile = cp_print_key_unit_nr(logger, mm_section, &
176  "PRINT%PROGRAM_RUN_INFO", extension=".mmLog")
177  IF (index_outfile > 0) WRITE (index_outfile, *) "Gau", -1.0_dp*exp(-bz*drji2*cosalpha*cosalpha &
178  - bxy*drji2*sinalpha*sinalpha)
179  IF (weight == 0 .AND. index_outfile > 0) WRITE (index_outfile, *) "Fermi 0"
180  IF (index_outfile > 0) WRITE (index_outfile, *) "expO", exp(-bo*sqrt(drji2))
181  IF (index_outfile > 0) WRITE (index_outfile, *) "cstpart", -(1.0 - exp(-bo*sqrt(drji2)) &
182  - bo*sqrt(drji2)*exp(-bo*sqrt(drji2)) &
183  - (((bo*sqrt(drji2))**2)/2)*exp(-bo*sqrt(drji2)) &
184  - (((bo*sqrt(drji2))**3)/6)*exp(-bo*sqrt(drji2)) &
185  - (((bo*sqrt(drji2))**4)/24)*exp(-bo*sqrt(drji2)) &
186  - (((bo*sqrt(drji2))**5)/120)*exp(-bo*sqrt(drji2)) &
187  - (((bo*sqrt(drji2))**6)/720)*exp(-bo*sqrt(drji2))) &
188  *gal21%c/(sqrt(drji2)**6)
189  IF (index_outfile > 0) WRITE (index_outfile, *) "params_lin_eps", gal21%epsilon1, gal21%epsilon2, gal21%epsilon3
190  IF (index_outfile > 0) WRITE (index_outfile, *) "params_lin_A0", ao
191  CALL cp_print_key_finished_output(index_outfile, logger, mm_section, &
192  "PRINT%PROGRAM_RUN_INFO")
193  END IF
194  !Compute the total energy
195  pot_loc = vgaussian + vang + vtt + vh
196 
197  END IF
198 
199  END SUBROUTINE gal21_energy
200 
201 ! **************************************************************************************************
202 !> \brief The idea is to build a vector normal to the local surface by using the symetry of the
203 !> surface that make the opposite vectors compensate themself. The vector is therefore in the
204 !>. direction of the missing atoms of a large coordination sphere
205 !> \param gal21 ...
206 !> \param r_last_update_pbc ...
207 !> \param jparticle ...
208 !> \param particle_set ...
209 !> \param cell ...
210 !> \return ...
211 !> \retval normale ...
212 ! **************************************************************************************************
213  FUNCTION normale(gal21, r_last_update_pbc, jparticle, particle_set, cell)
214  TYPE(gal21_pot_type), POINTER :: gal21
215  TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
216  INTEGER, INTENT(IN) :: jparticle
217  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
218  TYPE(cell_type), POINTER :: cell
219  REAL(kind=dp) :: normale(3)
220 
221  CHARACTER(LEN=2) :: element_symbol_k
222  INTEGER :: kparticle, natom
223  REAL(kind=dp) :: drjk, rjk(3)
224 
225  natom = SIZE(particle_set)
226  normale(:) = 0.0_dp
227 
228  DO kparticle = 1, natom !Loop on every atom of the system
229  IF (kparticle == jparticle) cycle !Avoid the principal Me atom (j) in the counting
230  CALL get_atomic_kind(atomic_kind=particle_set(kparticle)%atomic_kind, &
231  element_symbol=element_symbol_k)
232  IF (element_symbol_k /= gal21%met1 .AND. element_symbol_k /= gal21%met2) cycle !Keep only metals
233  rjk(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(kparticle)%r(:), cell)
234  drjk = sqrt(dot_product(rjk, rjk))
235  IF (drjk > gal21%rcutsq) cycle !Keep only those within square root of the force-field cutoff distance of the metallic atom of the evaluated pair
236  normale(:) = normale(:) - rjk(:)/(drjk*drjk*drjk*drjk*drjk) !Build the normal, vector by vector
237  END DO
238 
239  ! Normalisation of the vector
240  normale(:) = normale(:)/sqrt(dot_product(normale, normale))
241 
242  END FUNCTION normale
243 
244 ! **************************************************************************************************
245 !> \brief Scan all the Me atoms that have been counted in the O-Me paires and sum their exp. weights
246 !> \param gal21 ...
247 !> \param r_last_update_pbc ...
248 !> \param iparticle ...
249 !> \param particle_set ...
250 !> \param cell ...
251 !> \return ...
252 !> \retval somme ...
253 ! **************************************************************************************************
254  FUNCTION somme(gal21, r_last_update_pbc, iparticle, particle_set, cell)
255  TYPE(gal21_pot_type), POINTER :: gal21
256  TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
257  INTEGER, INTENT(IN) :: iparticle
258  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
259  TYPE(cell_type), POINTER :: cell
260  REAL(kind=dp) :: somme
261 
262  CHARACTER(LEN=2) :: element_symbol_k
263  INTEGER :: kparticle, natom
264  REAL(kind=dp) :: rki(3)
265 
266  natom = SIZE(particle_set)
267  somme = 0.0_dp
268 
269  DO kparticle = 1, natom !Loop on every atom of the system
270  CALL get_atomic_kind(atomic_kind=particle_set(kparticle)%atomic_kind, &
271  element_symbol=element_symbol_k)
272  IF (element_symbol_k /= gal21%met1 .AND. element_symbol_k /= gal21%met2) cycle !Keep only metals
273  rki(:) = pbc(r_last_update_pbc(kparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell)
274  IF (sqrt(dot_product(rki, rki)) > gal21%rcutsq) cycle !Keep only those within cutoff distance of the oxygen atom of the evaluated pair (the omega ensemble)
275  IF (element_symbol_k == gal21%met1) somme = somme + exp(-sqrt(dot_product(rki, rki))/gal21%r1) !Build the sum of the exponential weights
276  IF (element_symbol_k == gal21%met2) somme = somme + exp(-sqrt(dot_product(rki, rki))/gal21%r2) !Build the sum of the exponential weights
277  END DO
278 
279  END FUNCTION somme
280 
281 ! **************************************************************************************************
282 !> \brief Compute the angular dependance (on theta) of the forcefield
283 !> \param anglepart ...
284 !> \param VH ...
285 !> \param gal21 ...
286 !> \param r_last_update_pbc ...
287 !> \param iparticle ...
288 !> \param jparticle ...
289 !> \param cell ...
290 !> \param particle_set ...
291 !> \param nvec ...
292 !> \param energy ...
293 !> \param mm_section ...
294 !> \return ...
295 !> \retval angular ...
296 ! **************************************************************************************************
297  SUBROUTINE angular(anglepart, VH, gal21, r_last_update_pbc, iparticle, jparticle, cell, &
298  particle_set, nvec, energy, mm_section)
299  REAL(kind=dp) :: anglepart, vh
300  TYPE(gal21_pot_type), POINTER :: gal21
301  TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
302  INTEGER, INTENT(IN) :: iparticle, jparticle
303  TYPE(cell_type), POINTER :: cell
304  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
305  REAL(kind=dp), DIMENSION(3) :: nvec
306  LOGICAL :: energy
307  TYPE(section_vals_type), POINTER :: mm_section
308 
309  CHARACTER(LEN=2) :: element_symbol
310  INTEGER :: count_h, iatom, index_h1, index_h2, &
311  index_outfile, natom
312  REAL(kind=dp) :: a1, a2, a3, a4, bh, costheta, &
313  h_max_dist, rih(3), rih1(3), rih2(3), &
314  rix(3), rjh1(3), rjh2(3), theta
315  TYPE(cp_logger_type), POINTER :: logger
316 
317  count_h = 0
318  index_h1 = 0
319  index_h2 = 0
320  h_max_dist = 2.1_dp ! 1.1 angstrom
321  natom = SIZE(particle_set)
322 
323  DO iatom = 1, natom !Loop on every atom of the system
324  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
325  element_symbol=element_symbol)
326  IF (element_symbol /= "H") cycle !Kepp only hydrogen
327  rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell)
328  IF (sqrt(dot_product(rih, rih)) >= h_max_dist) cycle !Keep only hydrogen that are bounded to the considered O
329  count_h = count_h + 1
330  IF (count_h == 1) THEN
331  index_h1 = iatom
332  ELSEIF (count_h == 2) THEN
333  index_h2 = iatom
334  END IF
335  END DO
336 
337  ! Abort if the oxygen is not part of a water molecule (2 H)
338  IF (count_h /= 2) THEN
339  CALL cp_abort(__location__, &
340  " Error: Found "//cp_to_string(count_h)//" H atoms for O atom "//cp_to_string(iparticle))
341  END IF
342 
343  a1 = gal21%a11 + gal21%a12*gal21%gcn(jparticle) + gal21%a13*gal21%gcn(jparticle)*gal21%gcn(jparticle)
344  a2 = gal21%a21 + gal21%a22*gal21%gcn(jparticle) + gal21%a23*gal21%gcn(jparticle)*gal21%gcn(jparticle)
345  a3 = gal21%a31 + gal21%a32*gal21%gcn(jparticle) + gal21%a33*gal21%gcn(jparticle)*gal21%gcn(jparticle)
346  a4 = gal21%a41 + gal21%a42*gal21%gcn(jparticle) + gal21%a43*gal21%gcn(jparticle)*gal21%gcn(jparticle)
347 
348  rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
349  rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
350  rix(:) = rih1(:) + rih2(:) ! build the dipole vector rix of the H2O molecule
351  costheta = dot_product(rix, nvec)/sqrt(dot_product(rix, rix))
352  IF (costheta < -1.0_dp) costheta = -1.0_dp
353  IF (costheta > +1.0_dp) costheta = +1.0_dp
354  theta = acos(costheta) ! Theta is the angle between the normal to the surface and the dipole
355  anglepart = a1*costheta + a2*cos(2.0_dp*theta) + a3*cos(3.0_dp*theta) &
356  + a4*cos(4.0_dp*theta) ! build the fourier series
357 
358  bh = gal21%BH1 + gal21%gcn(jparticle)*gal21%BH2
359 
360  rjh1(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
361  rjh2(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
362  vh = (gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)*(exp(-bh*sqrt(dot_product(rjh1, rjh1))) + &
363  exp(-bh*sqrt(dot_product(rjh2, rjh2))))
364 
365  ! For fit purpose
366  IF (gal21%express .AND. energy) THEN
367  logger => cp_get_default_logger()
368  index_outfile = cp_print_key_unit_nr(logger, mm_section, &
369  "PRINT%PROGRAM_RUN_INFO", extension=".mmLog")
370 
371  IF (index_outfile > 0) WRITE (index_outfile, *) "Fourier", costheta, cos(2.0_dp*theta), cos(3.0_dp*theta), &
372  cos(4.0_dp*theta) !, theta
373  IF (index_outfile > 0) WRITE (index_outfile, *) "H_rep", exp(-bh*sqrt(dot_product(rjh1, rjh1))) + &
374  exp(-bh*sqrt(dot_product(rjh2, rjh2)))
375  !IF (index_outfile > 0) WRITE (index_outfile, *) "H_r6", -1/DOT_PRODUCT(rjh1,rjh1)**3 -1/DOT_PRODUCT(rjh2,rjh2)**3
376 
377  CALL cp_print_key_finished_output(index_outfile, logger, mm_section, &
378  "PRINT%PROGRAM_RUN_INFO")
379  END IF
380 
381  END SUBROUTINE
382 
383 ! **************************************************************************************************
384 !> \brief forces generated by the GAL2119 potential
385 !> \param gal21 all parameters of GAL2119
386 !> \param r_last_update_pbc position of every atoms on previous frame
387 !> \param iparticle first index of the atom of the evaluated pair
388 !> \param jparticle second index of the atom of the evaluated pair
389 !> \param f_nonbond all the forces applying on the system
390 !> \param pv_nonbond ...
391 !> \param use_virial request of usage of virial (for barostat)
392 !> \param cell dimension of the pbc cell
393 !> \param particle_set full list of atoms of the system
394 !> \author Clabaut Paul - 2019 - ENS de Lyon
395 ! **************************************************************************************************
396  SUBROUTINE gal21_forces(gal21, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, &
397  use_virial, cell, particle_set)
398  TYPE(gal21_pot_type), POINTER :: gal21
399  TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
400  INTEGER, INTENT(IN) :: iparticle, jparticle
401  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond
402  LOGICAL, INTENT(IN) :: use_virial
403  TYPE(cell_type), POINTER :: cell
404  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
405 
406  CHARACTER(LEN=2) :: element_symbol
407  REAL(kind=dp) :: anglepart, ao, bo, bxy, bz, cosalpha, dgauss(3), drji, drjicosalpha(3), &
408  drjisinalpha(3), dtt(3), dweight(3), eps, nvec(3), prefactor, rji(3), rji_hat(3), &
409  sinalpha, sum_weight, vgaussian, vh, weight
410  TYPE(section_vals_type), POINTER :: mm_section
411 
412  CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, &
413  element_symbol=element_symbol)
414 
415  IF (element_symbol == "O") THEN !To avoid counting two times each pair
416 
417  rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell)
418  drji = sqrt(dot_product(rji, rji))
419  rji_hat(:) = rji(:)/drji ! hat = pure directional component of a given vector
420 
421  IF (.NOT. ALLOCATED(gal21%n_vectors)) THEN !First calling of the forcefield only
422  ALLOCATE (gal21%n_vectors(3, SIZE(particle_set)))
423  gal21%n_vectors(:, :) = 0.0_dp
424  END IF
425 
426  !Build epsilon attraction and the a parameters of the Fourier serie as quadratic fucntion of gcn
427  eps = gal21%epsilon1*gal21%gcn(jparticle)*gal21%gcn(jparticle) + gal21%epsilon2*gal21%gcn(jparticle) + gal21%epsilon3
428  bxy = gal21%bxy1 + gal21%bxy2*gal21%gcn(jparticle)
429  bz = gal21%bz1 + gal21%bz2*gal21%gcn(jparticle)
430 
431  !Angular dependance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
432 
433  !Calculation of the normal vector centered on the Me atom of the pair, only the first time that an interaction with the metal atom of the pair is evaluated
434  IF (gal21%n_vectors(1, jparticle) == 0.0_dp .AND. &
435  gal21%n_vectors(2, jparticle) == 0.0_dp .AND. &
436  gal21%n_vectors(3, jparticle) == 0.0_dp) THEN
437  gal21%n_vectors(:, jparticle) = normale(gal21, r_last_update_pbc, jparticle, &
438  particle_set, cell)
439  END IF
440 
441  nvec(:) = gal21%n_vectors(:, jparticle) !Else, retrive it, should not have moved sinc metal is supposed to be frozen
442 
443  !Calculation of the sum of the expontial weights of each Me surrounding the principal one
444  sum_weight = somme(gal21, r_last_update_pbc, iparticle, particle_set, cell)
445 
446  !Exponential damping weight for angular dependance
447  weight = exp(-drji/gal21%r1)
448  dweight(:) = 1.0_dp/gal21%r1*weight*rji_hat(:) !Derivativ of it
449 
450  !Calculation of the truncated fourier series of the water-dipole/surface-normal angle
451  NULLIFY (mm_section)
452  anglepart = 0.0_dp
453  vh = 0.0_dp
454  CALL angular(anglepart, vh, gal21, r_last_update_pbc, iparticle, jparticle, cell, particle_set, nvec, &
455  .false., mm_section)
456 
457  !Build the average of the exponential weight while avoiding division by 0
458  IF (weight /= 0) THEN
459  ! Calculate the first component of the derivativ of the angular term
460  f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) + 2.0_dp*dweight(1:3)*weight* &
461  anglepart/sum_weight
462 
463  IF (use_virial) THEN
464  pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rji(1)*2.0_dp*dweight(1:3)*weight* &
465  anglepart/sum_weight
466  pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rji(2)*2.0_dp*dweight(1:3)*weight* &
467  anglepart/sum_weight
468  pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rji(3)*2.0_dp*dweight(1:3)*weight* &
469  anglepart/sum_weight
470  END IF
471 
472  ! Calculate the second component of the derivativ of the angular term
473  CALL somme_d(gal21, r_last_update_pbc, iparticle, jparticle, &
474  f_nonbond, pv_nonbond, use_virial, particle_set, cell, anglepart, sum_weight)
475 
476  prefactor = (-1.0_dp)*weight*weight/sum_weight ! Avoiding division by 0
477 
478  ! Calculate the third component of the derivativ of the angular term
479  CALL angular_d(gal21, r_last_update_pbc, iparticle, jparticle, &
480  f_nonbond, pv_nonbond, use_virial, prefactor, cell, particle_set, nvec)
481  END IF
482  !END Angular%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
483 
484  !Attractive Gaussian %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
485  !Alpha is the angle of the Me-O vector with the normale vector. Used for gaussian attaction
486  cosalpha = dot_product(rji, nvec)/drji
487  IF (cosalpha < -1.0_dp) cosalpha = -1.0_dp
488  IF (cosalpha > +1.0_dp) cosalpha = +1.0_dp
489  sinalpha = sin(acos(cosalpha))
490 
491  !Gaussian component of the energy
492  vgaussian = -1.0_dp*eps*exp(-bz*dot_product(rji, rji)*cosalpha*cosalpha &
493  - bxy*dot_product(rji, rji)*sinalpha*sinalpha)
494 
495  ! Calculation of partial derivativ of the gaussian components
496  drjicosalpha(:) = rji_hat(:)*cosalpha + nvec(:) - cosalpha*rji_hat(:)
497  drjisinalpha(:) = rji_hat(:)*sinalpha - (cosalpha/sinalpha)*(nvec(:) - cosalpha*rji_hat(:))
498  dgauss(:) = (-1.0_dp*bz*2*drji*cosalpha*drjicosalpha - &
499  1.0_dp*bxy*2*drji*sinalpha*drjisinalpha)*vgaussian*(-1.0_dp)
500 
501  ! Force due to gaussian term
502  f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) + dgauss(1:3)
503 
504  IF (use_virial) THEN
505  pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rji(1)*dgauss(1:3)
506  pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rji(2)*dgauss(1:3)
507  pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rji(3)*dgauss(1:3)
508  END IF
509  !END Gaussian%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
510 
511  ao = gal21%AO1 + gal21%AO2*gal21%gcn(jparticle)
512  bo = gal21%BO1 + gal21%BO2*gal21%gcn(jparticle)
513 
514  !Derivativ of the Tang and Toennies term
515  dtt(:) = (-(ao*bo + (bo**7)*gal21%c/720)*exp(-bo*drji) + 6*(gal21%c/drji**7)* &
516  (1.0 - exp(-bo*drji) &
517  - bo*drji*exp(-bo*drji) &
518  - (((bo*drji)**2)/2)*exp(-bo*drji) &
519  - (((bo*drji)**3)/6)*exp(-bo*drji) &
520  - (((bo*drji)**4)/24)*exp(-bo*drji) &
521  - (((bo*drji)**5)/120)*exp(-bo*drji) &
522  - (((bo*drji)**6)/720)*exp(-bo*drji)) &
523  )*rji_hat(:)
524 
525  ! Force of Tang & Toennies
526  f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) - dtt(1:3)
527 
528  IF (use_virial) THEN
529  pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) - rji(1)*dtt(1:3)
530  pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) - rji(2)*dtt(1:3)
531  pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) - rji(3)*dtt(1:3)
532  END IF
533 
534  END IF
535 
536  END SUBROUTINE gal21_forces
537 
538 ! **************************************************************************************************
539 !> \brief Derivativ of the second component of angular dependance
540 !> \param gal21 ...
541 !> \param r_last_update_pbc ...
542 !> \param iparticle ...
543 !> \param jparticle ...
544 !> \param f_nonbond ...
545 !> \param pv_nonbond ...
546 !> \param use_virial ...
547 !> \param particle_set ...
548 !> \param cell ...
549 !> \param anglepart ...
550 !> \param sum_weight ...
551 ! **************************************************************************************************
552  SUBROUTINE somme_d(gal21, r_last_update_pbc, iparticle, jparticle, &
553  f_nonbond, pv_nonbond, use_virial, particle_set, cell, anglepart, sum_weight)
554  TYPE(gal21_pot_type), POINTER :: gal21
555  TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
556  INTEGER, INTENT(IN) :: iparticle, jparticle
557  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond
558  LOGICAL, INTENT(IN) :: use_virial
559  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
560  TYPE(cell_type), POINTER :: cell
561  REAL(kind=dp), INTENT(IN) :: anglepart, sum_weight
562 
563  CHARACTER(LEN=2) :: element_symbol_k
564  INTEGER :: kparticle, natom
565  REAL(kind=dp) :: drki, dwdr(3), rji(3), rki(3), &
566  rki_hat(3), weight_rji
567 
568  rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell)
569  weight_rji = exp(-sqrt(dot_product(rji, rji))/gal21%r1)
570 
571  natom = SIZE(particle_set)
572  DO kparticle = 1, natom !Loop on every atom of the system
573  CALL get_atomic_kind(atomic_kind=particle_set(kparticle)%atomic_kind, &
574  element_symbol=element_symbol_k)
575  IF (element_symbol_k /= gal21%met1 .AND. element_symbol_k /= gal21%met2) cycle !Keep only metals
576  rki(:) = pbc(r_last_update_pbc(kparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell)
577  IF (sqrt(dot_product(rki, rki)) > gal21%rcutsq) cycle !Keep only those within cutoff distance of the oxygen atom of the evaluated pair (the omega ensemble)
578  drki = sqrt(dot_product(rki, rki))
579  rki_hat(:) = rki(:)/drki
580 
581  IF (element_symbol_k == gal21%met1) dwdr(:) = (-1.0_dp)*(1.0_dp/gal21%r1)*exp(-drki/gal21%r1)*rki_hat(:) !Build the sum of derivativs
582  IF (element_symbol_k == gal21%met2) dwdr(:) = (-1.0_dp)*(1.0_dp/gal21%r2)*exp(-drki/gal21%r2)*rki_hat(:)
583 
584  f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) + dwdr(1:3)*weight_rji &
585  *weight_rji*anglepart/(sum_weight**2)
586 
587  IF (use_virial) THEN
588  pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rki(1)*dwdr(1:3)*weight_rji &
589  *weight_rji*anglepart/(sum_weight**2)
590  pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rki(2)*dwdr(1:3)*weight_rji &
591  *weight_rji*anglepart/(sum_weight**2)
592  pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rki(3)*dwdr(1:3)*weight_rji &
593  *weight_rji*anglepart/(sum_weight**2)
594  END IF
595 
596  END DO
597 
598  END SUBROUTINE somme_d
599 
600 ! **************************************************************************************************
601 !> \brief Derivativ of the third component of angular term
602 !> \param gal21 ...
603 !> \param r_last_update_pbc ...
604 !> \param iparticle ...
605 !> \param jparticle ...
606 !> \param f_nonbond ...
607 !> \param pv_nonbond ...
608 !> \param use_virial ...
609 !> \param prefactor ...
610 !> \param cell ...
611 !> \param particle_set ...
612 !> \param nvec ...
613 ! **************************************************************************************************
614  SUBROUTINE angular_d(gal21, r_last_update_pbc, iparticle, jparticle, f_nonbond, &
615  pv_nonbond, use_virial, prefactor, cell, particle_set, nvec)
616  TYPE(gal21_pot_type), POINTER :: gal21
617  TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
618  INTEGER, INTENT(IN) :: iparticle, jparticle
619  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond
620  LOGICAL, INTENT(IN) :: use_virial
621  REAL(kind=dp), INTENT(IN) :: prefactor
622  TYPE(cell_type), POINTER :: cell
623  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
624  REAL(kind=dp), DIMENSION(3) :: nvec
625 
626  CHARACTER(LEN=2) :: element_symbol
627  INTEGER :: count_h, iatom, index_h1, index_h2, natom
628  REAL(kind=dp) :: a1, a2, a3, a4, bh, costheta, &
629  dsumdtheta, h_max_dist, theta
630  REAL(kind=dp), DIMENSION(3) :: dangular, dcostheta, rih, rih1, rih2, &
631  rix, rix_hat, rjh1, rjh2, rji, rji_hat
632 
633  count_h = 0
634  index_h1 = 0
635  index_h2 = 0
636  h_max_dist = 2.1_dp ! 1.1 angstrom
637  natom = SIZE(particle_set)
638 
639  DO iatom = 1, natom !Loop on every atom of the system
640  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
641  element_symbol=element_symbol)
642  IF (element_symbol /= "H") cycle !Kepp only hydrogen
643  rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell)
644  IF (sqrt(dot_product(rih, rih)) >= h_max_dist) cycle !Keep only hydrogen that are bounded to the considered O
645  count_h = count_h + 1
646  IF (count_h == 1) THEN
647  index_h1 = iatom
648  ELSEIF (count_h == 2) THEN
649  index_h2 = iatom
650  END IF
651  END DO
652 
653  ! Abort if the oxygen is not part of a water molecule (2 H)
654  IF (count_h /= 2) THEN
655  CALL cp_abort(__location__, &
656  " Error: Found "//cp_to_string(count_h)//" H atoms for O atom "//cp_to_string(iparticle))
657  END IF
658 
659  a1 = gal21%a11 + gal21%a12*gal21%gcn(jparticle) + gal21%a13*gal21%gcn(jparticle)*gal21%gcn(jparticle)
660  a2 = gal21%a21 + gal21%a22*gal21%gcn(jparticle) + gal21%a23*gal21%gcn(jparticle)*gal21%gcn(jparticle)
661  a3 = gal21%a31 + gal21%a32*gal21%gcn(jparticle) + gal21%a33*gal21%gcn(jparticle)*gal21%gcn(jparticle)
662  a4 = gal21%a41 + gal21%a42*gal21%gcn(jparticle) + gal21%a43*gal21%gcn(jparticle)*gal21%gcn(jparticle)
663 
664  rji(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(iparticle)%r(:), cell)
665  rji_hat(:) = rji(:)/sqrt(dot_product(rji, rji)) ! hat = pure directional component of a given vector
666 
667  !dipole vector rix of the H2O molecule
668  rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
669  rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
670  rix(:) = rih1(:) + rih2(:) ! build the dipole vector rix of the H2O molecule
671  rix_hat(:) = rix(:)/sqrt(dot_product(rix, rix)) ! hat = pure directional component of a given vector
672  costheta = dot_product(rix, nvec)/sqrt(dot_product(rix, rix)) ! Theta is the angle between the normal to the surface and the dipole
673  IF (costheta < -1.0_dp) costheta = -1.0_dp
674  IF (costheta > +1.0_dp) costheta = +1.0_dp
675  theta = acos(costheta) ! Theta is the angle between the normal to the surface and the dipole
676 
677  ! Calculation of partial derivativ of the angular components
678  dsumdtheta = -1.0_dp*a1*sin(theta) - a2*2.0_dp*sin(2.0_dp*theta) - &
679  a3*3.0_dp*sin(3.0_dp*theta) - a4*4.0_dp*sin(4.0_dp*theta)
680  dcostheta(:) = (1.0_dp/sqrt(dot_product(rix, rix)))*(nvec(:) - costheta*rix_hat(:))
681  dangular(:) = prefactor*dsumdtheta*(-1.0_dp/sin(theta))*dcostheta(:)
682 
683  !Force due to the third component of the derivativ of the angular term
684  f_nonbond(1:3, iparticle) = f_nonbond(1:3, iparticle) - dangular(1:3)*2.0_dp !(one per H)
685  f_nonbond(1:3, index_h1) = f_nonbond(1:3, index_h1) + dangular(1:3)
686  f_nonbond(1:3, index_h2) = f_nonbond(1:3, index_h2) + dangular(1:3)
687 
688  IF (use_virial) THEN
689  pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rix(1)*dangular(1:3)
690  pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rix(2)*dangular(1:3)
691  pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rix(3)*dangular(1:3)
692  END IF
693 
694  bh = gal21%BH1 + gal21%gcn(jparticle)*gal21%BH2
695 
696  rjh1(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell)
697  f_nonbond(1:3, index_h1) = f_nonbond(1:3, index_h1) + (gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
698  bh*exp(-bh*sqrt(dot_product(rjh1, rjh1)))*rjh1(:)/sqrt(dot_product(rjh1, rjh1))
699 
700  IF (use_virial) THEN
701  pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rjh1(1)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
702  bh*exp(-bh*sqrt(dot_product(rjh1, rjh1)))) &
703  *rjh1(:)/sqrt(dot_product(rjh1, rjh1))
704  pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rjh1(2)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
705  bh*exp(-bh*sqrt(dot_product(rjh1, rjh1)))) &
706  *rjh1(:)/sqrt(dot_product(rjh1, rjh1))
707  pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rjh1(3)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
708  bh*exp(-bh*sqrt(dot_product(rjh1, rjh1)))) &
709  *rjh1(:)/sqrt(dot_product(rjh1, rjh1))
710  END IF
711 
712  rjh2(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell)
713  f_nonbond(1:3, index_h2) = f_nonbond(1:3, index_h2) + ((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
714  bh*exp(-bh*sqrt(dot_product(rjh2, rjh2)))) &
715  *rjh2(:)/sqrt(dot_product(rjh2, rjh2))
716 
717  IF (use_virial) THEN
718  pv_nonbond(1, 1:3) = pv_nonbond(1, 1:3) + rjh2(1)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
719  bh*exp(-bh*sqrt(dot_product(rjh2, rjh2)))) &
720  *rjh2(:)/sqrt(dot_product(rjh2, rjh2))
721  pv_nonbond(2, 1:3) = pv_nonbond(2, 1:3) + rjh2(2)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
722  bh*exp(-bh*sqrt(dot_product(rjh2, rjh2)))) &
723  *rjh2(:)/sqrt(dot_product(rjh2, rjh2))
724  pv_nonbond(3, 1:3) = pv_nonbond(3, 1:3) + rjh2(3)*((gal21%AH2*gal21%gcn(jparticle) + gal21%AH1)* &
725  bh*exp(-bh*sqrt(dot_product(rjh2, rjh2)))) &
726  *rjh2(:)/sqrt(dot_product(rjh2, rjh2))
727  END IF
728 
729  END SUBROUTINE angular_d
730 
731 ! **************************************************************************************************
732 !> \brief ...
733 !> \param nonbonded ...
734 !> \param potparm ...
735 !> \param glob_loc_list ...
736 !> \param glob_cell_v ...
737 !> \param glob_loc_list_a ...
738 !> \param cell ...
739 !> \par History
740 ! **************************************************************************************************
741  SUBROUTINE setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, &
742  glob_loc_list_a, cell)
743  TYPE(fist_neighbor_type), POINTER :: nonbonded
744  TYPE(pair_potential_pp_type), POINTER :: potparm
745  INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list
746  REAL(kind=dp), DIMENSION(:, :), POINTER :: glob_cell_v
747  INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a
748  TYPE(cell_type), POINTER :: cell
749 
750  CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_gal21_arrays'
751 
752  INTEGER :: handle, i, iend, igrp, ikind, ilist, &
753  ipair, istart, jkind, nkinds, npairs, &
754  npairs_tot
755  INTEGER, DIMENSION(:), POINTER :: work_list, work_list2
756  INTEGER, DIMENSION(:, :), POINTER :: list
757  REAL(kind=dp), DIMENSION(3) :: cell_v, cvi
758  REAL(kind=dp), DIMENSION(:, :), POINTER :: rwork_list
759  TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
760  TYPE(pair_potential_single_type), POINTER :: pot
761 
762  cpassert(.NOT. ASSOCIATED(glob_loc_list))
763  cpassert(.NOT. ASSOCIATED(glob_loc_list_a))
764  cpassert(.NOT. ASSOCIATED(glob_cell_v))
765  CALL timeset(routinen, handle)
766  npairs_tot = 0
767  nkinds = SIZE(potparm%pot, 1)
768  DO ilist = 1, nonbonded%nlists
769  neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
770  npairs = neighbor_kind_pair%npairs
771  IF (npairs == 0) cycle
772  kind_group_loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
773  istart = neighbor_kind_pair%grp_kind_start(igrp)
774  iend = neighbor_kind_pair%grp_kind_end(igrp)
775  ikind = neighbor_kind_pair%ij_kind(1, igrp)
776  jkind = neighbor_kind_pair%ij_kind(2, igrp)
777  pot => potparm%pot(ikind, jkind)%pot
778  npairs = iend - istart + 1
779  IF (pot%no_mb) cycle
780  DO i = 1, SIZE(pot%type)
781  IF (pot%type(i) == gal21_type) npairs_tot = npairs_tot + npairs
782  END DO
783  END DO kind_group_loop1
784  END DO
785  ALLOCATE (work_list(npairs_tot))
786  ALLOCATE (work_list2(npairs_tot))
787  ALLOCATE (glob_loc_list(2, npairs_tot))
788  ALLOCATE (glob_cell_v(3, npairs_tot))
789  ! Fill arrays with data
790  npairs_tot = 0
791  DO ilist = 1, nonbonded%nlists
792  neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
793  npairs = neighbor_kind_pair%npairs
794  IF (npairs == 0) cycle
795  kind_group_loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
796  istart = neighbor_kind_pair%grp_kind_start(igrp)
797  iend = neighbor_kind_pair%grp_kind_end(igrp)
798  ikind = neighbor_kind_pair%ij_kind(1, igrp)
799  jkind = neighbor_kind_pair%ij_kind(2, igrp)
800  list => neighbor_kind_pair%list
801  cvi = neighbor_kind_pair%cell_vector
802  pot => potparm%pot(ikind, jkind)%pot
803  npairs = iend - istart + 1
804  IF (pot%no_mb) cycle
805  cell_v = matmul(cell%hmat, cvi)
806  DO i = 1, SIZE(pot%type)
807  ! gal21
808  IF (pot%type(i) == gal21_type) THEN
809  DO ipair = 1, npairs
810  glob_loc_list(:, npairs_tot + ipair) = list(:, istart - 1 + ipair)
811  glob_cell_v(1:3, npairs_tot + ipair) = cell_v(1:3)
812  END DO
813  npairs_tot = npairs_tot + npairs
814  END IF
815  END DO
816  END DO kind_group_loop2
817  END DO
818  ! Order the arrays w.r.t. the first index of glob_loc_list
819  CALL sort(glob_loc_list(1, :), npairs_tot, work_list)
820  DO ipair = 1, npairs_tot
821  work_list2(ipair) = glob_loc_list(2, work_list(ipair))
822  END DO
823  glob_loc_list(2, :) = work_list2
824  DEALLOCATE (work_list2)
825  ALLOCATE (rwork_list(3, npairs_tot))
826  DO ipair = 1, npairs_tot
827  rwork_list(:, ipair) = glob_cell_v(:, work_list(ipair))
828  END DO
829  glob_cell_v = rwork_list
830  DEALLOCATE (rwork_list)
831  DEALLOCATE (work_list)
832  ALLOCATE (glob_loc_list_a(npairs_tot))
833  glob_loc_list_a = glob_loc_list(1, :)
834  CALL timestop(handle)
835  END SUBROUTINE setup_gal21_arrays
836 
837 ! **************************************************************************************************
838 !> \brief ...
839 !> \param glob_loc_list ...
840 !> \param glob_cell_v ...
841 !> \param glob_loc_list_a ...
842 ! **************************************************************************************************
843  SUBROUTINE destroy_gal21_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
844  INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list
845  REAL(kind=dp), DIMENSION(:, :), POINTER :: glob_cell_v
846  INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a
847 
848  IF (ASSOCIATED(glob_loc_list)) THEN
849  DEALLOCATE (glob_loc_list)
850  END IF
851  IF (ASSOCIATED(glob_loc_list_a)) THEN
852  DEALLOCATE (glob_loc_list_a)
853  END IF
854  IF (ASSOCIATED(glob_cell_v)) THEN
855  DEALLOCATE (glob_cell_v)
856  END IF
857 
858  END SUBROUTINE destroy_gal21_arrays
859 
860 ! **************************************************************************************************
861 !> \brief prints the number of OH- ions or H3O+ ions near surface
862 !> \param nr_ions number of ions
863 !> \param mm_section ...
864 !> \param para_env ...
865 !> \param print_oh flag indicating if number OH- is printed
866 !> \param print_h3o flag indicating if number H3O+ is printed
867 !> \param print_o flag indicating if number O^(2-) is printed
868 ! **************************************************************************************************
869  SUBROUTINE print_nr_ions_gal21(nr_ions, mm_section, para_env, print_oh, &
870  print_h3o, print_o)
871  INTEGER, INTENT(INOUT) :: nr_ions
872  TYPE(section_vals_type), POINTER :: mm_section
873  TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
874  LOGICAL, INTENT(IN) :: print_oh, print_h3o, print_o
875 
876  INTEGER :: iw
877  TYPE(cp_logger_type), POINTER :: logger
878 
879  NULLIFY (logger)
880 
881  CALL para_env%sum(nr_ions)
882  logger => cp_get_default_logger()
883 
884  iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%PROGRAM_RUN_INFO", &
885  extension=".mmLog")
886 
887  IF (iw > 0 .AND. nr_ions > 0 .AND. print_oh) THEN
888  WRITE (iw, '(/,A,T71,I10,/)') " gal21: number of OH- ions at surface", nr_ions
889  END IF
890  IF (iw > 0 .AND. nr_ions > 0 .AND. print_h3o) THEN
891  WRITE (iw, '(/,A,T71,I10,/)') " gal21: number of H3O+ ions at surface", nr_ions
892  END IF
893  IF (iw > 0 .AND. nr_ions > 0 .AND. print_o) THEN
894  WRITE (iw, '(/,A,T71,I10,/)') " gal21: number of O^2- ions at surface", nr_ions
895  END IF
896 
897  CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%PROGRAM_RUN_INFO")
898 
899  END SUBROUTINE print_nr_ions_gal21
900 
901 END MODULE manybody_gal21
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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_gal21
Implementation of the GAL21 potential.
Definition: manybody_gal21.F:13
manybody_gal21::destroy_gal21_arrays
subroutine, public destroy_gal21_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
...
Definition: manybody_gal21.F:844
manybody_gal21::gal21_forces
subroutine, public gal21_forces(gal21, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, use_virial, cell, particle_set)
forces generated by the GAL2119 potential
Definition: manybody_gal21.F:398
manybody_gal21::setup_gal21_arrays
subroutine, public setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
...
Definition: manybody_gal21.F:743
manybody_gal21::print_nr_ions_gal21
subroutine, public print_nr_ions_gal21(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_gal21.F:871
manybody_gal21::gal21_energy
subroutine, public gal21_energy(pot_loc, gal21, r_last_update_pbc, iparticle, jparticle, cell, particle_set, mm_section)
Main part of the energy evaluation of GAL2119.
Definition: manybody_gal21.F:61
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::gal21_type
integer, parameter, public gal21_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