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 : USE input_section_types, ONLY: section_get_ivals,&
52 : section_get_lval,&
53 : section_vals_get_subs_vals,&
54 : section_vals_type
55 : USE kinds, ONLY: default_path_length,&
56 : default_string_length,&
57 : dp,&
58 : int_8
59 : USE mathconstants, ONLY: fourpi,&
60 : twopi
61 : USE message_passing, ONLY: mp_para_env_type
62 : USE particle_list_types, ONLY: particle_list_type
63 : USE pw_env_types, ONLY: pw_env_get,&
64 : pw_env_type
65 : USE pw_methods, ONLY: pw_axpy,&
66 : pw_copy,&
67 : pw_derive,&
68 : pw_integral_ab,&
69 : pw_integrate_function,&
70 : pw_transfer,&
71 : pw_zero
72 : USE pw_poisson_methods, ONLY: pw_poisson_solve
73 : USE pw_poisson_types, ONLY: pw_poisson_analytic,&
74 : pw_poisson_mt,&
75 : pw_poisson_type
76 : USE pw_pool_types, ONLY: pw_pool_p_type,&
77 : pw_pool_type
78 : USE pw_types, ONLY: pw_c1d_gs_type,&
79 : pw_r3d_rs_type
80 : USE qs_energy_types, ONLY: qs_energy_type
81 : USE qs_environment_types, ONLY: get_qs_env,&
82 : qs_environment_type
83 : USE qs_rho_types, ONLY: qs_rho_get,&
84 : qs_rho_type
85 : USE qs_scf_types, ONLY: qs_scf_env_type
86 : USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
87 : realspace_grid_input_type,&
88 : realspace_grid_type,&
89 : rs_grid_create,&
90 : rs_grid_create_descriptor,&
91 : rs_grid_release,&
92 : rs_grid_release_descriptor
93 : USE rs_methods, ONLY: derive_fdm_cd3,&
94 : derive_fdm_cd5,&
95 : derive_fdm_cd7
96 : #include "./base/base_uses.f90"
97 :
98 : IMPLICIT NONE
99 :
100 : PRIVATE
101 :
102 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_sccs'
103 :
104 : PUBLIC :: print_sccs_results, sccs
105 :
106 : CONTAINS
107 :
108 : ! **************************************************************************************************
109 : !> \brief Self-consistent continuum solvation (SCCS) model implementation
110 : !> \param qs_env ...
111 : !> \param rho_tot_gspace ...
112 : !> \param v_hartree_gspace ...
113 : !> \param v_sccs ...
114 : !> \param h_stress ...
115 : !> \par History:
116 : !> - Creation (10.10.2013,MK)
117 : !> \author Matthias Krack (MK)
118 : !> \version 1.0
119 : ! **************************************************************************************************
120 :
121 132 : SUBROUTINE sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs, h_stress)
122 :
123 : TYPE(qs_environment_type), POINTER :: qs_env
124 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_tot_gspace, v_hartree_gspace
125 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_sccs
126 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
127 : OPTIONAL :: h_stress
128 :
129 : CHARACTER(LEN=*), PARAMETER :: routineN = 'sccs'
130 : REAL(KIND=dp), PARAMETER :: epstol = 1.0E-8_dp
131 :
132 : CHARACTER(LEN=4*default_string_length) :: message, my_pos_cube
133 : CHARACTER(LEN=default_path_length) :: cube_path, filename, mpi_filename, &
134 : print_path
135 : INTEGER :: cube_unit, handle, i, ispin, iter, j, k, &
136 : nspin, output_unit, print_level
137 : INTEGER(KIND=int_8) :: ngpts
138 : INTEGER, DIMENSION(3) :: lb, ub
139 : LOGICAL :: append_cube, calculate_stress_tensor, &
140 : do_kpoints, mpi_io, should_output
141 : REAL(KIND=dp) :: cavity_surface, cavity_volume, cell_volume, dphi2, dvol, e_tot, &
142 : epsilon_solvent, f, polarisation_charge, rho_delta, rho_delta_avg, rho_delta_max, &
143 : rho_iter_new, tot_rho_elec, tot_rho_solute
144 : REAL(KIND=dp), DIMENSION(3) :: abc
145 : TYPE(cp_logger_type), POINTER :: logger
146 : TYPE(cp_subsys_type), POINTER :: cp_subsys
147 : TYPE(dft_control_type), POINTER :: dft_control
148 : TYPE(mp_para_env_type), POINTER :: para_env
149 : TYPE(particle_list_type), POINTER :: particles
150 : TYPE(pw_env_type), POINTER :: pw_env
151 : TYPE(pw_poisson_type), POINTER :: poisson_env
152 132 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
153 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
154 : TYPE(pw_r3d_rs_type) :: deps_elec, eps_elec
155 924 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: dln_eps_elec, dphi_tot
156 132 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_pw_r
157 : TYPE(pw_r3d_rs_type), POINTER :: rho_pw_r_sccs
158 : TYPE(qs_energy_type), POINTER :: energy
159 : TYPE(qs_rho_type), POINTER :: rho
160 : TYPE(qs_scf_env_type), POINTER :: scf_env
161 : TYPE(sccs_control_type), POINTER :: sccs_control
162 : TYPE(section_vals_type), POINTER :: input
163 :
164 132 : CALL timeset(routineN, handle)
165 :
166 132 : NULLIFY (auxbas_pw_pool)
167 132 : NULLIFY (cp_subsys)
168 132 : NULLIFY (dft_control)
169 132 : NULLIFY (energy)
170 132 : NULLIFY (input)
171 132 : NULLIFY (logger)
172 132 : NULLIFY (para_env)
173 132 : NULLIFY (particles)
174 132 : NULLIFY (poisson_env)
175 132 : NULLIFY (pw_env)
176 132 : NULLIFY (pw_pools)
177 132 : NULLIFY (rho)
178 132 : NULLIFY (sccs_control)
179 132 : NULLIFY (scf_env)
180 :
181 : ! Load data from Quickstep environment
182 : CALL get_qs_env(qs_env=qs_env, &
183 : cp_subsys=cp_subsys, &
184 : dft_control=dft_control, &
185 : do_kpoints=do_kpoints, &
186 : energy=energy, &
187 : input=input, &
188 : para_env=para_env, &
189 : pw_env=pw_env, &
190 : rho=rho, &
191 132 : scf_env=scf_env)
192 132 : CALL cp_subsys_get(cp_subsys, particles=particles)
193 :
194 132 : sccs_control => dft_control%sccs_control
195 :
196 132 : CPASSERT(ASSOCIATED(qs_env))
197 :
198 132 : IF (do_kpoints) THEN
199 0 : CPWARN("SCCS with k-points has not yet been fully validated")
200 : END IF
201 :
202 132 : IF (PRESENT(h_stress)) THEN
203 0 : calculate_stress_tensor = .TRUE.
204 0 : h_stress(:, :) = 0.0_dp
205 0 : CPWARN("The stress tensor for SCCS has not yet been fully validated")
206 : ELSE
207 : calculate_stress_tensor = .FALSE.
208 : END IF
209 :
210 : ! Get access to the PW grid pool
211 : CALL pw_env_get(pw_env, &
212 : auxbas_pw_pool=auxbas_pw_pool, &
213 : pw_pools=pw_pools, &
214 132 : poisson_env=poisson_env)
215 :
216 132 : CALL pw_zero(v_sccs)
217 :
218 : ! Calculate no SCCS contribution, if the requested SCF convergence threshold is not reached yet
219 132 : IF (.NOT. sccs_control%sccs_activated) THEN
220 50 : IF (sccs_control%eps_scf > 0.0_dp) THEN
221 48 : IF ((scf_env%iter_delta > sccs_control%eps_scf) .OR. &
222 : ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
223 : (qs_env%scf_env%iter_count <= 1))) THEN
224 40 : IF (calculate_stress_tensor) THEN
225 : ! Request also the calculation of the stress tensor contribution
226 : CALL pw_poisson_solve(poisson_env=poisson_env, &
227 : density=rho_tot_gspace, &
228 : ehartree=energy%hartree, &
229 : vhartree=v_hartree_gspace, &
230 0 : h_stress=h_stress)
231 : ELSE
232 : CALL pw_poisson_solve(poisson_env=poisson_env, &
233 : density=rho_tot_gspace, &
234 : ehartree=energy%hartree, &
235 40 : vhartree=v_hartree_gspace)
236 : END IF
237 40 : energy%sccs_pol = 0.0_dp
238 40 : energy%sccs_cav = 0.0_dp
239 40 : energy%sccs_dis = 0.0_dp
240 40 : energy%sccs_rep = 0.0_dp
241 40 : energy%sccs_sol = 0.0_dp
242 40 : energy%sccs_hartree = energy%hartree
243 40 : CALL timestop(handle)
244 40 : RETURN
245 : END IF
246 : END IF
247 10 : sccs_control%sccs_activated = .TRUE.
248 : END IF
249 :
250 92 : nspin = dft_control%nspins
251 :
252 : ! Manage print output control
253 92 : logger => cp_get_default_logger()
254 92 : print_level = logger%iter_info%print_level
255 92 : print_path = "DFT%PRINT%SCCS"
256 : should_output = (BTEST(cp_print_key_should_output(logger%iter_info, input, &
257 92 : TRIM(print_path)), cp_p_file))
258 : output_unit = cp_print_key_unit_nr(logger, input, TRIM(print_path), &
259 : extension=".sccs", &
260 : ignore_should_output=should_output, &
261 92 : log_filename=.FALSE.)
262 :
263 : ! Get rho
264 : CALL qs_rho_get(rho_struct=rho, &
265 : rho_r=rho_pw_r, &
266 92 : rho_r_sccs=rho_pw_r_sccs)
267 :
268 : ! Retrieve the last rho_iter from the previous SCCS cycle if available
269 92 : CPASSERT(ASSOCIATED(rho_pw_r_sccs))
270 :
271 : ! Retrieve the total electronic density in r-space
272 : BLOCK
273 : TYPE(pw_r3d_rs_type) :: rho_elec
274 92 : CALL auxbas_pw_pool%create_pw(rho_elec)
275 :
276 : ! Retrieve grid parameters
277 92 : ngpts = rho_elec%pw_grid%ngpts
278 92 : dvol = rho_elec%pw_grid%dvol
279 92 : cell_volume = rho_elec%pw_grid%vol
280 368 : abc(1:3) = REAL(rho_elec%pw_grid%npts(1:3), KIND=dp)*rho_elec%pw_grid%dr(1:3)
281 368 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
282 368 : ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
283 :
284 92 : CALL pw_copy(rho_pw_r(1), rho_elec)
285 92 : DO ispin = 2, nspin
286 92 : CALL pw_axpy(rho_pw_r(ispin), rho_elec)
287 : END DO
288 92 : tot_rho_elec = pw_integrate_function(rho_elec)
289 :
290 : ! Calculate the dielectric (smoothed) function of rho_elec in r-space
291 92 : CALL auxbas_pw_pool%create_pw(eps_elec)
292 92 : CALL auxbas_pw_pool%create_pw(deps_elec)
293 : ! Relative permittivity or dielectric constant of the solvent (medium)
294 92 : epsilon_solvent = sccs_control%epsilon_solvent
295 162 : SELECT CASE (sccs_control%method_id)
296 : CASE (sccs_andreussi)
297 : CALL andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%rho_max, &
298 70 : sccs_control%rho_min)
299 : CASE (sccs_fattebert_gygi)
300 : CALL fattebert_gygi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%beta, &
301 22 : sccs_control%rho_zero)
302 : CASE DEFAULT
303 92 : CPABORT("Invalid method specified for SCCS model")
304 : END SELECT
305 :
306 : ! Optional output of the dielectric function in cube file format
307 92 : filename = "DIELECTRIC_FUNCTION"
308 92 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
309 92 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), &
310 : cp_p_file)) THEN
311 4 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
312 4 : my_pos_cube = "REWIND"
313 4 : IF (append_cube) my_pos_cube = "APPEND"
314 4 : mpi_io = .TRUE.
315 : cube_unit = cp_print_key_unit_nr(logger, input, TRIM(cube_path), &
316 : extension=".cube", middle_name=TRIM(filename), &
317 : file_position=my_pos_cube, log_filename=.FALSE., &
318 4 : mpi_io=mpi_io, fout=mpi_filename)
319 4 : IF (output_unit > 0) THEN
320 2 : IF (.NOT. mpi_io) THEN
321 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
322 : ELSE
323 2 : filename = mpi_filename
324 : END IF
325 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
326 2 : "SCCS| The dielectric function is written in cube file format to the file:", &
327 4 : "SCCS| "//TRIM(filename)
328 : END IF
329 : CALL cp_pw_to_cube(eps_elec, cube_unit, TRIM(filename), particles=particles, &
330 : stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), &
331 4 : mpi_io=mpi_io)
332 4 : CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io)
333 : END IF
334 :
335 : ! Calculate the (quantum) volume and surface of the solute cavity
336 92 : cavity_surface = 0.0_dp
337 92 : cavity_volume = 0.0_dp
338 :
339 92 : IF (ABS(epsilon_solvent - 1.0_dp) > epstol) THEN
340 :
341 : BLOCK
342 : TYPE(pw_r3d_rs_type) :: theta, norm_drho_elec
343 368 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: drho_elec
344 92 : CALL auxbas_pw_pool%create_pw(theta)
345 92 : CALL pw_zero(theta)
346 :
347 : ! Calculate the (quantum) volume of the solute cavity
348 92 : f = 1.0_dp/(epsilon_solvent - 1.0_dp)
349 : !$OMP PARALLEL DO DEFAULT(NONE) &
350 : !$OMP PRIVATE(i,j,k) &
351 92 : !$OMP SHARED(epsilon_solvent,eps_elec,f,lb,theta,ub)
352 : DO k = lb(3), ub(3)
353 : DO j = lb(2), ub(2)
354 : DO i = lb(1), ub(1)
355 : theta%array(i, j, k) = f*(epsilon_solvent - eps_elec%array(i, j, k))
356 : END DO
357 : END DO
358 : END DO
359 : !$OMP END PARALLEL DO
360 92 : cavity_volume = pw_integrate_function(theta)
361 :
362 : ! Calculate the derivative of the electronic density in r-space
363 : ! TODO: Could be retrieved from the QS environment
364 368 : DO i = 1, 3
365 368 : CALL auxbas_pw_pool%create_pw(drho_elec(i))
366 : END DO
367 92 : CALL derive(rho_elec, drho_elec, sccs_derivative_fft, pw_env, input)
368 :
369 92 : CALL auxbas_pw_pool%create_pw(norm_drho_elec)
370 :
371 : ! Calculate the norm of the gradient of the electronic density in r-space
372 : !$OMP PARALLEL DO DEFAULT(NONE) &
373 : !$OMP PRIVATE(i,j,k) &
374 92 : !$OMP SHARED(drho_elec,lb,norm_drho_elec,ub)
375 : DO k = lb(3), ub(3)
376 : DO j = lb(2), ub(2)
377 : DO i = lb(1), ub(1)
378 : norm_drho_elec%array(i, j, k) = SQRT(drho_elec(1)%array(i, j, k)* &
379 : drho_elec(1)%array(i, j, k) + &
380 : drho_elec(2)%array(i, j, k)* &
381 : drho_elec(2)%array(i, j, k) + &
382 : drho_elec(3)%array(i, j, k)* &
383 : drho_elec(3)%array(i, j, k))
384 : END DO
385 : END DO
386 : END DO
387 : !$OMP END PARALLEL DO
388 :
389 : ! Optional output of the norm of the density gradient in cube file format
390 92 : filename = "DENSITY_GRADIENT"
391 92 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
392 92 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), &
393 : cp_p_file)) THEN
394 0 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
395 0 : my_pos_cube = "REWIND"
396 0 : IF (append_cube) my_pos_cube = "APPEND"
397 0 : mpi_io = .TRUE.
398 : cube_unit = cp_print_key_unit_nr(logger, input, TRIM(cube_path), &
399 : extension=".cube", middle_name=TRIM(filename), &
400 : file_position=my_pos_cube, log_filename=.FALSE., &
401 0 : mpi_io=mpi_io, fout=mpi_filename)
402 0 : IF (output_unit > 0) THEN
403 0 : IF (.NOT. mpi_io) THEN
404 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
405 : ELSE
406 0 : filename = mpi_filename
407 : END IF
408 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
409 0 : "SCCS| The norm of the density gradient is written in cube file format to the file:", &
410 0 : "SCCS| "//TRIM(filename)
411 : END IF
412 : CALL cp_pw_to_cube(norm_drho_elec, cube_unit, TRIM(filename), particles=particles, &
413 : stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), &
414 0 : mpi_io=mpi_io)
415 0 : CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io)
416 : END IF
417 :
418 : ! Calculate the (quantum) surface of the solute cavity
419 162 : SELECT CASE (sccs_control%method_id)
420 : CASE (sccs_andreussi)
421 : CALL surface_andreussi(rho_elec, norm_drho_elec, theta, epsilon_solvent, &
422 : sccs_control%rho_max, sccs_control%rho_min, &
423 70 : sccs_control%delta_rho)
424 : CASE (sccs_fattebert_gygi)
425 : CALL surface_fattebert_gygi(rho_elec, norm_drho_elec, theta, epsilon_solvent, &
426 : sccs_control%beta, sccs_control%rho_zero, &
427 22 : sccs_control%delta_rho)
428 : CASE DEFAULT
429 92 : CPABORT("Invalid method specified for SCCS model")
430 : END SELECT
431 92 : cavity_surface = pw_integrate_function(theta)
432 :
433 : ! Release storage
434 92 : CALL auxbas_pw_pool%give_back_pw(theta)
435 92 : CALL auxbas_pw_pool%give_back_pw(norm_drho_elec)
436 460 : DO i = 1, 3
437 368 : CALL auxbas_pw_pool%give_back_pw(drho_elec(i))
438 : END DO
439 : END BLOCK
440 :
441 : END IF ! epsilon_solvent > 1
442 :
443 92 : CALL auxbas_pw_pool%give_back_pw(rho_elec)
444 : END BLOCK
445 :
446 : BLOCK
447 : TYPE(pw_r3d_rs_type) :: rho_tot, phi_tot, rho_solute, rho_tot_zero
448 : ! Retrieve the total charge density (core + elec) of the solute in r-space
449 92 : CALL auxbas_pw_pool%create_pw(rho_solute)
450 92 : CALL pw_zero(rho_solute)
451 92 : CALL pw_transfer(rho_tot_gspace, rho_solute)
452 92 : tot_rho_solute = pw_integrate_function(rho_solute)
453 :
454 : ! Check total charge
455 92 : IF (ABS(tot_rho_solute) >= 1.0E-6_dp) THEN
456 92 : IF ((poisson_env%parameters%solver /= pw_poisson_analytic) .AND. &
457 : (poisson_env%parameters%solver /= pw_poisson_mt)) THEN
458 : WRITE (UNIT=message, FMT="(A,SP,F0.6,A)") &
459 92 : "The system (solute) has a non-negligible charge of ", -tot_rho_solute, &
460 : ". It is recommended to use non-periodic boundary conditions (PERIODIC none) "// &
461 184 : "combined with an appropriate Poisson solver (POISSON_SOLVER MT or analytic)"
462 92 : CPWARN(message)
463 : END IF
464 : END IF
465 :
466 : ! Reassign work storage to rho_tot_zero, because rho_elec is no longer needed
467 92 : CALL auxbas_pw_pool%create_pw(rho_tot_zero)
468 :
469 : ! Build the initial (rho_iter = 0) total charge density (solute plus polarisation) in r-space
470 : ! eps_elec <- ln(eps_elec)
471 : !$OMP PARALLEL DO DEFAULT(NONE) &
472 : !$OMP PRIVATE(i,j,k) &
473 : !$OMP SHARED(eps_elec,lb,message,output_unit,para_env,ub) &
474 92 : !$OMP SHARED(rho_solute,rho_tot_zero)
475 : DO k = lb(3), ub(3)
476 : DO j = lb(2), ub(2)
477 : DO i = lb(1), ub(1)
478 : IF (eps_elec%array(i, j, k) < 1.0_dp) THEN
479 : WRITE (UNIT=message, FMT="(A,ES12.3,A,3(I0,A))") &
480 : "SCCS| Invalid dielectric function value ", eps_elec%array(i, j, k), &
481 : " encountered at grid point (", i, ",", j, ",", k, ")"
482 : CPABORT(message)
483 : END IF
484 : rho_tot_zero%array(i, j, k) = rho_solute%array(i, j, k)/eps_elec%array(i, j, k)
485 : eps_elec%array(i, j, k) = LOG(eps_elec%array(i, j, k))
486 : END DO
487 : END DO
488 : END DO
489 : !$OMP END PARALLEL DO
490 :
491 : ! Build the derivative of LOG(eps_elec)
492 368 : DO i = 1, 3
493 276 : CALL auxbas_pw_pool%create_pw(dln_eps_elec(i))
494 368 : CALL pw_zero(dln_eps_elec(i))
495 : END DO
496 92 : CALL derive(eps_elec, dln_eps_elec, sccs_control%derivative_method, pw_env, input)
497 92 : CALL auxbas_pw_pool%give_back_pw(eps_elec)
498 :
499 : ! Print header for the SCCS cycle
500 92 : IF (should_output .AND. (output_unit > 0)) THEN
501 46 : IF (print_level > low_print_level) THEN
502 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.12)") &
503 46 : "SCCS| Total electronic charge density ", -tot_rho_elec, &
504 92 : "SCCS| Total charge density (solute) ", -tot_rho_solute
505 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.3)") &
506 46 : "SCCS| Volume of the cell [bohr^3]", cell_volume, &
507 46 : "SCCS| [angstrom^3]", &
508 92 : cp_unit_from_cp2k(cell_volume, "angstrom^3")
509 46 : IF (ABS(epsilon_solvent - 1.0_dp) > epstol) THEN
510 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.3)") &
511 46 : "SCCS| Volume of the solute cavity [bohr^3]", cavity_volume, &
512 46 : "SCCS| [angstrom^3]", &
513 46 : cp_unit_from_cp2k(cavity_volume, "angstrom^3"), &
514 46 : "SCCS| Surface of the solute cavity [bohr^2]", cavity_surface, &
515 46 : "SCCS| [angstrom^2]", &
516 92 : cp_unit_from_cp2k(cavity_surface, "angstrom^2")
517 : END IF
518 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
519 46 : "SCCS|", &
520 92 : "SCCS| Step Average residual Maximum residual E(Hartree) [a.u.]"
521 : END IF
522 : END IF
523 :
524 : ! Get storage for the derivative of the total potential (field) in r-space
525 368 : DO i = 1, 3
526 368 : CALL auxbas_pw_pool%create_pw(dphi_tot(i))
527 : END DO
528 :
529 : ! Initialise the total charge density in r-space rho_tot with rho_tot_zero + rho_iter_zero
530 92 : CALL auxbas_pw_pool%create_pw(rho_tot)
531 92 : CALL pw_copy(rho_tot_zero, rho_tot)
532 92 : CALL pw_axpy(rho_pw_r_sccs, rho_tot)
533 :
534 92 : CALL auxbas_pw_pool%create_pw(phi_tot)
535 92 : CALL pw_zero(phi_tot)
536 :
537 : ! Main SCCS iteration loop
538 92 : iter = 0
539 :
540 : iter_loop: DO
541 :
542 : ! Increment iteration counter
543 294 : iter = iter + 1
544 :
545 : ! Check if the requested maximum number of SCCS iterations is reached
546 294 : IF (iter > sccs_control%max_iter) THEN
547 0 : IF (output_unit > 0) THEN
548 : WRITE (UNIT=output_unit, FMT="(T3,A,/,T3,A,I0,A)") &
549 0 : "SCCS| Maximum number of SCCS iterations reached", &
550 0 : "SCCS| Iteration cycle did not converge in ", sccs_control%max_iter, " steps"
551 : ELSE
552 : WRITE (UNIT=message, FMT="(A,I0,A)") &
553 0 : "The SCCS iteration cycle did not converge in ", sccs_control%max_iter, " steps"
554 0 : CPWARN(message)
555 : END IF
556 : EXIT iter_loop
557 : END IF
558 :
559 : ! Calculate derivative of the current total potential in r-space
560 : CALL pw_poisson_solve(poisson_env=poisson_env, &
561 : density=rho_tot, &
562 : vhartree=phi_tot, &
563 294 : dvhartree=dphi_tot)
564 294 : energy%sccs_hartree = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
565 :
566 : ! Update total charge density (solute plus polarisation) in r-space
567 : ! based on the iterated polarisation charge density
568 294 : f = 1.0_dp/fourpi
569 294 : rho_delta_avg = 0.0_dp
570 294 : rho_delta_max = 0.0_dp
571 : !$OMP PARALLEL DO DEFAULT(NONE) &
572 : !$OMP PRIVATE(i,j,k,rho_delta,rho_iter_new) &
573 : !$OMP SHARED(dln_eps_elec,dphi_tot,f,lb,rho_pw_r_sccs,ub) &
574 : !$OMP SHARED(rho_tot,rho_tot_zero,sccs_control) &
575 : !$OMP REDUCTION(+:rho_delta_avg) &
576 294 : !$OMP REDUCTION(MAX:rho_delta_max)
577 : DO k = lb(3), ub(3)
578 : DO j = lb(2), ub(2)
579 : DO i = lb(1), ub(1)
580 : rho_iter_new = (dln_eps_elec(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + &
581 : dln_eps_elec(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + &
582 : dln_eps_elec(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k))*f
583 : rho_iter_new = rho_pw_r_sccs%array(i, j, k) + &
584 : sccs_control%mixing*(rho_iter_new - rho_pw_r_sccs%array(i, j, k))
585 : rho_delta = ABS(rho_iter_new - rho_pw_r_sccs%array(i, j, k))
586 : rho_delta_max = MAX(rho_delta, rho_delta_max)
587 : rho_delta_avg = rho_delta_avg + rho_delta
588 : rho_tot%array(i, j, k) = rho_tot_zero%array(i, j, k) + rho_iter_new
589 : rho_pw_r_sccs%array(i, j, k) = rho_iter_new
590 : END DO
591 : END DO
592 : END DO
593 : !$OMP END PARALLEL DO
594 :
595 294 : CALL para_env%sum(rho_delta_avg)
596 294 : rho_delta_avg = rho_delta_avg/REAL(ngpts, KIND=dp)
597 294 : CALL para_env%max(rho_delta_max)
598 :
599 294 : IF (should_output .AND. (output_unit > 0)) THEN
600 147 : IF (print_level > low_print_level) THEN
601 147 : IF ((ABS(rho_delta_avg) < 1.0E-8_dp) .OR. &
602 147 : (ABS(rho_delta_avg) >= 1.0E5_dp)) THEN
603 : WRITE (UNIT=output_unit, FMT="(T3,A,I6,4X,ES16.4,4X,ES16.4,1X,F25.12)") &
604 0 : "SCCS| ", iter, rho_delta_avg, rho_delta_max, energy%sccs_hartree
605 : ELSE
606 : WRITE (UNIT=output_unit, FMT="(T3,A,I6,4X,F16.8,4X,F16.8,1X,F25.12)") &
607 147 : "SCCS| ", iter, rho_delta_avg, rho_delta_max, energy%sccs_hartree
608 : END IF
609 : END IF
610 : END IF
611 :
612 : ! Check if the SCCS iteration cycle is converged to the requested tolerance
613 294 : IF (rho_delta_max <= sccs_control%eps_sccs) THEN
614 92 : IF (should_output .AND. (output_unit > 0)) THEN
615 : WRITE (UNIT=output_unit, FMT="(T3,A,I0,A)") &
616 46 : "SCCS| Iteration cycle converged in ", iter, " steps"
617 : END IF
618 : EXIT iter_loop
619 : END IF
620 :
621 : END DO iter_loop
622 :
623 : ! Release work storage which is no longer needed
624 92 : CALL auxbas_pw_pool%give_back_pw(rho_tot_zero)
625 368 : DO i = 1, 3
626 368 : CALL auxbas_pw_pool%give_back_pw(dln_eps_elec(i))
627 : END DO
628 :
629 : ! Optional output of the total charge density in cube file format
630 92 : filename = "TOTAL_CHARGE_DENSITY"
631 92 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
632 92 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), cp_p_file)) THEN
633 0 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
634 0 : my_pos_cube = "REWIND"
635 0 : IF (append_cube) my_pos_cube = "APPEND"
636 0 : mpi_io = .TRUE.
637 : cube_unit = cp_print_key_unit_nr(logger, input, TRIM(cube_path), &
638 : extension=".cube", middle_name=TRIM(filename), &
639 : file_position=my_pos_cube, log_filename=.FALSE., &
640 0 : mpi_io=mpi_io, fout=mpi_filename)
641 0 : IF (output_unit > 0) THEN
642 0 : IF (.NOT. mpi_io) THEN
643 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
644 : ELSE
645 0 : filename = mpi_filename
646 : END IF
647 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
648 0 : "SCCS| The total SCCS charge density is written in cube file format to the file:", &
649 0 : "SCCS| "//TRIM(filename)
650 : END IF
651 : CALL cp_pw_to_cube(rho_tot, cube_unit, TRIM(filename), particles=particles, &
652 0 : stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), mpi_io=mpi_io)
653 0 : CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io)
654 : END IF
655 :
656 : ! Calculate the total SCCS Hartree energy, potential, and its
657 : ! derivatives of the solute and the implicit solvent
658 92 : CALL pw_transfer(rho_tot, rho_tot_gspace)
659 92 : IF (calculate_stress_tensor) THEN
660 : ! Request also the calculation of the stress tensor contribution
661 : CALL pw_poisson_solve(poisson_env=poisson_env, &
662 : density=rho_tot_gspace, &
663 : ehartree=e_tot, &
664 : vhartree=v_hartree_gspace, &
665 : dvhartree=dphi_tot, &
666 0 : h_stress=h_stress)
667 : ELSE
668 : CALL pw_poisson_solve(poisson_env=poisson_env, &
669 : density=rho_tot_gspace, &
670 : ehartree=e_tot, &
671 : vhartree=v_hartree_gspace, &
672 92 : dvhartree=dphi_tot)
673 : END IF
674 92 : CALL pw_transfer(v_hartree_gspace, phi_tot)
675 92 : energy%sccs_hartree = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
676 :
677 : ! Calculate the Hartree energy and potential of the solute only
678 : BLOCK
679 : TYPE(pw_r3d_rs_type) :: phi_solute
680 92 : CALL auxbas_pw_pool%create_pw(phi_solute)
681 92 : CALL pw_zero(phi_solute)
682 : CALL pw_poisson_solve(poisson_env=poisson_env, &
683 : density=rho_solute, &
684 : ehartree=energy%hartree, &
685 92 : vhartree=phi_solute)
686 :
687 : ! Calculate the polarisation potential (store it in phi_tot)
688 : ! phi_pol = phi_tot - phi_solute
689 92 : CALL pw_axpy(phi_solute, phi_tot, alpha=-1.0_dp)
690 92 : CALL auxbas_pw_pool%give_back_pw(phi_solute)
691 : END BLOCK
692 :
693 : ! Optional output of the SCCS polarisation potential in cube file format
694 92 : filename = "POLARISATION_POTENTIAL"
695 92 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
696 92 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), &
697 : cp_p_file)) THEN
698 4 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
699 4 : my_pos_cube = "REWIND"
700 4 : IF (append_cube) my_pos_cube = "APPEND"
701 4 : mpi_io = .TRUE.
702 : cube_unit = cp_print_key_unit_nr(logger, input, TRIM(cube_path), &
703 : extension=".cube", middle_name=TRIM(filename), &
704 : file_position=my_pos_cube, log_filename=.FALSE., &
705 4 : mpi_io=mpi_io, fout=mpi_filename)
706 4 : IF (output_unit > 0) THEN
707 2 : IF (.NOT. mpi_io) THEN
708 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
709 : ELSE
710 2 : filename = mpi_filename
711 : END IF
712 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
713 2 : "SCCS| The SCCS polarisation potential is written in cube file format to the file:", &
714 4 : "SCCS| "//TRIM(filename)
715 : END IF
716 : CALL cp_pw_to_cube(phi_tot, cube_unit, TRIM(filename), particles=particles, &
717 4 : stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), mpi_io=mpi_io)
718 4 : CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io)
719 : END IF
720 :
721 : ! Calculate the polarisation charge (store it in rho_tot)
722 : ! rho_pol = rho_tot - rho_solute
723 92 : CALL pw_axpy(rho_solute, rho_tot, alpha=-1.0_dp)
724 92 : polarisation_charge = pw_integrate_function(rho_tot)
725 :
726 : ! Optional output of the SCCS polarisation charge in cube file format
727 92 : filename = "POLARISATION_CHARGE_DENSITY"
728 92 : cube_path = TRIM(print_path)//"%"//TRIM(filename)
729 92 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), &
730 : cp_p_file)) THEN
731 0 : append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND")
732 0 : my_pos_cube = "REWIND"
733 0 : IF (append_cube) my_pos_cube = "APPEND"
734 0 : mpi_io = .TRUE.
735 : cube_unit = cp_print_key_unit_nr(logger, input, TRIM(cube_path), &
736 : extension=".cube", middle_name=TRIM(filename), &
737 : file_position=my_pos_cube, log_filename=.FALSE., &
738 0 : mpi_io=mpi_io, fout=mpi_filename)
739 0 : IF (output_unit > 0) THEN
740 0 : IF (.NOT. mpi_io) THEN
741 0 : INQUIRE (UNIT=cube_unit, NAME=filename)
742 : ELSE
743 0 : filename = mpi_filename
744 : END IF
745 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
746 0 : "SCCS| The SCCS polarisation charge density is written in cube file format to the file:", &
747 0 : "SCCS| "//TRIM(filename)
748 : END IF
749 : CALL cp_pw_to_cube(rho_tot, cube_unit, TRIM(filename), particles=particles, &
750 0 : stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), mpi_io=mpi_io)
751 0 : CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io)
752 : END IF
753 :
754 : ! Calculate SCCS polarisation energy
755 92 : energy%sccs_pol = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
756 92 : CALL auxbas_pw_pool%give_back_pw(rho_solute)
757 92 : CALL auxbas_pw_pool%give_back_pw(phi_tot)
758 184 : CALL auxbas_pw_pool%give_back_pw(rho_tot)
759 : END BLOCK
760 :
761 : ! Calculate additional solvation terms
762 92 : energy%sccs_cav = sccs_control%gamma_solvent*cavity_surface
763 92 : energy%sccs_dis = sccs_control%beta_solvent*cavity_volume
764 92 : energy%sccs_rep = sccs_control%alpha_solvent*cavity_surface
765 : ! Calculate solvation free energy: \delta G^el + (alpha + gamma)*S + beta*V
766 92 : energy%sccs_sol = energy%sccs_pol + energy%sccs_rep + energy%sccs_cav + energy%sccs_dis
767 :
768 92 : IF (should_output .AND. (output_unit > 0)) THEN
769 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
770 46 : "SCCS|"
771 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.12)") &
772 46 : "SCCS| Polarisation charge", polarisation_charge
773 : !MK "SCCS| Total interaction energy [a.u.]", e_tot
774 : WRITE (UNIT=output_unit, FMT="(T3,A)") &
775 46 : "SCCS|"
776 46 : CALL print_sccs_results(energy, sccs_control, output_unit)
777 : END IF
778 :
779 : ! Calculate SCCS contribution to the Kohn-Sham potential
780 92 : f = -0.5_dp*dvol/fourpi
781 : !$OMP PARALLEL DO DEFAULT(NONE) &
782 : !$OMP PRIVATE(dphi2,i,j,k) &
783 : !$OMP SHARED(f,deps_elec,dphi_tot) &
784 92 : !$OMP SHARED(lb,ub,v_sccs)
785 : DO k = lb(3), ub(3)
786 : DO j = lb(2), ub(2)
787 : DO i = lb(1), ub(1)
788 : dphi2 = dphi_tot(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + &
789 : dphi_tot(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + &
790 : dphi_tot(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k)
791 : v_sccs%array(i, j, k) = v_sccs%array(i, j, k) + f*deps_elec%array(i, j, k)*dphi2
792 : END DO
793 : END DO
794 : END DO
795 : !$OMP END PARALLEL DO
796 :
797 92 : CALL auxbas_pw_pool%give_back_pw(deps_elec)
798 368 : DO i = 1, 3
799 368 : CALL auxbas_pw_pool%give_back_pw(dphi_tot(i))
800 : END DO
801 :
802 : ! Release the SCCS printout environment
803 : CALL cp_print_key_finished_output(output_unit, logger, input, TRIM(print_path), &
804 92 : ignore_should_output=should_output)
805 :
806 92 : CALL timestop(handle)
807 :
808 264 : END SUBROUTINE sccs
809 :
810 : ! **************************************************************************************************
811 : !> \brief Calculate the smoothed dielectric function of Andreussi et al.
812 : !> \param rho_elec ...
813 : !> \param eps_elec ...
814 : !> \param deps_elec ...
815 : !> \param epsilon_solvent ...
816 : !> \param rho_max ...
817 : !> \param rho_min ...
818 : !> \par History:
819 : !> - Creation (16.10.2013,MK)
820 : !> - Finite difference of isosurfaces implemented (21.12.2013,MK)
821 : !> \author Matthias Krack (MK)
822 : !> \version 1.1
823 : ! **************************************************************************************************
824 70 : SUBROUTINE andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, rho_max, &
825 : rho_min)
826 :
827 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec
828 : REAL(KIND=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min
829 :
830 : CHARACTER(LEN=*), PARAMETER :: routineN = 'andreussi'
831 : REAL(KIND=dp), PARAMETER :: rhotol = 1.0E-12_dp
832 :
833 : INTEGER :: handle, i, j, k
834 : INTEGER, DIMENSION(3) :: lb, ub
835 : REAL(KIND=dp) :: diff, dq, dt, f, ln_rho_max, ln_rho_min, &
836 : q, rho, t, x, y
837 :
838 70 : CALL timeset(routineN, handle)
839 :
840 70 : f = LOG(epsilon_solvent)/twopi
841 70 : diff = rho_max - rho_min
842 70 : IF (diff < SQRT(rhotol)) CPABORT("SCCS: Difference between rho(min) and rho(max) is too small")
843 70 : IF (rho_min >= rhotol) THEN
844 70 : ln_rho_max = LOG(rho_max)
845 70 : ln_rho_min = LOG(rho_min)
846 70 : q = twopi/(ln_rho_max - ln_rho_min)
847 70 : dq = -f*q
848 : END IF
849 :
850 280 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
851 280 : ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
852 :
853 : ! Calculate the dielectric function and its derivative
854 : !$OMP PARALLEL DO DEFAULT(NONE) &
855 : !$OMP PRIVATE(dt,i,j,k,rho,t,x,y) &
856 : !$OMP SHARED(deps_elec,dq,eps_elec,epsilon_solvent,f,lb,ub) &
857 70 : !$OMP SHARED(ln_rho_max,rho_elec,q,rho_max,rho_min)
858 : DO k = lb(3), ub(3)
859 : DO j = lb(2), ub(2)
860 : DO i = lb(1), ub(1)
861 : rho = rho_elec%array(i, j, k)
862 : IF (rho < rho_min) THEN
863 : eps_elec%array(i, j, k) = epsilon_solvent
864 : deps_elec%array(i, j, k) = 0.0_dp
865 : ELSE IF (rho <= rho_max) THEN
866 : x = LOG(rho)
867 : y = q*(ln_rho_max - x)
868 : t = f*(y - SIN(y))
869 : eps_elec%array(i, j, k) = EXP(t)
870 : dt = dq*(1.0_dp - COS(y))
871 : deps_elec%array(i, j, k) = eps_elec%array(i, j, k)*dt/rho
872 : ELSE
873 : eps_elec%array(i, j, k) = 1.0_dp
874 : deps_elec%array(i, j, k) = 0.0_dp
875 : END IF
876 : END DO
877 : END DO
878 : END DO
879 : !$OMP END PARALLEL DO
880 :
881 70 : CALL timestop(handle)
882 :
883 70 : END SUBROUTINE andreussi
884 :
885 : ! **************************************************************************************************
886 : !> \brief Calculate the smoothed dielectric function of Fattebert and Gygi
887 : !> \param rho_elec ...
888 : !> \param eps_elec ...
889 : !> \param deps_elec ...
890 : !> \param epsilon_solvent ...
891 : !> \param beta ...
892 : !> \param rho_zero ...
893 : !> \par History:
894 : !> - Creation (15.10.2013,MK)
895 : !> \author Matthias Krack (MK)
896 : !> \version 1.0
897 : ! **************************************************************************************************
898 22 : SUBROUTINE fattebert_gygi(rho_elec, eps_elec, deps_elec, epsilon_solvent, beta, &
899 : rho_zero)
900 :
901 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec
902 : REAL(KIND=dp), INTENT(IN) :: epsilon_solvent, beta, rho_zero
903 :
904 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fattebert_gygi'
905 : REAL(KIND=dp), PARAMETER :: rhotol = 1.0E-12_dp
906 :
907 : INTEGER :: handle, i, j, k
908 : INTEGER, DIMENSION(3) :: lb, ub
909 : REAL(KIND=dp) :: df, f, p, q, rho, s, t, twobeta
910 :
911 22 : CALL timeset(routineN, handle)
912 :
913 22 : df = (1.0_dp - epsilon_solvent)/rho_zero
914 22 : f = 0.5_dp*(epsilon_solvent - 1.0_dp)
915 22 : q = 1.0_dp/rho_zero
916 22 : twobeta = 2.0_dp*beta
917 :
918 88 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
919 88 : ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
920 :
921 : ! Calculate the smoothed dielectric function and its derivative
922 : !$OMP PARALLEL DO DEFAULT(NONE) &
923 : !$OMP PRIVATE(i,j,k,p,rho,s,t) &
924 : !$OMP SHARED(df,deps_elec,eps_elec,epsilon_solvent,f,lb,ub) &
925 22 : !$OMP SHARED(q,rho_elec,twobeta)
926 : DO k = lb(3), ub(3)
927 : DO j = lb(2), ub(2)
928 : DO i = lb(1), ub(1)
929 : rho = rho_elec%array(i, j, k)
930 : IF (rho < rhotol) THEN
931 : eps_elec%array(i, j, k) = epsilon_solvent
932 : deps_elec%array(i, j, k) = 0.0_dp
933 : ELSE
934 : s = rho*q
935 : p = s**twobeta
936 : t = 1.0_dp/(1.0_dp + p)
937 : eps_elec%array(i, j, k) = 1.0_dp + f*(1.0_dp + (1.0_dp - p)*t)
938 : deps_elec%array(i, j, k) = df*twobeta*t*t*p/s
939 : END IF
940 : END DO
941 : END DO
942 : END DO
943 : !$OMP END PARALLEL DO
944 :
945 22 : CALL timestop(handle)
946 :
947 22 : END SUBROUTINE fattebert_gygi
948 :
949 : ! **************************************************************************************************
950 : !> \brief Build the numerical derivative of a function on realspace grid
951 : !> \param f ...
952 : !> \param df ...
953 : !> \param method ...
954 : !> \param pw_env ...
955 : !> \param input ...
956 : !> \par History:
957 : !> - Creation (15.11.2013,MK)
958 : !> \author Matthias Krack (MK)
959 : !> \version 1.0
960 : ! **************************************************************************************************
961 184 : SUBROUTINE derive(f, df, method, pw_env, input)
962 :
963 : TYPE(pw_r3d_rs_type), INTENT(IN) :: f
964 : TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: df
965 : INTEGER, INTENT(IN) :: method
966 : TYPE(pw_env_type), POINTER :: pw_env
967 : TYPE(section_vals_type), POINTER :: input
968 :
969 : CHARACTER(LEN=*), PARAMETER :: routineN = 'derive'
970 :
971 : INTEGER :: border_points, handle, i
972 : INTEGER, DIMENSION(3) :: lb, n, ub
973 552 : TYPE(pw_c1d_gs_type), DIMENSION(2) :: work_g1d
974 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
975 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
976 : TYPE(realspace_grid_input_type) :: input_settings
977 : TYPE(realspace_grid_type), POINTER :: rs_grid
978 : TYPE(section_vals_type), POINTER :: rs_grid_section
979 :
980 184 : CALL timeset(routineN, handle)
981 :
982 184 : CPASSERT(ASSOCIATED(pw_env))
983 :
984 : ! Perform method specific setup
985 262 : SELECT CASE (method)
986 : CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7)
987 78 : NULLIFY (rs_desc)
988 78 : rs_grid_section => section_vals_get_subs_vals(input, "DFT%MGRID%RS_GRID")
989 0 : SELECT CASE (method)
990 : CASE (sccs_derivative_cd3)
991 0 : border_points = 1
992 : CASE (sccs_derivative_cd5)
993 78 : border_points = 2
994 : CASE (sccs_derivative_cd7)
995 78 : border_points = 3
996 : END SELECT
997 : CALL init_input_type(input_settings, 2*border_points + 1, rs_grid_section, &
998 78 : 1, [-1, -1, -1])
999 : CALL rs_grid_create_descriptor(rs_desc, f%pw_grid, input_settings, &
1000 78 : border_points=border_points)
1001 1638 : ALLOCATE (rs_grid)
1002 78 : CALL rs_grid_create(rs_grid, rs_desc)
1003 : !MK CALL rs_grid_print(rs_grid, 6)
1004 : CASE (sccs_derivative_fft)
1005 424 : lb(1:3) = f%pw_grid%bounds_local(1, 1:3)
1006 424 : ub(1:3) = f%pw_grid%bounds_local(2, 1:3)
1007 106 : NULLIFY (auxbas_pw_pool)
1008 106 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1009 : ! Get work storage for the 1d grids in g-space (derivative calculation)
1010 502 : DO i = 1, SIZE(work_g1d)
1011 318 : CALL auxbas_pw_pool%create_pw(work_g1d(i))
1012 : END DO
1013 : END SELECT
1014 :
1015 : ! Calculate the derivatives
1016 0 : SELECT CASE (method)
1017 : CASE (sccs_derivative_cd3)
1018 0 : CALL derive_fdm_cd3(f, df, rs_grid)
1019 : CASE (sccs_derivative_cd5)
1020 78 : CALL derive_fdm_cd5(f, df, rs_grid)
1021 : CASE (sccs_derivative_cd7)
1022 0 : CALL derive_fdm_cd7(f, df, rs_grid)
1023 : CASE (sccs_derivative_fft)
1024 : ! FFT
1025 106 : CALL pw_transfer(f, work_g1d(1))
1026 424 : DO i = 1, 3
1027 318 : n(:) = 0
1028 318 : n(i) = 1
1029 318 : CALL pw_copy(work_g1d(1), work_g1d(2))
1030 318 : CALL pw_derive(work_g1d(2), n(:))
1031 424 : CALL pw_transfer(work_g1d(2), df(i))
1032 : END DO
1033 : CASE DEFAULT
1034 184 : CPABORT("Invalid derivative method for SCCS specified")
1035 : END SELECT
1036 :
1037 : ! Perform method specific cleanup
1038 78 : SELECT CASE (method)
1039 : CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7)
1040 78 : CALL rs_grid_release(rs_grid)
1041 78 : DEALLOCATE (rs_grid)
1042 78 : CALL rs_grid_release_descriptor(rs_desc)
1043 : CASE (sccs_derivative_fft)
1044 502 : DO i = 1, SIZE(work_g1d)
1045 318 : CALL auxbas_pw_pool%give_back_pw(work_g1d(i))
1046 : END DO
1047 : END SELECT
1048 :
1049 184 : CALL timestop(handle)
1050 :
1051 736 : END SUBROUTINE derive
1052 :
1053 : ! **************************************************************************************************
1054 : !> \brief Calculate the finite difference between two isosurfaces of the
1055 : !> electronic density. The smoothed dielectric function of
1056 : !> Andreussi et al. is used as switching function eventually
1057 : !> defining the quantum volume and surface of the cavity.
1058 : !> \param rho_elec ...
1059 : !> \param norm_drho_elec ...
1060 : !> \param dtheta ...
1061 : !> \param epsilon_solvent ...
1062 : !> \param rho_max ...
1063 : !> \param rho_min ...
1064 : !> \param delta_rho ...
1065 : !> \par History:
1066 : !> - Creation (21.12.2013,MK)
1067 : !> \author Matthias Krack (MK)
1068 : !> \version 1.0
1069 : ! **************************************************************************************************
1070 70 : SUBROUTINE surface_andreussi(rho_elec, norm_drho_elec, dtheta, &
1071 : epsilon_solvent, rho_max, rho_min, delta_rho)
1072 :
1073 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, norm_drho_elec, dtheta
1074 : REAL(KIND=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min, &
1075 : delta_rho
1076 :
1077 : CHARACTER(LEN=*), PARAMETER :: routineN = 'surface_andreussi'
1078 : REAL(KIND=dp), PARAMETER :: rhotol = 1.0E-12_dp
1079 :
1080 : INTEGER :: handle, i, j, k, l
1081 : INTEGER, DIMENSION(3) :: lb, ub
1082 : REAL(KIND=dp) :: diff, e, eps_elec, f, ln_rho_max, &
1083 : ln_rho_min, q, rho, t, x, y
1084 : REAL(KIND=dp), DIMENSION(2) :: theta
1085 :
1086 70 : CALL timeset(routineN, handle)
1087 :
1088 70 : e = epsilon_solvent - 1.0_dp
1089 70 : f = LOG(epsilon_solvent)/twopi
1090 70 : diff = rho_max - rho_min
1091 70 : IF (diff < SQRT(rhotol)) CPABORT("SCCS: Difference between rho(min) and rho(max) is too small")
1092 70 : IF (rho_min >= rhotol) THEN
1093 70 : ln_rho_max = LOG(rho_max)
1094 70 : ln_rho_min = LOG(rho_min)
1095 70 : q = twopi/(ln_rho_max - ln_rho_min)
1096 : END IF
1097 :
1098 280 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1099 280 : ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1100 :
1101 : ! Calculate finite difference between two isosurfaces
1102 : !$OMP PARALLEL DO DEFAULT(NONE) &
1103 : !$OMP PRIVATE(eps_elec,i,j,k,l,rho,t,theta,x,y) &
1104 : !$OMP SHARED(delta_rho,dtheta,e,epsilon_solvent,f,lb) &
1105 70 : !$OMP SHARED(ln_rho_max,norm_drho_elec,rho_elec,q,rho_max,rho_min,ub)
1106 : DO k = lb(3), ub(3)
1107 : DO j = lb(2), ub(2)
1108 : DO i = lb(1), ub(1)
1109 : DO l = 1, 2
1110 : rho = rho_elec%array(i, j, k) + (REAL(l, KIND=dp) - 1.5_dp)*delta_rho
1111 : IF (rho < rho_min) THEN
1112 : eps_elec = epsilon_solvent
1113 : ELSE IF (rho <= rho_max) THEN
1114 : x = LOG(rho)
1115 : y = q*(ln_rho_max - x)
1116 : t = f*(y - SIN(y))
1117 : eps_elec = EXP(t)
1118 : ELSE
1119 : eps_elec = 1.0_dp
1120 : END IF
1121 : theta(l) = (epsilon_solvent - eps_elec)/e
1122 : END DO
1123 : dtheta%array(i, j, k) = (theta(2) - theta(1))*norm_drho_elec%array(i, j, k)/delta_rho
1124 : END DO
1125 : END DO
1126 : END DO
1127 : !$OMP END PARALLEL DO
1128 :
1129 70 : CALL timestop(handle)
1130 :
1131 70 : END SUBROUTINE surface_andreussi
1132 :
1133 : ! **************************************************************************************************
1134 : !> \brief Calculate the finite difference between two isosurfaces of the
1135 : !> the electronic density. The smoothed dielectric function of
1136 : !> Fattebert and Gygi is used as switching function eventually
1137 : !> defining the quantum volume and surface of the cavity.
1138 : !> \param rho_elec ...
1139 : !> \param norm_drho_elec ...
1140 : !> \param dtheta ...
1141 : !> \param epsilon_solvent ...
1142 : !> \param beta ...
1143 : !> \param rho_zero ...
1144 : !> \param delta_rho ...
1145 : !> \par History:
1146 : !> - Creation (21.12.2013,MK)
1147 : !> \author Matthias Krack (MK)
1148 : !> \version 1.0
1149 : ! **************************************************************************************************
1150 22 : SUBROUTINE surface_fattebert_gygi(rho_elec, norm_drho_elec, dtheta, &
1151 : epsilon_solvent, beta, rho_zero, delta_rho)
1152 :
1153 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, norm_drho_elec, dtheta
1154 : REAL(KIND=dp), INTENT(IN) :: epsilon_solvent, beta, rho_zero, &
1155 : delta_rho
1156 :
1157 : CHARACTER(LEN=*), PARAMETER :: routineN = 'surface_fattebert_gygi'
1158 : REAL(KIND=dp), PARAMETER :: rhotol = 1.0E-12_dp
1159 :
1160 : INTEGER :: handle, i, j, k, l
1161 : INTEGER, DIMENSION(3) :: lb, ub
1162 : REAL(KIND=dp) :: e, eps_elec, f, p, q, rho, s, t, twobeta
1163 : REAL(KIND=dp), DIMENSION(2) :: theta
1164 :
1165 22 : CALL timeset(routineN, handle)
1166 :
1167 22 : e = epsilon_solvent - 1.0_dp
1168 22 : f = 0.5_dp*e
1169 22 : q = 1.0_dp/rho_zero
1170 22 : twobeta = 2.0_dp*beta
1171 :
1172 88 : lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1173 88 : ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1174 :
1175 : ! Calculate finite difference between two isosurfaces
1176 : !$OMP PARALLEL DO DEFAULT(NONE) &
1177 : !$OMP PRIVATE(eps_elec,i,j,k,l,p,rho,s,t,theta) &
1178 : !$OMP SHARED(delta_rho,dtheta,e,epsilon_solvent,f,lb) &
1179 22 : !$OMP SHARED(norm_drho_elec,q,rho_elec,twobeta,ub)
1180 : DO k = lb(3), ub(3)
1181 : DO j = lb(2), ub(2)
1182 : DO i = lb(1), ub(1)
1183 : DO l = 1, 2
1184 : rho = rho_elec%array(i, j, k) + (REAL(l, KIND=dp) - 1.5_dp)*delta_rho
1185 : IF (rho < rhotol) THEN
1186 : eps_elec = epsilon_solvent
1187 : ELSE
1188 : s = rho*q
1189 : p = s**twobeta
1190 : t = 1.0_dp/(1.0_dp + p)
1191 : eps_elec = 1.0_dp + f*(1.0_dp + (1.0_dp - p)*t)
1192 : END IF
1193 : theta(l) = (epsilon_solvent - eps_elec)/e
1194 : END DO
1195 : dtheta%array(i, j, k) = (theta(2) - theta(1))*norm_drho_elec%array(i, j, k)/delta_rho
1196 : END DO
1197 : END DO
1198 : END DO
1199 : !$OMP END PARALLEL DO
1200 :
1201 22 : CALL timestop(handle)
1202 :
1203 22 : END SUBROUTINE surface_fattebert_gygi
1204 :
1205 : ! **************************************************************************************************
1206 : !> \brief Print SCCS results
1207 : !> \param energy ...
1208 : !> \param sccs_control ...
1209 : !> \param output_unit ...
1210 : !> \par History:
1211 : !> - Creation (11.11.2022,MK)
1212 : !> \author Matthias Krack (MK)
1213 : !> \version 1.0
1214 : ! **************************************************************************************************
1215 52 : SUBROUTINE print_sccs_results(energy, sccs_control, output_unit)
1216 :
1217 : TYPE(qs_energy_type), POINTER :: energy
1218 : TYPE(sccs_control_type), POINTER :: sccs_control
1219 : INTEGER, INTENT(IN) :: output_unit
1220 :
1221 52 : IF (output_unit > 0) THEN
1222 52 : CPASSERT(ASSOCIATED(energy))
1223 52 : CPASSERT(ASSOCIATED(sccs_control))
1224 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
1225 52 : "SCCS| Hartree energy of solute and solvent [Hartree]", energy%sccs_hartree, &
1226 104 : "SCCS| Hartree energy of the solute only [Hartree]", energy%hartree
1227 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14,/,T3,A,T61,F20.3)") &
1228 52 : "SCCS| Polarisation energy [Hartree]", energy%sccs_pol, &
1229 52 : "SCCS| [kcal/mol]", &
1230 52 : cp_unit_from_cp2k(energy%sccs_pol, "kcalmol"), &
1231 52 : "SCCS| Cavitation energy [Hartree]", energy%sccs_cav, &
1232 52 : "SCCS| [kcal/mol]", &
1233 52 : cp_unit_from_cp2k(energy%sccs_cav, "kcalmol"), &
1234 52 : "SCCS| Dispersion free energy [Hartree]", energy%sccs_dis, &
1235 52 : "SCCS| [kcal/mol]", &
1236 52 : cp_unit_from_cp2k(energy%sccs_dis, "kcalmol"), &
1237 52 : "SCCS| Repulsion free energy [Hartree]", energy%sccs_rep, &
1238 52 : "SCCS| [kcal/mol]", &
1239 52 : cp_unit_from_cp2k(energy%sccs_rep, "kcalmol"), &
1240 52 : "SCCS| Solvation free energy [Hartree]", energy%sccs_sol, &
1241 52 : "SCCS| [kcal/mol]", &
1242 104 : cp_unit_from_cp2k(energy%sccs_sol, "kcalmol")
1243 : END IF
1244 :
1245 52 : END SUBROUTINE print_sccs_results
1246 :
1247 : END MODULE qs_sccs
|
