Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par History
10 : !> CJM APRIL-30-2009: only uses fist_env
11 : !> Teodoro Laino [tlaino] - 05.2009 : Generalization to different Ewald
12 : !> methods (initial framework)
13 : !> \author CJM
14 : ! **************************************************************************************************
15 :
16 : MODULE fist_pol_scf
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind
19 : USE cell_types, ONLY: cell_type
20 : USE cp_log_handling, ONLY: cp_get_default_logger,&
21 : cp_logger_type
22 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
23 : cp_print_key_unit_nr
24 : USE distribution_1d_types, ONLY: distribution_1d_type
25 : USE ewald_environment_types, ONLY: ewald_env_get,&
26 : ewald_environment_type
27 : USE ewald_pw_types, ONLY: ewald_pw_type
28 : USE ewalds_multipole, ONLY: ewald_multipole_evaluate
29 : USE fist_energy_types, ONLY: fist_energy_type
30 : USE fist_nonbond_env_types, ONLY: fist_nonbond_env_type
31 : USE input_constants, ONLY: do_fist_pol_cg,&
32 : do_fist_pol_sc
33 : USE input_section_types, ONLY: section_vals_type
34 : USE kinds, ONLY: dp
35 : USE multipole_types, ONLY: multipole_type
36 : USE particle_types, ONLY: particle_type
37 : USE pw_poisson_types, ONLY: do_ewald_ewald,&
38 : do_ewald_pme,&
39 : do_ewald_spme
40 : #include "./base/base_uses.f90"
41 :
42 : IMPLICIT NONE
43 : PRIVATE
44 : LOGICAL, PRIVATE :: debug_this_module = .FALSE.
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_pol_scf'
46 :
47 : PUBLIC :: fist_pol_evaluate
48 :
49 : CONTAINS
50 :
51 : ! **************************************************************************************************
52 : !> \brief Main driver for evaluating energy/forces in a polarizable forcefield
53 : !> \param atomic_kind_set ...
54 : !> \param multipoles ...
55 : !> \param ewald_env ...
56 : !> \param ewald_pw ...
57 : !> \param fist_nonbond_env ...
58 : !> \param cell ...
59 : !> \param particle_set ...
60 : !> \param local_particles ...
61 : !> \param thermo ...
62 : !> \param vg_coulomb ...
63 : !> \param pot_nonbond ...
64 : !> \param f_nonbond ...
65 : !> \param fg_coulomb ...
66 : !> \param use_virial ...
67 : !> \param pv_g ...
68 : !> \param pv_nonbond ...
69 : !> \param mm_section ...
70 : !> \param do_ipol ...
71 : !> \author Toon.Verstraelen@gmail.com (2010-03-01)
72 : ! **************************************************************************************************
73 104 : SUBROUTINE fist_pol_evaluate(atomic_kind_set, multipoles, ewald_env, &
74 : ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, &
75 104 : thermo, vg_coulomb, pot_nonbond, f_nonbond, fg_coulomb, use_virial, &
76 : pv_g, pv_nonbond, mm_section, do_ipol)
77 :
78 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
79 : TYPE(multipole_type), POINTER :: multipoles
80 : TYPE(ewald_environment_type), POINTER :: ewald_env
81 : TYPE(ewald_pw_type), POINTER :: ewald_pw
82 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
83 : TYPE(cell_type), POINTER :: cell
84 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
85 : TYPE(distribution_1d_type), POINTER :: local_particles
86 : TYPE(fist_energy_type), POINTER :: thermo
87 : REAL(KIND=dp) :: vg_coulomb, pot_nonbond
88 : REAL(KIND=dp), DIMENSION(:, :) :: f_nonbond, fg_coulomb
89 : LOGICAL, INTENT(IN) :: use_virial
90 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_g, pv_nonbond
91 : TYPE(section_vals_type), POINTER :: mm_section
92 : INTEGER :: do_ipol
93 :
94 140 : SELECT CASE (do_ipol)
95 : CASE (do_fist_pol_sc)
96 : CALL fist_pol_evaluate_sc(atomic_kind_set, multipoles, ewald_env, &
97 : ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, &
98 : thermo, vg_coulomb, pot_nonbond, f_nonbond, fg_coulomb, use_virial, &
99 36 : pv_g, pv_nonbond, mm_section)
100 : CASE (do_fist_pol_cg)
101 : CALL fist_pol_evaluate_cg(atomic_kind_set, multipoles, ewald_env, &
102 : ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, &
103 : thermo, vg_coulomb, pot_nonbond, f_nonbond, fg_coulomb, use_virial, &
104 104 : pv_g, pv_nonbond, mm_section)
105 : END SELECT
106 :
107 104 : END SUBROUTINE fist_pol_evaluate
108 :
109 : ! **************************************************************************************************
110 : !> \brief Self-consistent solver for a polarizable force-field
111 : !> \param atomic_kind_set ...
112 : !> \param multipoles ...
113 : !> \param ewald_env ...
114 : !> \param ewald_pw ...
115 : !> \param fist_nonbond_env ...
116 : !> \param cell ...
117 : !> \param particle_set ...
118 : !> \param local_particles ...
119 : !> \param thermo ...
120 : !> \param vg_coulomb ...
121 : !> \param pot_nonbond ...
122 : !> \param f_nonbond ...
123 : !> \param fg_coulomb ...
124 : !> \param use_virial ...
125 : !> \param pv_g ...
126 : !> \param pv_nonbond ...
127 : !> \param mm_section ...
128 : !> \author Toon.Verstraelen@gmail.com (2010-03-01)
129 : !> \note
130 : !> Method: Given an initial guess of the induced dipoles, the electrostatic
131 : !> field is computed at each dipole. Then new induced dipoles are computed
132 : !> following p = alpha x E. This is repeated until a convergence criterion is
133 : !> met. The convergence is measured as the RSMD of the derivatives of the
134 : !> electrostatic energy (including dipole self-energy) towards the components
135 : !> of the dipoles.
136 : ! **************************************************************************************************
137 36 : SUBROUTINE fist_pol_evaluate_sc(atomic_kind_set, multipoles, ewald_env, ewald_pw, &
138 : fist_nonbond_env, cell, particle_set, local_particles, thermo, vg_coulomb, &
139 36 : pot_nonbond, f_nonbond, fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section)
140 :
141 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
142 : TYPE(multipole_type), POINTER :: multipoles
143 : TYPE(ewald_environment_type), POINTER :: ewald_env
144 : TYPE(ewald_pw_type), POINTER :: ewald_pw
145 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
146 : TYPE(cell_type), POINTER :: cell
147 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
148 : TYPE(distribution_1d_type), POINTER :: local_particles
149 : TYPE(fist_energy_type), POINTER :: thermo
150 : REAL(KIND=dp) :: vg_coulomb, pot_nonbond
151 : REAL(KIND=dp), DIMENSION(:, :) :: f_nonbond, fg_coulomb
152 : LOGICAL, INTENT(IN) :: use_virial
153 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_g, pv_nonbond
154 : TYPE(section_vals_type), POINTER :: mm_section
155 :
156 : CHARACTER(len=*), PARAMETER :: routineN = 'fist_pol_evaluate_sc'
157 :
158 : INTEGER :: ewald_type, handle, i, iatom, ii, ikind, &
159 : iter, iw, iw2, j, max_ipol_iter, &
160 : natom_of_kind, natoms, nkind, ntot
161 36 : INTEGER, DIMENSION(:), POINTER :: atom_list
162 : LOGICAL :: iwarn
163 : REAL(KIND=dp) :: apol, cpol, eps_pol, pot_nonbond_local, &
164 : rmsd, tmp_trace
165 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: efield1, efield2
166 : TYPE(atomic_kind_type), POINTER :: atomic_kind
167 : TYPE(cp_logger_type), POINTER :: logger
168 :
169 36 : CALL timeset(routineN, handle)
170 36 : NULLIFY (logger, atomic_kind)
171 36 : logger => cp_get_default_logger()
172 :
173 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%ITER_INFO", &
174 36 : extension=".mmLog")
175 : iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
176 36 : extension=".mmLog")
177 :
178 : CALL ewald_env_get(ewald_env, max_ipol_iter=max_ipol_iter, eps_pol=eps_pol, &
179 36 : ewald_type=ewald_type)
180 :
181 36 : natoms = SIZE(particle_set)
182 108 : ALLOCATE (efield1(3, natoms))
183 108 : ALLOCATE (efield2(9, natoms))
184 :
185 36 : nkind = SIZE(atomic_kind_set)
186 36 : IF (iw > 0) THEN
187 12 : WRITE (iw, FMT='(/,T2,"POL_SCF|" ,"Method: self-consistent")')
188 12 : WRITE (iw, FMT='(T2,"POL_SCF| "," Iteration",7X,"Conv.",T49,"Electrostatic & Induction Energy")')
189 : END IF
190 294 : pol_scf: DO iter = 1, max_ipol_iter
191 : ! Evaluate the electrostatic with Ewald schemes
192 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
193 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
194 : multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
195 : do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
196 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
197 294 : efield1=efield1, efield2=efield2)
198 294 : CALL logger%para_env%sum(pot_nonbond_local)
199 :
200 : ! compute the new dipoles, qudrupoles, and check for convergence
201 294 : ntot = 0
202 294 : rmsd = 0.0_dp
203 294 : thermo%e_induction = 0.0_dp
204 1198 : DO ikind = 1, nkind
205 904 : atomic_kind => atomic_kind_set(ikind)
206 904 : CALL get_atomic_kind(atomic_kind, apol=apol, cpol=cpol, natom=natom_of_kind, atom_list=atom_list)
207 : ! ignore atoms with dipole and quadrupole polarizability zero
208 904 : IF (apol == 0 .AND. cpol == 0) CYCLE
209 : ! increment counter correctly
210 562 : IF (apol /= 0) ntot = ntot + natom_of_kind
211 562 : IF (cpol /= 0) ntot = ntot + natom_of_kind
212 :
213 35924 : DO iatom = 1, natom_of_kind
214 34164 : ii = atom_list(iatom)
215 34164 : IF (apol /= 0) THEN
216 136656 : DO i = 1, 3
217 : ! the rmsd of the derivatives of the energy towards the
218 : ! components of the atomic dipole moments
219 136656 : rmsd = rmsd + (multipoles%dipoles(i, ii)/apol - efield1(i, ii))**2
220 : END DO
221 : END IF
222 34164 : IF (cpol /= 0) THEN
223 0 : rmsd = rmsd + (multipoles%quadrupoles(1, 1, ii)/cpol - efield2(1, ii))**2
224 0 : rmsd = rmsd + (multipoles%quadrupoles(2, 1, ii)/cpol - efield2(2, ii))**2
225 0 : rmsd = rmsd + (multipoles%quadrupoles(3, 1, ii)/cpol - efield2(3, ii))**2
226 0 : rmsd = rmsd + (multipoles%quadrupoles(1, 2, ii)/cpol - efield2(4, ii))**2
227 0 : rmsd = rmsd + (multipoles%quadrupoles(2, 2, ii)/cpol - efield2(5, ii))**2
228 0 : rmsd = rmsd + (multipoles%quadrupoles(3, 2, ii)/cpol - efield2(6, ii))**2
229 0 : rmsd = rmsd + (multipoles%quadrupoles(1, 3, ii)/cpol - efield2(7, ii))**2
230 0 : rmsd = rmsd + (multipoles%quadrupoles(2, 3, ii)/cpol - efield2(8, ii))**2
231 0 : rmsd = rmsd + (multipoles%quadrupoles(3, 3, ii)/cpol - efield2(9, ii))**2
232 : END IF
233 : ! compute dipole
234 136656 : multipoles%dipoles(:, ii) = apol*efield1(:, ii)
235 : ! compute quadrupole
236 34164 : IF (cpol /= 0) THEN
237 0 : multipoles%quadrupoles(1, 1, ii) = cpol*efield2(1, ii)
238 0 : multipoles%quadrupoles(2, 1, ii) = cpol*efield2(2, ii)
239 0 : multipoles%quadrupoles(3, 1, ii) = cpol*efield2(3, ii)
240 0 : multipoles%quadrupoles(1, 2, ii) = cpol*efield2(4, ii)
241 0 : multipoles%quadrupoles(2, 2, ii) = cpol*efield2(5, ii)
242 0 : multipoles%quadrupoles(3, 2, ii) = cpol*efield2(6, ii)
243 0 : multipoles%quadrupoles(1, 3, ii) = cpol*efield2(7, ii)
244 0 : multipoles%quadrupoles(2, 3, ii) = cpol*efield2(8, ii)
245 0 : multipoles%quadrupoles(3, 3, ii) = cpol*efield2(9, ii)
246 : END IF
247 : ! Compute the new induction term while we are here
248 34164 : IF (apol /= 0) THEN
249 : thermo%e_induction = thermo%e_induction + &
250 : DOT_PRODUCT(multipoles%dipoles(:, ii), &
251 136656 : multipoles%dipoles(:, ii))/apol/2.0_dp
252 : END IF
253 35068 : IF (cpol /= 0) THEN
254 : tmp_trace = 0._dp
255 0 : DO i = 1, 3
256 0 : DO j = 1, 3
257 : tmp_trace = tmp_trace + &
258 0 : multipoles%quadrupoles(i, j, ii)*multipoles%quadrupoles(i, j, ii)
259 : END DO
260 : END DO
261 0 : thermo%e_induction = thermo%e_induction + tmp_trace/cpol/6.0_dp
262 : END IF
263 : END DO
264 : END DO
265 294 : rmsd = SQRT(rmsd/REAL(ntot, KIND=dp))
266 294 : IF (iw > 0) THEN
267 : ! print the energy that is minimized (this is electrostatic + induction)
268 119 : WRITE (iw, FMT='(T5,"POL_SCF|",5X,I5,5X,E12.6,T61,F20.10)') iter, &
269 238 : rmsd, vg_coulomb + pot_nonbond_local + thermo%e_induction
270 : END IF
271 294 : IF (rmsd <= eps_pol) THEN
272 36 : IF (iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|",1X,"Self-consistent Polarization achieved.")')
273 : EXIT pol_scf
274 : END IF
275 :
276 258 : iwarn = ((rmsd > eps_pol) .AND. (iter == max_ipol_iter))
277 258 : IF (iwarn .AND. iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|",1X,"Self-consistent Polarization not converged!")')
278 258 : CPWARN_IF(iwarn, "Self-consistent Polarization not converged!")
279 : END DO pol_scf
280 :
281 : ! Now evaluate after convergence to obtain forces and converged energies
282 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
283 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
284 : multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., &
285 : do_stress=use_virial, do_efield=.FALSE., iw2=iw2, do_debug=.FALSE., &
286 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
287 : forces_local=fg_coulomb, forces_glob=f_nonbond, &
288 36 : pv_local=pv_g, pv_glob=pv_nonbond)
289 36 : pot_nonbond = pot_nonbond + pot_nonbond_local
290 36 : CALL logger%para_env%sum(pot_nonbond_local)
291 :
292 36 : IF (iw > 0) THEN
293 : ! print the energy that is minimized (this is electrostatic + induction)
294 : WRITE (iw, FMT='(T5,"POL_SCF|",5X,"Final",T61,F20.10,/)') &
295 12 : vg_coulomb + pot_nonbond_local + thermo%e_induction
296 : END IF
297 :
298 : ! Deallocate working arrays
299 36 : DEALLOCATE (efield1)
300 36 : DEALLOCATE (efield2)
301 : CALL cp_print_key_finished_output(iw2, logger, mm_section, &
302 36 : "PRINT%EWALD_INFO")
303 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
304 36 : "PRINT%ITER_INFO")
305 :
306 36 : CALL timestop(handle)
307 72 : END SUBROUTINE fist_pol_evaluate_sc
308 :
309 : ! **************************************************************************************************
310 : !> \brief Conjugate-gradient solver for a polarizable force-field
311 : !> \param atomic_kind_set ...
312 : !> \param multipoles ...
313 : !> \param ewald_env ...
314 : !> \param ewald_pw ...
315 : !> \param fist_nonbond_env ...
316 : !> \param cell ...
317 : !> \param particle_set ...
318 : !> \param local_particles ...
319 : !> \param thermo ...
320 : !> \param vg_coulomb ...
321 : !> \param pot_nonbond ...
322 : !> \param f_nonbond ...
323 : !> \param fg_coulomb ...
324 : !> \param use_virial ...
325 : !> \param pv_g ...
326 : !> \param pv_nonbond ...
327 : !> \param mm_section ...
328 : !> \author Toon.Verstraelen@gmail.com (2010-03-01)
329 : !> \note
330 : !> Method: The dipoles are found by minimizing the sum of the electrostatic
331 : !> and the induction energy directly using a conjugate gradient method. This
332 : !> routine assumes that the energy is a quadratic function of the dipoles.
333 : !> Finding the minimum is then done by solving a linear system. This will
334 : !> not work for polarizable force fields that include hyperpolarizability.
335 : !>
336 : !> The implementation of the conjugate gradient solver for linear systems
337 : !> is described in chapter 2.7 Sparse Linear Systems, under the section
338 : !> "Conjugate Gradient Method for a Sparse System". Although the inducible
339 : !> dipoles are the solution of a dense linear system, the same algorithm is
340 : !> still recommended for this situation. One does not have access to the
341 : !> entire hardness kernel to compute the solution with conventional linear
342 : !> algebra routines, but one only has a function that computes the dot
343 : !> product of the hardness kernel and a vector. (This is the routine that
344 : !> computes the electrostatic field at the atoms for a given vector of
345 : !> inducible dipoles.) Given such function, the conjugate gradient method
346 : !> is an efficient way to compute the solution of a linear system, and it
347 : !> scales well towards many degrees of freedom in terms of efficiency and
348 : !> memory usage.
349 : ! **************************************************************************************************
350 68 : SUBROUTINE fist_pol_evaluate_cg(atomic_kind_set, multipoles, ewald_env, ewald_pw, &
351 : fist_nonbond_env, cell, particle_set, local_particles, thermo, vg_coulomb, &
352 68 : pot_nonbond, f_nonbond, fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section)
353 :
354 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
355 : TYPE(multipole_type), POINTER :: multipoles
356 : TYPE(ewald_environment_type), POINTER :: ewald_env
357 : TYPE(ewald_pw_type), POINTER :: ewald_pw
358 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
359 : TYPE(cell_type), POINTER :: cell
360 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
361 : TYPE(distribution_1d_type), POINTER :: local_particles
362 : TYPE(fist_energy_type), POINTER :: thermo
363 : REAL(KIND=dp) :: vg_coulomb, pot_nonbond
364 : REAL(KIND=dp), DIMENSION(:, :) :: f_nonbond, fg_coulomb
365 : LOGICAL, INTENT(IN) :: use_virial
366 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_g, pv_nonbond
367 : TYPE(section_vals_type), POINTER :: mm_section
368 :
369 : CHARACTER(len=*), PARAMETER :: routineN = 'fist_pol_evaluate_cg'
370 :
371 : INTEGER :: ewald_type, handle, i, iatom, ii, ikind, &
372 : iter, iw, iw2, max_ipol_iter, &
373 : natom_of_kind, natoms, nkind, ntot
374 68 : INTEGER, DIMENSION(:), POINTER :: atom_list
375 : LOGICAL :: iwarn
376 : REAL(KIND=dp) :: alpha, apol, beta, denom, eps_pol, &
377 : pot_nonbond_local, rmsd
378 68 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: conjugate, conjugate_applied, efield1, &
379 68 : efield1_ext, residual, tmp_dipoles
380 : TYPE(atomic_kind_type), POINTER :: atomic_kind
381 : TYPE(cp_logger_type), POINTER :: logger
382 :
383 68 : CALL timeset(routineN, handle)
384 68 : NULLIFY (logger, atomic_kind)
385 68 : logger => cp_get_default_logger()
386 :
387 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%ITER_INFO", &
388 68 : extension=".mmLog")
389 : iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
390 68 : extension=".mmLog")
391 :
392 : CALL ewald_env_get(ewald_env, max_ipol_iter=max_ipol_iter, eps_pol=eps_pol, &
393 68 : ewald_type=ewald_type)
394 :
395 : ! allocate work arrays
396 68 : natoms = SIZE(particle_set)
397 204 : ALLOCATE (efield1(3, natoms))
398 136 : ALLOCATE (tmp_dipoles(3, natoms))
399 136 : ALLOCATE (residual(3, natoms))
400 136 : ALLOCATE (conjugate(3, natoms))
401 136 : ALLOCATE (conjugate_applied(3, natoms))
402 136 : ALLOCATE (efield1_ext(3, natoms))
403 :
404 : ! Compute the 'external' electrostatic field (all inducible dipoles
405 : ! equal to zero). this is required for the conjugate gradient solver.
406 : ! We assume that all dipoles are inducible dipoles.
407 31020 : tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles
408 31020 : multipoles%dipoles = 0.0_dp
409 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
410 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
411 : multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
412 : do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
413 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
414 68 : efield1=efield1_ext)
415 31020 : multipoles%dipoles = tmp_dipoles ! restore backup
416 :
417 : ! Compute the electric field with the initial guess of the dipoles.
418 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
419 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
420 : multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
421 : do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
422 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
423 68 : efield1=efield1)
424 :
425 : ! Compute the first residual explicitly.
426 68 : nkind = SIZE(atomic_kind_set)
427 68 : ntot = 0
428 31020 : residual = 0.0_dp
429 232 : DO ikind = 1, nkind
430 164 : atomic_kind => atomic_kind_set(ikind)
431 164 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
432 : ! ignore atoms with polarizability zero
433 164 : IF (apol == 0) CYCLE
434 78 : ntot = ntot + natom_of_kind
435 4466 : DO iatom = 1, natom_of_kind
436 4156 : ii = atom_list(iatom)
437 16788 : DO i = 1, 3
438 : ! residual = b - A x
439 16624 : residual(i, ii) = efield1(i, ii) - multipoles%dipoles(i, ii)/apol
440 : END DO
441 : END DO
442 : END DO
443 : ! The first conjugate residual is equal to the residual.
444 31020 : conjugate(:, :) = residual
445 :
446 68 : IF (iw > 0) THEN
447 34 : WRITE (iw, FMT="(/,T2,A,T63,A)") "POL_SCF| Method", "conjugate gradient"
448 34 : WRITE (iw, FMT="(T2,A,T26,A,T49,A)") "POL_SCF| Iteration", &
449 68 : "Convergence", "Electrostatic & Induction Energy"
450 : END IF
451 308 : pol_scf: DO iter = 1, max_ipol_iter
452 308 : IF (debug_this_module) THEN
453 : ! In principle the residual must not be computed explicitly any more. It
454 : ! is obtained in an indirect way below. When the DEBUG flag is set, the
455 : ! explicit residual is computed and compared with the implicitly derived
456 : ! residual as a double check.
457 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
458 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
459 : multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
460 : do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
461 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
462 0 : efield1=efield1)
463 : ! inapropriate use of denom to check the error on the residual
464 0 : denom = 0.0_dp
465 : END IF
466 308 : rmsd = 0.0_dp
467 : ! Compute the rmsd of the residual.
468 1090 : DO ikind = 1, nkind
469 782 : atomic_kind => atomic_kind_set(ikind)
470 782 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
471 : ! ignore atoms with polarizability zero
472 782 : IF (apol == 0) CYCLE
473 22890 : DO iatom = 1, natom_of_kind
474 21370 : ii = atom_list(iatom)
475 86262 : DO i = 1, 3
476 : ! residual = b - A x
477 64110 : rmsd = rmsd + residual(i, ii)**2
478 85480 : IF (debug_this_module) THEN
479 : denom = denom + (residual(i, ii) - (efield1(i, ii) - &
480 0 : multipoles%dipoles(i, ii)/apol))**2
481 : END IF
482 : END DO
483 : END DO
484 : END DO
485 308 : rmsd = SQRT(rmsd/ntot)
486 308 : IF (iw > 0) THEN
487 154 : WRITE (iw, FMT="(T2,A,T11,I9,T22,E15.6,T67,A)") "POL_SCF|", iter, rmsd, "(not computed)"
488 154 : IF (debug_this_module) THEN
489 0 : denom = SQRT(denom/ntot)
490 0 : WRITE (iw, FMT="(T2,A,T66,E15.6)") "POL_SCF| Error on implicit residual", denom
491 : END IF
492 : END IF
493 :
494 : ! Apply the hardness kernel to the conjugate residual.
495 : ! We assume that all non-zero dipoles are inducible dipoles.
496 153100 : tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles
497 153100 : multipoles%dipoles = conjugate
498 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
499 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
500 : multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
501 : do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
502 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
503 308 : efield1=conjugate_applied)
504 153100 : multipoles%dipoles = tmp_dipoles ! restore backup
505 153100 : conjugate_applied(:, :) = efield1_ext - conjugate_applied
506 :
507 : ! Finish conjugate_applied and compute alpha from the conjugate gradient algorithm.
508 308 : alpha = 0.0_dp
509 308 : denom = 0.0_dp
510 1090 : DO ikind = 1, nkind
511 782 : atomic_kind => atomic_kind_set(ikind)
512 782 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
513 : ! ignore atoms with polarizability zero
514 782 : IF (apol == 0) CYCLE
515 22890 : DO iatom = 1, natom_of_kind
516 21370 : ii = atom_list(iatom)
517 85480 : DO i = 1, 3
518 85480 : conjugate_applied(i, ii) = conjugate_applied(i, ii) + conjugate(i, ii)/apol
519 : END DO
520 85480 : alpha = alpha + DOT_PRODUCT(residual(:, ii), residual(:, ii))
521 86262 : denom = denom + DOT_PRODUCT(conjugate(:, ii), conjugate_applied(:, ii))
522 : END DO
523 : END DO
524 308 : alpha = alpha/denom
525 :
526 : ! Compute the new residual and beta from the conjugate gradient method.
527 308 : beta = 0.0_dp
528 308 : denom = 0.0_dp
529 1090 : DO ikind = 1, nkind
530 782 : atomic_kind => atomic_kind_set(ikind)
531 782 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
532 782 : IF (apol == 0) CYCLE
533 22890 : DO iatom = 1, natom_of_kind
534 21370 : ii = atom_list(iatom)
535 85480 : denom = denom + DOT_PRODUCT(residual(:, ii), residual(:, ii))
536 85480 : DO i = 1, 3
537 85480 : residual(i, ii) = residual(i, ii) - alpha*conjugate_applied(i, ii)
538 : END DO
539 86262 : beta = beta + DOT_PRODUCT(residual(:, ii), residual(:, ii))
540 : END DO
541 : END DO
542 308 : beta = beta/denom
543 :
544 : ! Compute the new dipoles, the new conjugate residual, and the induction
545 : ! energy.
546 308 : thermo%e_induction = 0.0_dp
547 1090 : DO ikind = 1, nkind
548 782 : atomic_kind => atomic_kind_set(ikind)
549 782 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
550 : ! ignore atoms with polarizability zero
551 782 : IF (apol == 0) CYCLE
552 22890 : DO iatom = 1, natom_of_kind
553 21370 : ii = atom_list(iatom)
554 86262 : DO i = 1, 3
555 64110 : multipoles%dipoles(i, ii) = multipoles%dipoles(i, ii) + alpha*conjugate(i, ii)
556 64110 : conjugate(i, ii) = residual(i, ii) + beta*conjugate(i, ii)
557 85480 : thermo%e_induction = thermo%e_induction + multipoles%dipoles(i, ii)**2/apol/2.0_dp
558 : END DO
559 : END DO
560 : END DO
561 :
562 : ! Quit if rmsd is low enough
563 308 : IF (rmsd <= eps_pol) THEN
564 68 : IF (iw > 0) WRITE (iw, FMT="(T2,A)") "POL_SCF| Self-consistent polarization converged"
565 : EXIT pol_scf
566 : END IF
567 :
568 : ! Print warning when not converged
569 240 : iwarn = ((rmsd > eps_pol) .AND. (iter >= max_ipol_iter))
570 0 : IF (iwarn) THEN
571 0 : IF (iw > 0) THEN
572 : WRITE (iw, FMT="(T2,A,I0,A,ES9.3)") &
573 0 : "POL_SCF| Self-consistent polarization not converged in ", max_ipol_iter, &
574 0 : " steps to ", eps_pol
575 : END IF
576 0 : CPWARN("Self-consistent Polarization not converged!")
577 : END IF
578 : END DO pol_scf
579 :
580 68 : IF (debug_this_module) THEN
581 : ! Now evaluate after convergence to obtain forces and converged energies
582 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
583 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
584 : multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., &
585 : do_stress=use_virial, do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
586 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
587 : forces_local=fg_coulomb, forces_glob=f_nonbond, &
588 0 : pv_local=pv_g, pv_glob=pv_nonbond, efield1=efield1)
589 :
590 : ! Do a final check on the convergence: compute the residual explicitely
591 0 : rmsd = 0.0_dp
592 0 : DO ikind = 1, nkind
593 0 : atomic_kind => atomic_kind_set(ikind)
594 0 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
595 : ! ignore atoms with polarizability zero
596 0 : IF (apol == 0) CYCLE
597 0 : DO iatom = 1, natom_of_kind
598 0 : ii = atom_list(iatom)
599 0 : DO i = 1, 3
600 : ! residual = b - A x
601 0 : rmsd = rmsd + (efield1(i, ii) - multipoles%dipoles(i, ii)/apol)**2
602 : END DO
603 : END DO
604 : END DO
605 0 : rmsd = SQRT(rmsd/ntot)
606 0 : IF (iw > 0) WRITE (iw, FMT="(T2,A,T66,E15.6)") "POL_SCF| Final RMSD of residual", rmsd
607 : ! Stop program when congergence is not reached after all
608 0 : IF (rmsd > eps_pol) THEN
609 0 : CPABORT("Error in the conjugate gradient method for self-consistent polarization!")
610 : END IF
611 : ELSE
612 : ! Now evaluate after convergence to obtain forces and converged energies
613 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
614 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
615 : multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., &
616 : do_stress=use_virial, do_efield=.FALSE., iw2=iw2, do_debug=.FALSE., &
617 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
618 : forces_local=fg_coulomb, forces_glob=f_nonbond, &
619 68 : pv_local=pv_g, pv_glob=pv_nonbond)
620 : END IF
621 68 : pot_nonbond = pot_nonbond + pot_nonbond_local
622 68 : CALL logger%para_env%sum(pot_nonbond_local)
623 :
624 102 : IF (iw > 0) WRITE (iw, FMT="(T2,A,T61,F20.10)") "POL_SCF| Final", &
625 68 : vg_coulomb + pot_nonbond_local + thermo%e_induction
626 :
627 : ! Deallocate working arrays
628 68 : DEALLOCATE (efield1)
629 68 : DEALLOCATE (tmp_dipoles)
630 68 : DEALLOCATE (residual)
631 68 : DEALLOCATE (conjugate)
632 68 : DEALLOCATE (conjugate_applied)
633 68 : DEALLOCATE (efield1_ext)
634 : CALL cp_print_key_finished_output(iw2, logger, mm_section, &
635 68 : "PRINT%EWALD_INFO")
636 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
637 68 : "PRINT%ITER_INFO")
638 :
639 68 : CALL timestop(handle)
640 :
641 136 : END SUBROUTINE fist_pol_evaluate_cg
642 :
643 : ! **************************************************************************************************
644 : !> \brief Main driver for evaluating electrostatic in polarible forcefields
645 : !> All the dependence on the Ewald method should go here!
646 : !> \param ewald_type ...
647 : !> \param ewald_env ...
648 : !> \param ewald_pw ...
649 : !> \param fist_nonbond_env ...
650 : !> \param cell ...
651 : !> \param particle_set ...
652 : !> \param local_particles ...
653 : !> \param vg_coulomb ...
654 : !> \param pot_nonbond ...
655 : !> \param thermo ...
656 : !> \param multipoles ...
657 : !> \param do_correction_bonded ...
658 : !> \param do_forces ...
659 : !> \param do_stress ...
660 : !> \param do_efield ...
661 : !> \param iw2 ...
662 : !> \param do_debug ...
663 : !> \param atomic_kind_set ...
664 : !> \param mm_section ...
665 : !> \param efield0 ...
666 : !> \param efield1 ...
667 : !> \param efield2 ...
668 : !> \param forces_local ...
669 : !> \param forces_glob ...
670 : !> \param pv_local ...
671 : !> \param pv_glob ...
672 : !> \author Teodoro Laino [tlaino] 05.2009
673 : ! **************************************************************************************************
674 1684 : SUBROUTINE eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, &
675 : cell, particle_set, local_particles, vg_coulomb, pot_nonbond, thermo, &
676 : multipoles, do_correction_bonded, do_forces, do_stress, do_efield, iw2, &
677 1684 : do_debug, atomic_kind_set, mm_section, efield0, efield1, efield2, forces_local, &
678 842 : forces_glob, pv_local, pv_glob)
679 :
680 : INTEGER, INTENT(IN) :: ewald_type
681 : TYPE(ewald_environment_type), POINTER :: ewald_env
682 : TYPE(ewald_pw_type), POINTER :: ewald_pw
683 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
684 : TYPE(cell_type), POINTER :: cell
685 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
686 : TYPE(distribution_1d_type), POINTER :: local_particles
687 : REAL(KIND=dp), INTENT(OUT) :: vg_coulomb, pot_nonbond
688 : TYPE(fist_energy_type), POINTER :: thermo
689 : TYPE(multipole_type), POINTER :: multipoles
690 : LOGICAL, INTENT(IN) :: do_correction_bonded, do_forces, &
691 : do_stress, do_efield
692 : INTEGER, INTENT(IN) :: iw2
693 : LOGICAL, INTENT(IN) :: do_debug
694 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
695 : TYPE(section_vals_type), POINTER :: mm_section
696 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: efield0
697 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: efield1, efield2
698 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
699 : OPTIONAL :: forces_local, forces_glob, pv_local, &
700 : pv_glob
701 :
702 : CHARACTER(len=*), PARAMETER :: routineN = 'eval_pol_ewald'
703 :
704 : INTEGER :: handle
705 :
706 842 : CALL timeset(routineN, handle)
707 :
708 842 : pot_nonbond = 0.0_dp ! Initialization..
709 842 : vg_coulomb = 0.0_dp ! Initialization..
710 1684 : SELECT CASE (ewald_type)
711 : CASE (do_ewald_ewald)
712 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, fist_nonbond_env, cell, &
713 : particle_set, local_particles, vg_coulomb, pot_nonbond, thermo%e_neut, &
714 : thermo%e_self, multipoles%task, do_correction_bonded=do_correction_bonded, &
715 : do_forces=do_forces, do_stress=do_stress, do_efield=do_efield, &
716 : radii=multipoles%radii, charges=multipoles%charges, &
717 : dipoles=multipoles%dipoles, quadrupoles=multipoles%quadrupoles, &
718 : forces_local=forces_local, forces_glob=forces_glob, pv_local=pv_local, &
719 : pv_glob=pv_glob, iw=iw2, do_debug=do_debug, atomic_kind_set=atomic_kind_set, &
720 5288 : mm_section=mm_section, efield0=efield0, efield1=efield1, efield2=efield2)
721 : CASE (do_ewald_pme)
722 0 : CPABORT("Multipole Ewald not yet implemented within a PME scheme!")
723 : CASE (do_ewald_spme)
724 842 : CPABORT("Multipole Ewald not yet implemented within a SPME scheme!")
725 : END SELECT
726 842 : CALL timestop(handle)
727 842 : END SUBROUTINE eval_pol_ewald
728 :
729 : END MODULE fist_pol_scf
|
