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