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