(git:374b731)
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-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
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 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
730END 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
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:55
ewald_pw_types::ewald_pw_type
Definition ewald_pw_types.F:64
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:126
multipole_types::multipole_type
Define multipole type.
Definition multipole_types.F:44
particle_types::particle_type
Definition particle_types.F:35