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