(git:3add494)
fist_pol_scf.F
Go to the documentation of this file.
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  SUBROUTINE fist_pol_evaluate(atomic_kind_set, multipoles, ewald_env, &
74  ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, &
75  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  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  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  pv_g, pv_nonbond, mm_section)
105  END SELECT
106 
107  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  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  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  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  CALL timeset(routinen, handle)
170  NULLIFY (logger, atomic_kind)
171  logger => cp_get_default_logger()
172 
173  iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%ITER_INFO", &
174  extension=".mmLog")
175  iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
176  extension=".mmLog")
177 
178  CALL ewald_env_get(ewald_env, max_ipol_iter=max_ipol_iter, eps_pol=eps_pol, &
179  ewald_type=ewald_type)
180 
181  natoms = SIZE(particle_set)
182  ALLOCATE (efield1(3, natoms))
183  ALLOCATE (efield2(9, natoms))
184 
185  nkind = SIZE(atomic_kind_set)
186  IF (iw > 0) THEN
187  WRITE (iw, fmt='(/,T2,"POL_SCF|" ,"Method: self-consistent")')
188  WRITE (iw, fmt='(T2,"POL_SCF| "," Iteration",7X,"Conv.",T49,"Electrostatic & Induction Energy")')
189  END IF
190  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  efield1=efield1, efield2=efield2)
198  CALL logger%para_env%sum(pot_nonbond_local)
199 
200  ! compute the new dipoles, qudrupoles, and check for convergence
201  ntot = 0
202  rmsd = 0.0_dp
203  thermo%e_induction = 0.0_dp
204  DO ikind = 1, nkind
205  atomic_kind => atomic_kind_set(ikind)
206  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  IF (apol == 0 .AND. cpol == 0) cycle
209  ! increment counter correctly
210  IF (apol /= 0) ntot = ntot + natom_of_kind
211  IF (cpol /= 0) ntot = ntot + natom_of_kind
212 
213  DO iatom = 1, natom_of_kind
214  ii = atom_list(iatom)
215  IF (apol /= 0) THEN
216  DO i = 1, 3
217  ! the rmsd of the derivatives of the energy towards the
218  ! components of the atomic dipole moments
219  rmsd = rmsd + (multipoles%dipoles(i, ii)/apol - efield1(i, ii))**2
220  END DO
221  END IF
222  IF (cpol /= 0) THEN
223  rmsd = rmsd + (multipoles%quadrupoles(1, 1, ii)/cpol - efield2(1, ii))**2
224  rmsd = rmsd + (multipoles%quadrupoles(2, 1, ii)/cpol - efield2(2, ii))**2
225  rmsd = rmsd + (multipoles%quadrupoles(3, 1, ii)/cpol - efield2(3, ii))**2
226  rmsd = rmsd + (multipoles%quadrupoles(1, 2, ii)/cpol - efield2(4, ii))**2
227  rmsd = rmsd + (multipoles%quadrupoles(2, 2, ii)/cpol - efield2(5, ii))**2
228  rmsd = rmsd + (multipoles%quadrupoles(3, 2, ii)/cpol - efield2(6, ii))**2
229  rmsd = rmsd + (multipoles%quadrupoles(1, 3, ii)/cpol - efield2(7, ii))**2
230  rmsd = rmsd + (multipoles%quadrupoles(2, 3, ii)/cpol - efield2(8, ii))**2
231  rmsd = rmsd + (multipoles%quadrupoles(3, 3, ii)/cpol - efield2(9, ii))**2
232  END IF
233 ! compute dipole
234  multipoles%dipoles(:, ii) = apol*efield1(:, ii)
235 ! compute quadrupole
236  IF (cpol /= 0) THEN
237  multipoles%quadrupoles(1, 1, ii) = cpol*efield2(1, ii)
238  multipoles%quadrupoles(2, 1, ii) = cpol*efield2(2, ii)
239  multipoles%quadrupoles(3, 1, ii) = cpol*efield2(3, ii)
240  multipoles%quadrupoles(1, 2, ii) = cpol*efield2(4, ii)
241  multipoles%quadrupoles(2, 2, ii) = cpol*efield2(5, ii)
242  multipoles%quadrupoles(3, 2, ii) = cpol*efield2(6, ii)
243  multipoles%quadrupoles(1, 3, ii) = cpol*efield2(7, ii)
244  multipoles%quadrupoles(2, 3, ii) = cpol*efield2(8, ii)
245  multipoles%quadrupoles(3, 3, ii) = cpol*efield2(9, ii)
246  END IF
247  ! Compute the new induction term while we are here
248  IF (apol /= 0) THEN
249  thermo%e_induction = thermo%e_induction + &
250  dot_product(multipoles%dipoles(:, ii), &
251  multipoles%dipoles(:, ii))/apol/2.0_dp
252  END IF
253  IF (cpol /= 0) THEN
254  tmp_trace = 0._dp
255  DO i = 1, 3
256  DO j = 1, 3
257  tmp_trace = tmp_trace + &
258  multipoles%quadrupoles(i, j, ii)*multipoles%quadrupoles(i, j, ii)
259  END DO
260  END DO
261  thermo%e_induction = thermo%e_induction + tmp_trace/cpol/6.0_dp
262  END IF
263  END DO
264  END DO
265  rmsd = sqrt(rmsd/real(ntot, kind=dp))
266  IF (iw > 0) THEN
267  ! print the energy that is minimized (this is electrostatic + induction)
268  WRITE (iw, fmt='(T5,"POL_SCF|",5X,I5,5X,E12.6,T61,F20.10)') iter, &
269  rmsd, vg_coulomb + pot_nonbond_local + thermo%e_induction
270  END IF
271  IF (rmsd <= eps_pol) THEN
272  IF (iw > 0) WRITE (iw, fmt='(T5,"POL_SCF|",1X,"Self-consistent Polarization achieved.")')
273  EXIT pol_scf
274  END IF
275 
276  iwarn = ((rmsd > eps_pol) .AND. (iter == max_ipol_iter))
277  IF (iwarn .AND. iw > 0) WRITE (iw, fmt='(T5,"POL_SCF|",1X,"Self-consistent Polarization not converged!")')
278  IF (iwarn) &
279  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  pv_local=pv_g, pv_glob=pv_nonbond)
290  pot_nonbond = pot_nonbond + pot_nonbond_local
291  CALL logger%para_env%sum(pot_nonbond_local)
292 
293  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  vg_coulomb + pot_nonbond_local + thermo%e_induction
297  END IF
298 
299  ! Deallocate working arrays
300  DEALLOCATE (efield1)
301  DEALLOCATE (efield2)
302  CALL cp_print_key_finished_output(iw2, logger, mm_section, &
303  "PRINT%EWALD_INFO")
304  CALL cp_print_key_finished_output(iw, logger, mm_section, &
305  "PRINT%ITER_INFO")
306 
307  CALL timestop(handle)
308  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  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  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  INTEGER, DIMENSION(:), POINTER :: atom_list
376  LOGICAL :: iwarn
377  REAL(kind=dp) :: alpha, apol, beta, denom, eps_pol, &
378  pot_nonbond_local, rmsd
379  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: conjugate, conjugate_applied, efield1, &
380  efield1_ext, residual, tmp_dipoles
381  TYPE(atomic_kind_type), POINTER :: atomic_kind
382  TYPE(cp_logger_type), POINTER :: logger
383 
384  CALL timeset(routinen, handle)
385  NULLIFY (logger, atomic_kind)
386  logger => cp_get_default_logger()
387 
388  iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%ITER_INFO", &
389  extension=".mmLog")
390  iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
391  extension=".mmLog")
392 
393  CALL ewald_env_get(ewald_env, max_ipol_iter=max_ipol_iter, eps_pol=eps_pol, &
394  ewald_type=ewald_type)
395 
396  ! allocate work arrays
397  natoms = SIZE(particle_set)
398  ALLOCATE (efield1(3, natoms))
399  ALLOCATE (tmp_dipoles(3, natoms))
400  ALLOCATE (residual(3, natoms))
401  ALLOCATE (conjugate(3, natoms))
402  ALLOCATE (conjugate_applied(3, natoms))
403  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  tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles
409  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  efield1=efield1_ext)
416  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  efield1=efield1)
425 
426  ! Compute the first residual explicitly.
427  nkind = SIZE(atomic_kind_set)
428  ntot = 0
429  residual = 0.0_dp
430  DO ikind = 1, nkind
431  atomic_kind => atomic_kind_set(ikind)
432  CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
433  ! ignore atoms with polarizability zero
434  IF (apol == 0) cycle
435  ntot = ntot + natom_of_kind
436  DO iatom = 1, natom_of_kind
437  ii = atom_list(iatom)
438  DO i = 1, 3
439  ! residual = b - A x
440  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  conjugate(:, :) = residual
446 
447  IF (iw > 0) THEN
448  WRITE (iw, fmt="(/,T2,A,T63,A)") "POL_SCF| Method", "conjugate gradient"
449  WRITE (iw, fmt="(T2,A,T26,A,T49,A)") "POL_SCF| Iteration", &
450  "Convergence", "Electrostatic & Induction Energy"
451  END IF
452  pol_scf: DO iter = 1, max_ipol_iter
453  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  efield1=efield1)
464  ! inapropriate use of denom to check the error on the residual
465  denom = 0.0_dp
466  END IF
467  rmsd = 0.0_dp
468  ! Compute the rmsd of the residual.
469  DO ikind = 1, nkind
470  atomic_kind => atomic_kind_set(ikind)
471  CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
472  ! ignore atoms with polarizability zero
473  IF (apol == 0) cycle
474  DO iatom = 1, natom_of_kind
475  ii = atom_list(iatom)
476  DO i = 1, 3
477  ! residual = b - A x
478  rmsd = rmsd + residual(i, ii)**2
479  IF (debug_this_module) THEN
480  denom = denom + (residual(i, ii) - (efield1(i, ii) - &
481  multipoles%dipoles(i, ii)/apol))**2
482  END IF
483  END DO
484  END DO
485  END DO
486  rmsd = sqrt(rmsd/ntot)
487  IF (iw > 0) THEN
488  WRITE (iw, fmt="(T2,A,T11,I9,T22,E15.6,T67,A)") "POL_SCF|", iter, rmsd, "(not computed)"
489  IF (debug_this_module) THEN
490  denom = sqrt(denom/ntot)
491  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  tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles
498  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  efield1=conjugate_applied)
505  multipoles%dipoles = tmp_dipoles ! restore backup
506  conjugate_applied(:, :) = efield1_ext - conjugate_applied
507 
508  ! Finish conjugate_applied and compute alpha from the conjugate gradient algorithm.
509  alpha = 0.0_dp
510  denom = 0.0_dp
511  DO ikind = 1, nkind
512  atomic_kind => atomic_kind_set(ikind)
513  CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
514  ! ignore atoms with polarizability zero
515  IF (apol == 0) cycle
516  DO iatom = 1, natom_of_kind
517  ii = atom_list(iatom)
518  DO i = 1, 3
519  conjugate_applied(i, ii) = conjugate_applied(i, ii) + conjugate(i, ii)/apol
520  END DO
521  alpha = alpha + dot_product(residual(:, ii), residual(:, ii))
522  denom = denom + dot_product(conjugate(:, ii), conjugate_applied(:, ii))
523  END DO
524  END DO
525  alpha = alpha/denom
526 
527  ! Compute the new residual and beta from the conjugate gradient method.
528  beta = 0.0_dp
529  denom = 0.0_dp
530  DO ikind = 1, nkind
531  atomic_kind => atomic_kind_set(ikind)
532  CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
533  IF (apol == 0) cycle
534  DO iatom = 1, natom_of_kind
535  ii = atom_list(iatom)
536  denom = denom + dot_product(residual(:, ii), residual(:, ii))
537  DO i = 1, 3
538  residual(i, ii) = residual(i, ii) - alpha*conjugate_applied(i, ii)
539  END DO
540  beta = beta + dot_product(residual(:, ii), residual(:, ii))
541  END DO
542  END DO
543  beta = beta/denom
544 
545  ! Compute the new dipoles, the new conjugate residual, and the induction
546  ! energy.
547  thermo%e_induction = 0.0_dp
548  DO ikind = 1, nkind
549  atomic_kind => atomic_kind_set(ikind)
550  CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
551  ! ignore atoms with polarizability zero
552  IF (apol == 0) cycle
553  DO iatom = 1, natom_of_kind
554  ii = atom_list(iatom)
555  DO i = 1, 3
556  multipoles%dipoles(i, ii) = multipoles%dipoles(i, ii) + alpha*conjugate(i, ii)
557  conjugate(i, ii) = residual(i, ii) + beta*conjugate(i, ii)
558  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  IF (rmsd <= eps_pol) THEN
565  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  iwarn = ((rmsd > eps_pol) .AND. (iter >= max_ipol_iter))
571  IF (iwarn) THEN
572  IF (iw > 0) THEN
573  WRITE (iw, fmt="(T2,A,I0,A,ES9.3)") &
574  "POL_SCF| Self-consistent polarization not converged in ", max_ipol_iter, &
575  " steps to ", eps_pol
576  END IF
577  cpwarn("Self-consistent Polarization not converged!")
578  END IF
579  END DO pol_scf
580 
581  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  pv_local=pv_g, pv_glob=pv_nonbond, efield1=efield1)
590 
591  ! Do a final check on the convergence: compute the residual explicitely
592  rmsd = 0.0_dp
593  DO ikind = 1, nkind
594  atomic_kind => atomic_kind_set(ikind)
595  CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
596  ! ignore atoms with polarizability zero
597  IF (apol == 0) cycle
598  DO iatom = 1, natom_of_kind
599  ii = atom_list(iatom)
600  DO i = 1, 3
601  ! residual = b - A x
602  rmsd = rmsd + (efield1(i, ii) - multipoles%dipoles(i, ii)/apol)**2
603  END DO
604  END DO
605  END DO
606  rmsd = sqrt(rmsd/ntot)
607  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  IF (rmsd > eps_pol) THEN
610  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  pv_local=pv_g, pv_glob=pv_nonbond)
621  END IF
622  pot_nonbond = pot_nonbond + pot_nonbond_local
623  CALL logger%para_env%sum(pot_nonbond_local)
624 
625  IF (iw > 0) WRITE (iw, fmt="(T2,A,T61,F20.10)") "POL_SCF| Final", &
626  vg_coulomb + pot_nonbond_local + thermo%e_induction
627 
628  ! Deallocate working arrays
629  DEALLOCATE (efield1)
630  DEALLOCATE (tmp_dipoles)
631  DEALLOCATE (residual)
632  DEALLOCATE (conjugate)
633  DEALLOCATE (conjugate_applied)
634  DEALLOCATE (efield1_ext)
635  CALL cp_print_key_finished_output(iw2, logger, mm_section, &
636  "PRINT%EWALD_INFO")
637  CALL cp_print_key_finished_output(iw, logger, mm_section, &
638  "PRINT%ITER_INFO")
639 
640  CALL timestop(handle)
641 
642  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  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  do_debug, atomic_kind_set, mm_section, efield0, efield1, efield2, forces_local, &
679  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  CALL timeset(routinen, handle)
708 
709  pot_nonbond = 0.0_dp ! Initialization..
710  vg_coulomb = 0.0_dp ! Initialization..
711  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  mm_section=mm_section, efield0=efield0, efield1=efield1, efield2=efield2)
722  CASE (do_ewald_pme)
723  cpabort("Multipole Ewald not yet implemented within a PME scheme!")
724  CASE (do_ewald_spme)
725  cpabort("Multipole Ewald not yet implemented within a SPME scheme!")
726  END SELECT
727  CALL timestop(handle)
728  END SUBROUTINE eval_pol_ewald
729 
730 END MODULE fist_pol_scf
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition: distribution_1d_types.F:24
ewald_environment_types
Definition: ewald_environment_types.F:14
ewald_environment_types::ewald_env_get
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Definition: ewald_environment_types.F:123
ewald_pw_types
pw_types
Definition: ewald_pw_types.F:12
ewalds_multipole
Treats the electrostatic for multipoles (up to quadrupoles)
Definition: ewalds_multipole.F:14
ewalds_multipole::ewald_multipole_evaluate
recursive subroutine, public ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task, do_correction_bonded, do_forces, do_stress, do_efield, radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, efield0, efield1, efield2, iw, do_debug, atomic_kind_set, mm_section)
Computes the potential and the force for a lattice sum of multipoles (up to quadrupole)
Definition: ewalds_multipole.F:122
fist_energy_types
Definition: fist_energy_types.F:13
fist_nonbond_env_types
Definition: fist_nonbond_env_types.F:13
fist_pol_scf
Definition: fist_pol_scf.F:16
fist_pol_scf::fist_pol_evaluate
subroutine, public fist_pol_evaluate(atomic_kind_set, multipoles, ewald_env, ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, thermo, vg_coulomb, pot_nonbond, f_nonbond, fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section, do_ipol)
Main driver for evaluating energy/forces in a polarizable forcefield.
Definition: fist_pol_scf.F:77
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_fist_pol_cg
integer, parameter, public do_fist_pol_cg
Definition: input_constants.F:941
input_constants::do_fist_pol_sc
integer, parameter, public do_fist_pol_sc
Definition: input_constants.F:941
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
multipole_types
Multipole structure: for multipole (fixed and induced) in FF based MD.
Definition: multipole_types.F:12
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
pw_poisson_types
functions related to the poisson solver on regular grids
Definition: pw_poisson_types.F:22
pw_poisson_types::do_ewald_pme
integer, parameter, public do_ewald_pme
Definition: pw_poisson_types.F:75
pw_poisson_types::do_ewald_ewald
integer, parameter, public do_ewald_ewald
Definition: pw_poisson_types.F:75
pw_poisson_types::do_ewald_spme
integer, parameter, public do_ewald_spme
Definition: pw_poisson_types.F:75
rmsd
Defines functions to perform rmsd in 3D.
Definition: rmsd.F:12