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