Line data Source code
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 :
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 : 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 :
110 : CONTAINS
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 160 : 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 160 : 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 1120 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: dln_eps_elec, dphi_tot
162 160 : 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 160 : CALL timeset(routineN, handle)
171 :
172 160 : NULLIFY (auxbas_pw_pool)
173 160 : NULLIFY (cp_subsys)
174 160 : NULLIFY (dft_control)
175 160 : NULLIFY (energy)
176 160 : NULLIFY (input)
177 160 : NULLIFY (logger)
178 160 : NULLIFY (para_env)
179 160 : NULLIFY (particles)
180 160 : NULLIFY (poisson_env)
181 160 : NULLIFY (pw_env)
182 160 : NULLIFY (pw_pools)
183 160 : NULLIFY (rho)
184 160 : NULLIFY (sccs_control)
185 160 : 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 160 : scf_env=scf_env)
198 160 : CALL cp_subsys_get(cp_subsys, particles=particles)
199 :
200 160 : sccs_control => dft_control%sccs_control
201 :
202 160 : CPASSERT(ASSOCIATED(qs_env))
203 :
204 160 : IF (do_kpoints) THEN
205 0 : CPWARN("SCCS with k-points has not yet been fully validated")
206 : END IF
207 :
208 160 : IF (PRESENT(h_stress)) THEN
209 0 : calculate_stress_tensor = .TRUE.
210 0 : h_stress(:, :) = 0.0_dp
211 0 : 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 160 : poisson_env=poisson_env)
221 :
222 160 : CALL pw_zero(v_sccs)
223 :
224 : ! Calculate no SCCS contribution, if the requested SCF convergence threshold is not reached yet
225 160 : IF (.NOT. sccs_control%sccs_activated) THEN
226 54 : IF (sccs_control%eps_scf > 0.0_dp) THEN
227 52 : 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 42 : 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 0 : 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 42 : vhartree=v_hartree_gspace)
242 : END IF
243 42 : energy%sccs_pol = 0.0_dp
244 42 : energy%sccs_cav = 0.0_dp
245 42 : energy%sccs_dis = 0.0_dp
246 42 : energy%sccs_rep = 0.0_dp
247 42 : energy%sccs_sol = 0.0_dp
248 42 : energy%sccs_hartree = energy%hartree
249 42 : CALL timestop(handle)
250 42 : RETURN
251 : END IF
252 : END IF
253 12 : sccs_control%sccs_activated = .TRUE.
254 : END IF
255 :
256 118 : nspin = dft_control%nspins
257 :
258 : ! Manage print output control
259 118 : logger => cp_get_default_logger()
260 118 : print_level = logger%iter_info%print_level
261 118 : print_path = "DFT%PRINT%SCCS"
262 : should_output = (BTEST(cp_print_key_should_output(logger%iter_info, input, &
263 118 : 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 118 : log_filename=.FALSE.)
268 :
269 : ! Get rho
270 : CALL qs_rho_get(rho_struct=rho, &
271 : rho_r=rho_pw_r, &
272 118 : rho_r_sccs=rho_pw_r_sccs)
273 :
274 : ! Retrieve the last rho_iter from the previous SCCS cycle if available
275 118 : 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 118 : CALL auxbas_pw_pool%create_pw(rho_elec)
281 :
282 : ! Retrieve grid parameters
283 118 : ngpts = rho_elec%pw_grid%ngpts
284 118 : dvol = rho_elec%pw_grid%dvol
285 118 : cell_volume = rho_elec%pw_grid%vol
286 472 : abc(1:3) = REAL(rho_elec%pw_grid%npts(1:3), KIND=dp)*rho_elec%pw_grid%dr(1:3)
287 472 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
288 472 : ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
289 :
290 118 : CALL pw_copy(rho_pw_r(1), rho_elec)
291 144 : DO ispin = 2, nspin
292 144 : CALL pw_axpy(rho_pw_r(ispin), rho_elec)
293 : END DO
294 118 : tot_rho_elec = pw_integrate_function(rho_elec)
295 :
296 : ! Calculate the dielectric (smoothed) function of rho_elec in r-space
297 118 : CALL auxbas_pw_pool%create_pw(eps_elec)
298 118 : CALL auxbas_pw_pool%create_pw(deps_elec)
299 118 : IF (sccs_control%method_id == sccs_saa_andreussi) THEN
300 0 : CALL auxbas_pw_pool%create_pw(s)
301 0 : CALL auxbas_pw_pool%create_pw(d_s_rhoel)
302 0 : CALL auxbas_pw_pool%create_pw(u)
303 0 : CALL auxbas_pw_pool%create_pw(ff)
304 0 : CALL auxbas_pw_pool%create_pw(t)
305 0 : CALL auxbas_pw_pool%create_pw(dtf)
306 0 : f0 = sccs_control%f0
307 0 : delta_eta = sccs_control%delta_eta
308 0 : alpha_zeta = sccs_control%alpha_zeta
309 0 : R_solv = sccs_control%R_solv
310 0 : delta_zeta = sccs_control%delta_zeta
311 : END IF
312 :
313 : ! Relative permittivity or dielectric constant of the solvent (medium)
314 118 : epsilon_solvent = sccs_control%epsilon_solvent
315 214 : 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 96 : 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 22 : 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 0 : sccs_control%rho_min, s, d_s_rhoel)
325 : CASE DEFAULT
326 118 : CPABORT("Invalid method specified for SCCS model")
327 : END SELECT
328 :
329 118 : IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
330 : ! Optional output of the dielectric function in cube file format
331 118 : filename = "DIELECTRIC_FUNCTION"
332 118 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
333 118 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), &
334 : cp_p_file)) THEN
335 4 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
336 4 : my_pos_cube = "REWIND"
337 4 : IF (append_cube) my_pos_cube = "APPEND"
338 4 : 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 4 : mpi_io=mpi_io, fout=mpi_filename)
343 4 : IF (output_unit > 0) THEN
344 2 : IF (.NOT. mpi_io) THEN
345 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
346 : ELSE
347 2 : filename = mpi_filename
348 : END IF
349 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
350 2 : "SCCS| The dielectric function is written in cube file format to the file:", &
351 4 : "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 4 : mpi_io=mpi_io)
356 4 : 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 118 : cavity_surface = 0.0_dp
362 118 : cavity_volume = 0.0_dp
363 :
364 118 : IF (ABS(epsilon_solvent - 1.0_dp) > epstol) THEN
365 :
366 : BLOCK
367 : TYPE(pw_r3d_rs_type) :: theta, norm_drho_elec
368 472 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: drho_elec
369 118 : CALL auxbas_pw_pool%create_pw(theta)
370 118 : CALL pw_zero(theta)
371 :
372 : ! Calculate the (quantum) volume of the solute cavity
373 118 : f = 1.0_dp/(epsilon_solvent - 1.0_dp)
374 : !$OMP PARALLEL DO DEFAULT(NONE) &
375 : !$OMP PRIVATE(i,j,k) &
376 118 : !$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 118 : 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 472 : DO i = 1, 3
390 472 : CALL auxbas_pw_pool%create_pw(drho_elec(i))
391 : END DO
392 118 : CALL derive(rho_elec, drho_elec, sccs_derivative_fft, pw_env, input)
393 :
394 118 : 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 118 : !$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 118 : filename = "DENSITY_GRADIENT"
416 118 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
417 118 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), &
418 : cp_p_file)) THEN
419 0 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
420 0 : my_pos_cube = "REWIND"
421 0 : IF (append_cube) my_pos_cube = "APPEND"
422 0 : 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 0 : mpi_io=mpi_io, fout=mpi_filename)
427 0 : IF (output_unit > 0) THEN
428 0 : IF (.NOT. mpi_io) THEN
429 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
430 : ELSE
431 0 : filename = mpi_filename
432 : END IF
433 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
434 0 : "SCCS| The norm of the density gradient is written in cube file format to the file:", &
435 0 : "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 0 : mpi_io=mpi_io)
440 0 : 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 214 : 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 96 : 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 22 : sccs_control%delta_rho)
453 : CASE (sccs_saa_andreussi)
454 : CASE DEFAULT
455 118 : CPABORT("Invalid method specified for SCCS model")
456 : END SELECT
457 118 : cavity_surface = pw_integrate_function(theta)
458 :
459 : ! Release storage
460 118 : CALL auxbas_pw_pool%give_back_pw(theta)
461 118 : CALL auxbas_pw_pool%give_back_pw(norm_drho_elec)
462 590 : DO i = 1, 3
463 472 : 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 118 : 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 118 : CALL auxbas_pw_pool%create_pw(rho_solute)
476 118 : CALL pw_zero(rho_solute)
477 118 : CALL pw_transfer(rho_tot_gspace, rho_solute)
478 118 : tot_rho_solute = pw_integrate_function(rho_solute)
479 :
480 : ! Check total charge
481 118 : IF (ABS(tot_rho_solute) >= 1.0E-6_dp) THEN
482 92 : 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 92 : "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 184 : "combined with an appropriate Poisson solver (POISSON_SOLVER MT or analytic)"
488 92 : 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 118 : CALL auxbas_pw_pool%create_pw(rho_tot_zero)
494 :
495 118 : 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 118 : !$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 0 : lxyz(:) = REAL(u%pw_grid%bounds(1, :), dp)
522 0 : 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 0 : !$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 0 : CALL pw_func_u_convolution(poisson_env=poisson_env, func=dtf, convolution=ff, u=u)
551 0 : CALL pw_scale(u, 1.0_dp/ff%array(0, 0, 0))
552 0 : 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 0 : !$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 0 : filename = "DIELECTRIC_FUNCTION"
574 0 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
575 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), &
576 : cp_p_file)) THEN
577 0 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
578 0 : my_pos_cube = "REWIND"
579 0 : IF (append_cube) my_pos_cube = "APPEND"
580 0 : 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 0 : mpi_io=mpi_io, fout=mpi_filename)
585 0 : IF (output_unit > 0) THEN
586 0 : IF (.NOT. mpi_io) THEN
587 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
588 : ELSE
589 0 : filename = mpi_filename
590 : END IF
591 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
592 0 : "SCCS| The dielectric function is written in cube file format to the file:", &
593 0 : "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 0 : mpi_io=mpi_io)
598 0 : 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 0 : !$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 472 : DO i = 1, 3
618 354 : CALL auxbas_pw_pool%create_pw(dln_eps_elec(i))
619 472 : CALL pw_zero(dln_eps_elec(i))
620 : END DO
621 118 : CALL derive(eps_elec, dln_eps_elec, sccs_control%derivative_method, pw_env, input)
622 118 : IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
623 118 : CALL auxbas_pw_pool%give_back_pw(eps_elec)
624 : END IF
625 :
626 : ! Print header for the SCCS cycle
627 118 : IF (should_output .AND. (output_unit > 0)) THEN
628 46 : IF (print_level > low_print_level) THEN
629 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.12)") &
630 46 : "SCCS| Total electronic charge density ", -tot_rho_elec, &
631 92 : "SCCS| Total charge density (solute) ", -tot_rho_solute
632 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.3)") &
633 46 : "SCCS| Volume of the cell [bohr^3]", cell_volume, &
634 46 : "SCCS| [angstrom^3]", &
635 92 : cp_unit_from_cp2k(cell_volume, "angstrom^3")
636 46 : IF (ABS(epsilon_solvent - 1.0_dp) > epstol) THEN
637 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.3)") &
638 46 : "SCCS| Volume of the solute cavity [bohr^3]", cavity_volume, &
639 46 : "SCCS| [angstrom^3]", &
640 46 : cp_unit_from_cp2k(cavity_volume, "angstrom^3"), &
641 46 : "SCCS| Surface of the solute cavity [bohr^2]", cavity_surface, &
642 46 : "SCCS| [angstrom^2]", &
643 92 : cp_unit_from_cp2k(cavity_surface, "angstrom^2")
644 : END IF
645 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
646 46 : "SCCS|", &
647 92 : "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 472 : DO i = 1, 3
653 472 : 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 118 : CALL auxbas_pw_pool%create_pw(rho_tot)
658 118 : CALL pw_copy(rho_tot_zero, rho_tot)
659 118 : CALL pw_axpy(rho_pw_r_sccs, rho_tot)
660 :
661 118 : CALL auxbas_pw_pool%create_pw(phi_tot)
662 118 : CALL pw_zero(phi_tot)
663 :
664 118 : IF (sccs_control%method_id == sccs_saa_andreussi) THEN
665 0 : CALL auxbas_pw_pool%create_pw(p_e_interface)
666 0 : CALL pw_zero(p_e_interface)
667 : END IF
668 :
669 : ! Main SCCS iteration loop
670 118 : iter = 0
671 :
672 : iter_loop: DO
673 :
674 : ! Increment iteration counter
675 692 : iter = iter + 1
676 :
677 : ! Check if the requested maximum number of SCCS iterations is reached
678 692 : IF (iter > sccs_control%max_iter) THEN
679 0 : IF (output_unit > 0) THEN
680 : WRITE (UNIT=output_unit, FMT="(T3,A,/,T3,A,I0,A)") &
681 0 : "SCCS| Maximum number of SCCS iterations reached", &
682 0 : "SCCS| Iteration cycle did not converge in ", sccs_control%max_iter, " steps"
683 : ELSE
684 : WRITE (UNIT=message, FMT="(A,I0,A)") &
685 0 : "The SCCS iteration cycle did not converge in ", sccs_control%max_iter, " steps"
686 0 : 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 692 : dvhartree=dphi_tot)
696 692 : 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 692 : f = 1.0_dp/fourpi
701 692 : rho_delta_avg = 0.0_dp
702 692 : 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 692 : !$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 692 : CALL para_env%sum(rho_delta_avg)
728 692 : rho_delta_avg = rho_delta_avg/REAL(ngpts, KIND=dp)
729 692 : CALL para_env%max(rho_delta_max)
730 :
731 692 : IF (should_output .AND. (output_unit > 0)) THEN
732 147 : IF (print_level > low_print_level) THEN
733 147 : IF ((ABS(rho_delta_avg) < 1.0E-8_dp) .OR. &
734 147 : (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 0 : "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 147 : "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 692 : IF (rho_delta_max <= sccs_control%eps_sccs) THEN
746 118 : IF (should_output .AND. (output_unit > 0)) THEN
747 : WRITE (UNIT=output_unit, FMT="(T3,A,I0,A)") &
748 46 : "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 118 : IF (sccs_control%method_id == sccs_saa_andreussi) THEN
756 0 : CALL auxbas_pw_pool%create_pw(tem_convolution)
757 0 : 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 0 : !$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 0 : 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 0 : !$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 118 : CALL auxbas_pw_pool%give_back_pw(rho_tot_zero)
797 472 : DO i = 1, 3
798 472 : 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 118 : filename = "TOTAL_CHARGE_DENSITY"
803 118 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
804 118 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), cp_p_file)) THEN
805 0 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
806 0 : my_pos_cube = "REWIND"
807 0 : IF (append_cube) my_pos_cube = "APPEND"
808 0 : 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 0 : mpi_io=mpi_io, fout=mpi_filename)
813 0 : IF (output_unit > 0) THEN
814 0 : IF (.NOT. mpi_io) THEN
815 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
816 : ELSE
817 0 : filename = mpi_filename
818 : END IF
819 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
820 0 : "SCCS| The total SCCS charge density is written in cube file format to the file:", &
821 0 : "SCCS| "//TRIM(filename)
822 : END IF
823 : CALL cp_pw_to_cube(rho_tot, cube_unit, TRIM(filename), particles=particles, &
824 0 : stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), mpi_io=mpi_io)
825 0 : 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 118 : CALL pw_transfer(rho_tot, rho_tot_gspace)
831 118 : 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 0 : 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 118 : dvhartree=dphi_tot)
845 : END IF
846 118 : CALL pw_transfer(v_hartree_gspace, phi_tot)
847 118 : 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 118 : CALL auxbas_pw_pool%create_pw(phi_solute)
853 118 : CALL pw_zero(phi_solute)
854 : CALL pw_poisson_solve(poisson_env=poisson_env, &
855 : density=rho_solute, &
856 : ehartree=energy%hartree, &
857 118 : vhartree=phi_solute)
858 :
859 : ! Calculate the polarisation potential (store it in phi_tot)
860 : ! phi_pol = phi_tot - phi_solute
861 118 : CALL pw_axpy(phi_solute, phi_tot, alpha=-1.0_dp)
862 118 : 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 118 : filename = "POLARISATION_POTENTIAL"
867 118 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
868 118 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), &
869 : cp_p_file)) THEN
870 4 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
871 4 : my_pos_cube = "REWIND"
872 4 : IF (append_cube) my_pos_cube = "APPEND"
873 4 : 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 4 : mpi_io=mpi_io, fout=mpi_filename)
878 4 : IF (output_unit > 0) THEN
879 2 : IF (.NOT. mpi_io) THEN
880 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
881 : ELSE
882 2 : filename = mpi_filename
883 : END IF
884 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
885 2 : "SCCS| The SCCS polarisation potential is written in cube file format to the file:", &
886 4 : "SCCS| "//TRIM(filename)
887 : END IF
888 : CALL cp_pw_to_cube(phi_tot, cube_unit, TRIM(filename), particles=particles, &
889 4 : stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), mpi_io=mpi_io)
890 4 : 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 118 : CALL pw_axpy(rho_solute, rho_tot, alpha=-1.0_dp)
896 118 : polarisation_charge = pw_integrate_function(rho_tot)
897 :
898 : ! Optional output of the SCCS polarisation charge in cube file format
899 118 : filename = "POLARISATION_CHARGE_DENSITY"
900 118 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
901 118 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), &
902 : cp_p_file)) THEN
903 0 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
904 0 : my_pos_cube = "REWIND"
905 0 : IF (append_cube) my_pos_cube = "APPEND"
906 0 : 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 0 : mpi_io=mpi_io, fout=mpi_filename)
911 0 : IF (output_unit > 0) THEN
912 0 : IF (.NOT. mpi_io) THEN
913 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
914 : ELSE
915 0 : filename = mpi_filename
916 : END IF
917 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
918 0 : "SCCS| The SCCS polarisation charge density is written in cube file format to the file:", &
919 0 : "SCCS| "//TRIM(filename)
920 : END IF
921 : CALL cp_pw_to_cube(rho_tot, cube_unit, TRIM(filename), particles=particles, &
922 0 : stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), mpi_io=mpi_io)
923 0 : 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 118 : energy%sccs_pol = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
928 118 : CALL auxbas_pw_pool%give_back_pw(rho_solute)
929 118 : CALL auxbas_pw_pool%give_back_pw(phi_tot)
930 118 : CALL auxbas_pw_pool%give_back_pw(rho_tot)
931 236 : IF (sccs_control%method_id == sccs_saa_andreussi) THEN
932 0 : CALL auxbas_pw_pool%give_back_pw(s)
933 0 : CALL auxbas_pw_pool%give_back_pw(u)
934 0 : CALL auxbas_pw_pool%give_back_pw(ff)
935 0 : CALL auxbas_pw_pool%give_back_pw(t)
936 0 : CALL auxbas_pw_pool%give_back_pw(dtf)
937 0 : CALL auxbas_pw_pool%give_back_pw(tem_convolution)
938 0 : CALL auxbas_pw_pool%give_back_pw(tem_func)
939 0 : CALL auxbas_pw_pool%give_back_pw(eps_elec)
940 : END IF
941 : END BLOCK
942 :
943 : ! Calculate additional solvation terms
944 118 : energy%sccs_cav = sccs_control%gamma_solvent*cavity_surface
945 118 : energy%sccs_dis = sccs_control%beta_solvent*cavity_volume
946 118 : energy%sccs_rep = sccs_control%alpha_solvent*cavity_surface
947 : ! Calculate solvation free energy: \delta G^el + (alpha + gamma)*S + beta*V
948 118 : energy%sccs_sol = energy%sccs_pol + energy%sccs_rep + energy%sccs_cav + energy%sccs_dis
949 :
950 118 : IF (should_output .AND. (output_unit > 0)) THEN
951 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
952 46 : "SCCS|"
953 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.12)") &
954 46 : "SCCS| Polarisation charge", polarisation_charge
955 : !MK "SCCS| Total interaction energy [a.u.]", e_tot
956 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
957 46 : "SCCS|"
958 46 : CALL print_sccs_results(energy, sccs_control, output_unit)
959 : END IF
960 :
961 : ! Calculate SCCS contribution to the Kohn-Sham potential
962 118 : IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
963 :
964 118 : 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 118 : !$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 0 : !$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 0 : CALL auxbas_pw_pool%give_back_pw(d_s_rhoel)
997 0 : CALL auxbas_pw_pool%give_back_pw(p_E_interface)
998 :
999 : END IF
1000 :
1001 118 : CALL auxbas_pw_pool%give_back_pw(deps_elec)
1002 472 : DO i = 1, 3
1003 472 : 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 118 : ignore_should_output=should_output)
1009 :
1010 118 : CALL timestop(handle)
1011 :
1012 320 : 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 96 : 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 96 : CALL timeset(routineN, handle)
1043 :
1044 96 : f = LOG(epsilon_solvent)/twopi
1045 96 : diff = rho_max - rho_min
1046 96 : IF (diff < SQRT(rhotol)) CPABORT("SCCS: Difference between rho(min) and rho(max) is too small")
1047 96 : IF (rho_min >= rhotol) THEN
1048 96 : ln_rho_max = LOG(rho_max)
1049 96 : ln_rho_min = LOG(rho_min)
1050 96 : q = twopi/(ln_rho_max - ln_rho_min)
1051 96 : dq = -f*q
1052 : END IF
1053 :
1054 384 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1055 384 : 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 96 : !$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 96 : CALL timestop(handle)
1086 :
1087 96 : 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 22 : 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 22 : CALL timeset(routineN, handle)
1116 :
1117 22 : df = (1.0_dp - epsilon_solvent)/rho_zero
1118 22 : f = 0.5_dp*(epsilon_solvent - 1.0_dp)
1119 22 : q = 1.0_dp/rho_zero
1120 22 : twobeta = 2.0_dp*beta
1121 :
1122 88 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1123 88 : 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 22 : !$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 22 : CALL timestop(handle)
1150 :
1151 22 : 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 0 : 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 0 : CALL timeset(routineN, handle)
1182 :
1183 0 : f = LOG(epsilon_solvent)/twopi
1184 0 : diff = rho_max - rho_min
1185 0 : IF (diff < SQRT(rhotol)) CPABORT("SCCS: Difference between rho(min) and rho(max) is too small")
1186 0 : IF (rho_min >= rhotol) THEN
1187 0 : ln_rho_max = LOG(rho_max)
1188 0 : ln_rho_min = LOG(rho_min)
1189 0 : q = twopi/(ln_rho_max - ln_rho_min)
1190 0 : dq = -f*q
1191 : END IF
1192 :
1193 0 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1194 0 : 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 0 : !$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 0 : CALL timestop(handle)
1231 :
1232 0 : 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 236 : 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 708 : 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 236 : CALL timeset(routineN, handle)
1266 :
1267 236 : CPASSERT(ASSOCIATED(pw_env))
1268 :
1269 : ! Perform method specific setup
1270 314 : SELECT CASE (method)
1271 : CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7)
1272 78 : NULLIFY (rs_desc)
1273 78 : rs_grid_section => section_vals_get_subs_vals(input, "DFT%MGRID%RS_GRID")
1274 0 : SELECT CASE (method)
1275 : CASE (sccs_derivative_cd3)
1276 0 : border_points = 1
1277 : CASE (sccs_derivative_cd5)
1278 78 : border_points = 2
1279 : CASE (sccs_derivative_cd7)
1280 78 : border_points = 3
1281 : END SELECT
1282 : CALL init_input_type(input_settings, 2*border_points + 1, rs_grid_section, &
1283 78 : 1, [-1, -1, -1])
1284 : CALL rs_grid_create_descriptor(rs_desc, f%pw_grid, input_settings, &
1285 78 : border_points=border_points)
1286 1638 : ALLOCATE (rs_grid)
1287 78 : CALL rs_grid_create(rs_grid, rs_desc)
1288 : !MK CALL rs_grid_print(rs_grid, 6)
1289 : CASE (sccs_derivative_fft)
1290 632 : lb(1:3) = f%pw_grid%bounds_local(1, 1:3)
1291 632 : ub(1:3) = f%pw_grid%bounds_local(2, 1:3)
1292 158 : NULLIFY (auxbas_pw_pool)
1293 158 : 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 710 : DO i = 1, SIZE(work_g1d)
1296 474 : CALL auxbas_pw_pool%create_pw(work_g1d(i))
1297 : END DO
1298 : END SELECT
1299 :
1300 : ! Calculate the derivatives
1301 0 : SELECT CASE (method)
1302 : CASE (sccs_derivative_cd3)
1303 0 : CALL derive_fdm_cd3(f, df, rs_grid)
1304 : CASE (sccs_derivative_cd5)
1305 78 : CALL derive_fdm_cd5(f, df, rs_grid)
1306 : CASE (sccs_derivative_cd7)
1307 0 : CALL derive_fdm_cd7(f, df, rs_grid)
1308 : CASE (sccs_derivative_fft)
1309 : ! FFT
1310 158 : CALL pw_transfer(f, work_g1d(1))
1311 632 : DO i = 1, 3
1312 474 : n(:) = 0
1313 474 : n(i) = 1
1314 474 : CALL pw_copy(work_g1d(1), work_g1d(2))
1315 474 : CALL pw_derive(work_g1d(2), n(:))
1316 632 : CALL pw_transfer(work_g1d(2), df(i))
1317 : END DO
1318 : CASE DEFAULT
1319 236 : CPABORT("Invalid derivative method for SCCS specified")
1320 : END SELECT
1321 :
1322 : ! Perform method specific cleanup
1323 78 : SELECT CASE (method)
1324 : CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7)
1325 78 : CALL rs_grid_release(rs_grid)
1326 78 : DEALLOCATE (rs_grid)
1327 78 : CALL rs_grid_release_descriptor(rs_desc)
1328 : CASE (sccs_derivative_fft)
1329 710 : DO i = 1, SIZE(work_g1d)
1330 474 : CALL auxbas_pw_pool%give_back_pw(work_g1d(i))
1331 : END DO
1332 : END SELECT
1333 :
1334 236 : CALL timestop(handle)
1335 :
1336 944 : 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 96 : 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 96 : CALL timeset(routineN, handle)
1372 :
1373 96 : e = epsilon_solvent - 1.0_dp
1374 96 : f = LOG(epsilon_solvent)/twopi
1375 96 : diff = rho_max - rho_min
1376 96 : IF (diff < SQRT(rhotol)) CPABORT("SCCS: Difference between rho(min) and rho(max) is too small")
1377 96 : IF (rho_min >= rhotol) THEN
1378 96 : ln_rho_max = LOG(rho_max)
1379 96 : ln_rho_min = LOG(rho_min)
1380 96 : q = twopi/(ln_rho_max - ln_rho_min)
1381 : END IF
1382 :
1383 384 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1384 384 : 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 96 : !$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 96 : CALL timestop(handle)
1415 :
1416 96 : 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 22 : 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 22 : CALL timeset(routineN, handle)
1451 :
1452 22 : e = epsilon_solvent - 1.0_dp
1453 22 : f = 0.5_dp*e
1454 22 : q = 1.0_dp/rho_zero
1455 22 : twobeta = 2.0_dp*beta
1456 :
1457 88 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1458 88 : 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 22 : !$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 22 : CALL timestop(handle)
1487 :
1488 22 : 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 53 : 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 53 : IF (output_unit > 0) THEN
1507 53 : CPASSERT(ASSOCIATED(energy))
1508 53 : CPASSERT(ASSOCIATED(sccs_control))
1509 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
1510 53 : "SCCS| Hartree energy of solute and solvent [Hartree]", energy%sccs_hartree, &
1511 106 : "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 53 : "SCCS| Polarisation energy [Hartree]", energy%sccs_pol, &
1514 53 : "SCCS| [kcal/mol]", &
1515 53 : cp_unit_from_cp2k(energy%sccs_pol, "kcalmol"), &
1516 53 : "SCCS| Cavitation energy [Hartree]", energy%sccs_cav, &
1517 53 : "SCCS| [kcal/mol]", &
1518 53 : cp_unit_from_cp2k(energy%sccs_cav, "kcalmol"), &
1519 53 : "SCCS| Dispersion free energy [Hartree]", energy%sccs_dis, &
1520 53 : "SCCS| [kcal/mol]", &
1521 53 : cp_unit_from_cp2k(energy%sccs_dis, "kcalmol"), &
1522 53 : "SCCS| Repulsion free energy [Hartree]", energy%sccs_rep, &
1523 53 : "SCCS| [kcal/mol]", &
1524 53 : cp_unit_from_cp2k(energy%sccs_rep, "kcalmol"), &
1525 53 : "SCCS| Solvation free energy [Hartree]", energy%sccs_sol, &
1526 53 : "SCCS| [kcal/mol]", &
1527 106 : cp_unit_from_cp2k(energy%sccs_sol, "kcalmol")
1528 : END IF
1529 :
1530 53 : END SUBROUTINE print_sccs_results
1531 :
1532 : END MODULE qs_sccs
|
