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