(git:21ef868)
Loading...
Searching...
No Matches
qs_sccs.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 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
29MODULE 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 sccs_saa_andreussi
52 USE input_section_types, ONLY: section_get_ivals,&
53 section_get_lval,&
54 section_vals_get_subs_vals,&
55 section_vals_type
56 USE kinds, ONLY: default_path_length,&
57 default_string_length,&
58 dp,&
59 int_8
60 USE mathconstants, ONLY: fourpi,&
61 pi,&
62 twopi
63 USE message_passing, ONLY: mp_para_env_type
64 USE particle_list_types, ONLY: particle_list_type
65 USE pw_env_types, ONLY: pw_env_get,&
66 pw_env_type
67 USE pw_methods, ONLY: pw_axpy,&
68 pw_copy,&
69 pw_derive,&
70 pw_integral_ab,&
71 pw_integrate_function,&
72 pw_scale,&
73 pw_transfer,&
74 pw_zero
75 USE pw_poisson_methods, ONLY: pw_func_u_convolution,&
76 pw_poisson_solve
77 USE pw_poisson_types, ONLY: pw_poisson_analytic,&
78 pw_poisson_mt,&
79 pw_poisson_type
80 USE pw_pool_types, ONLY: pw_pool_p_type,&
81 pw_pool_type
82 USE pw_types, ONLY: pw_c1d_gs_type,&
83 pw_r3d_rs_type
84 USE qs_energy_types, ONLY: qs_energy_type
85 USE qs_environment_types, ONLY: get_qs_env,&
86 qs_environment_type
87 USE qs_rho_types, ONLY: qs_rho_get,&
88 qs_rho_type
89 USE qs_scf_types, ONLY: qs_scf_env_type
90 USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
91 realspace_grid_input_type,&
92 realspace_grid_type,&
93 rs_grid_create,&
94 rs_grid_create_descriptor,&
95 rs_grid_release,&
96 rs_grid_release_descriptor
97 USE rs_methods, ONLY: derive_fdm_cd3,&
98 derive_fdm_cd5,&
99 derive_fdm_cd7
100#include "./base/base_uses.f90"
101
102 IMPLICIT NONE
103
104 PRIVATE
105
106 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_sccs'
107
108 PUBLIC :: print_sccs_results, sccs
109
110CONTAINS
111
112! **************************************************************************************************
113!> \brief Self-consistent continuum solvation (SCCS) model implementation
114!> \param qs_env ...
115!> \param rho_tot_gspace ...
116!> \param v_hartree_gspace ...
117!> \param v_sccs ...
118!> \param h_stress ...
119!> \par History:
120!> - Creation (10.10.2013,MK)
121!> \author Matthias Krack (MK)
122!> \version 1.0
123! **************************************************************************************************
124
125 SUBROUTINE sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs, h_stress)
126
127 TYPE(qs_environment_type), POINTER :: qs_env
128 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_tot_gspace, v_hartree_gspace
129 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_sccs
130 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
131 OPTIONAL :: h_stress
132
133 CHARACTER(LEN=*), PARAMETER :: routinen = 'sccs'
134 REAL(kind=dp), PARAMETER :: epstol = 1.0e-8_dp
135
136 CHARACTER(LEN=4*default_string_length) :: message, my_pos_cube
137 CHARACTER(LEN=default_path_length) :: cube_path, filename, mpi_filename, &
138 print_path
139 INTEGER :: cube_unit, handle, i, ispin, iter, j, k, &
140 nspin, output_unit, print_level
141 INTEGER(KIND=int_8) :: ngpts
142 INTEGER, DIMENSION(3) :: lb, ub
143 LOGICAL :: append_cube, calculate_stress_tensor, &
144 do_kpoints, mpi_io, should_output
145 REAL(kind=dp) :: alpha_zeta, cavity_surface, cavity_volume, cell_volume, delta_eta, &
146 delta_zeta, dphi2, dvol, e_tot, epsilon_solvent, f, f0, polarisation_charge, r, r_solv, &
147 rho_delta, rho_delta_avg, rho_delta_max, rho_iter_new, tot_rho_elec, tot_rho_solute
148 REAL(kind=dp), DIMENSION(3) :: abc, lxyz, uxyz
149 TYPE(cp_logger_type), POINTER :: logger
150 TYPE(cp_subsys_type), POINTER :: cp_subsys
151 TYPE(dft_control_type), POINTER :: dft_control
152 TYPE(mp_para_env_type), POINTER :: para_env
153 TYPE(particle_list_type), POINTER :: particles
154 TYPE(pw_env_type), POINTER :: pw_env
155 TYPE(pw_poisson_type), POINTER :: poisson_env
156 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
157 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
158 TYPE(pw_r3d_rs_type) :: d_s_rhoel, deps_elec, dtf, eps_elec, ff, &
159 p_e_interface, s, t, tem_convolution, &
160 tem_func, u
161 TYPE(pw_r3d_rs_type), DIMENSION(3) :: dln_eps_elec, dphi_tot
162 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_pw_r
163 TYPE(pw_r3d_rs_type), POINTER :: rho_pw_r_sccs
164 TYPE(qs_energy_type), POINTER :: energy
165 TYPE(qs_rho_type), POINTER :: rho
166 TYPE(qs_scf_env_type), POINTER :: scf_env
167 TYPE(sccs_control_type), POINTER :: sccs_control
168 TYPE(section_vals_type), POINTER :: input
169
170 CALL timeset(routinen, handle)
171
172 NULLIFY (auxbas_pw_pool)
173 NULLIFY (cp_subsys)
174 NULLIFY (dft_control)
175 NULLIFY (energy)
176 NULLIFY (input)
177 NULLIFY (logger)
178 NULLIFY (para_env)
179 NULLIFY (particles)
180 NULLIFY (poisson_env)
181 NULLIFY (pw_env)
182 NULLIFY (pw_pools)
183 NULLIFY (rho)
184 NULLIFY (sccs_control)
185 NULLIFY (scf_env)
186
187 ! Load data from Quickstep environment
188 CALL get_qs_env(qs_env=qs_env, &
189 cp_subsys=cp_subsys, &
190 do_kpoints=do_kpoints, &
191 dft_control=dft_control, &
192 energy=energy, &
193 input=input, &
194 para_env=para_env, &
195 pw_env=pw_env, &
196 rho=rho, &
197 scf_env=scf_env)
198 CALL cp_subsys_get(cp_subsys, particles=particles)
199
200 sccs_control => dft_control%sccs_control
201
202 cpassert(ASSOCIATED(qs_env))
203
204 IF (do_kpoints) THEN
205 cpwarn("SCCS with k-points has not yet been fully validated")
206 END IF
207
208 IF (PRESENT(h_stress)) THEN
209 calculate_stress_tensor = .true.
210 h_stress(:, :) = 0.0_dp
211 cpwarn("The stress tensor for SCCS has not yet been fully validated")
212 ELSE
213 calculate_stress_tensor = .false.
214 END IF
215
216 ! Get access to the PW grid pool
217 CALL pw_env_get(pw_env, &
218 auxbas_pw_pool=auxbas_pw_pool, &
219 pw_pools=pw_pools, &
220 poisson_env=poisson_env)
221
222 CALL pw_zero(v_sccs)
223
224 ! Calculate no SCCS contribution, if the requested SCF convergence threshold is not reached yet
225 IF (.NOT. sccs_control%sccs_activated) THEN
226 IF (sccs_control%eps_scf > 0.0_dp) THEN
227 IF ((scf_env%iter_delta > sccs_control%eps_scf) .OR. &
228 ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
229 (qs_env%scf_env%iter_count <= 1))) THEN
230 IF (calculate_stress_tensor) THEN
231 ! Request also the calculation of the stress tensor contribution
232 CALL pw_poisson_solve(poisson_env=poisson_env, &
233 density=rho_tot_gspace, &
234 ehartree=energy%hartree, &
235 vhartree=v_hartree_gspace, &
236 h_stress=h_stress)
237 ELSE
238 CALL pw_poisson_solve(poisson_env=poisson_env, &
239 density=rho_tot_gspace, &
240 ehartree=energy%hartree, &
241 vhartree=v_hartree_gspace)
242 END IF
243 energy%sccs_pol = 0.0_dp
244 energy%sccs_cav = 0.0_dp
245 energy%sccs_dis = 0.0_dp
246 energy%sccs_rep = 0.0_dp
247 energy%sccs_sol = 0.0_dp
248 energy%sccs_hartree = energy%hartree
249 CALL timestop(handle)
250 RETURN
251 END IF
252 END IF
253 sccs_control%sccs_activated = .true.
254 END IF
255
256 nspin = dft_control%nspins
257
258 ! Manage print output control
259 logger => cp_get_default_logger()
260 print_level = logger%iter_info%print_level
261 print_path = "DFT%PRINT%SCCS"
262 should_output = (btest(cp_print_key_should_output(logger%iter_info, input, &
263 trim(print_path)), cp_p_file))
264 output_unit = cp_print_key_unit_nr(logger, input, trim(print_path), &
265 extension=".sccs", &
266 ignore_should_output=should_output, &
267 log_filename=.false.)
268
269 ! Get rho
270 CALL qs_rho_get(rho_struct=rho, &
271 rho_r=rho_pw_r, &
272 rho_r_sccs=rho_pw_r_sccs)
273
274 ! Retrieve the last rho_iter from the previous SCCS cycle if available
275 cpassert(ASSOCIATED(rho_pw_r_sccs))
276
277 ! Retrieve the total electronic density in r-space
278 block
279 TYPE(pw_r3d_rs_type) :: rho_elec
280 CALL auxbas_pw_pool%create_pw(rho_elec)
281
282 ! Retrieve grid parameters
283 ngpts = rho_elec%pw_grid%ngpts
284 dvol = rho_elec%pw_grid%dvol
285 cell_volume = rho_elec%pw_grid%vol
286 abc(1:3) = real(rho_elec%pw_grid%npts(1:3), kind=dp)*rho_elec%pw_grid%dr(1:3)
287 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
288 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
289
290 CALL pw_copy(rho_pw_r(1), rho_elec)
291 DO ispin = 2, nspin
292 CALL pw_axpy(rho_pw_r(ispin), rho_elec)
293 END DO
294 tot_rho_elec = pw_integrate_function(rho_elec)
295
296 ! Calculate the dielectric (smoothed) function of rho_elec in r-space
297 CALL auxbas_pw_pool%create_pw(eps_elec)
298 CALL auxbas_pw_pool%create_pw(deps_elec)
299 IF (sccs_control%method_id == sccs_saa_andreussi) THEN
300 CALL auxbas_pw_pool%create_pw(s)
301 CALL auxbas_pw_pool%create_pw(d_s_rhoel)
302 CALL auxbas_pw_pool%create_pw(u)
303 CALL auxbas_pw_pool%create_pw(ff)
304 CALL auxbas_pw_pool%create_pw(t)
305 CALL auxbas_pw_pool%create_pw(dtf)
306 f0 = sccs_control%f0
307 delta_eta = sccs_control%delta_eta
308 alpha_zeta = sccs_control%alpha_zeta
309 r_solv = sccs_control%R_solv
310 delta_zeta = sccs_control%delta_zeta
311 END IF
312
313 ! Relative permittivity or dielectric constant of the solvent (medium)
314 epsilon_solvent = sccs_control%epsilon_solvent
315 SELECT CASE (sccs_control%method_id)
316 CASE (sccs_andreussi)
317 CALL andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%rho_max, &
318 sccs_control%rho_min)
319 CASE (sccs_fattebert_gygi)
320 CALL fattebert_gygi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%beta, &
321 sccs_control%rho_zero)
322 CASE (sccs_saa_andreussi)
323 CALL sa_andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%rho_max, &
324 sccs_control%rho_min, s, d_s_rhoel)
325 CASE DEFAULT
326 cpabort("Invalid method specified for SCCS model")
327 END SELECT
328
329 IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
330 ! Optional output of the dielectric function in cube file format
331 filename = "DIELECTRIC_FUNCTION"
332 cube_path = trim(print_path)//"%"//trim(filename)
333 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
334 cp_p_file)) THEN
335 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
336 my_pos_cube = "REWIND"
337 IF (append_cube) my_pos_cube = "APPEND"
338 mpi_io = .true.
339 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
340 extension=".cube", middle_name=trim(filename), &
341 file_position=my_pos_cube, log_filename=.false., &
342 mpi_io=mpi_io, fout=mpi_filename)
343 IF (output_unit > 0) THEN
344 IF (.NOT. mpi_io) THEN
345 INQUIRE (unit=cube_unit, name=filename)
346 ELSE
347 filename = mpi_filename
348 END IF
349 WRITE (unit=output_unit, fmt="(T3,A)") &
350 "SCCS| The dielectric function is written in cube file format to the file:", &
351 "SCCS| "//trim(filename)
352 END IF
353 CALL cp_pw_to_cube(eps_elec, cube_unit, trim(filename), particles=particles, &
354 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), &
355 mpi_io=mpi_io)
356 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
357 END IF
358 END IF
359
360 ! Calculate the (quantum) volume and surface of the solute cavity
361 cavity_surface = 0.0_dp
362 cavity_volume = 0.0_dp
363
364 IF (abs(epsilon_solvent - 1.0_dp) > epstol) THEN
365
366 block
367 TYPE(pw_r3d_rs_type) :: theta, norm_drho_elec
368 TYPE(pw_r3d_rs_type), DIMENSION(3) :: drho_elec
369 CALL auxbas_pw_pool%create_pw(theta)
370 CALL pw_zero(theta)
371
372 ! Calculate the (quantum) volume of the solute cavity
373 f = 1.0_dp/(epsilon_solvent - 1.0_dp)
374!$OMP PARALLEL DO DEFAULT(NONE) &
375!$OMP PRIVATE(i,j,k) &
376!$OMP SHARED(epsilon_solvent,eps_elec,f,lb,theta,ub)
377 DO k = lb(3), ub(3)
378 DO j = lb(2), ub(2)
379 DO i = lb(1), ub(1)
380 theta%array(i, j, k) = f*(epsilon_solvent - eps_elec%array(i, j, k))
381 END DO
382 END DO
383 END DO
384!$OMP END PARALLEL DO
385 cavity_volume = pw_integrate_function(theta)
386
387 ! Calculate the derivative of the electronic density in r-space
388 ! TODO: Could be retrieved from the QS environment
389 DO i = 1, 3
390 CALL auxbas_pw_pool%create_pw(drho_elec(i))
391 END DO
392 CALL derive(rho_elec, drho_elec, sccs_derivative_fft, pw_env, input)
393
394 CALL auxbas_pw_pool%create_pw(norm_drho_elec)
395
396 ! Calculate the norm of the gradient of the electronic density in r-space
397!$OMP PARALLEL DO DEFAULT(NONE) &
398!$OMP PRIVATE(i,j,k) &
399!$OMP SHARED(drho_elec,lb,norm_drho_elec,ub)
400 DO k = lb(3), ub(3)
401 DO j = lb(2), ub(2)
402 DO i = lb(1), ub(1)
403 norm_drho_elec%array(i, j, k) = sqrt(drho_elec(1)%array(i, j, k)* &
404 drho_elec(1)%array(i, j, k) + &
405 drho_elec(2)%array(i, j, k)* &
406 drho_elec(2)%array(i, j, k) + &
407 drho_elec(3)%array(i, j, k)* &
408 drho_elec(3)%array(i, j, k))
409 END DO
410 END DO
411 END DO
412!$OMP END PARALLEL DO
413
414 ! Optional output of the norm of the density gradient in cube file format
415 filename = "DENSITY_GRADIENT"
416 cube_path = trim(print_path)//"%"//trim(filename)
417 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
418 cp_p_file)) THEN
419 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
420 my_pos_cube = "REWIND"
421 IF (append_cube) my_pos_cube = "APPEND"
422 mpi_io = .true.
423 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
424 extension=".cube", middle_name=trim(filename), &
425 file_position=my_pos_cube, log_filename=.false., &
426 mpi_io=mpi_io, fout=mpi_filename)
427 IF (output_unit > 0) THEN
428 IF (.NOT. mpi_io) THEN
429 INQUIRE (unit=cube_unit, name=filename)
430 ELSE
431 filename = mpi_filename
432 END IF
433 WRITE (unit=output_unit, fmt="(T3,A)") &
434 "SCCS| The norm of the density gradient is written in cube file format to the file:", &
435 "SCCS| "//trim(filename)
436 END IF
437 CALL cp_pw_to_cube(norm_drho_elec, cube_unit, trim(filename), particles=particles, &
438 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), &
439 mpi_io=mpi_io)
440 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
441 END IF
442
443 ! Calculate the (quantum) surface of the solute cavity
444 SELECT CASE (sccs_control%method_id)
445 CASE (sccs_andreussi)
446 CALL surface_andreussi(rho_elec, norm_drho_elec, theta, epsilon_solvent, &
447 sccs_control%rho_max, sccs_control%rho_min, &
448 sccs_control%delta_rho)
449 CASE (sccs_fattebert_gygi)
450 CALL surface_fattebert_gygi(rho_elec, norm_drho_elec, theta, epsilon_solvent, &
451 sccs_control%beta, sccs_control%rho_zero, &
452 sccs_control%delta_rho)
453 CASE (sccs_saa_andreussi)
454 CASE DEFAULT
455 cpabort("Invalid method specified for SCCS model")
456 END SELECT
457 cavity_surface = pw_integrate_function(theta)
458
459 ! Release storage
460 CALL auxbas_pw_pool%give_back_pw(theta)
461 CALL auxbas_pw_pool%give_back_pw(norm_drho_elec)
462 DO i = 1, 3
463 CALL auxbas_pw_pool%give_back_pw(drho_elec(i))
464 END DO
465 END block
466
467 END IF ! epsilon_solvent > 1
468
469 CALL auxbas_pw_pool%give_back_pw(rho_elec)
470 END block
471
472 block
473 TYPE(pw_r3d_rs_type) :: rho_tot, phi_tot, rho_solute, rho_tot_zero
474 ! Retrieve the total charge density (core + elec) of the solute in r-space
475 CALL auxbas_pw_pool%create_pw(rho_solute)
476 CALL pw_zero(rho_solute)
477 CALL pw_transfer(rho_tot_gspace, rho_solute)
478 tot_rho_solute = pw_integrate_function(rho_solute)
479
480 ! Check total charge
481 IF (abs(tot_rho_solute) >= 1.0e-6_dp) THEN
482 IF ((poisson_env%parameters%solver /= pw_poisson_analytic) .AND. &
483 (poisson_env%parameters%solver /= pw_poisson_mt)) THEN
484 WRITE (unit=message, fmt="(A,SP,F0.6,A)") &
485 "The system (solute) has a non-negligible charge of ", -tot_rho_solute, &
486 ". It is recommended to use non-periodic boundary conditions (PERIODIC none) "// &
487 "combined with an appropriate Poisson solver (POISSON_SOLVER MT or analytic)"
488 cpwarn(message)
489 END IF
490 END IF
491
492 ! Reassign work storage to rho_tot_zero, because rho_elec is no longer needed
493 CALL auxbas_pw_pool%create_pw(rho_tot_zero)
494
495 IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
496
497 ! Build the initial (rho_iter = 0) total charge density (solute plus polarisation) in r-space
498 ! eps_elec <- ln(eps_elec)
499!$OMP PARALLEL DO DEFAULT(NONE) &
500!$OMP PRIVATE(i,j,k) &
501!$OMP SHARED(eps_elec,lb,message,output_unit,para_env,ub) &
502!$OMP SHARED(rho_solute,rho_tot_zero)
503 DO k = lb(3), ub(3)
504 DO j = lb(2), ub(2)
505 DO i = lb(1), ub(1)
506 IF (eps_elec%array(i, j, k) < 1.0_dp) THEN
507 WRITE (unit=message, fmt="(A,ES12.3,A,3(I0,A))") &
508 "SCCS| Invalid dielectric function value ", eps_elec%array(i, j, k), &
509 " encountered at grid point (", i, ",", j, ",", k, ")"
510 cpabort(message)
511 END IF
512 rho_tot_zero%array(i, j, k) = rho_solute%array(i, j, k)/eps_elec%array(i, j, k)
513 eps_elec%array(i, j, k) = log(eps_elec%array(i, j, k))
514 END DO
515 END DO
516 END DO
517!$OMP END PARALLEL DO
518
519 ELSE
520
521 lxyz(:) = real(u%pw_grid%bounds(1, :), dp)
522 uxyz(:) = real(u%pw_grid%bounds(2, :), dp)
523!$OMP PARALLEL DO DEFAULT(NONE) &
524!$OMP PRIVATE(i,j,k,r) &
525!$OMP SHARED(lb,message,output_unit,para_env,ub,dtf,uxyz,lxyz) &
526!$OMP SHARED(u, alpha_zeta, R_solv, delta_zeta, eps_elec)
527 DO k = lb(3), ub(3)
528 DO j = lb(2), ub(2)
529 DO i = lb(1), ub(1)
530 IF (eps_elec%array(i, j, k) < 1.0_dp) THEN
531 WRITE (unit=message, fmt="(A,ES12.3,A,3(I0,A))") &
532 "SCCS| Invalid dielectric function value ", eps_elec%array(i, j, k), &
533 " encountered at grid point (", i, ",", j, ",", k, ")"
534 cpabort(message)
535 END IF
536 IF (i > -1) r = ((real(uxyz(1), dp) - real(i, dp) + 0.5_dp)*u%pw_grid%dr(1))**2.0_dp
537 IF (i < 0) r = ((real(i, dp) - real(lxyz(1), dp) + 0.5_dp)*u%pw_grid%dr(1))**2.0_dp
538 IF (j > -1) r = ((real(uxyz(2), dp) - real(j, dp) + 0.5_dp)*u%pw_grid%dr(2))**2.0_dp + r
539 IF (j < 0) r = ((real(j, dp) - real(lxyz(2), dp) + 0.5_dp)*u%pw_grid%dr(2))**2.0_dp + r
540 IF (k > -1) r = ((real(uxyz(3), dp) - real(k, dp) + 0.5_dp)*u%pw_grid%dr(3))**2.0_dp + r
541 IF (k < 0) r = ((real(k, dp) - real(lxyz(3), dp) + 0.5_dp)*u%pw_grid%dr(3))**2.0_dp + r
542 r = sqrt(r)
543 u%array(i, j, k) = (1.0_dp/2.0_dp)*erfc((r - alpha_zeta*r_solv)/delta_zeta)
544 dtf%array(i, j, k) = 1.0_dp
545 END DO
546 END DO
547 END DO
548!$OMP END PARALLEL DO
549
550 CALL pw_func_u_convolution(poisson_env=poisson_env, func=dtf, convolution=ff, u=u)
551 CALL pw_scale(u, 1.0_dp/ff%array(0, 0, 0))
552 CALL pw_func_u_convolution(poisson_env=poisson_env, func=s, convolution=ff, u=u)
553
554!$OMP PARALLEL DO DEFAULT(NONE) &
555!$OMP PRIVATE(i,j,k) &
556!$OMP SHARED(lb,para_env,ub, delta_eta, f0, epsilon_solvent) &
557!$OMP SHARED(ff, t, dtf, eps_elec, rho_tot_zero, rho_solute, s)
558 DO k = lb(3), ub(3)
559 DO j = lb(2), ub(2)
560 DO i = lb(1), ub(1)
561 t%array(i, j, k) = 0.5_dp*(1.0_dp + erf((ff%array(i, j, k) - f0)/delta_eta))
562 dtf%array(i, j, k) = (1.0_dp/sqrt(pi))*exp(-(ff%array(i, j, k) - f0)**2.0/(delta_eta**2.0))
563 eps_elec%array(i, j, k) = exp(log(epsilon_solvent)*(1 - (s%array(i, j, k) &
564 + (1 - s%array(i, j, k))*t%array(i, j, k))))
565 rho_tot_zero%array(i, j, k) = rho_solute%array(i, j, k)/eps_elec%array(i, j, k)
566! eps_elec%array(i, j, k) = LOG(eps_elec%array(i, j, k))
567 END DO
568 END DO
569 END DO
570!$OMP END PARALLEL DO
571
572 ! Optional output of the dielectric function in cube file format (solvent aware algorithm case)
573 filename = "DIELECTRIC_FUNCTION"
574 cube_path = trim(print_path)//"%"//trim(filename)
575 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
576 cp_p_file)) THEN
577 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
578 my_pos_cube = "REWIND"
579 IF (append_cube) my_pos_cube = "APPEND"
580 mpi_io = .true.
581 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
582 extension=".cube", middle_name=trim(filename), &
583 file_position=my_pos_cube, log_filename=.false., &
584 mpi_io=mpi_io, fout=mpi_filename)
585 IF (output_unit > 0) THEN
586 IF (.NOT. mpi_io) THEN
587 INQUIRE (unit=cube_unit, name=filename)
588 ELSE
589 filename = mpi_filename
590 END IF
591 WRITE (unit=output_unit, fmt="(T3,A)") &
592 "SCCS| The dielectric function is written in cube file format to the file:", &
593 "SCCS| "//trim(filename)
594 END IF
595 CALL cp_pw_to_cube(eps_elec, cube_unit, trim(filename), particles=particles, &
596 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), &
597 mpi_io=mpi_io)
598 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
599 END IF
600
601!$OMP PARALLEL DO DEFAULT(NONE) &
602!$OMP PRIVATE(i,j,k) &
603!$OMP SHARED(lb,para_env,ub, delta_eta, f0, epsilon_solvent) &
604!$OMP SHARED(ff, t, dtf, eps_elec, rho_tot_zero, rho_solute, s)
605 DO k = lb(3), ub(3)
606 DO j = lb(2), ub(2)
607 DO i = lb(1), ub(1)
608 eps_elec%array(i, j, k) = log(eps_elec%array(i, j, k))
609 END DO
610 END DO
611 END DO
612!$OMP END PARALLEL DO
613
614 END IF
615
616 ! Build the derivative of LOG(eps_elec)
617 DO i = 1, 3
618 CALL auxbas_pw_pool%create_pw(dln_eps_elec(i))
619 CALL pw_zero(dln_eps_elec(i))
620 END DO
621 CALL derive(eps_elec, dln_eps_elec, sccs_control%derivative_method, pw_env, input)
622 IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
623 CALL auxbas_pw_pool%give_back_pw(eps_elec)
624 END IF
625
626 ! Print header for the SCCS cycle
627 IF (should_output .AND. (output_unit > 0)) THEN
628 IF (print_level > low_print_level) THEN
629 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.12)") &
630 "SCCS| Total electronic charge density ", -tot_rho_elec, &
631 "SCCS| Total charge density (solute) ", -tot_rho_solute
632 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.3)") &
633 "SCCS| Volume of the cell [bohr^3]", cell_volume, &
634 "SCCS| [angstrom^3]", &
635 cp_unit_from_cp2k(cell_volume, "angstrom^3")
636 IF (abs(epsilon_solvent - 1.0_dp) > epstol) THEN
637 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.3)") &
638 "SCCS| Volume of the solute cavity [bohr^3]", cavity_volume, &
639 "SCCS| [angstrom^3]", &
640 cp_unit_from_cp2k(cavity_volume, "angstrom^3"), &
641 "SCCS| Surface of the solute cavity [bohr^2]", cavity_surface, &
642 "SCCS| [angstrom^2]", &
643 cp_unit_from_cp2k(cavity_surface, "angstrom^2")
644 END IF
645 WRITE (unit=output_unit, fmt="(T3,A)") &
646 "SCCS|", &
647 "SCCS| Step Average residual Maximum residual E(Hartree) [a.u.]"
648 END IF
649 END IF
650
651 ! Get storage for the derivative of the total potential (field) in r-space
652 DO i = 1, 3
653 CALL auxbas_pw_pool%create_pw(dphi_tot(i))
654 END DO
655
656 ! Initialise the total charge density in r-space rho_tot with rho_tot_zero + rho_iter_zero
657 CALL auxbas_pw_pool%create_pw(rho_tot)
658 CALL pw_copy(rho_tot_zero, rho_tot)
659 CALL pw_axpy(rho_pw_r_sccs, rho_tot)
660
661 CALL auxbas_pw_pool%create_pw(phi_tot)
662 CALL pw_zero(phi_tot)
663
664 IF (sccs_control%method_id == sccs_saa_andreussi) THEN
665 CALL auxbas_pw_pool%create_pw(p_e_interface)
666 CALL pw_zero(p_e_interface)
667 END IF
668
669 ! Main SCCS iteration loop
670 iter = 0
671
672 iter_loop: DO
673
674 ! Increment iteration counter
675 iter = iter + 1
676
677 ! Check if the requested maximum number of SCCS iterations is reached
678 IF (iter > sccs_control%max_iter) THEN
679 IF (output_unit > 0) THEN
680 WRITE (unit=output_unit, fmt="(T3,A,/,T3,A,I0,A)") &
681 "SCCS| Maximum number of SCCS iterations reached", &
682 "SCCS| Iteration cycle did not converge in ", sccs_control%max_iter, " steps"
683 ELSE
684 WRITE (unit=message, fmt="(A,I0,A)") &
685 "The SCCS iteration cycle did not converge in ", sccs_control%max_iter, " steps"
686 cpwarn(message)
687 END IF
688 EXIT iter_loop
689 END IF
690
691 ! Calculate derivative of the current total potential in r-space
692 CALL pw_poisson_solve(poisson_env=poisson_env, &
693 density=rho_tot, &
694 vhartree=phi_tot, &
695 dvhartree=dphi_tot)
696 energy%sccs_hartree = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
697
698 ! Update total charge density (solute plus polarisation) in r-space
699 ! based on the iterated polarisation charge density
700 f = 1.0_dp/fourpi
701 rho_delta_avg = 0.0_dp
702 rho_delta_max = 0.0_dp
703!$OMP PARALLEL DO DEFAULT(NONE) &
704!$OMP PRIVATE(i,j,k,rho_delta,rho_iter_new) &
705!$OMP SHARED(dln_eps_elec,dphi_tot,f,lb,rho_pw_r_sccs,ub) &
706!$OMP SHARED(rho_tot,rho_tot_zero,sccs_control) &
707!$OMP REDUCTION(+:rho_delta_avg) &
708!$OMP REDUCTION(MAX:rho_delta_max)
709 DO k = lb(3), ub(3)
710 DO j = lb(2), ub(2)
711 DO i = lb(1), ub(1)
712 rho_iter_new = (dln_eps_elec(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + &
713 dln_eps_elec(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + &
714 dln_eps_elec(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k))*f
715 rho_iter_new = rho_pw_r_sccs%array(i, j, k) + &
716 sccs_control%mixing*(rho_iter_new - rho_pw_r_sccs%array(i, j, k))
717 rho_delta = abs(rho_iter_new - rho_pw_r_sccs%array(i, j, k))
718 rho_delta_max = max(rho_delta, rho_delta_max)
719 rho_delta_avg = rho_delta_avg + rho_delta
720 rho_tot%array(i, j, k) = rho_tot_zero%array(i, j, k) + rho_iter_new
721 rho_pw_r_sccs%array(i, j, k) = rho_iter_new
722 END DO
723 END DO
724 END DO
725!$OMP END PARALLEL DO
726
727 CALL para_env%sum(rho_delta_avg)
728 rho_delta_avg = rho_delta_avg/real(ngpts, kind=dp)
729 CALL para_env%max(rho_delta_max)
730
731 IF (should_output .AND. (output_unit > 0)) THEN
732 IF (print_level > low_print_level) THEN
733 IF ((abs(rho_delta_avg) < 1.0e-8_dp) .OR. &
734 (abs(rho_delta_avg) >= 1.0e5_dp)) THEN
735 WRITE (unit=output_unit, fmt="(T3,A,I6,4X,ES16.4,4X,ES16.4,1X,F25.12)") &
736 "SCCS| ", iter, rho_delta_avg, rho_delta_max, energy%sccs_hartree
737 ELSE
738 WRITE (unit=output_unit, fmt="(T3,A,I6,4X,F16.8,4X,F16.8,1X,F25.12)") &
739 "SCCS| ", iter, rho_delta_avg, rho_delta_max, energy%sccs_hartree
740 END IF
741 END IF
742 END IF
743
744 ! Check if the SCCS iteration cycle is converged to the requested tolerance
745 IF (rho_delta_max <= sccs_control%eps_sccs) THEN
746 IF (should_output .AND. (output_unit > 0)) THEN
747 WRITE (unit=output_unit, fmt="(T3,A,I0,A)") &
748 "SCCS| Iteration cycle converged in ", iter, " steps"
749 END IF
750 EXIT iter_loop
751 END IF
752
753 END DO iter_loop
754
755 IF (sccs_control%method_id == sccs_saa_andreussi) THEN
756 CALL auxbas_pw_pool%create_pw(tem_convolution)
757 CALL auxbas_pw_pool%create_pw(tem_func)
758
759!$OMP PARALLEL DO DEFAULT(NONE) &
760!$OMP PRIVATE(i,j,k,dphi2) &
761!$OMP SHARED(lb,ub, tem_func, s, dtf) &
762!$OMP SHARED(dphi_tot, p_e_interface, eps_elec, epsilon_solvent)
763 DO k = lb(3), ub(3)
764 DO j = lb(2), ub(2)
765 DO i = lb(1), ub(1)
766 p_e_interface%array(i, j, k) = 1.0_dp/(4.0_dp*twopi)*(exp(eps_elec%array(i, j, k))*log(epsilon_solvent))
767 dphi2 = dphi_tot(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + &
768 dphi_tot(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + &
769 dphi_tot(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k)
770 p_e_interface%array(i, j, k) = p_e_interface%array(i, j, k)*dphi2
771 tem_func%array(i, j, k) = (1.0_dp - s%array(i, j, k))*p_e_interface%array(i, j, k)*dtf%array(i, j, k)
772 END DO
773 END DO
774 END DO
775!$OMP END PARALLEL DO
776
777 CALL pw_func_u_convolution(poisson_env=poisson_env, func=tem_func, convolution=tem_convolution, u=u)
778
779!$OMP PARALLEL DO DEFAULT(NONE) &
780!$OMP PRIVATE(i,j,k) &
781!$OMP SHARED(lb,ub) &
782!$OMP SHARED(p_e_interface, t, tem_convolution, phi_tot, d_s_rhoel)
783 DO k = lb(3), ub(3)
784 DO j = lb(2), ub(2)
785 DO i = lb(1), ub(1)
786 p_e_interface%array(i, j, k) = p_e_interface%array(i, j, k)*(1.0_dp - t%array(i, j, k)) &
787 + tem_convolution%array(i, j, k)
788 END DO
789 END DO
790 END DO
791!$OMP END PARALLEL DO
792
793 END IF
794
795 ! Release work storage which is no longer needed
796 CALL auxbas_pw_pool%give_back_pw(rho_tot_zero)
797 DO i = 1, 3
798 CALL auxbas_pw_pool%give_back_pw(dln_eps_elec(i))
799 END DO
800
801 ! Optional output of the total charge density in cube file format
802 filename = "TOTAL_CHARGE_DENSITY"
803 cube_path = trim(print_path)//"%"//trim(filename)
804 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), cp_p_file)) THEN
805 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
806 my_pos_cube = "REWIND"
807 IF (append_cube) my_pos_cube = "APPEND"
808 mpi_io = .true.
809 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
810 extension=".cube", middle_name=trim(filename), &
811 file_position=my_pos_cube, log_filename=.false., &
812 mpi_io=mpi_io, fout=mpi_filename)
813 IF (output_unit > 0) THEN
814 IF (.NOT. mpi_io) THEN
815 INQUIRE (unit=cube_unit, name=filename)
816 ELSE
817 filename = mpi_filename
818 END IF
819 WRITE (unit=output_unit, fmt="(T3,A)") &
820 "SCCS| The total SCCS charge density is written in cube file format to the file:", &
821 "SCCS| "//trim(filename)
822 END IF
823 CALL cp_pw_to_cube(rho_tot, cube_unit, trim(filename), particles=particles, &
824 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), mpi_io=mpi_io)
825 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
826 END IF
827
828 ! Calculate the total SCCS Hartree energy, potential, and its
829 ! derivatives of the solute and the implicit solvent
830 CALL pw_transfer(rho_tot, rho_tot_gspace)
831 IF (calculate_stress_tensor) THEN
832 ! Request also the calculation of the stress tensor contribution
833 CALL pw_poisson_solve(poisson_env=poisson_env, &
834 density=rho_tot_gspace, &
835 ehartree=e_tot, &
836 vhartree=v_hartree_gspace, &
837 dvhartree=dphi_tot, &
838 h_stress=h_stress)
839 ELSE
840 CALL pw_poisson_solve(poisson_env=poisson_env, &
841 density=rho_tot_gspace, &
842 ehartree=e_tot, &
843 vhartree=v_hartree_gspace, &
844 dvhartree=dphi_tot)
845 END IF
846 CALL pw_transfer(v_hartree_gspace, phi_tot)
847 energy%sccs_hartree = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
848
849 ! Calculate the Hartree energy and potential of the solute only
850 block
851 TYPE(pw_r3d_rs_type) :: phi_solute
852 CALL auxbas_pw_pool%create_pw(phi_solute)
853 CALL pw_zero(phi_solute)
854 CALL pw_poisson_solve(poisson_env=poisson_env, &
855 density=rho_solute, &
856 ehartree=energy%hartree, &
857 vhartree=phi_solute)
858
859 ! Calculate the polarisation potential (store it in phi_tot)
860 ! phi_pol = phi_tot - phi_solute
861 CALL pw_axpy(phi_solute, phi_tot, alpha=-1.0_dp)
862 CALL auxbas_pw_pool%give_back_pw(phi_solute)
863 END block
864
865 ! Optional output of the SCCS polarisation potential in cube file format
866 filename = "POLARISATION_POTENTIAL"
867 cube_path = trim(print_path)//"%"//trim(filename)
868 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
869 cp_p_file)) THEN
870 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
871 my_pos_cube = "REWIND"
872 IF (append_cube) my_pos_cube = "APPEND"
873 mpi_io = .true.
874 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
875 extension=".cube", middle_name=trim(filename), &
876 file_position=my_pos_cube, log_filename=.false., &
877 mpi_io=mpi_io, fout=mpi_filename)
878 IF (output_unit > 0) THEN
879 IF (.NOT. mpi_io) THEN
880 INQUIRE (unit=cube_unit, name=filename)
881 ELSE
882 filename = mpi_filename
883 END IF
884 WRITE (unit=output_unit, fmt="(T3,A)") &
885 "SCCS| The SCCS polarisation potential is written in cube file format to the file:", &
886 "SCCS| "//trim(filename)
887 END IF
888 CALL cp_pw_to_cube(phi_tot, cube_unit, trim(filename), particles=particles, &
889 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), mpi_io=mpi_io)
890 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
891 END IF
892
893 ! Calculate the polarisation charge (store it in rho_tot)
894 ! rho_pol = rho_tot - rho_solute
895 CALL pw_axpy(rho_solute, rho_tot, alpha=-1.0_dp)
896 polarisation_charge = pw_integrate_function(rho_tot)
897
898 ! Optional output of the SCCS polarisation charge in cube file format
899 filename = "POLARISATION_CHARGE_DENSITY"
900 cube_path = trim(print_path)//"%"//trim(filename)
901 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
902 cp_p_file)) THEN
903 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
904 my_pos_cube = "REWIND"
905 IF (append_cube) my_pos_cube = "APPEND"
906 mpi_io = .true.
907 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
908 extension=".cube", middle_name=trim(filename), &
909 file_position=my_pos_cube, log_filename=.false., &
910 mpi_io=mpi_io, fout=mpi_filename)
911 IF (output_unit > 0) THEN
912 IF (.NOT. mpi_io) THEN
913 INQUIRE (unit=cube_unit, name=filename)
914 ELSE
915 filename = mpi_filename
916 END IF
917 WRITE (unit=output_unit, fmt="(T3,A)") &
918 "SCCS| The SCCS polarisation charge density is written in cube file format to the file:", &
919 "SCCS| "//trim(filename)
920 END IF
921 CALL cp_pw_to_cube(rho_tot, cube_unit, trim(filename), particles=particles, &
922 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), mpi_io=mpi_io)
923 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
924 END IF
925
926 ! Calculate SCCS polarisation energy
927 energy%sccs_pol = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
928 CALL auxbas_pw_pool%give_back_pw(rho_solute)
929 CALL auxbas_pw_pool%give_back_pw(phi_tot)
930 CALL auxbas_pw_pool%give_back_pw(rho_tot)
931 IF (sccs_control%method_id == sccs_saa_andreussi) THEN
932 CALL auxbas_pw_pool%give_back_pw(s)
933 CALL auxbas_pw_pool%give_back_pw(u)
934 CALL auxbas_pw_pool%give_back_pw(ff)
935 CALL auxbas_pw_pool%give_back_pw(t)
936 CALL auxbas_pw_pool%give_back_pw(dtf)
937 CALL auxbas_pw_pool%give_back_pw(tem_convolution)
938 CALL auxbas_pw_pool%give_back_pw(tem_func)
939 CALL auxbas_pw_pool%give_back_pw(eps_elec)
940 END IF
941 END block
942
943 ! Calculate additional solvation terms
944 energy%sccs_cav = sccs_control%gamma_solvent*cavity_surface
945 energy%sccs_dis = sccs_control%beta_solvent*cavity_volume
946 energy%sccs_rep = sccs_control%alpha_solvent*cavity_surface
947 ! Calculate solvation free energy: \delta G^el + (alpha + gamma)*S + beta*V
948 energy%sccs_sol = energy%sccs_pol + energy%sccs_rep + energy%sccs_cav + energy%sccs_dis
949
950 IF (should_output .AND. (output_unit > 0)) THEN
951 WRITE (unit=output_unit, fmt="(T3,A)") &
952 "SCCS|"
953 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.12)") &
954 "SCCS| Polarisation charge", polarisation_charge
955 !MK "SCCS| Total interaction energy [a.u.]", e_tot
956 WRITE (unit=output_unit, fmt="(T3,A)") &
957 "SCCS|"
958 CALL print_sccs_results(energy, sccs_control, output_unit)
959 END IF
960
961 ! Calculate SCCS contribution to the Kohn-Sham potential
962 IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
963
964 f = -0.5_dp*dvol/fourpi
965!$OMP PARALLEL DO DEFAULT(NONE) &
966!$OMP PRIVATE(dphi2,i,j,k) &
967!$OMP SHARED(f,deps_elec,dphi_tot) &
968!$OMP SHARED(lb,ub,v_sccs)
969 DO k = lb(3), ub(3)
970 DO j = lb(2), ub(2)
971 DO i = lb(1), ub(1)
972 dphi2 = dphi_tot(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + &
973 dphi_tot(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + &
974 dphi_tot(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k)
975 v_sccs%array(i, j, k) = v_sccs%array(i, j, k) + f*deps_elec%array(i, j, k)*dphi2
976 END DO
977 END DO
978 END DO
979!$OMP END PARALLEL DO
980
981 ELSE
982
983!$OMP PARALLEL DO DEFAULT(NONE) &
984!$OMP PRIVATE(i,j,k) &
985!$OMP SHARED(lb,ub, dvol) &
986!$OMP SHARED(p_E_interface, v_sccs, d_s_rhoel)
987 DO k = lb(3), ub(3)
988 DO j = lb(2), ub(2)
989 DO i = lb(1), ub(1)
990 v_sccs%array(i, j, k) = dvol*d_s_rhoel%array(i, j, k)*p_e_interface%array(i, j, k)
991 END DO
992 END DO
993 END DO
994!$OMP END PARALLEL DO
995
996 CALL auxbas_pw_pool%give_back_pw(d_s_rhoel)
997 CALL auxbas_pw_pool%give_back_pw(p_e_interface)
998
999 END IF
1000
1001 CALL auxbas_pw_pool%give_back_pw(deps_elec)
1002 DO i = 1, 3
1003 CALL auxbas_pw_pool%give_back_pw(dphi_tot(i))
1004 END DO
1005
1006 ! Release the SCCS printout environment
1007 CALL cp_print_key_finished_output(output_unit, logger, input, trim(print_path), &
1008 ignore_should_output=should_output)
1009
1010 CALL timestop(handle)
1011
1012 END SUBROUTINE sccs
1013
1014! **************************************************************************************************
1015!> \brief Calculate the smoothed dielectric function of Andreussi et al.
1016!> \param rho_elec ...
1017!> \param eps_elec ...
1018!> \param deps_elec ...
1019!> \param epsilon_solvent ...
1020!> \param rho_max ...
1021!> \param rho_min ...
1022!> \par History:
1023!> - Creation (16.10.2013,MK)
1024!> - Finite difference of isosurfaces implemented (21.12.2013,MK)
1025!> \author Matthias Krack (MK)
1026!> \version 1.1
1027! **************************************************************************************************
1028 SUBROUTINE andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, rho_max, &
1029 rho_min)
1030
1031 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec
1032 REAL(kind=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min
1033
1034 CHARACTER(LEN=*), PARAMETER :: routinen = 'andreussi'
1035 REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1036
1037 INTEGER :: handle, i, j, k
1038 INTEGER, DIMENSION(3) :: lb, ub
1039 REAL(kind=dp) :: diff, dq, dt, f, ln_rho_max, ln_rho_min, &
1040 q, rho, t, x, y
1041
1042 CALL timeset(routinen, handle)
1043
1044 f = log(epsilon_solvent)/twopi
1045 diff = rho_max - rho_min
1046 IF (diff < sqrt(rhotol)) cpabort("SCCS: Difference between rho(min) and rho(max) is too small")
1047 IF (rho_min >= rhotol) THEN
1048 ln_rho_max = log(rho_max)
1049 ln_rho_min = log(rho_min)
1050 q = twopi/(ln_rho_max - ln_rho_min)
1051 dq = -f*q
1052 END IF
1053
1054 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1055 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1056
1057 ! Calculate the dielectric function and its derivative
1058!$OMP PARALLEL DO DEFAULT(NONE) &
1059!$OMP PRIVATE(dt,i,j,k,rho,t,x,y) &
1060!$OMP SHARED(deps_elec,dq,eps_elec,epsilon_solvent,f,lb,ub) &
1061!$OMP SHARED(ln_rho_max,rho_elec,q,rho_max,rho_min)
1062 DO k = lb(3), ub(3)
1063 DO j = lb(2), ub(2)
1064 DO i = lb(1), ub(1)
1065 rho = rho_elec%array(i, j, k)
1066 IF (rho < rho_min) THEN
1067 eps_elec%array(i, j, k) = epsilon_solvent
1068 deps_elec%array(i, j, k) = 0.0_dp
1069 ELSE IF (rho <= rho_max) THEN
1070 x = log(rho)
1071 y = q*(ln_rho_max - x)
1072 t = f*(y - sin(y))
1073 eps_elec%array(i, j, k) = exp(t)
1074 dt = dq*(1.0_dp - cos(y))
1075 deps_elec%array(i, j, k) = eps_elec%array(i, j, k)*dt/rho
1076 ELSE
1077 eps_elec%array(i, j, k) = 1.0_dp
1078 deps_elec%array(i, j, k) = 0.0_dp
1079 END IF
1080 END DO
1081 END DO
1082 END DO
1083!$OMP END PARALLEL DO
1084
1085 CALL timestop(handle)
1086
1087 END SUBROUTINE andreussi
1088
1089! **************************************************************************************************
1090!> \brief Calculate the smoothed dielectric function of Fattebert and Gygi
1091!> \param rho_elec ...
1092!> \param eps_elec ...
1093!> \param deps_elec ...
1094!> \param epsilon_solvent ...
1095!> \param beta ...
1096!> \param rho_zero ...
1097!> \par History:
1098!> - Creation (15.10.2013,MK)
1099!> \author Matthias Krack (MK)
1100!> \version 1.0
1101! **************************************************************************************************
1102 SUBROUTINE fattebert_gygi(rho_elec, eps_elec, deps_elec, epsilon_solvent, beta, &
1103 rho_zero)
1104
1105 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec
1106 REAL(kind=dp), INTENT(IN) :: epsilon_solvent, beta, rho_zero
1107
1108 CHARACTER(LEN=*), PARAMETER :: routinen = 'fattebert_gygi'
1109 REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1110
1111 INTEGER :: handle, i, j, k
1112 INTEGER, DIMENSION(3) :: lb, ub
1113 REAL(kind=dp) :: df, f, p, q, rho, s, t, twobeta
1114
1115 CALL timeset(routinen, handle)
1116
1117 df = (1.0_dp - epsilon_solvent)/rho_zero
1118 f = 0.5_dp*(epsilon_solvent - 1.0_dp)
1119 q = 1.0_dp/rho_zero
1120 twobeta = 2.0_dp*beta
1121
1122 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1123 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1124
1125 ! Calculate the smoothed dielectric function and its derivative
1126!$OMP PARALLEL DO DEFAULT(NONE) &
1127!$OMP PRIVATE(i,j,k,p,rho,s,t) &
1128!$OMP SHARED(df,deps_elec,eps_elec,epsilon_solvent,f,lb,ub) &
1129!$OMP SHARED(q,rho_elec,twobeta)
1130 DO k = lb(3), ub(3)
1131 DO j = lb(2), ub(2)
1132 DO i = lb(1), ub(1)
1133 rho = rho_elec%array(i, j, k)
1134 IF (rho < rhotol) THEN
1135 eps_elec%array(i, j, k) = epsilon_solvent
1136 deps_elec%array(i, j, k) = 0.0_dp
1137 ELSE
1138 s = rho*q
1139 p = s**twobeta
1140 t = 1.0_dp/(1.0_dp + p)
1141 eps_elec%array(i, j, k) = 1.0_dp + f*(1.0_dp + (1.0_dp - p)*t)
1142 deps_elec%array(i, j, k) = df*twobeta*t*t*p/s
1143 END IF
1144 END DO
1145 END DO
1146 END DO
1147!$OMP END PARALLEL DO
1148
1149 CALL timestop(handle)
1150
1151 END SUBROUTINE fattebert_gygi
1152
1153! **************************************************************************************************
1154!> \brief Calculate the smoothed dielectric function of solvent aware algorithm
1155!> \param rho_elec ...
1156!> \param eps_elec ...
1157!> \param deps_elec ...
1158!> \param epsilon_solvent ...
1159!> \param rho_max ...
1160!> \param rho_min ...
1161!> \param s ...
1162!> \param d_s_rhoel ...
1163!> \author Ziwei Chai
1164!> \version 1.0
1165! **************************************************************************************************
1166 SUBROUTINE sa_andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, rho_max, &
1167 rho_min, s, d_s_rhoel)
1168
1169 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec
1170 REAL(kind=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min
1171 TYPE(pw_r3d_rs_type), INTENT(IN) :: s, d_s_rhoel
1172
1173 CHARACTER(LEN=*), PARAMETER :: routinen = 'sa_andreussi'
1174 REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1175
1176 INTEGER :: handle, i, j, k
1177 INTEGER, DIMENSION(3) :: lb, ub
1178 REAL(kind=dp) :: diff, dq, dt, f, ln_rho_max, ln_rho_min, &
1179 q, rho, t, x, y
1180
1181 CALL timeset(routinen, handle)
1182
1183 f = log(epsilon_solvent)/twopi
1184 diff = rho_max - rho_min
1185 IF (diff < sqrt(rhotol)) cpabort("SCCS: Difference between rho(min) and rho(max) is too small")
1186 IF (rho_min >= rhotol) THEN
1187 ln_rho_max = log(rho_max)
1188 ln_rho_min = log(rho_min)
1189 q = twopi/(ln_rho_max - ln_rho_min)
1190 dq = -f*q
1191 END IF
1192
1193 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1194 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1195
1196 ! Calculate the dielectric function and its derivative
1197!$OMP PARALLEL DO DEFAULT(NONE) &
1198!$OMP PRIVATE(dt,i,j,k,rho,t,x,y) &
1199!$OMP SHARED(deps_elec,dq,eps_elec,epsilon_solvent,f,lb,ub) &
1200!$OMP SHARED(ln_rho_max,rho_elec,q,rho_max,rho_min, s, d_s_rhoel)
1201 DO k = lb(3), ub(3)
1202 DO j = lb(2), ub(2)
1203 DO i = lb(1), ub(1)
1204 rho = rho_elec%array(i, j, k)
1205 IF (rho < rho_min) THEN
1206 eps_elec%array(i, j, k) = epsilon_solvent
1207 s%array(i, j, k) = 0.0_dp
1208 deps_elec%array(i, j, k) = 0.0_dp
1209 d_s_rhoel%array(i, j, k) = 0.0_dp
1210 ELSE IF (rho <= rho_max) THEN
1211 x = log(rho)
1212 y = q*(ln_rho_max - x)
1213 t = f*(y - sin(y))
1214 s%array(i, j, k) = 1 - t/(f*twopi)
1215 eps_elec%array(i, j, k) = exp(t)
1216 dt = dq*(1.0_dp - cos(y))
1217 deps_elec%array(i, j, k) = eps_elec%array(i, j, k)*dt/rho
1218 d_s_rhoel%array(i, j, k) = deps_elec%array(i, j, k)/(-eps_elec%array(i, j, k)*log(epsilon_solvent))
1219 ELSE
1220 eps_elec%array(i, j, k) = 1.0_dp
1221 deps_elec%array(i, j, k) = 0.0_dp
1222 s%array(i, j, k) = 1.0_dp
1223 d_s_rhoel%array(i, j, k) = 0.0_dp
1224 END IF
1225 END DO
1226 END DO
1227 END DO
1228!$OMP END PARALLEL DO
1229
1230 CALL timestop(handle)
1231
1232 END SUBROUTINE sa_andreussi
1233
1234! **************************************************************************************************
1235!> \brief Build the numerical derivative of a function on realspace grid
1236!> \param f ...
1237!> \param df ...
1238!> \param method ...
1239!> \param pw_env ...
1240!> \param input ...
1241!> \par History:
1242!> - Creation (15.11.2013,MK)
1243!> \author Matthias Krack (MK)
1244!> \version 1.0
1245! **************************************************************************************************
1246 SUBROUTINE derive(f, df, method, pw_env, input)
1247
1248 TYPE(pw_r3d_rs_type), INTENT(IN) :: f
1249 TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: df
1250 INTEGER, INTENT(IN) :: method
1251 TYPE(pw_env_type), POINTER :: pw_env
1252 TYPE(section_vals_type), POINTER :: input
1253
1254 CHARACTER(LEN=*), PARAMETER :: routinen = 'derive'
1255
1256 INTEGER :: border_points, handle, i
1257 INTEGER, DIMENSION(3) :: lb, n, ub
1258 TYPE(pw_c1d_gs_type), DIMENSION(2) :: work_g1d
1259 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1260 TYPE(realspace_grid_desc_type), POINTER :: rs_desc
1261 TYPE(realspace_grid_input_type) :: input_settings
1262 TYPE(realspace_grid_type), POINTER :: rs_grid
1263 TYPE(section_vals_type), POINTER :: rs_grid_section
1264
1265 CALL timeset(routinen, handle)
1266
1267 cpassert(ASSOCIATED(pw_env))
1268
1269 ! Perform method specific setup
1270 SELECT CASE (method)
1271 CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7)
1272 NULLIFY (rs_desc)
1273 rs_grid_section => section_vals_get_subs_vals(input, "DFT%MGRID%RS_GRID")
1274 SELECT CASE (method)
1275 CASE (sccs_derivative_cd3)
1276 border_points = 1
1277 CASE (sccs_derivative_cd5)
1278 border_points = 2
1279 CASE (sccs_derivative_cd7)
1280 border_points = 3
1281 END SELECT
1282 CALL init_input_type(input_settings, 2*border_points + 1, rs_grid_section, &
1283 1, [-1, -1, -1])
1284 CALL rs_grid_create_descriptor(rs_desc, f%pw_grid, input_settings, &
1285 border_points=border_points)
1286 ALLOCATE (rs_grid)
1287 CALL rs_grid_create(rs_grid, rs_desc)
1288!MK CALL rs_grid_print(rs_grid, 6)
1289 CASE (sccs_derivative_fft)
1290 lb(1:3) = f%pw_grid%bounds_local(1, 1:3)
1291 ub(1:3) = f%pw_grid%bounds_local(2, 1:3)
1292 NULLIFY (auxbas_pw_pool)
1293 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1294 ! Get work storage for the 1d grids in g-space (derivative calculation)
1295 DO i = 1, SIZE(work_g1d)
1296 CALL auxbas_pw_pool%create_pw(work_g1d(i))
1297 END DO
1298 END SELECT
1299
1300 ! Calculate the derivatives
1301 SELECT CASE (method)
1302 CASE (sccs_derivative_cd3)
1303 CALL derive_fdm_cd3(f, df, rs_grid)
1304 CASE (sccs_derivative_cd5)
1305 CALL derive_fdm_cd5(f, df, rs_grid)
1306 CASE (sccs_derivative_cd7)
1307 CALL derive_fdm_cd7(f, df, rs_grid)
1308 CASE (sccs_derivative_fft)
1309 ! FFT
1310 CALL pw_transfer(f, work_g1d(1))
1311 DO i = 1, 3
1312 n(:) = 0
1313 n(i) = 1
1314 CALL pw_copy(work_g1d(1), work_g1d(2))
1315 CALL pw_derive(work_g1d(2), n(:))
1316 CALL pw_transfer(work_g1d(2), df(i))
1317 END DO
1318 CASE DEFAULT
1319 cpabort("Invalid derivative method for SCCS specified")
1320 END SELECT
1321
1322 ! Perform method specific cleanup
1323 SELECT CASE (method)
1324 CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7)
1325 CALL rs_grid_release(rs_grid)
1326 DEALLOCATE (rs_grid)
1327 CALL rs_grid_release_descriptor(rs_desc)
1328 CASE (sccs_derivative_fft)
1329 DO i = 1, SIZE(work_g1d)
1330 CALL auxbas_pw_pool%give_back_pw(work_g1d(i))
1331 END DO
1332 END SELECT
1333
1334 CALL timestop(handle)
1335
1336 END SUBROUTINE derive
1337
1338! **************************************************************************************************
1339!> \brief Calculate the finite difference between two isosurfaces of the
1340!> electronic density. The smoothed dielectric function of
1341!> Andreussi et al. is used as switching function eventually
1342!> defining the quantum volume and surface of the cavity.
1343!> \param rho_elec ...
1344!> \param norm_drho_elec ...
1345!> \param dtheta ...
1346!> \param epsilon_solvent ...
1347!> \param rho_max ...
1348!> \param rho_min ...
1349!> \param delta_rho ...
1350!> \par History:
1351!> - Creation (21.12.2013,MK)
1352!> \author Matthias Krack (MK)
1353!> \version 1.0
1354! **************************************************************************************************
1355 SUBROUTINE surface_andreussi(rho_elec, norm_drho_elec, dtheta, &
1356 epsilon_solvent, rho_max, rho_min, delta_rho)
1357
1358 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, norm_drho_elec, dtheta
1359 REAL(kind=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min, &
1360 delta_rho
1361
1362 CHARACTER(LEN=*), PARAMETER :: routinen = 'surface_andreussi'
1363 REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1364
1365 INTEGER :: handle, i, j, k, l
1366 INTEGER, DIMENSION(3) :: lb, ub
1367 REAL(kind=dp) :: diff, e, eps_elec, f, ln_rho_max, &
1368 ln_rho_min, q, rho, t, x, y
1369 REAL(kind=dp), DIMENSION(2) :: theta
1370
1371 CALL timeset(routinen, handle)
1372
1373 e = epsilon_solvent - 1.0_dp
1374 f = log(epsilon_solvent)/twopi
1375 diff = rho_max - rho_min
1376 IF (diff < sqrt(rhotol)) cpabort("SCCS: Difference between rho(min) and rho(max) is too small")
1377 IF (rho_min >= rhotol) THEN
1378 ln_rho_max = log(rho_max)
1379 ln_rho_min = log(rho_min)
1380 q = twopi/(ln_rho_max - ln_rho_min)
1381 END IF
1382
1383 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1384 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1385
1386 ! Calculate finite difference between two isosurfaces
1387!$OMP PARALLEL DO DEFAULT(NONE) &
1388!$OMP PRIVATE(eps_elec,i,j,k,l,rho,t,theta,x,y) &
1389!$OMP SHARED(delta_rho,dtheta,e,epsilon_solvent,f,lb) &
1390!$OMP SHARED(ln_rho_max,norm_drho_elec,rho_elec,q,rho_max,rho_min,ub)
1391 DO k = lb(3), ub(3)
1392 DO j = lb(2), ub(2)
1393 DO i = lb(1), ub(1)
1394 DO l = 1, 2
1395 rho = rho_elec%array(i, j, k) + (real(l, kind=dp) - 1.5_dp)*delta_rho
1396 IF (rho < rho_min) THEN
1397 eps_elec = epsilon_solvent
1398 ELSE IF (rho <= rho_max) THEN
1399 x = log(rho)
1400 y = q*(ln_rho_max - x)
1401 t = f*(y - sin(y))
1402 eps_elec = exp(t)
1403 ELSE
1404 eps_elec = 1.0_dp
1405 END IF
1406 theta(l) = (epsilon_solvent - eps_elec)/e
1407 END DO
1408 dtheta%array(i, j, k) = (theta(2) - theta(1))*norm_drho_elec%array(i, j, k)/delta_rho
1409 END DO
1410 END DO
1411 END DO
1412!$OMP END PARALLEL DO
1413
1414 CALL timestop(handle)
1415
1416 END SUBROUTINE surface_andreussi
1417
1418! **************************************************************************************************
1419!> \brief Calculate the finite difference between two isosurfaces of the
1420!> the electronic density. The smoothed dielectric function of
1421!> Fattebert and Gygi is used as switching function eventually
1422!> defining the quantum volume and surface of the cavity.
1423!> \param rho_elec ...
1424!> \param norm_drho_elec ...
1425!> \param dtheta ...
1426!> \param epsilon_solvent ...
1427!> \param beta ...
1428!> \param rho_zero ...
1429!> \param delta_rho ...
1430!> \par History:
1431!> - Creation (21.12.2013,MK)
1432!> \author Matthias Krack (MK)
1433!> \version 1.0
1434! **************************************************************************************************
1435 SUBROUTINE surface_fattebert_gygi(rho_elec, norm_drho_elec, dtheta, &
1436 epsilon_solvent, beta, rho_zero, delta_rho)
1437
1438 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, norm_drho_elec, dtheta
1439 REAL(kind=dp), INTENT(IN) :: epsilon_solvent, beta, rho_zero, &
1440 delta_rho
1441
1442 CHARACTER(LEN=*), PARAMETER :: routinen = 'surface_fattebert_gygi'
1443 REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1444
1445 INTEGER :: handle, i, j, k, l
1446 INTEGER, DIMENSION(3) :: lb, ub
1447 REAL(kind=dp) :: e, eps_elec, f, p, q, rho, s, t, twobeta
1448 REAL(kind=dp), DIMENSION(2) :: theta
1449
1450 CALL timeset(routinen, handle)
1451
1452 e = epsilon_solvent - 1.0_dp
1453 f = 0.5_dp*e
1454 q = 1.0_dp/rho_zero
1455 twobeta = 2.0_dp*beta
1456
1457 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1458 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1459
1460 ! Calculate finite difference between two isosurfaces
1461!$OMP PARALLEL DO DEFAULT(NONE) &
1462!$OMP PRIVATE(eps_elec,i,j,k,l,p,rho,s,t,theta) &
1463!$OMP SHARED(delta_rho,dtheta,e,epsilon_solvent,f,lb) &
1464!$OMP SHARED(norm_drho_elec,q,rho_elec,twobeta,ub)
1465 DO k = lb(3), ub(3)
1466 DO j = lb(2), ub(2)
1467 DO i = lb(1), ub(1)
1468 DO l = 1, 2
1469 rho = rho_elec%array(i, j, k) + (real(l, kind=dp) - 1.5_dp)*delta_rho
1470 IF (rho < rhotol) THEN
1471 eps_elec = epsilon_solvent
1472 ELSE
1473 s = rho*q
1474 p = s**twobeta
1475 t = 1.0_dp/(1.0_dp + p)
1476 eps_elec = 1.0_dp + f*(1.0_dp + (1.0_dp - p)*t)
1477 END IF
1478 theta(l) = (epsilon_solvent - eps_elec)/e
1479 END DO
1480 dtheta%array(i, j, k) = (theta(2) - theta(1))*norm_drho_elec%array(i, j, k)/delta_rho
1481 END DO
1482 END DO
1483 END DO
1484!$OMP END PARALLEL DO
1485
1486 CALL timestop(handle)
1487
1488 END SUBROUTINE surface_fattebert_gygi
1489
1490! **************************************************************************************************
1491!> \brief Print SCCS results
1492!> \param energy ...
1493!> \param sccs_control ...
1494!> \param output_unit ...
1495!> \par History:
1496!> - Creation (11.11.2022,MK)
1497!> \author Matthias Krack (MK)
1498!> \version 1.0
1499! **************************************************************************************************
1500 SUBROUTINE print_sccs_results(energy, sccs_control, output_unit)
1501
1502 TYPE(qs_energy_type), POINTER :: energy
1503 TYPE(sccs_control_type), POINTER :: sccs_control
1504 INTEGER, INTENT(IN) :: output_unit
1505
1506 IF (output_unit > 0) THEN
1507 cpassert(ASSOCIATED(energy))
1508 cpassert(ASSOCIATED(sccs_control))
1509 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1510 "SCCS| Hartree energy of solute and solvent [Hartree]", energy%sccs_hartree, &
1511 "SCCS| Hartree energy of the solute only [Hartree]", energy%hartree
1512 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14,/,T3,A,T61,F20.3)") &
1513 "SCCS| Polarisation energy [Hartree]", energy%sccs_pol, &
1514 "SCCS| [kcal/mol]", &
1515 cp_unit_from_cp2k(energy%sccs_pol, "kcalmol"), &
1516 "SCCS| Cavitation energy [Hartree]", energy%sccs_cav, &
1517 "SCCS| [kcal/mol]", &
1518 cp_unit_from_cp2k(energy%sccs_cav, "kcalmol"), &
1519 "SCCS| Dispersion free energy [Hartree]", energy%sccs_dis, &
1520 "SCCS| [kcal/mol]", &
1521 cp_unit_from_cp2k(energy%sccs_dis, "kcalmol"), &
1522 "SCCS| Repulsion free energy [Hartree]", energy%sccs_rep, &
1523 "SCCS| [kcal/mol]", &
1524 cp_unit_from_cp2k(energy%sccs_rep, "kcalmol"), &
1525 "SCCS| Solvation free energy [Hartree]", energy%sccs_sol, &
1526 "SCCS| [kcal/mol]", &
1527 cp_unit_from_cp2k(energy%sccs_sol, "kcalmol")
1528 END IF
1529
1530 END SUBROUTINE print_sccs_results
1531
1532END MODULE qs_sccs
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_integrate_function
Definition pw_methods.F:116
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
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:853
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:1073
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, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
...
Definition cp_realspace_grid_cube.F:48
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, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
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:1251
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::sccs_derivative_cd5
integer, parameter, public sccs_derivative_cd5
Definition input_constants.F:1240
input_constants::sccs_saa_andreussi
integer, parameter, public sccs_saa_andreussi
Definition input_constants.F:1235
input_constants::sccs_fattebert_gygi
integer, parameter, public sccs_fattebert_gygi
Definition input_constants.F:1235
input_constants::sccs_derivative_cd7
integer, parameter, public sccs_derivative_cd7
Definition input_constants.F:1240
input_constants::sccs_derivative_fft
integer, parameter, public sccs_derivative_fft
Definition input_constants.F:1240
input_constants::sccs_derivative_cd3
integer, parameter, public sccs_derivative_cd3
Definition input_constants.F:1240
input_constants::sccs_andreussi
integer, parameter, public sccs_andreussi
Definition input_constants.F:1235
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:1001
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:735
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition input_section_types.F:1017
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::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
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:10257
pw_poisson_methods
Definition pw_poisson_methods.F:13
pw_poisson_methods::pw_func_u_convolution
subroutine, public pw_func_u_convolution(poisson_env, func, convolution, u)
computes the convolution f * u using reciprocal-space multiplication
Definition pw_poisson_methods.F:4097
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:67
pw_poisson_types::pw_poisson_analytic
integer, parameter, public pw_poisson_analytic
Definition pw_poisson_types.F:67
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, mimic, 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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
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:126
qs_sccs::print_sccs_results
subroutine, public print_sccs_results(energy, sccs_control, output_unit)
Print SCCS results.
Definition qs_sccs.F:1501
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:478
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:202
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:2034
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:2014
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
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_control_types::sccs_control_type
Definition cp_control_types.F:562
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
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:123
pw_pool_types::pw_pool_p_type
to create arrays of pools
Definition pw_pool_types.F:135
pw_pool_types::pw_pool_type
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:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:97
realspace_grid_types::realspace_grid_desc_type
Definition realspace_grid_types.F:90
realspace_grid_types::realspace_grid_input_type
Definition realspace_grid_types.F:80
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136