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 : !> \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 : IF (iwarn) &
279 0 : CPWARN("Self-consistent Polarization not converged! ")
280 : END DO pol_scf
281 :
282 : ! Now evaluate after convergence to obtain forces and converged energies
283 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
284 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
285 : multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., &
286 : do_stress=use_virial, do_efield=.FALSE., iw2=iw2, do_debug=.FALSE., &
287 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
288 : forces_local=fg_coulomb, forces_glob=f_nonbond, &
289 36 : pv_local=pv_g, pv_glob=pv_nonbond)
290 36 : pot_nonbond = pot_nonbond + pot_nonbond_local
291 36 : CALL logger%para_env%sum(pot_nonbond_local)
292 :
293 36 : IF (iw > 0) THEN
294 : ! print the energy that is minimized (this is electrostatic + induction)
295 : WRITE (iw, FMT='(T5,"POL_SCF|",5X,"Final",T61,F20.10,/)') &
296 12 : vg_coulomb + pot_nonbond_local + thermo%e_induction
297 : END IF
298 :
299 : ! Deallocate working arrays
300 36 : DEALLOCATE (efield1)
301 36 : DEALLOCATE (efield2)
302 : CALL cp_print_key_finished_output(iw2, logger, mm_section, &
303 36 : "PRINT%EWALD_INFO")
304 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
305 36 : "PRINT%ITER_INFO")
306 :
307 36 : CALL timestop(handle)
308 72 : END SUBROUTINE fist_pol_evaluate_sc
309 :
310 : ! **************************************************************************************************
311 : !> \brief Conjugate-gradient solver for a polarizable force-field
312 : !> \param atomic_kind_set ...
313 : !> \param multipoles ...
314 : !> \param ewald_env ...
315 : !> \param ewald_pw ...
316 : !> \param fist_nonbond_env ...
317 : !> \param cell ...
318 : !> \param particle_set ...
319 : !> \param local_particles ...
320 : !> \param thermo ...
321 : !> \param vg_coulomb ...
322 : !> \param pot_nonbond ...
323 : !> \param f_nonbond ...
324 : !> \param fg_coulomb ...
325 : !> \param use_virial ...
326 : !> \param pv_g ...
327 : !> \param pv_nonbond ...
328 : !> \param mm_section ...
329 : !> \author Toon.Verstraelen@gmail.com (2010-03-01)
330 : !> \note
331 : !> Method: The dipoles are found by minimizing the sum of the electrostatic
332 : !> and the induction energy directly using a conjugate gradient method. This
333 : !> routine assumes that the energy is a quadratic function of the dipoles.
334 : !> Finding the minimum is then done by solving a linear system. This will
335 : !> not work for polarizable force fields that include hyperpolarizability.
336 : !>
337 : !> The implementation of the conjugate gradient solver for linear systems
338 : !> is described in chapter 2.7 Sparse Linear Systems, under the section
339 : !> "Conjugate Gradient Method for a Sparse System". Although the inducible
340 : !> dipoles are the solution of a dense linear system, the same algorithm is
341 : !> still recommended for this situation. One does not have access to the
342 : !> entire hardness kernel to compute the solution with conventional linear
343 : !> algebra routines, but one only has a function that computes the dot
344 : !> product of the hardness kernel and a vector. (This is the routine that
345 : !> computes the electrostatic field at the atoms for a given vector of
346 : !> inducible dipoles.) Given such function, the conjugate gradient method
347 : !> is an efficient way to compute the solution of a linear system, and it
348 : !> scales well towards many degrees of freedom in terms of efficiency and
349 : !> memory usage.
350 : ! **************************************************************************************************
351 68 : SUBROUTINE fist_pol_evaluate_cg(atomic_kind_set, multipoles, ewald_env, ewald_pw, &
352 : fist_nonbond_env, cell, particle_set, local_particles, thermo, vg_coulomb, &
353 68 : pot_nonbond, f_nonbond, fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section)
354 :
355 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
356 : TYPE(multipole_type), POINTER :: multipoles
357 : TYPE(ewald_environment_type), POINTER :: ewald_env
358 : TYPE(ewald_pw_type), POINTER :: ewald_pw
359 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
360 : TYPE(cell_type), POINTER :: cell
361 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
362 : TYPE(distribution_1d_type), POINTER :: local_particles
363 : TYPE(fist_energy_type), POINTER :: thermo
364 : REAL(KIND=dp) :: vg_coulomb, pot_nonbond
365 : REAL(KIND=dp), DIMENSION(:, :) :: f_nonbond, fg_coulomb
366 : LOGICAL, INTENT(IN) :: use_virial
367 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_g, pv_nonbond
368 : TYPE(section_vals_type), POINTER :: mm_section
369 :
370 : CHARACTER(len=*), PARAMETER :: routineN = 'fist_pol_evaluate_cg'
371 :
372 : INTEGER :: ewald_type, handle, i, iatom, ii, ikind, &
373 : iter, iw, iw2, max_ipol_iter, &
374 : natom_of_kind, natoms, nkind, ntot
375 68 : INTEGER, DIMENSION(:), POINTER :: atom_list
376 : LOGICAL :: iwarn
377 : REAL(KIND=dp) :: alpha, apol, beta, denom, eps_pol, &
378 : pot_nonbond_local, rmsd
379 68 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: conjugate, conjugate_applied, efield1, &
380 68 : efield1_ext, residual, tmp_dipoles
381 : TYPE(atomic_kind_type), POINTER :: atomic_kind
382 : TYPE(cp_logger_type), POINTER :: logger
383 :
384 68 : CALL timeset(routineN, handle)
385 68 : NULLIFY (logger, atomic_kind)
386 68 : logger => cp_get_default_logger()
387 :
388 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%ITER_INFO", &
389 68 : extension=".mmLog")
390 : iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
391 68 : extension=".mmLog")
392 :
393 : CALL ewald_env_get(ewald_env, max_ipol_iter=max_ipol_iter, eps_pol=eps_pol, &
394 68 : ewald_type=ewald_type)
395 :
396 : ! allocate work arrays
397 68 : natoms = SIZE(particle_set)
398 204 : ALLOCATE (efield1(3, natoms))
399 204 : ALLOCATE (tmp_dipoles(3, natoms))
400 204 : ALLOCATE (residual(3, natoms))
401 204 : ALLOCATE (conjugate(3, natoms))
402 204 : ALLOCATE (conjugate_applied(3, natoms))
403 204 : ALLOCATE (efield1_ext(3, natoms))
404 :
405 : ! Compute the 'external' electrostatic field (all inducible dipoles
406 : ! equal to zero). this is required for the conjugate gradient solver.
407 : ! We assume that all dipoles are inducible dipoles.
408 31020 : tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles
409 31020 : multipoles%dipoles = 0.0_dp
410 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
411 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
412 : multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
413 : do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
414 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
415 68 : efield1=efield1_ext)
416 31020 : multipoles%dipoles = tmp_dipoles ! restore backup
417 :
418 : ! Compute the electric field with the initial guess of the dipoles.
419 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
420 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
421 : multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
422 : do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
423 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
424 68 : efield1=efield1)
425 :
426 : ! Compute the first residual explicitly.
427 68 : nkind = SIZE(atomic_kind_set)
428 68 : ntot = 0
429 31020 : residual = 0.0_dp
430 232 : DO ikind = 1, nkind
431 164 : atomic_kind => atomic_kind_set(ikind)
432 164 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
433 : ! ignore atoms with polarizability zero
434 164 : IF (apol == 0) CYCLE
435 78 : ntot = ntot + natom_of_kind
436 4466 : DO iatom = 1, natom_of_kind
437 4156 : ii = atom_list(iatom)
438 16788 : DO i = 1, 3
439 : ! residual = b - A x
440 16624 : residual(i, ii) = efield1(i, ii) - multipoles%dipoles(i, ii)/apol
441 : END DO
442 : END DO
443 : END DO
444 : ! The first conjugate residual is equal to the residual.
445 31020 : conjugate(:, :) = residual
446 :
447 68 : IF (iw > 0) THEN
448 34 : WRITE (iw, FMT="(/,T2,A,T63,A)") "POL_SCF| Method", "conjugate gradient"
449 34 : WRITE (iw, FMT="(T2,A,T26,A,T49,A)") "POL_SCF| Iteration", &
450 68 : "Convergence", "Electrostatic & Induction Energy"
451 : END IF
452 308 : pol_scf: DO iter = 1, max_ipol_iter
453 308 : IF (debug_this_module) THEN
454 : ! In principle the residual must not be computed explicitly any more. It
455 : ! is obtained in an indirect way below. When the DEBUG flag is set, the
456 : ! explicit residual is computed and compared with the implicitly derived
457 : ! residual as a double check.
458 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
459 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
460 : multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
461 : do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
462 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
463 0 : efield1=efield1)
464 : ! inapropriate use of denom to check the error on the residual
465 0 : denom = 0.0_dp
466 : END IF
467 308 : rmsd = 0.0_dp
468 : ! Compute the rmsd of the residual.
469 1090 : DO ikind = 1, nkind
470 782 : atomic_kind => atomic_kind_set(ikind)
471 782 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
472 : ! ignore atoms with polarizability zero
473 782 : IF (apol == 0) CYCLE
474 22890 : DO iatom = 1, natom_of_kind
475 21370 : ii = atom_list(iatom)
476 86262 : DO i = 1, 3
477 : ! residual = b - A x
478 64110 : rmsd = rmsd + residual(i, ii)**2
479 85480 : IF (debug_this_module) THEN
480 : denom = denom + (residual(i, ii) - (efield1(i, ii) - &
481 0 : multipoles%dipoles(i, ii)/apol))**2
482 : END IF
483 : END DO
484 : END DO
485 : END DO
486 308 : rmsd = SQRT(rmsd/ntot)
487 308 : IF (iw > 0) THEN
488 154 : WRITE (iw, FMT="(T2,A,T11,I9,T22,E15.6,T67,A)") "POL_SCF|", iter, rmsd, "(not computed)"
489 154 : IF (debug_this_module) THEN
490 0 : denom = SQRT(denom/ntot)
491 0 : WRITE (iw, FMT="(T2,A,T66,E15.6)") "POL_SCF| Error on implicit residual", denom
492 : END IF
493 : END IF
494 :
495 : ! Apply the hardness kernel to the conjugate residual.
496 : ! We assume that all non-zero dipoles are inducible dipoles.
497 153100 : tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles
498 153100 : multipoles%dipoles = conjugate
499 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
500 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
501 : multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., &
502 : do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
503 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
504 308 : efield1=conjugate_applied)
505 153100 : multipoles%dipoles = tmp_dipoles ! restore backup
506 153100 : conjugate_applied(:, :) = efield1_ext - conjugate_applied
507 :
508 : ! Finish conjugate_applied and compute alpha from the conjugate gradient algorithm.
509 308 : alpha = 0.0_dp
510 308 : denom = 0.0_dp
511 1090 : DO ikind = 1, nkind
512 782 : atomic_kind => atomic_kind_set(ikind)
513 782 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
514 : ! ignore atoms with polarizability zero
515 782 : IF (apol == 0) CYCLE
516 22890 : DO iatom = 1, natom_of_kind
517 21370 : ii = atom_list(iatom)
518 85480 : DO i = 1, 3
519 85480 : conjugate_applied(i, ii) = conjugate_applied(i, ii) + conjugate(i, ii)/apol
520 : END DO
521 85480 : alpha = alpha + DOT_PRODUCT(residual(:, ii), residual(:, ii))
522 86262 : denom = denom + DOT_PRODUCT(conjugate(:, ii), conjugate_applied(:, ii))
523 : END DO
524 : END DO
525 308 : alpha = alpha/denom
526 :
527 : ! Compute the new residual and beta from the conjugate gradient method.
528 308 : beta = 0.0_dp
529 308 : denom = 0.0_dp
530 1090 : DO ikind = 1, nkind
531 782 : atomic_kind => atomic_kind_set(ikind)
532 782 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
533 782 : IF (apol == 0) CYCLE
534 22890 : DO iatom = 1, natom_of_kind
535 21370 : ii = atom_list(iatom)
536 85480 : denom = denom + DOT_PRODUCT(residual(:, ii), residual(:, ii))
537 85480 : DO i = 1, 3
538 85480 : residual(i, ii) = residual(i, ii) - alpha*conjugate_applied(i, ii)
539 : END DO
540 86262 : beta = beta + DOT_PRODUCT(residual(:, ii), residual(:, ii))
541 : END DO
542 : END DO
543 308 : beta = beta/denom
544 :
545 : ! Compute the new dipoles, the new conjugate residual, and the induction
546 : ! energy.
547 308 : thermo%e_induction = 0.0_dp
548 1090 : DO ikind = 1, nkind
549 782 : atomic_kind => atomic_kind_set(ikind)
550 782 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
551 : ! ignore atoms with polarizability zero
552 782 : IF (apol == 0) CYCLE
553 22890 : DO iatom = 1, natom_of_kind
554 21370 : ii = atom_list(iatom)
555 86262 : DO i = 1, 3
556 64110 : multipoles%dipoles(i, ii) = multipoles%dipoles(i, ii) + alpha*conjugate(i, ii)
557 64110 : conjugate(i, ii) = residual(i, ii) + beta*conjugate(i, ii)
558 85480 : thermo%e_induction = thermo%e_induction + multipoles%dipoles(i, ii)**2/apol/2.0_dp
559 : END DO
560 : END DO
561 : END DO
562 :
563 : ! Quit if rmsd is low enough
564 308 : IF (rmsd <= eps_pol) THEN
565 68 : IF (iw > 0) WRITE (iw, FMT="(T2,A)") "POL_SCF| Self-consistent polarization converged"
566 : EXIT pol_scf
567 : END IF
568 :
569 : ! Print warning when not converged
570 240 : iwarn = ((rmsd > eps_pol) .AND. (iter >= max_ipol_iter))
571 0 : IF (iwarn) THEN
572 0 : IF (iw > 0) THEN
573 : WRITE (iw, FMT="(T2,A,I0,A,ES9.3)") &
574 0 : "POL_SCF| Self-consistent polarization not converged in ", max_ipol_iter, &
575 0 : " steps to ", eps_pol
576 : END IF
577 0 : CPWARN("Self-consistent Polarization not converged!")
578 : END IF
579 : END DO pol_scf
580 :
581 68 : IF (debug_this_module) THEN
582 : ! Now evaluate after convergence to obtain forces and converged energies
583 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
584 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
585 : multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., &
586 : do_stress=use_virial, do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., &
587 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
588 : forces_local=fg_coulomb, forces_glob=f_nonbond, &
589 0 : pv_local=pv_g, pv_glob=pv_nonbond, efield1=efield1)
590 :
591 : ! Do a final check on the convergence: compute the residual explicitely
592 0 : rmsd = 0.0_dp
593 0 : DO ikind = 1, nkind
594 0 : atomic_kind => atomic_kind_set(ikind)
595 0 : CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
596 : ! ignore atoms with polarizability zero
597 0 : IF (apol == 0) CYCLE
598 0 : DO iatom = 1, natom_of_kind
599 0 : ii = atom_list(iatom)
600 0 : DO i = 1, 3
601 : ! residual = b - A x
602 0 : rmsd = rmsd + (efield1(i, ii) - multipoles%dipoles(i, ii)/apol)**2
603 : END DO
604 : END DO
605 : END DO
606 0 : rmsd = SQRT(rmsd/ntot)
607 0 : IF (iw > 0) WRITE (iw, FMT="(T2,A,T66,E15.6)") "POL_SCF| Final RMSD of residual", rmsd
608 : ! Stop program when congergence is not reached after all
609 0 : IF (rmsd > eps_pol) THEN
610 0 : CPABORT("Error in the conjugate gradient method for self-consistent polarization!")
611 : END IF
612 : ELSE
613 : ! Now evaluate after convergence to obtain forces and converged energies
614 : CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, &
615 : particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, &
616 : multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., &
617 : do_stress=use_virial, do_efield=.FALSE., iw2=iw2, do_debug=.FALSE., &
618 : atomic_kind_set=atomic_kind_set, mm_section=mm_section, &
619 : forces_local=fg_coulomb, forces_glob=f_nonbond, &
620 68 : pv_local=pv_g, pv_glob=pv_nonbond)
621 : END IF
622 68 : pot_nonbond = pot_nonbond + pot_nonbond_local
623 68 : CALL logger%para_env%sum(pot_nonbond_local)
624 :
625 102 : IF (iw > 0) WRITE (iw, FMT="(T2,A,T61,F20.10)") "POL_SCF| Final", &
626 68 : vg_coulomb + pot_nonbond_local + thermo%e_induction
627 :
628 : ! Deallocate working arrays
629 68 : DEALLOCATE (efield1)
630 68 : DEALLOCATE (tmp_dipoles)
631 68 : DEALLOCATE (residual)
632 68 : DEALLOCATE (conjugate)
633 68 : DEALLOCATE (conjugate_applied)
634 68 : DEALLOCATE (efield1_ext)
635 : CALL cp_print_key_finished_output(iw2, logger, mm_section, &
636 68 : "PRINT%EWALD_INFO")
637 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
638 68 : "PRINT%ITER_INFO")
639 :
640 68 : CALL timestop(handle)
641 :
642 136 : END SUBROUTINE fist_pol_evaluate_cg
643 :
644 : ! **************************************************************************************************
645 : !> \brief Main driver for evaluating electrostatic in polarible forcefields
646 : !> All the dependence on the Ewald method should go here!
647 : !> \param ewald_type ...
648 : !> \param ewald_env ...
649 : !> \param ewald_pw ...
650 : !> \param fist_nonbond_env ...
651 : !> \param cell ...
652 : !> \param particle_set ...
653 : !> \param local_particles ...
654 : !> \param vg_coulomb ...
655 : !> \param pot_nonbond ...
656 : !> \param thermo ...
657 : !> \param multipoles ...
658 : !> \param do_correction_bonded ...
659 : !> \param do_forces ...
660 : !> \param do_stress ...
661 : !> \param do_efield ...
662 : !> \param iw2 ...
663 : !> \param do_debug ...
664 : !> \param atomic_kind_set ...
665 : !> \param mm_section ...
666 : !> \param efield0 ...
667 : !> \param efield1 ...
668 : !> \param efield2 ...
669 : !> \param forces_local ...
670 : !> \param forces_glob ...
671 : !> \param pv_local ...
672 : !> \param pv_glob ...
673 : !> \author Teodoro Laino [tlaino] 05.2009
674 : ! **************************************************************************************************
675 1684 : SUBROUTINE eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, &
676 : cell, particle_set, local_particles, vg_coulomb, pot_nonbond, thermo, &
677 : multipoles, do_correction_bonded, do_forces, do_stress, do_efield, iw2, &
678 1684 : do_debug, atomic_kind_set, mm_section, efield0, efield1, efield2, forces_local, &
679 842 : forces_glob, pv_local, pv_glob)
680 :
681 : INTEGER, INTENT(IN) :: ewald_type
682 : TYPE(ewald_environment_type), POINTER :: ewald_env
683 : TYPE(ewald_pw_type), POINTER :: ewald_pw
684 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
685 : TYPE(cell_type), POINTER :: cell
686 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
687 : TYPE(distribution_1d_type), POINTER :: local_particles
688 : REAL(KIND=dp), INTENT(OUT) :: vg_coulomb, pot_nonbond
689 : TYPE(fist_energy_type), POINTER :: thermo
690 : TYPE(multipole_type), POINTER :: multipoles
691 : LOGICAL, INTENT(IN) :: do_correction_bonded, do_forces, &
692 : do_stress, do_efield
693 : INTEGER, INTENT(IN) :: iw2
694 : LOGICAL, INTENT(IN) :: do_debug
695 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
696 : TYPE(section_vals_type), POINTER :: mm_section
697 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: efield0
698 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: efield1, efield2
699 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
700 : OPTIONAL :: forces_local, forces_glob, pv_local, &
701 : pv_glob
702 :
703 : CHARACTER(len=*), PARAMETER :: routineN = 'eval_pol_ewald'
704 :
705 : INTEGER :: handle
706 :
707 842 : CALL timeset(routineN, handle)
708 :
709 842 : pot_nonbond = 0.0_dp ! Initialization..
710 842 : vg_coulomb = 0.0_dp ! Initialization..
711 1684 : SELECT CASE (ewald_type)
712 : CASE (do_ewald_ewald)
713 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, fist_nonbond_env, cell, &
714 : particle_set, local_particles, vg_coulomb, pot_nonbond, thermo%e_neut, &
715 : thermo%e_self, multipoles%task, do_correction_bonded=do_correction_bonded, &
716 : do_forces=do_forces, do_stress=do_stress, do_efield=do_efield, &
717 : radii=multipoles%radii, charges=multipoles%charges, &
718 : dipoles=multipoles%dipoles, quadrupoles=multipoles%quadrupoles, &
719 : forces_local=forces_local, forces_glob=forces_glob, pv_local=pv_local, &
720 : pv_glob=pv_glob, iw=iw2, do_debug=do_debug, atomic_kind_set=atomic_kind_set, &
721 5288 : mm_section=mm_section, efield0=efield0, efield1=efield1, efield2=efield2)
722 : CASE (do_ewald_pme)
723 0 : CPABORT("Multipole Ewald not yet implemented within a PME scheme!")
724 : CASE (do_ewald_spme)
725 842 : CPABORT("Multipole Ewald not yet implemented within a SPME scheme!")
726 : END SELECT
727 842 : CALL timestop(handle)
728 842 : END SUBROUTINE eval_pol_ewald
729 :
730 : END MODULE fist_pol_scf
|
