(git:e7e05ae)
qs_sccs.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 !> \brief Self-consistent continuum solvation (SCCS) model implementation
10 !> \author Matthias Krack (MK)
11 !> \version 1.0
12 !> \par Literature:
13 !> - J.-L. Fattebert and F. Gygi,
14 !> Density functional theory for efficient ab initio molecular dynamics
15 !> simulations in solution, J. Comput. Chem. 23, 662-666 (2002)
16 !> - O. Andreussi, I. Dabo, and N. Marzari,
17 !> Revised self-consistent continuum solvation in electronic-structure
18 !> calculations, J. Chem. Phys. 136, 064102-20 (2012)
19 !> \par History:
20 !> - Creation (10.10.2013,MK)
21 !> - Derivatives using finite differences (26.11.2013,MK)
22 !> - Cube file dump of the dielectric function (19.12.2013,MK)
23 !> - Cube file dump of the polarisation potential (20.12.2013,MK)
24 !> - Calculation of volume and surface of the cavity (21.12.2013,MK)
25 !> - Functional derivative of the cavitation energy (28.12.2013,MK)
26 !> - Update printout (11.11.2022,MK)
27 ! **************************************************************************************************
28 
29 MODULE qs_sccs
30 
31  USE cp_control_types, ONLY: dft_control_type,&
32  sccs_control_type
33  USE cp_log_handling, ONLY: cp_get_default_logger,&
34  cp_logger_type
35  USE cp_output_handling, ONLY: cp_p_file,&
36  cp_print_key_finished_output,&
37  cp_print_key_should_output,&
38  cp_print_key_unit_nr,&
39  low_print_level
40  USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
41  USE cp_realspace_grid_init, ONLY: init_input_type
42  USE cp_subsys_types, ONLY: cp_subsys_get,&
43  cp_subsys_type
44  USE cp_units, ONLY: cp_unit_from_cp2k
45  USE input_constants, ONLY: sccs_andreussi,&
46  sccs_derivative_cd3,&
47  sccs_derivative_cd5,&
48  sccs_derivative_cd7,&
49  sccs_derivative_fft,&
50  sccs_fattebert_gygi
51  USE input_section_types, ONLY: section_get_ivals,&
52  section_get_lval,&
53  section_vals_get_subs_vals,&
54  section_vals_type
55  USE kinds, ONLY: default_path_length,&
56  default_string_length,&
57  dp,&
58  int_8
59  USE mathconstants, ONLY: fourpi,&
60  twopi
61  USE message_passing, ONLY: mp_para_env_type
62  USE particle_list_types, ONLY: particle_list_type
63  USE pw_env_types, ONLY: pw_env_get,&
64  pw_env_type
65  USE pw_methods, ONLY: pw_axpy,&
66  pw_copy,&
67  pw_derive,&
68  pw_integral_ab,&
69  pw_integrate_function,&
70  pw_transfer,&
71  pw_zero
72  USE pw_poisson_methods, ONLY: pw_poisson_solve
73  USE pw_poisson_types, ONLY: pw_poisson_analytic,&
74  pw_poisson_mt,&
75  pw_poisson_type
76  USE pw_pool_types, ONLY: pw_pool_p_type,&
77  pw_pool_type
78  USE pw_types, ONLY: pw_c1d_gs_type,&
79  pw_r3d_rs_type
80  USE qs_energy_types, ONLY: qs_energy_type
81  USE qs_environment_types, ONLY: get_qs_env,&
82  qs_environment_type
83  USE qs_rho_types, ONLY: qs_rho_get,&
84  qs_rho_type
85  USE qs_scf_types, ONLY: qs_scf_env_type
86  USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
87  realspace_grid_input_type,&
88  realspace_grid_type,&
89  rs_grid_create,&
90  rs_grid_create_descriptor,&
91  rs_grid_release,&
92  rs_grid_release_descriptor
93  USE rs_methods, ONLY: derive_fdm_cd3,&
94  derive_fdm_cd5,&
95  derive_fdm_cd7
96 #include "./base/base_uses.f90"
97 
98  IMPLICIT NONE
99 
100  PRIVATE
101 
102  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_sccs'
103 
104  PUBLIC :: print_sccs_results, sccs
105 
106 CONTAINS
107 
108 ! **************************************************************************************************
109 !> \brief Self-consistent continuum solvation (SCCS) model implementation
110 !> \param qs_env ...
111 !> \param rho_tot_gspace ...
112 !> \param v_hartree_gspace ...
113 !> \param v_sccs ...
114 !> \param h_stress ...
115 !> \par History:
116 !> - Creation (10.10.2013,MK)
117 !> \author Matthias Krack (MK)
118 !> \version 1.0
119 ! **************************************************************************************************
120 
121  SUBROUTINE sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs, h_stress)
122 
123  TYPE(qs_environment_type), POINTER :: qs_env
124  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_tot_gspace, v_hartree_gspace
125  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_sccs
126  REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
127  OPTIONAL :: h_stress
128 
129  CHARACTER(LEN=*), PARAMETER :: routinen = 'sccs'
130  REAL(kind=dp), PARAMETER :: epstol = 1.0e-8_dp
131 
132  CHARACTER(LEN=4*default_string_length) :: message
133  CHARACTER(LEN=default_path_length) :: mpi_filename
134  CHARACTER(LEN=default_string_length) :: cube_path, filename, my_pos_cube, &
135  print_path
136  INTEGER :: cube_unit, handle, i, ispin, iter, j, k, &
137  nspin, output_unit, print_level
138  INTEGER(KIND=int_8) :: ngpts
139  INTEGER, DIMENSION(3) :: lb, ub
140  LOGICAL :: append_cube, calculate_stress_tensor, &
141  mpi_io, should_output
142  REAL(kind=dp) :: cavity_surface, cavity_volume, cell_volume, dphi2, dvol, e_tot, &
143  epsilon_solvent, f, polarisation_charge, rho_delta, rho_delta_avg, rho_delta_max, &
144  rho_iter_new, tot_rho_elec, tot_rho_solute
145  REAL(kind=dp), DIMENSION(3) :: abc
146  TYPE(cp_logger_type), POINTER :: logger
147  TYPE(cp_subsys_type), POINTER :: cp_subsys
148  TYPE(dft_control_type), POINTER :: dft_control
149  TYPE(mp_para_env_type), POINTER :: para_env
150  TYPE(particle_list_type), POINTER :: particles
151  TYPE(pw_env_type), POINTER :: pw_env
152  TYPE(pw_poisson_type), POINTER :: poisson_env
153  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
154  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
155  TYPE(pw_r3d_rs_type) :: deps_elec, eps_elec
156  TYPE(pw_r3d_rs_type), DIMENSION(3) :: dln_eps_elec, dphi_tot
157  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_pw_r
158  TYPE(pw_r3d_rs_type), POINTER :: rho_pw_r_sccs
159  TYPE(qs_energy_type), POINTER :: energy
160  TYPE(qs_rho_type), POINTER :: rho
161  TYPE(qs_scf_env_type), POINTER :: scf_env
162  TYPE(sccs_control_type), POINTER :: sccs_control
163  TYPE(section_vals_type), POINTER :: input
164 
165  CALL timeset(routinen, handle)
166 
167  NULLIFY (auxbas_pw_pool)
168  NULLIFY (cp_subsys)
169  NULLIFY (dft_control)
170  NULLIFY (energy)
171  NULLIFY (input)
172  NULLIFY (logger)
173  NULLIFY (para_env)
174  NULLIFY (particles)
175  NULLIFY (poisson_env)
176  NULLIFY (pw_env)
177  NULLIFY (pw_pools)
178  NULLIFY (rho)
179  NULLIFY (sccs_control)
180  NULLIFY (scf_env)
181 
182  ! Load data from Quickstep environment
183  CALL get_qs_env(qs_env=qs_env, &
184  cp_subsys=cp_subsys, &
185  dft_control=dft_control, &
186  energy=energy, &
187  input=input, &
188  para_env=para_env, &
189  pw_env=pw_env, &
190  rho=rho, &
191  scf_env=scf_env)
192  CALL cp_subsys_get(cp_subsys, particles=particles)
193 
194  sccs_control => dft_control%sccs_control
195 
196  cpassert(ASSOCIATED(qs_env))
197 
198  IF (PRESENT(h_stress)) THEN
199  calculate_stress_tensor = .true.
200  h_stress(:, :) = 0.0_dp
201  cpwarn("The stress tensor for SCCS has not yet been fully validated")
202  ELSE
203  calculate_stress_tensor = .false.
204  END IF
205 
206  ! Get access to the PW grid pool
207  CALL pw_env_get(pw_env, &
208  auxbas_pw_pool=auxbas_pw_pool, &
209  pw_pools=pw_pools, &
210  poisson_env=poisson_env)
211 
212  CALL pw_zero(v_sccs)
213 
214  ! Calculate no SCCS contribution, if the requested SCF convergence threshold is not reached yet
215  IF (.NOT. sccs_control%sccs_activated) THEN
216  IF (sccs_control%eps_scf > 0.0_dp) THEN
217  IF ((scf_env%iter_delta > sccs_control%eps_scf) .OR. &
218  ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
219  (qs_env%scf_env%iter_count <= 1))) THEN
220  IF (calculate_stress_tensor) THEN
221  ! Request also the calculation of the stress tensor contribution
222  CALL pw_poisson_solve(poisson_env=poisson_env, &
223  density=rho_tot_gspace, &
224  ehartree=energy%hartree, &
225  vhartree=v_hartree_gspace, &
226  h_stress=h_stress)
227  ELSE
228  CALL pw_poisson_solve(poisson_env=poisson_env, &
229  density=rho_tot_gspace, &
230  ehartree=energy%hartree, &
231  vhartree=v_hartree_gspace)
232  END IF
233  energy%sccs_pol = 0.0_dp
234  energy%sccs_cav = 0.0_dp
235  energy%sccs_dis = 0.0_dp
236  energy%sccs_rep = 0.0_dp
237  energy%sccs_sol = 0.0_dp
238  energy%sccs_hartree = energy%hartree
239  CALL timestop(handle)
240  RETURN
241  END IF
242  END IF
243  sccs_control%sccs_activated = .true.
244  END IF
245 
246  nspin = dft_control%nspins
247 
248  ! Manage print output control
249  logger => cp_get_default_logger()
250  print_level = logger%iter_info%print_level
251  print_path = "DFT%PRINT%SCCS"
252  should_output = (btest(cp_print_key_should_output(logger%iter_info, input, &
253  trim(print_path)), cp_p_file))
254  output_unit = cp_print_key_unit_nr(logger, input, trim(print_path), &
255  extension=".sccs", &
256  ignore_should_output=should_output, &
257  log_filename=.false.)
258 
259  ! Get rho
260  CALL qs_rho_get(rho_struct=rho, &
261  rho_r=rho_pw_r, &
262  rho_r_sccs=rho_pw_r_sccs)
263 
264  ! Retrieve the last rho_iter from the previous SCCS cycle if available
265  cpassert(ASSOCIATED(rho_pw_r_sccs))
266 
267  ! Retrieve the total electronic density in r-space
268  block
269  TYPE(pw_r3d_rs_type) :: rho_elec
270  CALL auxbas_pw_pool%create_pw(rho_elec)
271 
272  ! Retrieve grid parameters
273  ngpts = rho_elec%pw_grid%ngpts
274  dvol = rho_elec%pw_grid%dvol
275  cell_volume = rho_elec%pw_grid%vol
276  abc(1:3) = real(rho_elec%pw_grid%npts(1:3), kind=dp)*rho_elec%pw_grid%dr(1:3)
277  lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
278  ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
279 
280  CALL pw_copy(rho_pw_r(1), rho_elec)
281  DO ispin = 2, nspin
282  CALL pw_axpy(rho_pw_r(ispin), rho_elec)
283  END DO
284  tot_rho_elec = pw_integrate_function(rho_elec)
285 
286  ! Calculate the dielectric (smoothed) function of rho_elec in r-space
287  CALL auxbas_pw_pool%create_pw(eps_elec)
288  CALL auxbas_pw_pool%create_pw(deps_elec)
289  ! Relative permittivity or dielectric constant of the solvent (medium)
290  epsilon_solvent = sccs_control%epsilon_solvent
291  SELECT CASE (sccs_control%method_id)
292  CASE (sccs_andreussi)
293  CALL andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%rho_max, &
294  sccs_control%rho_min)
295  CASE (sccs_fattebert_gygi)
296  CALL fattebert_gygi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%beta, &
297  sccs_control%rho_zero)
298  CASE DEFAULT
299  cpabort("Invalid method specified for SCCS model")
300  END SELECT
301 
302  ! Optional output of the dielectric function in cube file format
303  filename = "DIELECTRIC_FUNCTION"
304  cube_path = trim(print_path)//"%"//trim(filename)
305  IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
306  cp_p_file)) THEN
307  append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
308  my_pos_cube = "REWIND"
309  IF (append_cube) my_pos_cube = "APPEND"
310  mpi_io = .true.
311  cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
312  extension=".cube", middle_name=trim(filename), &
313  file_position=my_pos_cube, log_filename=.false., &
314  mpi_io=mpi_io, fout=mpi_filename)
315  IF (output_unit > 0) THEN
316  IF (.NOT. mpi_io) THEN
317  INQUIRE (unit=cube_unit, name=filename)
318  ELSE
319  filename = mpi_filename
320  END IF
321  WRITE (unit=output_unit, fmt="(T3,A)") &
322  "SCCS| The dielectric function is written in cube file format to the file:", &
323  "SCCS| "//trim(filename)
324  END IF
325  CALL cp_pw_to_cube(eps_elec, cube_unit, trim(filename), particles=particles, &
326  stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), &
327  mpi_io=mpi_io)
328  CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
329  END IF
330 
331  ! Calculate the (quantum) volume and surface of the solute cavity
332  cavity_surface = 0.0_dp
333  cavity_volume = 0.0_dp
334 
335  IF (abs(epsilon_solvent - 1.0_dp) > epstol) THEN
336 
337  block
338  TYPE(pw_r3d_rs_type) :: theta, norm_drho_elec
339  TYPE(pw_r3d_rs_type), DIMENSION(3) :: drho_elec
340  CALL auxbas_pw_pool%create_pw(theta)
341  CALL pw_zero(theta)
342 
343  ! Calculate the (quantum) volume of the solute cavity
344  f = 1.0_dp/(epsilon_solvent - 1.0_dp)
345 !$OMP PARALLEL DO DEFAULT(NONE) &
346 !$OMP PRIVATE(i,j,k) &
347 !$OMP SHARED(epsilon_solvent,eps_elec,f,lb,theta,ub)
348  DO k = lb(3), ub(3)
349  DO j = lb(2), ub(2)
350  DO i = lb(1), ub(1)
351  theta%array(i, j, k) = f*(epsilon_solvent - eps_elec%array(i, j, k))
352  END DO
353  END DO
354  END DO
355 !$OMP END PARALLEL DO
356  cavity_volume = pw_integrate_function(theta)
357 
358  ! Calculate the derivative of the electronic density in r-space
359  ! TODO: Could be retrieved from the QS environment
360  DO i = 1, 3
361  CALL auxbas_pw_pool%create_pw(drho_elec(i))
362  END DO
363  CALL derive(rho_elec, drho_elec, sccs_derivative_fft, pw_env, input)
364 
365  CALL auxbas_pw_pool%create_pw(norm_drho_elec)
366 
367  ! Calculate the norm of the gradient of the electronic density in r-space
368 !$OMP PARALLEL DO DEFAULT(NONE) &
369 !$OMP PRIVATE(i,j,k) &
370 !$OMP SHARED(drho_elec,lb,norm_drho_elec,ub)
371  DO k = lb(3), ub(3)
372  DO j = lb(2), ub(2)
373  DO i = lb(1), ub(1)
374  norm_drho_elec%array(i, j, k) = sqrt(drho_elec(1)%array(i, j, k)* &
375  drho_elec(1)%array(i, j, k) + &
376  drho_elec(2)%array(i, j, k)* &
377  drho_elec(2)%array(i, j, k) + &
378  drho_elec(3)%array(i, j, k)* &
379  drho_elec(3)%array(i, j, k))
380  END DO
381  END DO
382  END DO
383 !$OMP END PARALLEL DO
384 
385  ! Optional output of the norm of the density gradient in cube file format
386  filename = "DENSITY_GRADIENT"
387  cube_path = trim(print_path)//"%"//trim(filename)
388  IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
389  cp_p_file)) THEN
390  append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
391  my_pos_cube = "REWIND"
392  IF (append_cube) my_pos_cube = "APPEND"
393  mpi_io = .true.
394  cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
395  extension=".cube", middle_name=trim(filename), &
396  file_position=my_pos_cube, log_filename=.false., &
397  mpi_io=mpi_io, fout=mpi_filename)
398  IF (output_unit > 0) THEN
399  IF (.NOT. mpi_io) THEN
400  INQUIRE (unit=cube_unit, name=filename)
401  ELSE
402  filename = mpi_filename
403  END IF
404  WRITE (unit=output_unit, fmt="(T3,A)") &
405  "SCCS| The norm of the density gradient is written in cube file format to the file:", &
406  "SCCS| "//trim(filename)
407  END IF
408  CALL cp_pw_to_cube(norm_drho_elec, cube_unit, trim(filename), particles=particles, &
409  stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), &
410  mpi_io=mpi_io)
411  CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
412  END IF
413 
414  ! Calculate the (quantum) surface of the solute cavity
415  SELECT CASE (sccs_control%method_id)
416  CASE (sccs_andreussi)
417  CALL surface_andreussi(rho_elec, norm_drho_elec, theta, epsilon_solvent, &
418  sccs_control%rho_max, sccs_control%rho_min, &
419  sccs_control%delta_rho)
420  CASE (sccs_fattebert_gygi)
421  CALL surface_fattebert_gygi(rho_elec, norm_drho_elec, theta, epsilon_solvent, &
422  sccs_control%beta, sccs_control%rho_zero, &
423  sccs_control%delta_rho)
424  CASE DEFAULT
425  cpabort("Invalid method specified for SCCS model")
426  END SELECT
427  cavity_surface = pw_integrate_function(theta)
428 
429  ! Release storage
430  CALL auxbas_pw_pool%give_back_pw(theta)
431  CALL auxbas_pw_pool%give_back_pw(norm_drho_elec)
432  DO i = 1, 3
433  CALL auxbas_pw_pool%give_back_pw(drho_elec(i))
434  END DO
435  END block
436 
437  END IF ! epsilon_solvent > 1
438 
439  CALL auxbas_pw_pool%give_back_pw(rho_elec)
440  END block
441 
442  block
443  TYPE(pw_r3d_rs_type) :: rho_tot, phi_tot, rho_solute, rho_tot_zero
444  ! Retrieve the total charge density (core + elec) of the solute in r-space
445  CALL auxbas_pw_pool%create_pw(rho_solute)
446  CALL pw_zero(rho_solute)
447  CALL pw_transfer(rho_tot_gspace, rho_solute)
448  tot_rho_solute = pw_integrate_function(rho_solute)
449 
450  ! Check total charge
451  IF (abs(tot_rho_solute) >= 1.0e-6_dp) THEN
452  IF ((poisson_env%parameters%solver /= pw_poisson_analytic) .AND. &
453  (poisson_env%parameters%solver /= pw_poisson_mt)) THEN
454  WRITE (unit=message, fmt="(A,SP,F0.6,A)") &
455  "The system (solute) has a non-negligible charge of ", -tot_rho_solute, &
456  ". It is recommended to use non-periodic boundary conditions (PERIODIC none) "// &
457  "combined with an appropriate Poisson solver (POISSON_SOLVER MT or analytic)"
458  cpwarn(message)
459  END IF
460  END IF
461 
462  ! Reassign work storage to rho_tot_zero, because rho_elec is no longer needed
463  CALL auxbas_pw_pool%create_pw(rho_tot_zero)
464 
465  ! Build the initial (rho_iter = 0) total charge density (solute plus polarisation) in r-space
466  ! eps_elec <- ln(eps_elec)
467 !$OMP PARALLEL DO DEFAULT(NONE) &
468 !$OMP PRIVATE(i,j,k) &
469 !$OMP SHARED(eps_elec,lb,message,output_unit,para_env,ub) &
470 !$OMP SHARED(rho_solute,rho_tot_zero)
471  DO k = lb(3), ub(3)
472  DO j = lb(2), ub(2)
473  DO i = lb(1), ub(1)
474  IF (eps_elec%array(i, j, k) < 1.0_dp) THEN
475  WRITE (unit=message, fmt="(A,ES12.3,A,3(I0,A))") &
476  "SCCS| Invalid dielectric function value ", eps_elec%array(i, j, k), &
477  " encountered at grid point (", i, ",", j, ",", k, ")"
478  cpabort(message)
479  END IF
480  rho_tot_zero%array(i, j, k) = rho_solute%array(i, j, k)/eps_elec%array(i, j, k)
481  eps_elec%array(i, j, k) = log(eps_elec%array(i, j, k))
482  END DO
483  END DO
484  END DO
485 !$OMP END PARALLEL DO
486 
487  ! Build the derivative of LOG(eps_elec)
488  DO i = 1, 3
489  CALL auxbas_pw_pool%create_pw(dln_eps_elec(i))
490  CALL pw_zero(dln_eps_elec(i))
491  END DO
492  CALL derive(eps_elec, dln_eps_elec, sccs_control%derivative_method, pw_env, input)
493  CALL auxbas_pw_pool%give_back_pw(eps_elec)
494 
495  ! Print header for the SCCS cycle
496  IF (should_output .AND. (output_unit > 0)) THEN
497  IF (print_level > low_print_level) THEN
498  WRITE (unit=output_unit, fmt="(T3,A,T56,F25.12)") &
499  "SCCS| Total electronic charge density ", -tot_rho_elec, &
500  "SCCS| Total charge density (solute) ", -tot_rho_solute
501  WRITE (unit=output_unit, fmt="(T3,A,T56,F25.3)") &
502  "SCCS| Volume of the cell [bohr^3]", cell_volume, &
503  "SCCS| [angstrom^3]", &
504  cp_unit_from_cp2k(cell_volume, "angstrom^3")
505  IF (abs(epsilon_solvent - 1.0_dp) > epstol) THEN
506  WRITE (unit=output_unit, fmt="(T3,A,T56,F25.3)") &
507  "SCCS| Volume of the solute cavity [bohr^3]", cavity_volume, &
508  "SCCS| [angstrom^3]", &
509  cp_unit_from_cp2k(cavity_volume, "angstrom^3"), &
510  "SCCS| Surface of the solute cavity [bohr^2]", cavity_surface, &
511  "SCCS| [angstrom^2]", &
512  cp_unit_from_cp2k(cavity_surface, "angstrom^2")
513  END IF
514  WRITE (unit=output_unit, fmt="(T3,A)") &
515  "SCCS|", &
516  "SCCS| Step Average residual Maximum residual E(Hartree) [a.u.]"
517  END IF
518  END IF
519 
520  ! Get storage for the derivative of the total potential (field) in r-space
521  DO i = 1, 3
522  CALL auxbas_pw_pool%create_pw(dphi_tot(i))
523  END DO
524 
525  ! Initialise the total charge density in r-space rho_tot with rho_tot_zero + rho_iter_zero
526  CALL auxbas_pw_pool%create_pw(rho_tot)
527  CALL pw_copy(rho_tot_zero, rho_tot)
528  CALL pw_axpy(rho_pw_r_sccs, rho_tot)
529 
530  CALL auxbas_pw_pool%create_pw(phi_tot)
531  CALL pw_zero(phi_tot)
532 
533  ! Main SCCS iteration loop
534  iter = 0
535 
536  iter_loop: DO
537 
538  ! Increment iteration counter
539  iter = iter + 1
540 
541  ! Check if the requested maximum number of SCCS iterations is reached
542  IF (iter > sccs_control%max_iter) THEN
543  IF (output_unit > 0) THEN
544  WRITE (unit=output_unit, fmt="(T3,A,/,T3,A,I0,A)") &
545  "SCCS| Maximum number of SCCS iterations reached", &
546  "SCCS| Iteration cycle did not converge in ", sccs_control%max_iter, " steps"
547  ELSE
548  WRITE (unit=message, fmt="(A,I0,A)") &
549  "The SCCS iteration cycle did not converge in ", sccs_control%max_iter, " steps"
550  cpwarn(message)
551  END IF
552  EXIT iter_loop
553  END IF
554 
555  ! Calculate derivative of the current total potential in r-space
556  CALL pw_poisson_solve(poisson_env=poisson_env, &
557  density=rho_tot, &
558  vhartree=phi_tot, &
559  dvhartree=dphi_tot)
560  energy%sccs_hartree = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
561 
562  ! Update total charge density (solute plus polarisation) in r-space
563  ! based on the iterated polarisation charge density
564  f = 1.0_dp/fourpi
565  rho_delta_avg = 0.0_dp
566  rho_delta_max = 0.0_dp
567 !$OMP PARALLEL DO DEFAULT(NONE) &
568 !$OMP PRIVATE(i,j,k,rho_delta,rho_iter_new) &
569 !$OMP SHARED(dln_eps_elec,dphi_tot,f,lb,rho_pw_r_sccs,ub) &
570 !$OMP SHARED(rho_tot,rho_tot_zero,sccs_control) &
571 !$OMP REDUCTION(+:rho_delta_avg) &
572 !$OMP REDUCTION(MAX:rho_delta_max)
573  DO k = lb(3), ub(3)
574  DO j = lb(2), ub(2)
575  DO i = lb(1), ub(1)
576  rho_iter_new = (dln_eps_elec(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + &
577  dln_eps_elec(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + &
578  dln_eps_elec(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k))*f
579  rho_iter_new = rho_pw_r_sccs%array(i, j, k) + &
580  sccs_control%mixing*(rho_iter_new - rho_pw_r_sccs%array(i, j, k))
581  rho_delta = abs(rho_iter_new - rho_pw_r_sccs%array(i, j, k))
582  rho_delta_max = max(rho_delta, rho_delta_max)
583  rho_delta_avg = rho_delta_avg + rho_delta
584  rho_tot%array(i, j, k) = rho_tot_zero%array(i, j, k) + rho_iter_new
585  rho_pw_r_sccs%array(i, j, k) = rho_iter_new
586  END DO
587  END DO
588  END DO
589 !$OMP END PARALLEL DO
590 
591  CALL para_env%sum(rho_delta_avg)
592  rho_delta_avg = rho_delta_avg/real(ngpts, kind=dp)
593  CALL para_env%max(rho_delta_max)
594 
595  IF (should_output .AND. (output_unit > 0)) THEN
596  IF (print_level > low_print_level) THEN
597  IF ((abs(rho_delta_avg) < 1.0e-8_dp) .OR. &
598  (abs(rho_delta_avg) >= 1.0e5_dp)) THEN
599  WRITE (unit=output_unit, fmt="(T3,A,I6,4X,ES16.4,4X,ES16.4,1X,F25.12)") &
600  "SCCS| ", iter, rho_delta_avg, rho_delta_max, energy%sccs_hartree
601  ELSE
602  WRITE (unit=output_unit, fmt="(T3,A,I6,4X,F16.8,4X,F16.8,1X,F25.12)") &
603  "SCCS| ", iter, rho_delta_avg, rho_delta_max, energy%sccs_hartree
604  END IF
605  END IF
606  END IF
607 
608  ! Check if the SCCS iteration cycle is converged to the requested tolerance
609  IF (rho_delta_max <= sccs_control%eps_sccs) THEN
610  IF (should_output .AND. (output_unit > 0)) THEN
611  WRITE (unit=output_unit, fmt="(T3,A,I0,A)") &
612  "SCCS| Iteration cycle converged in ", iter, " steps"
613  END IF
614  EXIT iter_loop
615  END IF
616 
617  END DO iter_loop
618 
619  ! Release work storage which is no longer needed
620  CALL auxbas_pw_pool%give_back_pw(rho_tot_zero)
621  DO i = 1, 3
622  CALL auxbas_pw_pool%give_back_pw(dln_eps_elec(i))
623  END DO
624 
625  ! Optional output of the total charge density in cube file format
626  filename = "TOTAL_CHARGE_DENSITY"
627  cube_path = trim(print_path)//"%"//trim(filename)
628  IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), cp_p_file)) THEN
629  append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
630  my_pos_cube = "REWIND"
631  IF (append_cube) my_pos_cube = "APPEND"
632  mpi_io = .true.
633  cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
634  extension=".cube", middle_name=trim(filename), &
635  file_position=my_pos_cube, log_filename=.false., &
636  mpi_io=mpi_io, fout=mpi_filename)
637  IF (output_unit > 0) THEN
638  IF (.NOT. mpi_io) THEN
639  INQUIRE (unit=cube_unit, name=filename)
640  ELSE
641  filename = mpi_filename
642  END IF
643  WRITE (unit=output_unit, fmt="(T3,A)") &
644  "SCCS| The total SCCS charge density is written in cube file format to the file:", &
645  "SCCS| "//trim(filename)
646  END IF
647  CALL cp_pw_to_cube(rho_tot, cube_unit, trim(filename), particles=particles, &
648  stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), mpi_io=mpi_io)
649  CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
650  END IF
651 
652  ! Calculate the total SCCS Hartree energy, potential, and its
653  ! derivatives of the solute and the implicit solvent
654  CALL pw_transfer(rho_tot, rho_tot_gspace)
655  IF (calculate_stress_tensor) THEN
656  ! Request also the calculation of the stress tensor contribution
657  CALL pw_poisson_solve(poisson_env=poisson_env, &
658  density=rho_tot_gspace, &
659  ehartree=e_tot, &
660  vhartree=v_hartree_gspace, &
661  dvhartree=dphi_tot, &
662  h_stress=h_stress)
663  ELSE
664  CALL pw_poisson_solve(poisson_env=poisson_env, &
665  density=rho_tot_gspace, &
666  ehartree=e_tot, &
667  vhartree=v_hartree_gspace, &
668  dvhartree=dphi_tot)
669  END IF
670  CALL pw_transfer(v_hartree_gspace, phi_tot)
671  energy%sccs_hartree = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
672 
673  ! Calculate the Hartree energy and potential of the solute only
674  block
675  TYPE(pw_r3d_rs_type) :: phi_solute
676  CALL auxbas_pw_pool%create_pw(phi_solute)
677  CALL pw_zero(phi_solute)
678  CALL pw_poisson_solve(poisson_env=poisson_env, &
679  density=rho_solute, &
680  ehartree=energy%hartree, &
681  vhartree=phi_solute)
682 
683  ! Calculate the polarisation potential (store it in phi_tot)
684  ! phi_pol = phi_tot - phi_solute
685  CALL pw_axpy(phi_solute, phi_tot, alpha=-1.0_dp)
686  CALL auxbas_pw_pool%give_back_pw(phi_solute)
687  END block
688 
689  ! Optional output of the SCCS polarisation potential in cube file format
690  filename = "POLARISATION_POTENTIAL"
691  cube_path = trim(print_path)//"%"//trim(filename)
692  IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
693  cp_p_file)) THEN
694  append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
695  my_pos_cube = "REWIND"
696  IF (append_cube) my_pos_cube = "APPEND"
697  mpi_io = .true.
698  cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
699  extension=".cube", middle_name=trim(filename), &
700  file_position=my_pos_cube, log_filename=.false., &
701  mpi_io=mpi_io, fout=mpi_filename)
702  IF (output_unit > 0) THEN
703  IF (.NOT. mpi_io) THEN
704  INQUIRE (unit=cube_unit, name=filename)
705  ELSE
706  filename = mpi_filename
707  END IF
708  WRITE (unit=output_unit, fmt="(T3,A)") &
709  "SCCS| The SCCS polarisation potential is written in cube file format to the file:", &
710  "SCCS| "//trim(filename)
711  END IF
712  CALL cp_pw_to_cube(phi_tot, cube_unit, trim(filename), particles=particles, &
713  stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), mpi_io=mpi_io)
714  CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
715  END IF
716 
717  ! Calculate the polarisation charge (store it in rho_tot)
718  ! rho_pol = rho_tot - rho_solute
719  CALL pw_axpy(rho_solute, rho_tot, alpha=-1.0_dp)
720  polarisation_charge = pw_integrate_function(rho_tot)
721 
722  ! Optional output of the SCCS polarisation charge in cube file format
723  filename = "POLARISATION_CHARGE_DENSITY"
724  cube_path = trim(print_path)//"%"//trim(filename)
725  IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
726  cp_p_file)) THEN
727  append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
728  my_pos_cube = "REWIND"
729  IF (append_cube) my_pos_cube = "APPEND"
730  mpi_io = .true.
731  cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
732  extension=".cube", middle_name=trim(filename), &
733  file_position=my_pos_cube, log_filename=.false., &
734  mpi_io=mpi_io, fout=mpi_filename)
735  IF (output_unit > 0) THEN
736  IF (.NOT. mpi_io) THEN
737  INQUIRE (unit=cube_unit, name=filename)
738  ELSE
739  filename = mpi_filename
740  END IF
741  WRITE (unit=output_unit, fmt="(T3,A)") &
742  "SCCS| The SCCS polarisation charge density is written in cube file format to the file:", &
743  "SCCS| "//trim(filename)
744  END IF
745  CALL cp_pw_to_cube(rho_tot, cube_unit, trim(filename), particles=particles, &
746  stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), mpi_io=mpi_io)
747  CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
748  END IF
749 
750  ! Calculate SCCS polarisation energy
751  energy%sccs_pol = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
752  CALL auxbas_pw_pool%give_back_pw(rho_solute)
753  CALL auxbas_pw_pool%give_back_pw(phi_tot)
754  CALL auxbas_pw_pool%give_back_pw(rho_tot)
755  END block
756 
757  ! Calculate additional solvation terms
758  energy%sccs_cav = sccs_control%gamma_solvent*cavity_surface
759  energy%sccs_dis = sccs_control%beta_solvent*cavity_volume
760  energy%sccs_rep = sccs_control%alpha_solvent*cavity_surface
761  ! Calculate solvation free energy: \delta G^el + (alpha + gamma)*S + beta*V
762  energy%sccs_sol = energy%sccs_pol + energy%sccs_rep + energy%sccs_cav + energy%sccs_dis
763 
764  IF (should_output .AND. (output_unit > 0)) THEN
765  WRITE (unit=output_unit, fmt="(T3,A)") &
766  "SCCS|"
767  WRITE (unit=output_unit, fmt="(T3,A,T56,F25.12)") &
768  "SCCS| Polarisation charge", polarisation_charge
769  !MK "SCCS| Total interaction energy [a.u.]", e_tot
770  WRITE (unit=output_unit, fmt="(T3,A)") &
771  "SCCS|"
772  CALL print_sccs_results(energy, sccs_control, output_unit)
773  END IF
774 
775  ! Calculate SCCS contribution to the Kohn-Sham potential
776  f = -0.5_dp*dvol/fourpi
777 !$OMP PARALLEL DO DEFAULT(NONE) &
778 !$OMP PRIVATE(dphi2,i,j,k) &
779 !$OMP SHARED(f,deps_elec,dphi_tot) &
780 !$OMP SHARED(lb,ub,v_sccs)
781  DO k = lb(3), ub(3)
782  DO j = lb(2), ub(2)
783  DO i = lb(1), ub(1)
784  dphi2 = dphi_tot(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + &
785  dphi_tot(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + &
786  dphi_tot(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k)
787  v_sccs%array(i, j, k) = v_sccs%array(i, j, k) + f*deps_elec%array(i, j, k)*dphi2
788  END DO
789  END DO
790  END DO
791 !$OMP END PARALLEL DO
792 
793  CALL auxbas_pw_pool%give_back_pw(deps_elec)
794  DO i = 1, 3
795  CALL auxbas_pw_pool%give_back_pw(dphi_tot(i))
796  END DO
797 
798  ! Release the SCCS printout environment
799  CALL cp_print_key_finished_output(output_unit, logger, input, trim(print_path), &
800  ignore_should_output=should_output)
801 
802  CALL timestop(handle)
803 
804  END SUBROUTINE sccs
805 
806 ! **************************************************************************************************
807 !> \brief Calculate the smoothed dielectric function of Andreussi et al.
808 !> \param rho_elec ...
809 !> \param eps_elec ...
810 !> \param deps_elec ...
811 !> \param epsilon_solvent ...
812 !> \param rho_max ...
813 !> \param rho_min ...
814 !> \par History:
815 !> - Creation (16.10.2013,MK)
816 !> - Finite difference of isosurfaces implemented (21.12.2013,MK)
817 !> \author Matthias Krack (MK)
818 !> \version 1.1
819 ! **************************************************************************************************
820  SUBROUTINE andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, rho_max, &
821  rho_min)
822 
823  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec
824  REAL(kind=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min
825 
826  CHARACTER(LEN=*), PARAMETER :: routinen = 'andreussi'
827  REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
828 
829  INTEGER :: handle, i, j, k
830  INTEGER, DIMENSION(3) :: lb, ub
831  REAL(kind=dp) :: diff, dq, dt, f, ln_rho_max, ln_rho_min, &
832  q, rho, t, x, y
833 
834  CALL timeset(routinen, handle)
835 
836  f = log(epsilon_solvent)/twopi
837  diff = rho_max - rho_min
838  IF (diff < sqrt(rhotol)) cpabort("SCCS: Difference between rho(min) and rho(max) is too small")
839  IF (rho_min >= rhotol) THEN
840  ln_rho_max = log(rho_max)
841  ln_rho_min = log(rho_min)
842  q = twopi/(ln_rho_max - ln_rho_min)
843  dq = -f*q
844  END IF
845 
846  lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
847  ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
848 
849  ! Calculate the dielectric function and its derivative
850 !$OMP PARALLEL DO DEFAULT(NONE) &
851 !$OMP PRIVATE(dt,i,j,k,rho,t,x,y) &
852 !$OMP SHARED(deps_elec,dq,eps_elec,epsilon_solvent,f,lb,ub) &
853 !$OMP SHARED(ln_rho_max,rho_elec,q,rho_max,rho_min)
854  DO k = lb(3), ub(3)
855  DO j = lb(2), ub(2)
856  DO i = lb(1), ub(1)
857  rho = rho_elec%array(i, j, k)
858  IF (rho < rho_min) THEN
859  eps_elec%array(i, j, k) = epsilon_solvent
860  deps_elec%array(i, j, k) = 0.0_dp
861  ELSE IF (rho <= rho_max) THEN
862  x = log(rho)
863  y = q*(ln_rho_max - x)
864  t = f*(y - sin(y))
865  eps_elec%array(i, j, k) = exp(t)
866  dt = dq*(1.0_dp - cos(y))
867  deps_elec%array(i, j, k) = eps_elec%array(i, j, k)*dt/rho
868  ELSE
869  eps_elec%array(i, j, k) = 1.0_dp
870  deps_elec%array(i, j, k) = 0.0_dp
871  END IF
872  END DO
873  END DO
874  END DO
875 !$OMP END PARALLEL DO
876 
877  CALL timestop(handle)
878 
879  END SUBROUTINE andreussi
880 
881 ! **************************************************************************************************
882 !> \brief Calculate the smoothed dielectric function of Fattebert and Gygi
883 !> \param rho_elec ...
884 !> \param eps_elec ...
885 !> \param deps_elec ...
886 !> \param epsilon_solvent ...
887 !> \param beta ...
888 !> \param rho_zero ...
889 !> \par History:
890 !> - Creation (15.10.2013,MK)
891 !> \author Matthias Krack (MK)
892 !> \version 1.0
893 ! **************************************************************************************************
894  SUBROUTINE fattebert_gygi(rho_elec, eps_elec, deps_elec, epsilon_solvent, beta, &
895  rho_zero)
896 
897  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec
898  REAL(kind=dp), INTENT(IN) :: epsilon_solvent, beta, rho_zero
899 
900  CHARACTER(LEN=*), PARAMETER :: routinen = 'fattebert_gygi'
901  REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
902 
903  INTEGER :: handle, i, j, k
904  INTEGER, DIMENSION(3) :: lb, ub
905  REAL(kind=dp) :: df, f, p, q, rho, s, t, twobeta
906 
907  CALL timeset(routinen, handle)
908 
909  df = (1.0_dp - epsilon_solvent)/rho_zero
910  f = 0.5_dp*(epsilon_solvent - 1.0_dp)
911  q = 1.0_dp/rho_zero
912  twobeta = 2.0_dp*beta
913 
914  lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
915  ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
916 
917  ! Calculate the smoothed dielectric function and its derivative
918 !$OMP PARALLEL DO DEFAULT(NONE) &
919 !$OMP PRIVATE(i,j,k,p,rho,s,t) &
920 !$OMP SHARED(df,deps_elec,eps_elec,epsilon_solvent,f,lb,ub) &
921 !$OMP SHARED(q,rho_elec,twobeta)
922  DO k = lb(3), ub(3)
923  DO j = lb(2), ub(2)
924  DO i = lb(1), ub(1)
925  rho = rho_elec%array(i, j, k)
926  IF (rho < rhotol) THEN
927  eps_elec%array(i, j, k) = epsilon_solvent
928  deps_elec%array(i, j, k) = 0.0_dp
929  ELSE
930  s = rho*q
931  p = s**twobeta
932  t = 1.0_dp/(1.0_dp + p)
933  eps_elec%array(i, j, k) = 1.0_dp + f*(1.0_dp + (1.0_dp - p)*t)
934  deps_elec%array(i, j, k) = df*twobeta*t*t*p/s
935  END IF
936  END DO
937  END DO
938  END DO
939 !$OMP END PARALLEL DO
940 
941  CALL timestop(handle)
942 
943  END SUBROUTINE fattebert_gygi
944 
945 ! **************************************************************************************************
946 !> \brief Build the numerical derivative of a function on realspace grid
947 !> \param f ...
948 !> \param df ...
949 !> \param method ...
950 !> \param pw_env ...
951 !> \param input ...
952 !> \par History:
953 !> - Creation (15.11.2013,MK)
954 !> \author Matthias Krack (MK)
955 !> \version 1.0
956 ! **************************************************************************************************
957  SUBROUTINE derive(f, df, method, pw_env, input)
958 
959  TYPE(pw_r3d_rs_type), INTENT(IN) :: f
960  TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: df
961  INTEGER, INTENT(IN) :: method
962  TYPE(pw_env_type), POINTER :: pw_env
963  TYPE(section_vals_type), POINTER :: input
964 
965  CHARACTER(LEN=*), PARAMETER :: routinen = 'derive'
966 
967  INTEGER :: border_points, handle, i
968  INTEGER, DIMENSION(3) :: lb, n, ub
969  TYPE(pw_c1d_gs_type), DIMENSION(2) :: work_g1d
970  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
971  TYPE(realspace_grid_desc_type), POINTER :: rs_desc
972  TYPE(realspace_grid_input_type) :: input_settings
973  TYPE(realspace_grid_type), POINTER :: rs_grid
974  TYPE(section_vals_type), POINTER :: rs_grid_section
975 
976  CALL timeset(routinen, handle)
977 
978  cpassert(ASSOCIATED(pw_env))
979 
980  ! Perform method specific setup
981  SELECT CASE (method)
982  CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7)
983  NULLIFY (rs_desc)
984  rs_grid_section => section_vals_get_subs_vals(input, "DFT%MGRID%RS_GRID")
985  SELECT CASE (method)
986  CASE (sccs_derivative_cd3)
987  border_points = 1
988  CASE (sccs_derivative_cd5)
989  border_points = 2
990  CASE (sccs_derivative_cd7)
991  border_points = 3
992  END SELECT
993  CALL init_input_type(input_settings, 2*border_points + 1, rs_grid_section, &
994  1, (/-1, -1, -1/))
995  CALL rs_grid_create_descriptor(rs_desc, f%pw_grid, input_settings, &
996  border_points=border_points)
997  ALLOCATE (rs_grid)
998  CALL rs_grid_create(rs_grid, rs_desc)
999 !MK CALL rs_grid_print(rs_grid, 6)
1000  CASE (sccs_derivative_fft)
1001  lb(1:3) = f%pw_grid%bounds_local(1, 1:3)
1002  ub(1:3) = f%pw_grid%bounds_local(2, 1:3)
1003  NULLIFY (auxbas_pw_pool)
1004  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1005  ! Get work storage for the 1d grids in g-space (derivative calculation)
1006  DO i = 1, SIZE(work_g1d)
1007  CALL auxbas_pw_pool%create_pw(work_g1d(i))
1008  END DO
1009  END SELECT
1010 
1011  ! Calculate the derivatives
1012  SELECT CASE (method)
1013  CASE (sccs_derivative_cd3)
1014  CALL derive_fdm_cd3(f, df, rs_grid)
1015  CASE (sccs_derivative_cd5)
1016  CALL derive_fdm_cd5(f, df, rs_grid)
1017  CASE (sccs_derivative_cd7)
1018  CALL derive_fdm_cd7(f, df, rs_grid)
1019  CASE (sccs_derivative_fft)
1020  ! FFT
1021  CALL pw_transfer(f, work_g1d(1))
1022  DO i = 1, 3
1023  n(:) = 0
1024  n(i) = 1
1025  CALL pw_copy(work_g1d(1), work_g1d(2))
1026  CALL pw_derive(work_g1d(2), n(:))
1027  CALL pw_transfer(work_g1d(2), df(i))
1028  END DO
1029  CASE DEFAULT
1030  cpabort("Invalid derivative method for SCCS specified")
1031  END SELECT
1032 
1033  ! Perform method specific cleanup
1034  SELECT CASE (method)
1035  CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7)
1036  CALL rs_grid_release(rs_grid)
1037  DEALLOCATE (rs_grid)
1038  CALL rs_grid_release_descriptor(rs_desc)
1039  CASE (sccs_derivative_fft)
1040  DO i = 1, SIZE(work_g1d)
1041  CALL auxbas_pw_pool%give_back_pw(work_g1d(i))
1042  END DO
1043  END SELECT
1044 
1045  CALL timestop(handle)
1046 
1047  END SUBROUTINE derive
1048 
1049 ! **************************************************************************************************
1050 !> \brief Calculate the finite difference between two isosurfaces of the
1051 !> electronic density. The smoothed dielectric function of
1052 !> Andreussi et al. is used as switching function eventually
1053 !> defining the quantum volume and surface of the cavity.
1054 !> \param rho_elec ...
1055 !> \param norm_drho_elec ...
1056 !> \param dtheta ...
1057 !> \param epsilon_solvent ...
1058 !> \param rho_max ...
1059 !> \param rho_min ...
1060 !> \param delta_rho ...
1061 !> \par History:
1062 !> - Creation (21.12.2013,MK)
1063 !> \author Matthias Krack (MK)
1064 !> \version 1.0
1065 ! **************************************************************************************************
1066  SUBROUTINE surface_andreussi(rho_elec, norm_drho_elec, dtheta, &
1067  epsilon_solvent, rho_max, rho_min, delta_rho)
1068 
1069  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, norm_drho_elec, dtheta
1070  REAL(kind=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min, &
1071  delta_rho
1072 
1073  CHARACTER(LEN=*), PARAMETER :: routinen = 'surface_andreussi'
1074  REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1075 
1076  INTEGER :: handle, i, j, k, l
1077  INTEGER, DIMENSION(3) :: lb, ub
1078  REAL(kind=dp) :: diff, e, eps_elec, f, ln_rho_max, &
1079  ln_rho_min, q, rho, t, x, y
1080  REAL(kind=dp), DIMENSION(2) :: theta
1081 
1082  CALL timeset(routinen, handle)
1083 
1084  e = epsilon_solvent - 1.0_dp
1085  f = log(epsilon_solvent)/twopi
1086  diff = rho_max - rho_min
1087  IF (diff < sqrt(rhotol)) cpabort("SCCS: Difference between rho(min) and rho(max) is too small")
1088  IF (rho_min >= rhotol) THEN
1089  ln_rho_max = log(rho_max)
1090  ln_rho_min = log(rho_min)
1091  q = twopi/(ln_rho_max - ln_rho_min)
1092  END IF
1093 
1094  lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1095  ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1096 
1097  ! Calculate finite difference between two isosurfaces
1098 !$OMP PARALLEL DO DEFAULT(NONE) &
1099 !$OMP PRIVATE(eps_elec,i,j,k,l,rho,t,theta,x,y) &
1100 !$OMP SHARED(delta_rho,dtheta,e,epsilon_solvent,f,lb) &
1101 !$OMP SHARED(ln_rho_max,norm_drho_elec,rho_elec,q,rho_max,rho_min,ub)
1102  DO k = lb(3), ub(3)
1103  DO j = lb(2), ub(2)
1104  DO i = lb(1), ub(1)
1105  DO l = 1, 2
1106  rho = rho_elec%array(i, j, k) + (real(l, kind=dp) - 1.5_dp)*delta_rho
1107  IF (rho < rho_min) THEN
1108  eps_elec = epsilon_solvent
1109  ELSE IF (rho <= rho_max) THEN
1110  x = log(rho)
1111  y = q*(ln_rho_max - x)
1112  t = f*(y - sin(y))
1113  eps_elec = exp(t)
1114  ELSE
1115  eps_elec = 1.0_dp
1116  END IF
1117  theta(l) = (epsilon_solvent - eps_elec)/e
1118  END DO
1119  dtheta%array(i, j, k) = (theta(2) - theta(1))*norm_drho_elec%array(i, j, k)/delta_rho
1120  END DO
1121  END DO
1122  END DO
1123 !$OMP END PARALLEL DO
1124 
1125  CALL timestop(handle)
1126 
1127  END SUBROUTINE surface_andreussi
1128 
1129 ! **************************************************************************************************
1130 !> \brief Calculate the finite difference between two isosurfaces of the
1131 !> the electronic density. The smoothed dielectric function of
1132 !> Fattebert and Gygi is used as switching function eventually
1133 !> defining the quantum volume and surface of the cavity.
1134 !> \param rho_elec ...
1135 !> \param norm_drho_elec ...
1136 !> \param dtheta ...
1137 !> \param epsilon_solvent ...
1138 !> \param beta ...
1139 !> \param rho_zero ...
1140 !> \param delta_rho ...
1141 !> \par History:
1142 !> - Creation (21.12.2013,MK)
1143 !> \author Matthias Krack (MK)
1144 !> \version 1.0
1145 ! **************************************************************************************************
1146  SUBROUTINE surface_fattebert_gygi(rho_elec, norm_drho_elec, dtheta, &
1147  epsilon_solvent, beta, rho_zero, delta_rho)
1148 
1149  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, norm_drho_elec, dtheta
1150  REAL(kind=dp), INTENT(IN) :: epsilon_solvent, beta, rho_zero, &
1151  delta_rho
1152 
1153  CHARACTER(LEN=*), PARAMETER :: routinen = 'surface_fattebert_gygi'
1154  REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1155 
1156  INTEGER :: handle, i, j, k, l
1157  INTEGER, DIMENSION(3) :: lb, ub
1158  REAL(kind=dp) :: e, eps_elec, f, p, q, rho, s, t, twobeta
1159  REAL(kind=dp), DIMENSION(2) :: theta
1160 
1161  CALL timeset(routinen, handle)
1162 
1163  e = epsilon_solvent - 1.0_dp
1164  f = 0.5_dp*e
1165  q = 1.0_dp/rho_zero
1166  twobeta = 2.0_dp*beta
1167 
1168  lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1169  ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1170 
1171  ! Calculate finite difference between two isosurfaces
1172 !$OMP PARALLEL DO DEFAULT(NONE) &
1173 !$OMP PRIVATE(eps_elec,i,j,k,l,p,rho,s,t,theta) &
1174 !$OMP SHARED(delta_rho,dtheta,e,epsilon_solvent,f,lb) &
1175 !$OMP SHARED(norm_drho_elec,q,rho_elec,twobeta,ub)
1176  DO k = lb(3), ub(3)
1177  DO j = lb(2), ub(2)
1178  DO i = lb(1), ub(1)
1179  DO l = 1, 2
1180  rho = rho_elec%array(i, j, k) + (real(l, kind=dp) - 1.5_dp)*delta_rho
1181  IF (rho < rhotol) THEN
1182  eps_elec = epsilon_solvent
1183  ELSE
1184  s = rho*q
1185  p = s**twobeta
1186  t = 1.0_dp/(1.0_dp + p)
1187  eps_elec = 1.0_dp + f*(1.0_dp + (1.0_dp - p)*t)
1188  END IF
1189  theta(l) = (epsilon_solvent - eps_elec)/e
1190  END DO
1191  dtheta%array(i, j, k) = (theta(2) - theta(1))*norm_drho_elec%array(i, j, k)/delta_rho
1192  END DO
1193  END DO
1194  END DO
1195 !$OMP END PARALLEL DO
1196 
1197  CALL timestop(handle)
1198 
1199  END SUBROUTINE surface_fattebert_gygi
1200 
1201 ! **************************************************************************************************
1202 !> \brief Print SCCS results
1203 !> \param energy ...
1204 !> \param sccs_control ...
1205 !> \param output_unit ...
1206 !> \par History:
1207 !> - Creation (11.11.2022,MK)
1208 !> \author Matthias Krack (MK)
1209 !> \version 1.0
1210 ! **************************************************************************************************
1211  SUBROUTINE print_sccs_results(energy, sccs_control, output_unit)
1212 
1213  TYPE(qs_energy_type), POINTER :: energy
1214  TYPE(sccs_control_type), POINTER :: sccs_control
1215  INTEGER, INTENT(IN) :: output_unit
1216 
1217  IF (output_unit > 0) THEN
1218  cpassert(ASSOCIATED(energy))
1219  cpassert(ASSOCIATED(sccs_control))
1220  WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1221  "SCCS| Hartree energy of solute and solvent [Hartree]", energy%sccs_hartree, &
1222  "SCCS| Hartree energy of the solute only [Hartree]", energy%hartree
1223  WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14,/,T3,A,T61,F20.3)") &
1224  "SCCS| Polarisation energy [Hartree]", energy%sccs_pol, &
1225  "SCCS| [kcal/mol]", &
1226  cp_unit_from_cp2k(energy%sccs_pol, "kcalmol"), &
1227  "SCCS| Cavitation energy [Hartree]", energy%sccs_cav, &
1228  "SCCS| [kcal/mol]", &
1229  cp_unit_from_cp2k(energy%sccs_cav, "kcalmol"), &
1230  "SCCS| Dispersion free energy [Hartree]", energy%sccs_dis, &
1231  "SCCS| [kcal/mol]", &
1232  cp_unit_from_cp2k(energy%sccs_dis, "kcalmol"), &
1233  "SCCS| Repulsion free energy [Hartree]", energy%sccs_rep, &
1234  "SCCS| [kcal/mol]", &
1235  cp_unit_from_cp2k(energy%sccs_rep, "kcalmol"), &
1236  "SCCS| Solvation free energy [Hartree]", energy%sccs_sol, &
1237  "SCCS| [kcal/mol]", &
1238  cp_unit_from_cp2k(energy%sccs_sol, "kcalmol")
1239  END IF
1240 
1241  END SUBROUTINE print_sccs_results
1242 
1243 END MODULE qs_sccs
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
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::low_print_level
integer, parameter, public low_print_level
Definition: cp_output_handling.F:73
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
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition: cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
Definition: cp_realspace_grid_cube.F:45
cp_realspace_grid_init
Definition: cp_realspace_grid_init.F:15
cp_realspace_grid_init::init_input_type
subroutine, public init_input_type(input_settings, nsmax, rs_grid_section, ilevel, higher_grid_layout)
parses an input section to assign the proper values to the input type
Definition: cp_realspace_grid_init.F:48
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition: cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition: cp_subsys_types.F:345
cp_units
unit conversion facility
Definition: cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition: cp_units.F:1179
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::sccs_derivative_cd5
integer, parameter, public sccs_derivative_cd5
Definition: input_constants.F:1169
input_constants::sccs_fattebert_gygi
integer, parameter, public sccs_fattebert_gygi
Definition: input_constants.F:1165
input_constants::sccs_derivative_cd7
integer, parameter, public sccs_derivative_cd7
Definition: input_constants.F:1169
input_constants::sccs_derivative_fft
integer, parameter, public sccs_derivative_fft
Definition: input_constants.F:1169
input_constants::sccs_derivative_cd3
integer, parameter, public sccs_derivative_cd3
Definition: input_constants.F:1169
input_constants::sccs_andreussi
integer, parameter, public sccs_andreussi
Definition: input_constants.F:1165
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
input_section_types::section_get_ivals
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
Definition: input_section_types.F:989
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition: input_section_types.F:1005
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition: kinds.F:54
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition: mathconstants.F:113
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_methods
Definition: pw_methods.F:25
pw_methods::pw_derive
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
Definition: pw_methods.F:10106
pw_poisson_methods
Definition: pw_poisson_methods.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition: pw_poisson_types.F:22
pw_poisson_types::pw_poisson_mt
integer, parameter, public pw_poisson_mt
Definition: pw_poisson_types.F:66
pw_poisson_types::pw_poisson_analytic
integer, parameter, public pw_poisson_analytic
Definition: pw_poisson_types.F:66
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
qs_sccs
Self-consistent continuum solvation (SCCS) model implementation.
Definition: qs_sccs.F:29
qs_sccs::sccs
subroutine, public sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs, h_stress)
Self-consistent continuum solvation (SCCS) model implementation.
Definition: qs_sccs.F:122
qs_sccs::print_sccs_results
subroutine, public print_sccs_results(energy, sccs_control, output_unit)
Print SCCS results.
Definition: qs_sccs.F:1212
qs_scf_types
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
realspace_grid_types
Definition: realspace_grid_types.F:18
realspace_grid_types::rs_grid_create
subroutine, public rs_grid_create(rs, desc)
...
Definition: realspace_grid_types.F:482
realspace_grid_types::rs_grid_create_descriptor
subroutine, public rs_grid_create_descriptor(desc, pw_grid, input_settings, border_points)
Determine the setup of real space grids - this is divided up into the creation of a descriptor and th...
Definition: realspace_grid_types.F:203
realspace_grid_types::rs_grid_release_descriptor
subroutine, public rs_grid_release_descriptor(rs_desc)
releases the given rs grid descriptor (see doc/ReferenceCounting.html)
Definition: realspace_grid_types.F:2036
realspace_grid_types::rs_grid_release
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
Definition: realspace_grid_types.F:2016
rs_methods
numerical operations on real-space grid
Definition: rs_methods.F:14
rs_methods::derive_fdm_cd7
subroutine, public derive_fdm_cd7(f, df, rs_grid)
6th order finite difference derivative of a function on realspace grid
Definition: rs_methods.F:198
rs_methods::derive_fdm_cd5
subroutine, public derive_fdm_cd5(f, df, rs_grid)
4th order finite difference derivative of a function on realspace grid
Definition: rs_methods.F:129
rs_methods::derive_fdm_cd3
subroutine, public derive_fdm_cd3(f, df, rs_grid)
2nd order finite difference derivative of a function on realspace grid
Definition: rs_methods.F:60