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