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