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