(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
13MODULE 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
45CONTAINS
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
801END MODULE manybody_gal
cell_types::pbc
Definition cell_types.F:92
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
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
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
fist_neighbor_list_types::fist_neighbor_type
Definition fist_neighbor_list_types.F:48
fist_neighbor_list_types::neighbor_kind_pairs_type
Definition fist_neighbor_list_types.F:26
fist_nonbond_env_types::pos_type
Definition fist_nonbond_env_types.F:44
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
pair_potential_types::gal_pot_type
Definition pair_potential_types.F:282
pair_potential_types::pair_potential_pp_type
Definition pair_potential_types.F:412
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:379
particle_types::particle_type
Definition particle_types.F:35