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