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