(git:ed6f26b)
Loading...
Searching...
No Matches
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-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
16MODULE 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
49CONTAINS
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 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 pv_local=pv_g, pv_glob=pv_nonbond)
289 pot_nonbond = pot_nonbond + pot_nonbond_local
290 CALL logger%para_env%sum(pot_nonbond_local)
291
292 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 vg_coulomb + pot_nonbond_local + thermo%e_induction
296 END IF
297
298 ! Deallocate working arrays
299 DEALLOCATE (efield1)
300 DEALLOCATE (efield2)
301 CALL cp_print_key_finished_output(iw2, logger, mm_section, &
302 "PRINT%EWALD_INFO")
303 CALL cp_print_key_finished_output(iw, logger, mm_section, &
304 "PRINT%ITER_INFO")
305
306 CALL timestop(handle)
307 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 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 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 INTEGER, DIMENSION(:), POINTER :: atom_list
375 LOGICAL :: iwarn
376 REAL(kind=dp) :: alpha, apol, beta, denom, eps_pol, &
377 pot_nonbond_local, rmsd
378 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: conjugate, conjugate_applied, efield1, &
379 efield1_ext, residual, tmp_dipoles
380 TYPE(atomic_kind_type), POINTER :: atomic_kind
381 TYPE(cp_logger_type), POINTER :: logger
382
383 CALL timeset(routinen, handle)
384 NULLIFY (logger, atomic_kind)
385 logger => cp_get_default_logger()
386
387 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%ITER_INFO", &
388 extension=".mmLog")
389 iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
390 extension=".mmLog")
391
392 CALL ewald_env_get(ewald_env, max_ipol_iter=max_ipol_iter, eps_pol=eps_pol, &
393 ewald_type=ewald_type)
394
395 ! allocate work arrays
396 natoms = SIZE(particle_set)
397 ALLOCATE (efield1(3, natoms))
398 ALLOCATE (tmp_dipoles(3, natoms))
399 ALLOCATE (residual(3, natoms))
400 ALLOCATE (conjugate(3, natoms))
401 ALLOCATE (conjugate_applied(3, natoms))
402 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 tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles
408 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 efield1=efield1_ext)
415 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 efield1=efield1)
424
425 ! Compute the first residual explicitly.
426 nkind = SIZE(atomic_kind_set)
427 ntot = 0
428 residual = 0.0_dp
429 DO ikind = 1, nkind
430 atomic_kind => atomic_kind_set(ikind)
431 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
432 ! ignore atoms with polarizability zero
433 IF (apol == 0) cycle
434 ntot = ntot + natom_of_kind
435 DO iatom = 1, natom_of_kind
436 ii = atom_list(iatom)
437 DO i = 1, 3
438 ! residual = b - A x
439 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 conjugate(:, :) = residual
445
446 IF (iw > 0) THEN
447 WRITE (iw, fmt="(/,T2,A,T63,A)") "POL_SCF| Method", "conjugate gradient"
448 WRITE (iw, fmt="(T2,A,T26,A,T49,A)") "POL_SCF| Iteration", &
449 "Convergence", "Electrostatic & Induction Energy"
450 END IF
451 pol_scf: DO iter = 1, max_ipol_iter
452 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 efield1=efield1)
463 ! inapropriate use of denom to check the error on the residual
464 denom = 0.0_dp
465 END IF
466 rmsd = 0.0_dp
467 ! Compute the rmsd of the residual.
468 DO ikind = 1, nkind
469 atomic_kind => atomic_kind_set(ikind)
470 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
471 ! ignore atoms with polarizability zero
472 IF (apol == 0) cycle
473 DO iatom = 1, natom_of_kind
474 ii = atom_list(iatom)
475 DO i = 1, 3
476 ! residual = b - A x
477 rmsd = rmsd + residual(i, ii)**2
478 IF (debug_this_module) THEN
479 denom = denom + (residual(i, ii) - (efield1(i, ii) - &
480 multipoles%dipoles(i, ii)/apol))**2
481 END IF
482 END DO
483 END DO
484 END DO
485 rmsd = sqrt(rmsd/ntot)
486 IF (iw > 0) THEN
487 WRITE (iw, fmt="(T2,A,T11,I9,T22,E15.6,T67,A)") "POL_SCF|", iter, rmsd, "(not computed)"
488 IF (debug_this_module) THEN
489 denom = sqrt(denom/ntot)
490 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 tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles
497 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 efield1=conjugate_applied)
504 multipoles%dipoles = tmp_dipoles ! restore backup
505 conjugate_applied(:, :) = efield1_ext - conjugate_applied
506
507 ! Finish conjugate_applied and compute alpha from the conjugate gradient algorithm.
508 alpha = 0.0_dp
509 denom = 0.0_dp
510 DO ikind = 1, nkind
511 atomic_kind => atomic_kind_set(ikind)
512 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
513 ! ignore atoms with polarizability zero
514 IF (apol == 0) cycle
515 DO iatom = 1, natom_of_kind
516 ii = atom_list(iatom)
517 DO i = 1, 3
518 conjugate_applied(i, ii) = conjugate_applied(i, ii) + conjugate(i, ii)/apol
519 END DO
520 alpha = alpha + dot_product(residual(:, ii), residual(:, ii))
521 denom = denom + dot_product(conjugate(:, ii), conjugate_applied(:, ii))
522 END DO
523 END DO
524 alpha = alpha/denom
525
526 ! Compute the new residual and beta from the conjugate gradient method.
527 beta = 0.0_dp
528 denom = 0.0_dp
529 DO ikind = 1, nkind
530 atomic_kind => atomic_kind_set(ikind)
531 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
532 IF (apol == 0) cycle
533 DO iatom = 1, natom_of_kind
534 ii = atom_list(iatom)
535 denom = denom + dot_product(residual(:, ii), residual(:, ii))
536 DO i = 1, 3
537 residual(i, ii) = residual(i, ii) - alpha*conjugate_applied(i, ii)
538 END DO
539 beta = beta + dot_product(residual(:, ii), residual(:, ii))
540 END DO
541 END DO
542 beta = beta/denom
543
544 ! Compute the new dipoles, the new conjugate residual, and the induction
545 ! energy.
546 thermo%e_induction = 0.0_dp
547 DO ikind = 1, nkind
548 atomic_kind => atomic_kind_set(ikind)
549 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
550 ! ignore atoms with polarizability zero
551 IF (apol == 0) cycle
552 DO iatom = 1, natom_of_kind
553 ii = atom_list(iatom)
554 DO i = 1, 3
555 multipoles%dipoles(i, ii) = multipoles%dipoles(i, ii) + alpha*conjugate(i, ii)
556 conjugate(i, ii) = residual(i, ii) + beta*conjugate(i, ii)
557 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 IF (rmsd <= eps_pol) THEN
564 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 iwarn = ((rmsd > eps_pol) .AND. (iter >= max_ipol_iter))
570 IF (iwarn) THEN
571 IF (iw > 0) THEN
572 WRITE (iw, fmt="(T2,A,I0,A,ES9.3)") &
573 "POL_SCF| Self-consistent polarization not converged in ", max_ipol_iter, &
574 " steps to ", eps_pol
575 END IF
576 cpwarn("Self-consistent Polarization not converged!")
577 END IF
578 END DO pol_scf
579
580 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 pv_local=pv_g, pv_glob=pv_nonbond, efield1=efield1)
589
590 ! Do a final check on the convergence: compute the residual explicitely
591 rmsd = 0.0_dp
592 DO ikind = 1, nkind
593 atomic_kind => atomic_kind_set(ikind)
594 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list)
595 ! ignore atoms with polarizability zero
596 IF (apol == 0) cycle
597 DO iatom = 1, natom_of_kind
598 ii = atom_list(iatom)
599 DO i = 1, 3
600 ! residual = b - A x
601 rmsd = rmsd + (efield1(i, ii) - multipoles%dipoles(i, ii)/apol)**2
602 END DO
603 END DO
604 END DO
605 rmsd = sqrt(rmsd/ntot)
606 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 IF (rmsd > eps_pol) THEN
609 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 pv_local=pv_g, pv_glob=pv_nonbond)
620 END IF
621 pot_nonbond = pot_nonbond + pot_nonbond_local
622 CALL logger%para_env%sum(pot_nonbond_local)
623
624 IF (iw > 0) WRITE (iw, fmt="(T2,A,T61,F20.10)") "POL_SCF| Final", &
625 vg_coulomb + pot_nonbond_local + thermo%e_induction
626
627 ! Deallocate working arrays
628 DEALLOCATE (efield1)
629 DEALLOCATE (tmp_dipoles)
630 DEALLOCATE (residual)
631 DEALLOCATE (conjugate)
632 DEALLOCATE (conjugate_applied)
633 DEALLOCATE (efield1_ext)
634 CALL cp_print_key_finished_output(iw2, logger, mm_section, &
635 "PRINT%EWALD_INFO")
636 CALL cp_print_key_finished_output(iw, logger, mm_section, &
637 "PRINT%ITER_INFO")
638
639 CALL timestop(handle)
640
641 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 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 do_debug, atomic_kind_set, mm_section, efield0, efield1, efield2, forces_local, &
678 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 CALL timeset(routinen, handle)
707
708 pot_nonbond = 0.0_dp ! Initialization..
709 vg_coulomb = 0.0_dp ! Initialization..
710 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 mm_section=mm_section, efield0=efield0, efield1=efield1, efield2=efield2)
721 CASE (do_ewald_pme)
722 cpabort("Multipole Ewald not yet implemented within a PME scheme!")
723 CASE (do_ewald_spme)
724 cpabort("Multipole Ewald not yet implemented within a SPME scheme!")
725 END SELECT
726 CALL timestop(handle)
727 END SUBROUTINE eval_pol_ewald
728
729END 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:853
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:1063
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:124
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:954
input_constants::do_fist_pol_sc
integer, parameter, public do_fist_pol_sc
Definition input_constants.F:954
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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:56
ewald_pw_types::ewald_pw_type
Definition ewald_pw_types.F:63
fist_energy_types::fist_energy_type
Definition fist_energy_types.F:24
fist_nonbond_env_types::fist_nonbond_env_type
Definition fist_nonbond_env_types.F:81
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
multipole_types::multipole_type
Define multipole type.
Definition multipole_types.F:44
particle_types::particle_type
Definition particle_types.F:35