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