(git:0761fd1)
Loading...
Searching...
No Matches
qs_sccs.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Self-consistent continuum solvation (SCCS) model implementation
10!> \author Matthias Krack (MK)
11!> \version 1.0
12!> \par Literature:
13!> - J.-L. Fattebert and F. Gygi,
14!> Density functional theory for efficient ab initio molecular dynamics
15!> simulations in solution, J. Comput. Chem. 23, 662-666 (2002)
16!> - O. Andreussi, I. Dabo, and N. Marzari,
17!> Revised self-consistent continuum solvation in electronic-structure
18!> calculations, J. Chem. Phys. 136, 064102-20 (2012)
19!> \par History:
20!> - Creation (10.10.2013,MK)
21!> - Derivatives using finite differences (26.11.2013,MK)
22!> - Cube file dump of the dielectric function (19.12.2013,MK)
23!> - Cube file dump of the polarisation potential (20.12.2013,MK)
24!> - Calculation of volume and surface of the cavity (21.12.2013,MK)
25!> - Functional derivative of the cavitation energy (28.12.2013,MK)
26!> - Update printout (11.11.2022,MK)
27! **************************************************************************************************
28
29MODULE qs_sccs
30
31 USE atomic_kind_types, ONLY: atomic_kind_type,&
32 get_atomic_kind
33 USE cp_control_types, ONLY: dft_control_type,&
34 sccs_control_type
35 USE cp_log_handling, ONLY: cp_get_default_logger,&
36 cp_logger_get_default_io_unit,&
37 cp_logger_type
38 USE cp_output_handling, ONLY: cp_p_file,&
39 cp_print_key_finished_output,&
40 cp_print_key_should_output,&
41 cp_print_key_unit_nr,&
42 low_print_level
43 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
44 USE cp_realspace_grid_init, ONLY: init_input_type
45 USE cp_subsys_types, ONLY: cp_subsys_get,&
46 cp_subsys_type
47 USE cp_units, ONLY: cp_unit_from_cp2k
48 USE input_constants, ONLY: sccs_andreussi,&
49 sccs_derivative_cd3,&
50 sccs_derivative_cd5,&
51 sccs_derivative_cd7,&
52 sccs_derivative_fft,&
53 sccs_fattebert_gygi,&
54 sccs_saa_andreussi
55 USE input_section_types, ONLY: section_get_ivals,&
56 section_get_lval,&
57 section_vals_get_subs_vals,&
58 section_vals_type
59 USE kinds, ONLY: default_path_length,&
60 default_string_length,&
61 dp,&
62 int_8
63 USE mathconstants, ONLY: fourpi,&
64 oorootpi,&
65 pi,&
66 twopi
67 USE message_passing, ONLY: mp_para_env_type
68 USE particle_list_types, ONLY: particle_list_type
69 USE pw_env_types, ONLY: pw_env_get,&
70 pw_env_type
71 USE pw_methods, ONLY: pw_axpy,&
72 pw_copy,&
73 pw_derive,&
74 pw_integral_ab,&
75 pw_integrate_function,&
76 pw_scale,&
77 pw_transfer,&
78 pw_zero
79 USE pw_poisson_methods, ONLY: pw_func_u_convolution,&
80 pw_poisson_solve
81 USE pw_poisson_types, ONLY: pw_poisson_analytic,&
82 pw_poisson_mt,&
83 pw_poisson_type
84 USE pw_pool_types, ONLY: pw_pool_p_type,&
85 pw_pool_type
86 USE pw_types, ONLY: pw_c1d_gs_type,&
87 pw_r3d_rs_type
88 USE qs_collocate_density, ONLY: calculate_rho_core
89 USE qs_energy_types, ONLY: qs_energy_type
90 USE qs_environment_types, ONLY: get_qs_env,&
91 qs_environment_type
92 USE qs_kind_types, ONLY: get_qs_kind,&
93 qs_kind_type
94 USE qs_rho_types, ONLY: qs_rho_get,&
95 qs_rho_type
96 USE qs_scf_types, ONLY: qs_scf_env_type
97 USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
98 realspace_grid_input_type,&
99 realspace_grid_type,&
100 rs_grid_create,&
101 rs_grid_create_descriptor,&
102 rs_grid_release,&
103 rs_grid_release_descriptor
104 USE rs_methods, ONLY: derive_fdm_cd3,&
105 derive_fdm_cd5,&
106 derive_fdm_cd7
107 USE string_utilities, ONLY: uppercase
108#include "./base/base_uses.f90"
109
110 IMPLICIT NONE
111
112 PRIVATE
113
114 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_sccs'
115
116 PUBLIC :: print_sccs_results, sccs, sccs_core_fill_coefficients
117
118CONTAINS
119
120! **************************************************************************************************
121!> \brief Self-consistent continuum solvation (SCCS) model implementation
122!> \param qs_env ...
123!> \param rho_tot_gspace ...
124!> \param v_hartree_gspace ...
125!> \param v_sccs ...
126!> \param h_stress ...
127!> \par History:
128!> - Creation (10.10.2013,MK)
129!> \author Matthias Krack (MK)
130!> \version 1.0
131! **************************************************************************************************
132
133 SUBROUTINE sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs, h_stress)
134
135 TYPE(qs_environment_type), POINTER :: qs_env
136 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_tot_gspace, v_hartree_gspace
137 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_sccs
138 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
139 OPTIONAL :: h_stress
140
141 CHARACTER(LEN=*), PARAMETER :: routinen = 'sccs'
142 REAL(kind=dp), PARAMETER :: epstol = 1.0e-8_dp
143
144 CHARACTER(LEN=4*default_string_length) :: message, my_pos_cube
145 CHARACTER(LEN=default_path_length) :: cube_path, filename, mpi_filename, &
146 print_path
147 INTEGER :: cube_unit, handle, i, ispin, iter, j, k, &
148 nspin, output_unit, print_level
149 INTEGER(KIND=int_8) :: ngpts
150 INTEGER, DIMENSION(3) :: lb, ub
151 LOGICAL :: append_cube, calculate_stress_tensor, &
152 do_kpoints, mpi_io, should_output
153 REAL(kind=dp) :: alpha_zeta, cavity_surface, cavity_volume, cell_volume, delta_eta, &
154 delta_zeta, dphi2, dvol, e_tot, epsilon_solvent, f, f0, polarisation_charge, r, r_solv, &
155 rho_delta, rho_delta_avg, rho_delta_max, rho_iter_new, tot_rho_elec, tot_rho_solute
156 REAL(kind=dp), DIMENSION(3) :: abc, lxyz, uxyz
157 TYPE(cp_logger_type), POINTER :: logger
158 TYPE(cp_subsys_type), POINTER :: cp_subsys
159 TYPE(dft_control_type), POINTER :: dft_control
160 TYPE(mp_para_env_type), POINTER :: para_env
161 TYPE(particle_list_type), POINTER :: particles
162 TYPE(pw_env_type), POINTER :: pw_env
163 TYPE(pw_poisson_type), POINTER :: poisson_env
164 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
165 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
166 TYPE(pw_r3d_rs_type) :: d_s_rhoel, deps_elec, dtf, eps_elec, ff, &
167 p_e_interface, s, t, tem_convolution, &
168 tem_func, u
169 TYPE(pw_r3d_rs_type), DIMENSION(3) :: dln_eps_elec, dphi_tot
170 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_pw_r
171 TYPE(pw_r3d_rs_type), POINTER :: rho_pw_r_sccs
172 TYPE(qs_energy_type), POINTER :: energy
173 TYPE(qs_rho_type), POINTER :: rho
174 TYPE(qs_scf_env_type), POINTER :: scf_env
175 TYPE(sccs_control_type), POINTER :: sccs_control
176 TYPE(section_vals_type), POINTER :: input
177
178 CALL timeset(routinen, handle)
179
180 NULLIFY (auxbas_pw_pool)
181 NULLIFY (cp_subsys)
182 NULLIFY (dft_control)
183 NULLIFY (energy)
184 NULLIFY (input)
185 NULLIFY (logger)
186 NULLIFY (para_env)
187 NULLIFY (particles)
188 NULLIFY (poisson_env)
189 NULLIFY (pw_env)
190 NULLIFY (pw_pools)
191 NULLIFY (rho)
192 NULLIFY (sccs_control)
193 NULLIFY (scf_env)
194
195 ! Load data from Quickstep environment
196 CALL get_qs_env(qs_env=qs_env, &
197 cp_subsys=cp_subsys, &
198 do_kpoints=do_kpoints, &
199 dft_control=dft_control, &
200 energy=energy, &
201 input=input, &
202 para_env=para_env, &
203 pw_env=pw_env, &
204 rho=rho, &
205 scf_env=scf_env)
206 CALL cp_subsys_get(cp_subsys, particles=particles)
207
208 sccs_control => dft_control%sccs_control
209
210 cpassert(ASSOCIATED(qs_env))
211
212 IF (do_kpoints) THEN
213 cpwarn("SCCS with k-points has not yet been fully validated")
214 END IF
215
216 IF (PRESENT(h_stress)) THEN
217 calculate_stress_tensor = .true.
218 h_stress(:, :) = 0.0_dp
219 cpwarn("The stress tensor for SCCS has not yet been fully validated")
220 ELSE
221 calculate_stress_tensor = .false.
222 END IF
223
224 ! Get access to the PW grid pool
225 CALL pw_env_get(pw_env, &
226 auxbas_pw_pool=auxbas_pw_pool, &
227 pw_pools=pw_pools, &
228 poisson_env=poisson_env)
229
230 CALL pw_zero(v_sccs)
231
232 ! Calculate no SCCS contribution, if the requested SCF convergence threshold is not reached yet
233 IF (.NOT. sccs_control%sccs_activated) THEN
234 IF (sccs_control%eps_scf > 0.0_dp) THEN
235 IF ((scf_env%iter_delta > sccs_control%eps_scf) .OR. &
236 ((qs_env%scf_env%outer_scf%iter_count == 0) .AND. &
237 (qs_env%scf_env%iter_count <= 1))) THEN
238 IF (calculate_stress_tensor) THEN
239 ! Request also the calculation of the stress tensor contribution
240 CALL pw_poisson_solve(poisson_env=poisson_env, &
241 density=rho_tot_gspace, &
242 ehartree=energy%hartree, &
243 vhartree=v_hartree_gspace, &
244 h_stress=h_stress)
245 ELSE
246 CALL pw_poisson_solve(poisson_env=poisson_env, &
247 density=rho_tot_gspace, &
248 ehartree=energy%hartree, &
249 vhartree=v_hartree_gspace)
250 END IF
251 energy%sccs_pol = 0.0_dp
252 energy%sccs_cav = 0.0_dp
253 energy%sccs_dis = 0.0_dp
254 energy%sccs_rep = 0.0_dp
255 energy%sccs_sol = 0.0_dp
256 energy%sccs_hartree = energy%hartree
257 CALL timestop(handle)
258 RETURN
259 END IF
260 END IF
261 sccs_control%sccs_activated = .true.
262 END IF
263
264 nspin = dft_control%nspins
265
266 ! Manage print output control
267 logger => cp_get_default_logger()
268 print_level = logger%iter_info%print_level
269 print_path = "DFT%PRINT%SCCS"
270 should_output = (btest(cp_print_key_should_output(logger%iter_info, input, &
271 trim(print_path)), cp_p_file))
272 output_unit = cp_print_key_unit_nr(logger, input, trim(print_path), &
273 extension=".sccs", &
274 ignore_should_output=should_output, &
275 log_filename=.false.)
276
277 ! Get rho
278 CALL qs_rho_get(rho_struct=rho, &
279 rho_r=rho_pw_r, &
280 rho_r_sccs=rho_pw_r_sccs)
281
282 ! Retrieve the last rho_iter from the previous SCCS cycle if available
283 cpassert(ASSOCIATED(rho_pw_r_sccs))
284
285 ! Retrieve the total electronic density in r-space
286 block
287 TYPE(pw_r3d_rs_type) :: rho_elec
288 CALL auxbas_pw_pool%create_pw(rho_elec)
289
290 ! Retrieve grid parameters
291 ngpts = rho_elec%pw_grid%ngpts
292 dvol = rho_elec%pw_grid%dvol
293 cell_volume = rho_elec%pw_grid%vol
294 abc(1:3) = real(rho_elec%pw_grid%npts(1:3), kind=dp)*rho_elec%pw_grid%dr(1:3)
295 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
296 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
297
298 CALL pw_copy(rho_pw_r(1), rho_elec)
299 DO ispin = 2, nspin
300 CALL pw_axpy(rho_pw_r(ispin), rho_elec)
301 END DO
302 tot_rho_elec = pw_integrate_function(rho_elec)
303
304 ! Optionally fill the low valence density region at the nuclei of pseudopotential
305 ! atoms with normalised Gaussians Z_val/(sqrt(pi)*sigma)**3*exp(-r**2/sigma**2)
306 ! (Environ "full" boundary). The filled density is only used for the definition
307 ! of the dielectric function (cavity)
308 IF (sccs_control%core_fill) THEN
309 block
310 REAL(kind=dp) :: tot_rho_fill
311 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: calpha, ccore
312 TYPE(pw_r3d_rs_type) :: rho_fill
313 CALL sccs_core_fill_coefficients(qs_env, calpha, ccore, cp_logger_get_default_io_unit(logger))
314 CALL auxbas_pw_pool%create_pw(rho_fill)
315 ! calculate_rho_core collocates -ccore*exp(-calpha*r**2), i.e. the negated fill
316 CALL calculate_rho_core(rho_fill, tot_rho_fill, qs_env, calpha=calpha, ccore=ccore)
317 CALL pw_axpy(rho_fill, rho_elec, alpha=-1.0_dp)
318 CALL auxbas_pw_pool%give_back_pw(rho_fill)
319 DEALLOCATE (calpha, ccore)
320 END block
321 END IF
322
323 ! Calculate the dielectric (smoothed) function of rho_elec in r-space
324 CALL auxbas_pw_pool%create_pw(eps_elec)
325 CALL auxbas_pw_pool%create_pw(deps_elec)
326 IF (sccs_control%method_id == sccs_saa_andreussi) THEN
327 CALL auxbas_pw_pool%create_pw(s)
328 CALL auxbas_pw_pool%create_pw(d_s_rhoel)
329 CALL auxbas_pw_pool%create_pw(u)
330 CALL auxbas_pw_pool%create_pw(ff)
331 CALL auxbas_pw_pool%create_pw(t)
332 CALL auxbas_pw_pool%create_pw(dtf)
333 f0 = sccs_control%f0
334 delta_eta = sccs_control%delta_eta
335 alpha_zeta = sccs_control%alpha_zeta
336 r_solv = sccs_control%R_solv
337 delta_zeta = sccs_control%delta_zeta
338 END IF
339
340 ! Relative permittivity or dielectric constant of the solvent (medium)
341 epsilon_solvent = sccs_control%epsilon_solvent
342 SELECT CASE (sccs_control%method_id)
343 CASE (sccs_andreussi)
344 CALL andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%rho_max, &
345 sccs_control%rho_min)
346 CASE (sccs_fattebert_gygi)
347 CALL fattebert_gygi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%beta, &
348 sccs_control%rho_zero)
349 CASE (sccs_saa_andreussi)
350 CALL sa_andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%rho_max, &
351 sccs_control%rho_min, s, d_s_rhoel)
352 CASE DEFAULT
353 cpabort("Invalid method specified for SCCS model")
354 END SELECT
355
356 IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
357 ! Optional output of the dielectric function in cube file format
358 filename = "DIELECTRIC_FUNCTION"
359 cube_path = trim(print_path)//"%"//trim(filename)
360 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
361 cp_p_file)) THEN
362 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
363 my_pos_cube = "REWIND"
364 IF (append_cube) my_pos_cube = "APPEND"
365 mpi_io = .true.
366 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
367 extension=".cube", middle_name=trim(filename), &
368 file_position=my_pos_cube, log_filename=.false., &
369 mpi_io=mpi_io, fout=mpi_filename)
370 IF (output_unit > 0) THEN
371 IF (.NOT. mpi_io) THEN
372 INQUIRE (unit=cube_unit, name=filename)
373 ELSE
374 filename = mpi_filename
375 END IF
376 WRITE (unit=output_unit, fmt="(T3,A)") &
377 "SCCS| The dielectric function is written in cube file format to the file:", &
378 "SCCS| "//trim(filename)
379 END IF
380 CALL cp_pw_to_cube(eps_elec, cube_unit, trim(filename), particles=particles, &
381 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), &
382 mpi_io=mpi_io)
383 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
384 END IF
385 END IF
386
387 ! Calculate the (quantum) volume and surface of the solute cavity
388 cavity_surface = 0.0_dp
389 cavity_volume = 0.0_dp
390
391 IF (abs(epsilon_solvent - 1.0_dp) > epstol) THEN
392
393 block
394 TYPE(pw_r3d_rs_type) :: theta, norm_drho_elec
395 TYPE(pw_r3d_rs_type), DIMENSION(3) :: drho_elec
396 CALL auxbas_pw_pool%create_pw(theta)
397 CALL pw_zero(theta)
398
399 ! Calculate the (quantum) volume of the solute cavity
400 f = 1.0_dp/(epsilon_solvent - 1.0_dp)
401!$OMP PARALLEL DO DEFAULT(NONE) &
402!$OMP PRIVATE(i,j,k) &
403!$OMP SHARED(epsilon_solvent,eps_elec,f,lb,theta,ub)
404 DO k = lb(3), ub(3)
405 DO j = lb(2), ub(2)
406 DO i = lb(1), ub(1)
407 theta%array(i, j, k) = f*(epsilon_solvent - eps_elec%array(i, j, k))
408 END DO
409 END DO
410 END DO
411!$OMP END PARALLEL DO
412 cavity_volume = pw_integrate_function(theta)
413
414 ! Calculate the derivative of the electronic density in r-space
415 ! TODO: Could be retrieved from the QS environment
416 DO i = 1, 3
417 CALL auxbas_pw_pool%create_pw(drho_elec(i))
418 END DO
419 CALL derive(rho_elec, drho_elec, sccs_derivative_fft, pw_env, input)
420
421 CALL auxbas_pw_pool%create_pw(norm_drho_elec)
422
423 ! Calculate the norm of the gradient of the electronic density in r-space
424!$OMP PARALLEL DO DEFAULT(NONE) &
425!$OMP PRIVATE(i,j,k) &
426!$OMP SHARED(drho_elec,lb,norm_drho_elec,ub)
427 DO k = lb(3), ub(3)
428 DO j = lb(2), ub(2)
429 DO i = lb(1), ub(1)
430 norm_drho_elec%array(i, j, k) = sqrt(drho_elec(1)%array(i, j, k)* &
431 drho_elec(1)%array(i, j, k) + &
432 drho_elec(2)%array(i, j, k)* &
433 drho_elec(2)%array(i, j, k) + &
434 drho_elec(3)%array(i, j, k)* &
435 drho_elec(3)%array(i, j, k))
436 END DO
437 END DO
438 END DO
439!$OMP END PARALLEL DO
440
441 ! Optional output of the norm of the density gradient in cube file format
442 filename = "DENSITY_GRADIENT"
443 cube_path = trim(print_path)//"%"//trim(filename)
444 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
445 cp_p_file)) THEN
446 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
447 my_pos_cube = "REWIND"
448 IF (append_cube) my_pos_cube = "APPEND"
449 mpi_io = .true.
450 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
451 extension=".cube", middle_name=trim(filename), &
452 file_position=my_pos_cube, log_filename=.false., &
453 mpi_io=mpi_io, fout=mpi_filename)
454 IF (output_unit > 0) THEN
455 IF (.NOT. mpi_io) THEN
456 INQUIRE (unit=cube_unit, name=filename)
457 ELSE
458 filename = mpi_filename
459 END IF
460 WRITE (unit=output_unit, fmt="(T3,A)") &
461 "SCCS| The norm of the density gradient is written in cube file format to the file:", &
462 "SCCS| "//trim(filename)
463 END IF
464 CALL cp_pw_to_cube(norm_drho_elec, cube_unit, trim(filename), particles=particles, &
465 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), &
466 mpi_io=mpi_io)
467 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
468 END IF
469
470 ! Calculate the (quantum) surface of the solute cavity
471 SELECT CASE (sccs_control%method_id)
472 CASE (sccs_andreussi)
473 CALL surface_andreussi(rho_elec, norm_drho_elec, theta, epsilon_solvent, &
474 sccs_control%rho_max, sccs_control%rho_min, &
475 sccs_control%delta_rho)
476 CASE (sccs_fattebert_gygi)
477 CALL surface_fattebert_gygi(rho_elec, norm_drho_elec, theta, epsilon_solvent, &
478 sccs_control%beta, sccs_control%rho_zero, &
479 sccs_control%delta_rho)
480 CASE (sccs_saa_andreussi)
481 CASE DEFAULT
482 cpabort("Invalid method specified for SCCS model")
483 END SELECT
484 cavity_surface = pw_integrate_function(theta)
485
486 ! Release storage
487 CALL auxbas_pw_pool%give_back_pw(theta)
488 CALL auxbas_pw_pool%give_back_pw(norm_drho_elec)
489 DO i = 1, 3
490 CALL auxbas_pw_pool%give_back_pw(drho_elec(i))
491 END DO
492 END block
493
494 END IF ! epsilon_solvent > 1
495
496 CALL auxbas_pw_pool%give_back_pw(rho_elec)
497 END block
498
499 block
500 TYPE(pw_r3d_rs_type) :: rho_tot, phi_tot, rho_solute, rho_tot_zero
501 ! Retrieve the total charge density (core + elec) of the solute in r-space
502 CALL auxbas_pw_pool%create_pw(rho_solute)
503 CALL pw_zero(rho_solute)
504 CALL pw_transfer(rho_tot_gspace, rho_solute)
505 tot_rho_solute = pw_integrate_function(rho_solute)
506
507 ! Check total charge
508 IF (abs(tot_rho_solute) >= 1.0e-6_dp) THEN
509 IF ((poisson_env%parameters%solver /= pw_poisson_analytic) .AND. &
510 (poisson_env%parameters%solver /= pw_poisson_mt)) THEN
511 WRITE (unit=message, fmt="(A,SP,F0.6,A)") &
512 "The system (solute) has a non-negligible charge of ", -tot_rho_solute, &
513 ". It is recommended to use non-periodic boundary conditions (PERIODIC none) "// &
514 "combined with an appropriate Poisson solver (POISSON_SOLVER MT or analytic)"
515 cpwarn(message)
516 END IF
517 END IF
518
519 ! Reassign work storage to rho_tot_zero, because rho_elec is no longer needed
520 CALL auxbas_pw_pool%create_pw(rho_tot_zero)
521
522 IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
523
524 ! Build the initial (rho_iter = 0) total charge density (solute plus polarisation) in r-space
525 ! eps_elec <- ln(eps_elec)
526!$OMP PARALLEL DO DEFAULT(NONE) &
527!$OMP PRIVATE(i,j,k) &
528!$OMP SHARED(eps_elec,lb,message,output_unit,para_env,ub) &
529!$OMP SHARED(rho_solute,rho_tot_zero)
530 DO k = lb(3), ub(3)
531 DO j = lb(2), ub(2)
532 DO i = lb(1), ub(1)
533 IF (eps_elec%array(i, j, k) < 1.0_dp) THEN
534 WRITE (unit=message, fmt="(A,ES12.3,A,3(I0,A))") &
535 "SCCS| Invalid dielectric function value ", eps_elec%array(i, j, k), &
536 " encountered at grid point (", i, ",", j, ",", k, ")"
537 cpabort(message)
538 END IF
539 rho_tot_zero%array(i, j, k) = rho_solute%array(i, j, k)/eps_elec%array(i, j, k)
540 eps_elec%array(i, j, k) = log(eps_elec%array(i, j, k))
541 END DO
542 END DO
543 END DO
544!$OMP END PARALLEL DO
545
546 ELSE
547
548 lxyz(:) = real(u%pw_grid%bounds(1, :), dp)
549 uxyz(:) = real(u%pw_grid%bounds(2, :), dp)
550!$OMP PARALLEL DO DEFAULT(NONE) &
551!$OMP PRIVATE(i,j,k,r) &
552!$OMP SHARED(lb,message,output_unit,para_env,ub,dtf,uxyz,lxyz) &
553!$OMP SHARED(u, alpha_zeta, R_solv, delta_zeta, eps_elec)
554 DO k = lb(3), ub(3)
555 DO j = lb(2), ub(2)
556 DO i = lb(1), ub(1)
557 IF (eps_elec%array(i, j, k) < 1.0_dp) THEN
558 WRITE (unit=message, fmt="(A,ES12.3,A,3(I0,A))") &
559 "SCCS| Invalid dielectric function value ", eps_elec%array(i, j, k), &
560 " encountered at grid point (", i, ",", j, ",", k, ")"
561 cpabort(message)
562 END IF
563 IF (i > -1) r = ((real(uxyz(1), dp) - real(i, dp) + 0.5_dp)*u%pw_grid%dr(1))**2.0_dp
564 IF (i < 0) r = ((real(i, dp) - real(lxyz(1), dp) + 0.5_dp)*u%pw_grid%dr(1))**2.0_dp
565 IF (j > -1) r = ((real(uxyz(2), dp) - real(j, dp) + 0.5_dp)*u%pw_grid%dr(2))**2.0_dp + r
566 IF (j < 0) r = ((real(j, dp) - real(lxyz(2), dp) + 0.5_dp)*u%pw_grid%dr(2))**2.0_dp + r
567 IF (k > -1) r = ((real(uxyz(3), dp) - real(k, dp) + 0.5_dp)*u%pw_grid%dr(3))**2.0_dp + r
568 IF (k < 0) r = ((real(k, dp) - real(lxyz(3), dp) + 0.5_dp)*u%pw_grid%dr(3))**2.0_dp + r
569 r = sqrt(r)
570 u%array(i, j, k) = (1.0_dp/2.0_dp)*erfc((r - alpha_zeta*r_solv)/delta_zeta)
571 dtf%array(i, j, k) = 1.0_dp
572 END DO
573 END DO
574 END DO
575!$OMP END PARALLEL DO
576
577 CALL pw_func_u_convolution(poisson_env=poisson_env, func=dtf, convolution=ff, u=u)
578 CALL pw_scale(u, 1.0_dp/ff%array(0, 0, 0))
579 CALL pw_func_u_convolution(poisson_env=poisson_env, func=s, convolution=ff, u=u)
580
581!$OMP PARALLEL DO DEFAULT(NONE) &
582!$OMP PRIVATE(i,j,k) &
583!$OMP SHARED(lb,para_env,ub, delta_eta, f0, epsilon_solvent) &
584!$OMP SHARED(ff, t, dtf, eps_elec, rho_tot_zero, rho_solute, s)
585 DO k = lb(3), ub(3)
586 DO j = lb(2), ub(2)
587 DO i = lb(1), ub(1)
588 t%array(i, j, k) = 0.5_dp*(1.0_dp + erf((ff%array(i, j, k) - f0)/delta_eta))
589 dtf%array(i, j, k) = oorootpi*exp(-((ff%array(i, j, k) - f0)/delta_eta)**2)
590 eps_elec%array(i, j, k) = exp(log(epsilon_solvent)*(1 - (s%array(i, j, k) &
591 + (1 - s%array(i, j, k))*t%array(i, j, k))))
592 rho_tot_zero%array(i, j, k) = rho_solute%array(i, j, k)/eps_elec%array(i, j, k)
593! eps_elec%array(i, j, k) = LOG(eps_elec%array(i, j, k))
594 END DO
595 END DO
596 END DO
597!$OMP END PARALLEL DO
598
599 ! Optional output of the dielectric function in cube file format (solvent aware algorithm case)
600 filename = "DIELECTRIC_FUNCTION"
601 cube_path = trim(print_path)//"%"//trim(filename)
602 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
603 cp_p_file)) THEN
604 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
605 my_pos_cube = "REWIND"
606 IF (append_cube) my_pos_cube = "APPEND"
607 mpi_io = .true.
608 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
609 extension=".cube", middle_name=trim(filename), &
610 file_position=my_pos_cube, log_filename=.false., &
611 mpi_io=mpi_io, fout=mpi_filename)
612 IF (output_unit > 0) THEN
613 IF (.NOT. mpi_io) THEN
614 INQUIRE (unit=cube_unit, name=filename)
615 ELSE
616 filename = mpi_filename
617 END IF
618 WRITE (unit=output_unit, fmt="(T3,A)") &
619 "SCCS| The dielectric function is written in cube file format to the file:", &
620 "SCCS| "//trim(filename)
621 END IF
622 CALL cp_pw_to_cube(eps_elec, cube_unit, trim(filename), particles=particles, &
623 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), &
624 mpi_io=mpi_io)
625 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
626 END IF
627
628!$OMP PARALLEL DO DEFAULT(NONE) &
629!$OMP PRIVATE(i,j,k) &
630!$OMP SHARED(lb,para_env,ub, delta_eta, f0, epsilon_solvent) &
631!$OMP SHARED(ff, t, dtf, eps_elec, rho_tot_zero, rho_solute, s)
632 DO k = lb(3), ub(3)
633 DO j = lb(2), ub(2)
634 DO i = lb(1), ub(1)
635 eps_elec%array(i, j, k) = log(eps_elec%array(i, j, k))
636 END DO
637 END DO
638 END DO
639!$OMP END PARALLEL DO
640
641 END IF
642
643 ! Build the derivative of LOG(eps_elec)
644 DO i = 1, 3
645 CALL auxbas_pw_pool%create_pw(dln_eps_elec(i))
646 CALL pw_zero(dln_eps_elec(i))
647 END DO
648 CALL derive(eps_elec, dln_eps_elec, sccs_control%derivative_method, pw_env, input)
649 IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
650 CALL auxbas_pw_pool%give_back_pw(eps_elec)
651 END IF
652
653 ! Print header for the SCCS cycle
654 IF (should_output .AND. (output_unit > 0)) THEN
655 IF (print_level > low_print_level) THEN
656 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.12)") &
657 "SCCS| Total electronic charge density ", -tot_rho_elec, &
658 "SCCS| Total charge density (solute) ", -tot_rho_solute
659 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.3)") &
660 "SCCS| Volume of the cell [bohr^3]", cell_volume, &
661 "SCCS| [angstrom^3]", &
662 cp_unit_from_cp2k(cell_volume, "angstrom^3")
663 IF (abs(epsilon_solvent - 1.0_dp) > epstol) THEN
664 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.3)") &
665 "SCCS| Volume of the solute cavity [bohr^3]", cavity_volume, &
666 "SCCS| [angstrom^3]", &
667 cp_unit_from_cp2k(cavity_volume, "angstrom^3"), &
668 "SCCS| Surface of the solute cavity [bohr^2]", cavity_surface, &
669 "SCCS| [angstrom^2]", &
670 cp_unit_from_cp2k(cavity_surface, "angstrom^2")
671 END IF
672 WRITE (unit=output_unit, fmt="(T3,A)") &
673 "SCCS|", &
674 "SCCS| Step Average residual Maximum residual E(Hartree) [a.u.]"
675 END IF
676 END IF
677
678 ! Get storage for the derivative of the total potential (field) in r-space
679 DO i = 1, 3
680 CALL auxbas_pw_pool%create_pw(dphi_tot(i))
681 END DO
682
683 ! Initialise the total charge density in r-space rho_tot with rho_tot_zero + rho_iter_zero
684 CALL auxbas_pw_pool%create_pw(rho_tot)
685 CALL pw_copy(rho_tot_zero, rho_tot)
686 CALL pw_axpy(rho_pw_r_sccs, rho_tot)
687
688 CALL auxbas_pw_pool%create_pw(phi_tot)
689 CALL pw_zero(phi_tot)
690
691 IF (sccs_control%method_id == sccs_saa_andreussi) THEN
692 CALL auxbas_pw_pool%create_pw(p_e_interface)
693 CALL pw_zero(p_e_interface)
694 END IF
695
696 ! Main SCCS iteration loop
697 iter = 0
698
699 iter_loop: DO
700
701 ! Increment iteration counter
702 iter = iter + 1
703
704 ! Check if the requested maximum number of SCCS iterations is reached
705 IF (iter > sccs_control%max_iter) THEN
706 IF (output_unit > 0) THEN
707 WRITE (unit=output_unit, fmt="(T3,A,/,T3,A,I0,A)") &
708 "SCCS| Maximum number of SCCS iterations reached", &
709 "SCCS| Iteration cycle did not converge in ", sccs_control%max_iter, " steps"
710 ELSE
711 WRITE (unit=message, fmt="(A,I0,A)") &
712 "The SCCS iteration cycle did not converge in ", sccs_control%max_iter, " steps"
713 cpwarn(message)
714 END IF
715 EXIT iter_loop
716 END IF
717
718 ! Calculate derivative of the current total potential in r-space
719 CALL pw_poisson_solve(poisson_env=poisson_env, &
720 density=rho_tot, &
721 vhartree=phi_tot, &
722 dvhartree=dphi_tot)
723 energy%sccs_hartree = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
724
725 ! Update total charge density (solute plus polarisation) in r-space
726 ! based on the iterated polarisation charge density
727 f = 1.0_dp/fourpi
728 rho_delta_avg = 0.0_dp
729 rho_delta_max = 0.0_dp
730!$OMP PARALLEL DO DEFAULT(NONE) &
731!$OMP PRIVATE(i,j,k,rho_delta,rho_iter_new) &
732!$OMP SHARED(dln_eps_elec,dphi_tot,f,lb,rho_pw_r_sccs,ub) &
733!$OMP SHARED(rho_tot,rho_tot_zero,sccs_control) &
734!$OMP REDUCTION(+:rho_delta_avg) &
735!$OMP REDUCTION(MAX:rho_delta_max)
736 DO k = lb(3), ub(3)
737 DO j = lb(2), ub(2)
738 DO i = lb(1), ub(1)
739 rho_iter_new = (dln_eps_elec(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + &
740 dln_eps_elec(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + &
741 dln_eps_elec(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k))*f
742 rho_iter_new = rho_pw_r_sccs%array(i, j, k) + &
743 sccs_control%mixing*(rho_iter_new - rho_pw_r_sccs%array(i, j, k))
744 rho_delta = abs(rho_iter_new - rho_pw_r_sccs%array(i, j, k))
745 rho_delta_max = max(rho_delta, rho_delta_max)
746 rho_delta_avg = rho_delta_avg + rho_delta
747 rho_tot%array(i, j, k) = rho_tot_zero%array(i, j, k) + rho_iter_new
748 rho_pw_r_sccs%array(i, j, k) = rho_iter_new
749 END DO
750 END DO
751 END DO
752!$OMP END PARALLEL DO
753
754 CALL para_env%sum(rho_delta_avg)
755 rho_delta_avg = rho_delta_avg/real(ngpts, kind=dp)
756 CALL para_env%max(rho_delta_max)
757
758 IF (should_output .AND. (output_unit > 0)) THEN
759 IF (print_level > low_print_level) THEN
760 IF ((abs(rho_delta_avg) < 1.0e-8_dp) .OR. &
761 (abs(rho_delta_avg) >= 1.0e5_dp)) THEN
762 WRITE (unit=output_unit, fmt="(T3,A,I6,4X,ES16.4,4X,ES16.4,1X,F25.12)") &
763 "SCCS| ", iter, rho_delta_avg, rho_delta_max, energy%sccs_hartree
764 ELSE
765 WRITE (unit=output_unit, fmt="(T3,A,I6,4X,F16.8,4X,F16.8,1X,F25.12)") &
766 "SCCS| ", iter, rho_delta_avg, rho_delta_max, energy%sccs_hartree
767 END IF
768 END IF
769 END IF
770
771 ! Check if the SCCS iteration cycle is converged to the requested tolerance
772 IF (rho_delta_max <= sccs_control%eps_sccs) THEN
773 IF (should_output .AND. (output_unit > 0)) THEN
774 WRITE (unit=output_unit, fmt="(T3,A,I0,A)") &
775 "SCCS| Iteration cycle converged in ", iter, " steps"
776 END IF
777 EXIT iter_loop
778 END IF
779
780 END DO iter_loop
781
782 IF (sccs_control%method_id == sccs_saa_andreussi) THEN
783 CALL auxbas_pw_pool%create_pw(tem_convolution)
784 CALL auxbas_pw_pool%create_pw(tem_func)
785
786!$OMP PARALLEL DO DEFAULT(NONE) &
787!$OMP PRIVATE(i,j,k,dphi2) &
788!$OMP SHARED(lb,ub, tem_func, s, dtf) &
789!$OMP SHARED(dphi_tot, p_e_interface, eps_elec, epsilon_solvent)
790 DO k = lb(3), ub(3)
791 DO j = lb(2), ub(2)
792 DO i = lb(1), ub(1)
793 p_e_interface%array(i, j, k) = 1.0_dp/(4.0_dp*twopi)*(exp(eps_elec%array(i, j, k))*log(epsilon_solvent))
794 dphi2 = dphi_tot(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + &
795 dphi_tot(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + &
796 dphi_tot(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k)
797 p_e_interface%array(i, j, k) = p_e_interface%array(i, j, k)*dphi2
798 tem_func%array(i, j, k) = (1.0_dp - s%array(i, j, k))*p_e_interface%array(i, j, k)*dtf%array(i, j, k)
799 END DO
800 END DO
801 END DO
802!$OMP END PARALLEL DO
803
804 CALL pw_func_u_convolution(poisson_env=poisson_env, func=tem_func, convolution=tem_convolution, u=u)
805
806!$OMP PARALLEL DO DEFAULT(NONE) &
807!$OMP PRIVATE(i,j,k) &
808!$OMP SHARED(lb,ub) &
809!$OMP SHARED(p_e_interface, t, tem_convolution, phi_tot, d_s_rhoel)
810 DO k = lb(3), ub(3)
811 DO j = lb(2), ub(2)
812 DO i = lb(1), ub(1)
813 p_e_interface%array(i, j, k) = p_e_interface%array(i, j, k)*(1.0_dp - t%array(i, j, k)) &
814 + tem_convolution%array(i, j, k)
815 END DO
816 END DO
817 END DO
818!$OMP END PARALLEL DO
819
820 END IF
821
822 ! Release work storage which is no longer needed
823 CALL auxbas_pw_pool%give_back_pw(rho_tot_zero)
824 DO i = 1, 3
825 CALL auxbas_pw_pool%give_back_pw(dln_eps_elec(i))
826 END DO
827
828 ! Optional output of the total charge density in cube file format
829 filename = "TOTAL_CHARGE_DENSITY"
830 cube_path = trim(print_path)//"%"//trim(filename)
831 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), cp_p_file)) THEN
832 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
833 my_pos_cube = "REWIND"
834 IF (append_cube) my_pos_cube = "APPEND"
835 mpi_io = .true.
836 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
837 extension=".cube", middle_name=trim(filename), &
838 file_position=my_pos_cube, log_filename=.false., &
839 mpi_io=mpi_io, fout=mpi_filename)
840 IF (output_unit > 0) THEN
841 IF (.NOT. mpi_io) THEN
842 INQUIRE (unit=cube_unit, name=filename)
843 ELSE
844 filename = mpi_filename
845 END IF
846 WRITE (unit=output_unit, fmt="(T3,A)") &
847 "SCCS| The total SCCS charge density is written in cube file format to the file:", &
848 "SCCS| "//trim(filename)
849 END IF
850 CALL cp_pw_to_cube(rho_tot, cube_unit, trim(filename), particles=particles, &
851 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), mpi_io=mpi_io)
852 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
853 END IF
854
855 ! Calculate the total SCCS Hartree energy, potential, and its
856 ! derivatives of the solute and the implicit solvent
857 CALL pw_transfer(rho_tot, rho_tot_gspace)
858 IF (calculate_stress_tensor) THEN
859 ! Request also the calculation of the stress tensor contribution
860 CALL pw_poisson_solve(poisson_env=poisson_env, &
861 density=rho_tot_gspace, &
862 ehartree=e_tot, &
863 vhartree=v_hartree_gspace, &
864 dvhartree=dphi_tot, &
865 h_stress=h_stress)
866 ELSE
867 CALL pw_poisson_solve(poisson_env=poisson_env, &
868 density=rho_tot_gspace, &
869 ehartree=e_tot, &
870 vhartree=v_hartree_gspace, &
871 dvhartree=dphi_tot)
872 END IF
873 CALL pw_transfer(v_hartree_gspace, phi_tot)
874 energy%sccs_hartree = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
875
876 ! Calculate the Hartree energy and potential of the solute only
877 block
878 TYPE(pw_r3d_rs_type) :: phi_solute
879 CALL auxbas_pw_pool%create_pw(phi_solute)
880 CALL pw_zero(phi_solute)
881 CALL pw_poisson_solve(poisson_env=poisson_env, &
882 density=rho_solute, &
883 ehartree=energy%hartree, &
884 vhartree=phi_solute)
885
886 ! Calculate the polarisation potential (store it in phi_tot)
887 ! phi_pol = phi_tot - phi_solute
888 CALL pw_axpy(phi_solute, phi_tot, alpha=-1.0_dp)
889 CALL auxbas_pw_pool%give_back_pw(phi_solute)
890 END block
891
892 ! Optional output of the SCCS polarisation potential in cube file format
893 filename = "POLARISATION_POTENTIAL"
894 cube_path = trim(print_path)//"%"//trim(filename)
895 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
896 cp_p_file)) THEN
897 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
898 my_pos_cube = "REWIND"
899 IF (append_cube) my_pos_cube = "APPEND"
900 mpi_io = .true.
901 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
902 extension=".cube", middle_name=trim(filename), &
903 file_position=my_pos_cube, log_filename=.false., &
904 mpi_io=mpi_io, fout=mpi_filename)
905 IF (output_unit > 0) THEN
906 IF (.NOT. mpi_io) THEN
907 INQUIRE (unit=cube_unit, name=filename)
908 ELSE
909 filename = mpi_filename
910 END IF
911 WRITE (unit=output_unit, fmt="(T3,A)") &
912 "SCCS| The SCCS polarisation potential is written in cube file format to the file:", &
913 "SCCS| "//trim(filename)
914 END IF
915 CALL cp_pw_to_cube(phi_tot, cube_unit, trim(filename), particles=particles, &
916 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), mpi_io=mpi_io)
917 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
918 END IF
919
920 ! Calculate the polarisation charge (store it in rho_tot)
921 ! rho_pol = rho_tot - rho_solute
922 CALL pw_axpy(rho_solute, rho_tot, alpha=-1.0_dp)
923 polarisation_charge = pw_integrate_function(rho_tot)
924
925 ! Optional output of the SCCS polarisation charge in cube file format
926 filename = "POLARISATION_CHARGE_DENSITY"
927 cube_path = trim(print_path)//"%"//trim(filename)
928 IF (btest(cp_print_key_should_output(logger%iter_info, input, trim(cube_path)), &
929 cp_p_file)) THEN
930 append_cube = section_get_lval(input, trim(cube_path)//"%APPEND")
931 my_pos_cube = "REWIND"
932 IF (append_cube) my_pos_cube = "APPEND"
933 mpi_io = .true.
934 cube_unit = cp_print_key_unit_nr(logger, input, trim(cube_path), &
935 extension=".cube", middle_name=trim(filename), &
936 file_position=my_pos_cube, log_filename=.false., &
937 mpi_io=mpi_io, fout=mpi_filename)
938 IF (output_unit > 0) THEN
939 IF (.NOT. mpi_io) THEN
940 INQUIRE (unit=cube_unit, name=filename)
941 ELSE
942 filename = mpi_filename
943 END IF
944 WRITE (unit=output_unit, fmt="(T3,A)") &
945 "SCCS| The SCCS polarisation charge density is written in cube file format to the file:", &
946 "SCCS| "//trim(filename)
947 END IF
948 CALL cp_pw_to_cube(rho_tot, cube_unit, trim(filename), particles=particles, &
949 stride=section_get_ivals(input, trim(cube_path)//"%STRIDE"), mpi_io=mpi_io)
950 CALL cp_print_key_finished_output(cube_unit, logger, input, trim(cube_path), mpi_io=mpi_io)
951 END IF
952
953 ! Calculate SCCS polarisation energy
954 energy%sccs_pol = 0.5_dp*pw_integral_ab(rho_solute, phi_tot)
955 CALL auxbas_pw_pool%give_back_pw(rho_solute)
956 CALL auxbas_pw_pool%give_back_pw(phi_tot)
957 CALL auxbas_pw_pool%give_back_pw(rho_tot)
958 IF (sccs_control%method_id == sccs_saa_andreussi) THEN
959 CALL auxbas_pw_pool%give_back_pw(s)
960 CALL auxbas_pw_pool%give_back_pw(u)
961 CALL auxbas_pw_pool%give_back_pw(ff)
962 CALL auxbas_pw_pool%give_back_pw(t)
963 CALL auxbas_pw_pool%give_back_pw(dtf)
964 CALL auxbas_pw_pool%give_back_pw(tem_convolution)
965 CALL auxbas_pw_pool%give_back_pw(tem_func)
966 CALL auxbas_pw_pool%give_back_pw(eps_elec)
967 END IF
968 END block
969
970 ! Calculate additional solvation terms
971 energy%sccs_cav = sccs_control%gamma_solvent*cavity_surface
972 energy%sccs_dis = sccs_control%beta_solvent*cavity_volume
973 energy%sccs_rep = sccs_control%alpha_solvent*cavity_surface
974 ! Calculate solvation free energy: \delta G^el + (alpha + gamma)*S + beta*V
975 energy%sccs_sol = energy%sccs_pol + energy%sccs_rep + energy%sccs_cav + energy%sccs_dis
976
977 IF (should_output .AND. (output_unit > 0)) THEN
978 WRITE (unit=output_unit, fmt="(T3,A)") &
979 "SCCS|"
980 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.12)") &
981 "SCCS| Polarisation charge", polarisation_charge
982 !MK "SCCS| Total interaction energy [a.u.]", e_tot
983 WRITE (unit=output_unit, fmt="(T3,A)") &
984 "SCCS|"
985 CALL print_sccs_results(energy, sccs_control, output_unit)
986 END IF
987
988 ! Calculate SCCS contribution to the Kohn-Sham potential
989 IF (sccs_control%method_id /= sccs_saa_andreussi) THEN
990
991 f = -0.5_dp*dvol/fourpi
992!$OMP PARALLEL DO DEFAULT(NONE) &
993!$OMP PRIVATE(dphi2,i,j,k) &
994!$OMP SHARED(f,deps_elec,dphi_tot) &
995!$OMP SHARED(lb,ub,v_sccs)
996 DO k = lb(3), ub(3)
997 DO j = lb(2), ub(2)
998 DO i = lb(1), ub(1)
999 dphi2 = dphi_tot(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + &
1000 dphi_tot(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + &
1001 dphi_tot(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k)
1002 v_sccs%array(i, j, k) = v_sccs%array(i, j, k) + f*deps_elec%array(i, j, k)*dphi2
1003 END DO
1004 END DO
1005 END DO
1006!$OMP END PARALLEL DO
1007
1008 ELSE
1009
1010!$OMP PARALLEL DO DEFAULT(NONE) &
1011!$OMP PRIVATE(i,j,k) &
1012!$OMP SHARED(lb,ub, dvol) &
1013!$OMP SHARED(p_E_interface, v_sccs, d_s_rhoel)
1014 DO k = lb(3), ub(3)
1015 DO j = lb(2), ub(2)
1016 DO i = lb(1), ub(1)
1017 v_sccs%array(i, j, k) = dvol*d_s_rhoel%array(i, j, k)*p_e_interface%array(i, j, k)
1018 END DO
1019 END DO
1020 END DO
1021!$OMP END PARALLEL DO
1022
1023 CALL auxbas_pw_pool%give_back_pw(d_s_rhoel)
1024 CALL auxbas_pw_pool%give_back_pw(p_e_interface)
1025
1026 END IF
1027
1028 CALL auxbas_pw_pool%give_back_pw(deps_elec)
1029 DO i = 1, 3
1030 CALL auxbas_pw_pool%give_back_pw(dphi_tot(i))
1031 END DO
1032
1033 ! Release the SCCS printout environment
1034 CALL cp_print_key_finished_output(output_unit, logger, input, trim(print_path), &
1035 ignore_should_output=should_output)
1036
1037 CALL timestop(handle)
1038
1039 END SUBROUTINE sccs
1040
1041! **************************************************************************************************
1042!> \brief Calculate the smoothed dielectric function of Andreussi et al.
1043!> \param rho_elec ...
1044!> \param eps_elec ...
1045!> \param deps_elec ...
1046!> \param epsilon_solvent ...
1047!> \param rho_max ...
1048!> \param rho_min ...
1049!> \par History:
1050!> - Creation (16.10.2013,MK)
1051!> - Finite difference of isosurfaces implemented (21.12.2013,MK)
1052!> \author Matthias Krack (MK)
1053!> \version 1.1
1054! **************************************************************************************************
1055 SUBROUTINE andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, rho_max, &
1056 rho_min)
1057
1058 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec
1059 REAL(kind=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min
1060
1061 CHARACTER(LEN=*), PARAMETER :: routinen = 'andreussi'
1062 REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1063
1064 INTEGER :: handle, i, j, k
1065 INTEGER, DIMENSION(3) :: lb, ub
1066 REAL(kind=dp) :: diff, dq, dt, f, ln_rho_max, ln_rho_min, &
1067 q, rho, t, x, y
1068
1069 CALL timeset(routinen, handle)
1070
1071 f = log(epsilon_solvent)/twopi
1072 diff = rho_max - rho_min
1073 IF (diff < sqrt(rhotol)) cpabort("SCCS: Difference between rho(min) and rho(max) is too small")
1074 IF (rho_min >= rhotol) THEN
1075 ln_rho_max = log(rho_max)
1076 ln_rho_min = log(rho_min)
1077 q = twopi/(ln_rho_max - ln_rho_min)
1078 dq = -f*q
1079 END IF
1080
1081 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1082 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1083
1084 ! Calculate the dielectric function and its derivative
1085!$OMP PARALLEL DO DEFAULT(NONE) &
1086!$OMP PRIVATE(dt,i,j,k,rho,t,x,y) &
1087!$OMP SHARED(deps_elec,dq,eps_elec,epsilon_solvent,f,lb,ub) &
1088!$OMP SHARED(ln_rho_max,rho_elec,q,rho_max,rho_min)
1089 DO k = lb(3), ub(3)
1090 DO j = lb(2), ub(2)
1091 DO i = lb(1), ub(1)
1092 rho = rho_elec%array(i, j, k)
1093 IF (rho < rho_min) THEN
1094 eps_elec%array(i, j, k) = epsilon_solvent
1095 deps_elec%array(i, j, k) = 0.0_dp
1096 ELSE IF (rho <= rho_max) THEN
1097 x = log(rho)
1098 y = q*(ln_rho_max - x)
1099 t = f*(y - sin(y))
1100 eps_elec%array(i, j, k) = exp(t)
1101 dt = dq*(1.0_dp - cos(y))
1102 deps_elec%array(i, j, k) = eps_elec%array(i, j, k)*dt/rho
1103 ELSE
1104 eps_elec%array(i, j, k) = 1.0_dp
1105 deps_elec%array(i, j, k) = 0.0_dp
1106 END IF
1107 END DO
1108 END DO
1109 END DO
1110!$OMP END PARALLEL DO
1111
1112 CALL timestop(handle)
1113
1114 END SUBROUTINE andreussi
1115
1116! **************************************************************************************************
1117!> \brief Calculate the smoothed dielectric function of Fattebert and Gygi
1118!> \param rho_elec ...
1119!> \param eps_elec ...
1120!> \param deps_elec ...
1121!> \param epsilon_solvent ...
1122!> \param beta ...
1123!> \param rho_zero ...
1124!> \par History:
1125!> - Creation (15.10.2013,MK)
1126!> \author Matthias Krack (MK)
1127!> \version 1.0
1128! **************************************************************************************************
1129 SUBROUTINE fattebert_gygi(rho_elec, eps_elec, deps_elec, epsilon_solvent, beta, &
1130 rho_zero)
1131
1132 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec
1133 REAL(kind=dp), INTENT(IN) :: epsilon_solvent, beta, rho_zero
1134
1135 CHARACTER(LEN=*), PARAMETER :: routinen = 'fattebert_gygi'
1136 REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1137
1138 INTEGER :: handle, i, j, k
1139 INTEGER, DIMENSION(3) :: lb, ub
1140 REAL(kind=dp) :: df, f, p, q, rho, s, t, twobeta
1141
1142 CALL timeset(routinen, handle)
1143
1144 df = (1.0_dp - epsilon_solvent)/rho_zero
1145 f = 0.5_dp*(epsilon_solvent - 1.0_dp)
1146 q = 1.0_dp/rho_zero
1147 twobeta = 2.0_dp*beta
1148
1149 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1150 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1151
1152 ! Calculate the smoothed dielectric function and its derivative
1153!$OMP PARALLEL DO DEFAULT(NONE) &
1154!$OMP PRIVATE(i,j,k,p,rho,s,t) &
1155!$OMP SHARED(df,deps_elec,eps_elec,epsilon_solvent,f,lb,ub) &
1156!$OMP SHARED(q,rho_elec,twobeta)
1157 DO k = lb(3), ub(3)
1158 DO j = lb(2), ub(2)
1159 DO i = lb(1), ub(1)
1160 rho = rho_elec%array(i, j, k)
1161 IF (rho < rhotol) THEN
1162 eps_elec%array(i, j, k) = epsilon_solvent
1163 deps_elec%array(i, j, k) = 0.0_dp
1164 ELSE
1165 s = rho*q
1166 p = s**twobeta
1167 t = 1.0_dp/(1.0_dp + p)
1168 eps_elec%array(i, j, k) = 1.0_dp + f*(1.0_dp + (1.0_dp - p)*t)
1169 deps_elec%array(i, j, k) = df*twobeta*t*t*p/s
1170 END IF
1171 END DO
1172 END DO
1173 END DO
1174!$OMP END PARALLEL DO
1175
1176 CALL timestop(handle)
1177
1178 END SUBROUTINE fattebert_gygi
1179
1180! **************************************************************************************************
1181!> \brief Calculate the smoothed dielectric function of solvent aware algorithm
1182!> \param rho_elec ...
1183!> \param eps_elec ...
1184!> \param deps_elec ...
1185!> \param epsilon_solvent ...
1186!> \param rho_max ...
1187!> \param rho_min ...
1188!> \param s ...
1189!> \param d_s_rhoel ...
1190!> \author Ziwei Chai
1191!> \version 1.0
1192! **************************************************************************************************
1193 SUBROUTINE sa_andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, rho_max, &
1194 rho_min, s, d_s_rhoel)
1195
1196 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec
1197 REAL(kind=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min
1198 TYPE(pw_r3d_rs_type), INTENT(IN) :: s, d_s_rhoel
1199
1200 CHARACTER(LEN=*), PARAMETER :: routinen = 'sa_andreussi'
1201 REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1202
1203 INTEGER :: handle, i, j, k
1204 INTEGER, DIMENSION(3) :: lb, ub
1205 REAL(kind=dp) :: diff, dq, dt, f, ln_rho_max, ln_rho_min, &
1206 q, rho, t, x, y
1207
1208 CALL timeset(routinen, handle)
1209
1210 f = log(epsilon_solvent)/twopi
1211 diff = rho_max - rho_min
1212 IF (diff < sqrt(rhotol)) cpabort("SCCS: Difference between rho(min) and rho(max) is too small")
1213 IF (rho_min >= rhotol) THEN
1214 ln_rho_max = log(rho_max)
1215 ln_rho_min = log(rho_min)
1216 q = twopi/(ln_rho_max - ln_rho_min)
1217 dq = -f*q
1218 END IF
1219
1220 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1221 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1222
1223 ! Calculate the dielectric function and its derivative
1224!$OMP PARALLEL DO DEFAULT(NONE) &
1225!$OMP PRIVATE(dt,i,j,k,rho,t,x,y) &
1226!$OMP SHARED(deps_elec,dq,eps_elec,epsilon_solvent,f,lb,ub) &
1227!$OMP SHARED(ln_rho_max,rho_elec,q,rho_max,rho_min, s, d_s_rhoel)
1228 DO k = lb(3), ub(3)
1229 DO j = lb(2), ub(2)
1230 DO i = lb(1), ub(1)
1231 rho = rho_elec%array(i, j, k)
1232 IF (rho < rho_min) THEN
1233 eps_elec%array(i, j, k) = epsilon_solvent
1234 s%array(i, j, k) = 0.0_dp
1235 deps_elec%array(i, j, k) = 0.0_dp
1236 d_s_rhoel%array(i, j, k) = 0.0_dp
1237 ELSE IF (rho <= rho_max) THEN
1238 x = log(rho)
1239 y = q*(ln_rho_max - x)
1240 t = f*(y - sin(y))
1241 s%array(i, j, k) = 1 - t/(f*twopi)
1242 eps_elec%array(i, j, k) = exp(t)
1243 dt = dq*(1.0_dp - cos(y))
1244 deps_elec%array(i, j, k) = eps_elec%array(i, j, k)*dt/rho
1245 d_s_rhoel%array(i, j, k) = deps_elec%array(i, j, k)/(-eps_elec%array(i, j, k)*log(epsilon_solvent))
1246 ELSE
1247 eps_elec%array(i, j, k) = 1.0_dp
1248 deps_elec%array(i, j, k) = 0.0_dp
1249 s%array(i, j, k) = 1.0_dp
1250 d_s_rhoel%array(i, j, k) = 0.0_dp
1251 END IF
1252 END DO
1253 END DO
1254 END DO
1255!$OMP END PARALLEL DO
1256
1257 CALL timestop(handle)
1258
1259 END SUBROUTINE sa_andreussi
1260
1261! **************************************************************************************************
1262!> \brief Build the numerical derivative of a function on realspace grid
1263!> \param f ...
1264!> \param df ...
1265!> \param method ...
1266!> \param pw_env ...
1267!> \param input ...
1268!> \par History:
1269!> - Creation (15.11.2013,MK)
1270!> \author Matthias Krack (MK)
1271!> \version 1.0
1272! **************************************************************************************************
1273 SUBROUTINE derive(f, df, method, pw_env, input)
1274
1275 TYPE(pw_r3d_rs_type), INTENT(IN) :: f
1276 TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: df
1277 INTEGER, INTENT(IN) :: method
1278 TYPE(pw_env_type), POINTER :: pw_env
1279 TYPE(section_vals_type), POINTER :: input
1280
1281 CHARACTER(LEN=*), PARAMETER :: routinen = 'derive'
1282
1283 INTEGER :: border_points, handle, i
1284 INTEGER, DIMENSION(3) :: lb, n, ub
1285 TYPE(pw_c1d_gs_type), DIMENSION(2) :: work_g1d
1286 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1287 TYPE(realspace_grid_desc_type), POINTER :: rs_desc
1288 TYPE(realspace_grid_input_type) :: input_settings
1289 TYPE(realspace_grid_type), POINTER :: rs_grid
1290 TYPE(section_vals_type), POINTER :: rs_grid_section
1291
1292 CALL timeset(routinen, handle)
1293
1294 cpassert(ASSOCIATED(pw_env))
1295
1296 ! Perform method specific setup
1297 SELECT CASE (method)
1298 CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7)
1299 NULLIFY (rs_desc)
1300 rs_grid_section => section_vals_get_subs_vals(input, "DFT%MGRID%RS_GRID")
1301 SELECT CASE (method)
1302 CASE (sccs_derivative_cd3)
1303 border_points = 1
1304 CASE (sccs_derivative_cd5)
1305 border_points = 2
1306 CASE (sccs_derivative_cd7)
1307 border_points = 3
1308 END SELECT
1309 CALL init_input_type(input_settings, 2*border_points + 1, rs_grid_section, &
1310 1, [-1, -1, -1])
1311 CALL rs_grid_create_descriptor(rs_desc, f%pw_grid, input_settings, &
1312 border_points=border_points)
1313 ALLOCATE (rs_grid)
1314 CALL rs_grid_create(rs_grid, rs_desc)
1315!MK CALL rs_grid_print(rs_grid, 6)
1316 CASE (sccs_derivative_fft)
1317 lb(1:3) = f%pw_grid%bounds_local(1, 1:3)
1318 ub(1:3) = f%pw_grid%bounds_local(2, 1:3)
1319 NULLIFY (auxbas_pw_pool)
1320 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1321 ! Get work storage for the 1d grids in g-space (derivative calculation)
1322 DO i = 1, SIZE(work_g1d)
1323 CALL auxbas_pw_pool%create_pw(work_g1d(i))
1324 END DO
1325 END SELECT
1326
1327 ! Calculate the derivatives
1328 SELECT CASE (method)
1329 CASE (sccs_derivative_cd3)
1330 CALL derive_fdm_cd3(f, df, rs_grid)
1331 CASE (sccs_derivative_cd5)
1332 CALL derive_fdm_cd5(f, df, rs_grid)
1333 CASE (sccs_derivative_cd7)
1334 CALL derive_fdm_cd7(f, df, rs_grid)
1335 CASE (sccs_derivative_fft)
1336 ! FFT
1337 CALL pw_transfer(f, work_g1d(1))
1338 DO i = 1, 3
1339 n(:) = 0
1340 n(i) = 1
1341 CALL pw_copy(work_g1d(1), work_g1d(2))
1342 CALL pw_derive(work_g1d(2), n(:))
1343 CALL pw_transfer(work_g1d(2), df(i))
1344 END DO
1345 CASE DEFAULT
1346 cpabort("Invalid derivative method for SCCS specified")
1347 END SELECT
1348
1349 ! Perform method specific cleanup
1350 SELECT CASE (method)
1351 CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7)
1352 CALL rs_grid_release(rs_grid)
1353 DEALLOCATE (rs_grid)
1354 CALL rs_grid_release_descriptor(rs_desc)
1355 CASE (sccs_derivative_fft)
1356 DO i = 1, SIZE(work_g1d)
1357 CALL auxbas_pw_pool%give_back_pw(work_g1d(i))
1358 END DO
1359 END SELECT
1360
1361 CALL timestop(handle)
1362
1363 END SUBROUTINE derive
1364
1365! **************************************************************************************************
1366!> \brief Calculate the finite difference between two isosurfaces of the
1367!> electronic density. The smoothed dielectric function of
1368!> Andreussi et al. is used as switching function eventually
1369!> defining the quantum volume and surface of the cavity.
1370!> \param rho_elec ...
1371!> \param norm_drho_elec ...
1372!> \param dtheta ...
1373!> \param epsilon_solvent ...
1374!> \param rho_max ...
1375!> \param rho_min ...
1376!> \param delta_rho ...
1377!> \par History:
1378!> - Creation (21.12.2013,MK)
1379!> \author Matthias Krack (MK)
1380!> \version 1.0
1381! **************************************************************************************************
1382 SUBROUTINE surface_andreussi(rho_elec, norm_drho_elec, dtheta, &
1383 epsilon_solvent, rho_max, rho_min, delta_rho)
1384
1385 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, norm_drho_elec, dtheta
1386 REAL(kind=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min, &
1387 delta_rho
1388
1389 CHARACTER(LEN=*), PARAMETER :: routinen = 'surface_andreussi'
1390 REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1391
1392 INTEGER :: handle, i, j, k, l
1393 INTEGER, DIMENSION(3) :: lb, ub
1394 REAL(kind=dp) :: diff, e, eps_elec, f, ln_rho_max, &
1395 ln_rho_min, q, rho, t, x, y
1396 REAL(kind=dp), DIMENSION(2) :: theta
1397
1398 CALL timeset(routinen, handle)
1399
1400 e = epsilon_solvent - 1.0_dp
1401 f = log(epsilon_solvent)/twopi
1402 diff = rho_max - rho_min
1403 IF (diff < sqrt(rhotol)) cpabort("SCCS: Difference between rho(min) and rho(max) is too small")
1404 IF (rho_min >= rhotol) THEN
1405 ln_rho_max = log(rho_max)
1406 ln_rho_min = log(rho_min)
1407 q = twopi/(ln_rho_max - ln_rho_min)
1408 END IF
1409
1410 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1411 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1412
1413 ! Calculate finite difference between two isosurfaces
1414!$OMP PARALLEL DO DEFAULT(NONE) &
1415!$OMP PRIVATE(eps_elec,i,j,k,l,rho,t,theta,x,y) &
1416!$OMP SHARED(delta_rho,dtheta,e,epsilon_solvent,f,lb) &
1417!$OMP SHARED(ln_rho_max,norm_drho_elec,rho_elec,q,rho_max,rho_min,ub)
1418 DO k = lb(3), ub(3)
1419 DO j = lb(2), ub(2)
1420 DO i = lb(1), ub(1)
1421 DO l = 1, 2
1422 rho = rho_elec%array(i, j, k) + (real(l, kind=dp) - 1.5_dp)*delta_rho
1423 IF (rho < rho_min) THEN
1424 eps_elec = epsilon_solvent
1425 ELSE IF (rho <= rho_max) THEN
1426 x = log(rho)
1427 y = q*(ln_rho_max - x)
1428 t = f*(y - sin(y))
1429 eps_elec = exp(t)
1430 ELSE
1431 eps_elec = 1.0_dp
1432 END IF
1433 theta(l) = (epsilon_solvent - eps_elec)/e
1434 END DO
1435 dtheta%array(i, j, k) = (theta(2) - theta(1))*norm_drho_elec%array(i, j, k)/delta_rho
1436 END DO
1437 END DO
1438 END DO
1439!$OMP END PARALLEL DO
1440
1441 CALL timestop(handle)
1442
1443 END SUBROUTINE surface_andreussi
1444
1445! **************************************************************************************************
1446!> \brief Calculate the finite difference between two isosurfaces of the
1447!> the electronic density. The smoothed dielectric function of
1448!> Fattebert and Gygi is used as switching function eventually
1449!> defining the quantum volume and surface of the cavity.
1450!> \param rho_elec ...
1451!> \param norm_drho_elec ...
1452!> \param dtheta ...
1453!> \param epsilon_solvent ...
1454!> \param beta ...
1455!> \param rho_zero ...
1456!> \param delta_rho ...
1457!> \par History:
1458!> - Creation (21.12.2013,MK)
1459!> \author Matthias Krack (MK)
1460!> \version 1.0
1461! **************************************************************************************************
1462 SUBROUTINE surface_fattebert_gygi(rho_elec, norm_drho_elec, dtheta, &
1463 epsilon_solvent, beta, rho_zero, delta_rho)
1464
1465 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, norm_drho_elec, dtheta
1466 REAL(kind=dp), INTENT(IN) :: epsilon_solvent, beta, rho_zero, &
1467 delta_rho
1468
1469 CHARACTER(LEN=*), PARAMETER :: routinen = 'surface_fattebert_gygi'
1470 REAL(kind=dp), PARAMETER :: rhotol = 1.0e-12_dp
1471
1472 INTEGER :: handle, i, j, k, l
1473 INTEGER, DIMENSION(3) :: lb, ub
1474 REAL(kind=dp) :: e, eps_elec, f, p, q, rho, s, t, twobeta
1475 REAL(kind=dp), DIMENSION(2) :: theta
1476
1477 CALL timeset(routinen, handle)
1478
1479 e = epsilon_solvent - 1.0_dp
1480 f = 0.5_dp*e
1481 q = 1.0_dp/rho_zero
1482 twobeta = 2.0_dp*beta
1483
1484 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3)
1485 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3)
1486
1487 ! Calculate finite difference between two isosurfaces
1488!$OMP PARALLEL DO DEFAULT(NONE) &
1489!$OMP PRIVATE(eps_elec,i,j,k,l,p,rho,s,t,theta) &
1490!$OMP SHARED(delta_rho,dtheta,e,epsilon_solvent,f,lb) &
1491!$OMP SHARED(norm_drho_elec,q,rho_elec,twobeta,ub)
1492 DO k = lb(3), ub(3)
1493 DO j = lb(2), ub(2)
1494 DO i = lb(1), ub(1)
1495 DO l = 1, 2
1496 rho = rho_elec%array(i, j, k) + (real(l, kind=dp) - 1.5_dp)*delta_rho
1497 IF (rho < rhotol) THEN
1498 eps_elec = epsilon_solvent
1499 ELSE
1500 s = rho*q
1501 p = s**twobeta
1502 t = 1.0_dp/(1.0_dp + p)
1503 eps_elec = 1.0_dp + f*(1.0_dp + (1.0_dp - p)*t)
1504 END IF
1505 theta(l) = (epsilon_solvent - eps_elec)/e
1506 END DO
1507 dtheta%array(i, j, k) = (theta(2) - theta(1))*norm_drho_elec%array(i, j, k)/delta_rho
1508 END DO
1509 END DO
1510 END DO
1511!$OMP END PARALLEL DO
1512
1513 CALL timestop(handle)
1514
1515 END SUBROUTINE surface_fattebert_gygi
1516
1517! **************************************************************************************************
1518!> \brief Print SCCS results
1519!> \param energy ...
1520!> \param sccs_control ...
1521!> \param output_unit ...
1522!> \par History:
1523!> - Creation (11.11.2022,MK)
1524!> \author Matthias Krack (MK)
1525!> \version 1.0
1526! **************************************************************************************************
1527 SUBROUTINE print_sccs_results(energy, sccs_control, output_unit)
1528
1529 TYPE(qs_energy_type), POINTER :: energy
1530 TYPE(sccs_control_type), POINTER :: sccs_control
1531 INTEGER, INTENT(IN) :: output_unit
1532
1533 IF (output_unit > 0) THEN
1534 cpassert(ASSOCIATED(energy))
1535 cpassert(ASSOCIATED(sccs_control))
1536 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1537 "SCCS| Hartree energy of solute and solvent [Hartree]", energy%sccs_hartree, &
1538 "SCCS| Hartree energy of the solute only [Hartree]", energy%hartree
1539 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14,/,T3,A,T61,F20.3)") &
1540 "SCCS| Polarisation energy [Hartree]", energy%sccs_pol, &
1541 "SCCS| [kcal/mol]", &
1542 cp_unit_from_cp2k(energy%sccs_pol, "kcalmol"), &
1543 "SCCS| Cavitation energy [Hartree]", energy%sccs_cav, &
1544 "SCCS| [kcal/mol]", &
1545 cp_unit_from_cp2k(energy%sccs_cav, "kcalmol"), &
1546 "SCCS| Dispersion free energy [Hartree]", energy%sccs_dis, &
1547 "SCCS| [kcal/mol]", &
1548 cp_unit_from_cp2k(energy%sccs_dis, "kcalmol"), &
1549 "SCCS| Repulsion free energy [Hartree]", energy%sccs_rep, &
1550 "SCCS| [kcal/mol]", &
1551 cp_unit_from_cp2k(energy%sccs_rep, "kcalmol"), &
1552 "SCCS| Solvation free energy [Hartree]", energy%sccs_sol, &
1553 "SCCS| [kcal/mol]", &
1554 cp_unit_from_cp2k(energy%sccs_sol, "kcalmol")
1555 END IF
1556
1557 END SUBROUTINE print_sccs_results
1558
1559! **************************************************************************************************
1560!> \brief Exponents and coefficients of the SCCS core fill Gaussians per atomic kind, such that
1561!> the fill density of an atom is ccore*exp(-calpha*r**2) = Z_val/(sqrt(pi)*sigma)**3*
1562!> exp(-r**2/sigma**2). A zero coefficient means that the kind is not filled.
1563!> \param qs_env ...
1564!> \param calpha Gaussian exponents 1/sigma**2 per kind
1565!> \param ccore Gaussian coefficients per kind
1566!> \param output_unit if > 0, the fill parameters are printed once
1567! **************************************************************************************************
1568 SUBROUTINE sccs_core_fill_coefficients(qs_env, calpha, ccore, output_unit)
1569
1570 TYPE(qs_environment_type), POINTER :: qs_env
1571 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
1572 INTENT(OUT) :: calpha, ccore
1573 INTEGER, INTENT(IN), OPTIONAL :: output_unit
1574
1575 CHARACTER(LEN=default_string_length) :: kind_name
1576 INTEGER :: ikind, iover, nkind, z
1577 LOGICAL :: found
1578 LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1579 REAL(kind=dp) :: r_thr, rho_thr, sigma, zeff
1580 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1581 TYPE(dft_control_type), POINTER :: dft_control
1582 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1583 TYPE(sccs_control_type), POINTER :: sccs_control
1584
1585 NULLIFY (atomic_kind_set, dft_control, qs_kind_set)
1586 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1587 dft_control=dft_control)
1588 sccs_control => dft_control%sccs_control
1589 nkind = SIZE(atomic_kind_set)
1590 ALLOCATE (calpha(nkind), ccore(nkind), used(SIZE(sccs_control%core_fill_kind_label)))
1591 used(:) = .false.
1592
1593 ! Density threshold used only for the diagnostic radius of the fill:
1594 ! the onset of the dielectric transition of the chosen method
1595 SELECT CASE (sccs_control%method_id)
1596 CASE (sccs_andreussi, sccs_saa_andreussi)
1597 rho_thr = sccs_control%rho_min
1598 CASE (sccs_fattebert_gygi)
1599 rho_thr = sccs_control%rho_zero
1600 CASE DEFAULT
1601 rho_thr = 0.0_dp
1602 END SELECT
1603
1604 IF (PRESENT(output_unit)) THEN
1605 IF ((output_unit > 0) .AND. (.NOT. sccs_control%core_fill_printed)) THEN
1606 WRITE (unit=output_unit, fmt="(/,T3,A,/,T3,A,ES10.3,/,T3,A)") &
1607 "SCCS| Core fill: Z_val/(sqrt(pi)*sigma)**3*exp(-r**2/sigma**2) added for the cavity", &
1608 "SCCS| r_thr: radius at which the fill density drops to rho_thr =", rho_thr, &
1609 "SCCS| Kind Z Z_val sigma [bohr] r_thr [Angstrom]"
1610 END IF
1611 END IF
1612
1613 DO ikind = 1, nkind
1614 CALL get_atomic_kind(atomic_kind_set(ikind), name=kind_name, z=z)
1615 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
1616 CALL uppercase(kind_name)
1617 ! Default spread. Kinds without removed core electrons (Z_val = Z, e.g. H-q1, He-q2,
1618 ! Li-q3, Be-q4) have no core hole in the valence density and are not filled
1619 ! (generalises the hydrogen exception of Environ)
1620 sigma = sccs_control%core_fill_spread
1621 IF (zeff >= real(z, kind=dp) - 0.5_dp) sigma = 0.0_dp
1622 found = .false.
1623 DO iover = 1, SIZE(sccs_control%core_fill_kind_label)
1624 IF (trim(sccs_control%core_fill_kind_label(iover)) == trim(kind_name)) THEN
1625 sigma = sccs_control%core_fill_kind_spread(iover)
1626 used(iover) = .true.
1627 found = .true.
1628 END IF
1629 END DO
1630 IF ((sigma > 0.0_dp) .AND. (zeff > 0.0_dp)) THEN
1631 calpha(ikind) = 1.0_dp/(sigma*sigma)
1632 ccore(ikind) = zeff*(calpha(ikind)/pi)**1.5_dp
1633 ELSE
1634 calpha(ikind) = 0.0_dp
1635 ccore(ikind) = 0.0_dp
1636 END IF
1637 IF (PRESENT(output_unit)) THEN
1638 IF ((output_unit > 0) .AND. (.NOT. sccs_control%core_fill_printed)) THEN
1639 IF ((ccore(ikind) > 0.0_dp) .AND. (rho_thr > 0.0_dp)) THEN
1640 r_thr = 0.0_dp
1641 IF (ccore(ikind) > rho_thr) THEN
1642 r_thr = sqrt(log(ccore(ikind)/rho_thr)/calpha(ikind))
1643 END IF
1644 WRITE (unit=output_unit, fmt="(T3,A,T10,A,T18,I3,T24,F6.2,T34,F8.4,T50,F8.4,A)") &
1645 "SCCS|", trim(kind_name), z, zeff, sigma, cp_unit_from_cp2k(r_thr, "angstrom"), &
1646 merge(" (input)", " ", found)
1647 ELSE IF (ccore(ikind) > 0.0_dp) THEN
1648 WRITE (unit=output_unit, fmt="(T3,A,T10,A,T18,I3,T24,F6.2,T34,F8.4,T50,A)") &
1649 "SCCS|", trim(kind_name), z, zeff, sigma, "-"//merge(" (input)", " ", found)
1650 ELSE
1651 WRITE (unit=output_unit, fmt="(T3,A,T10,A,T18,I3,T24,F6.2,T34,A)") &
1652 "SCCS|", trim(kind_name), z, zeff, "no fill"//merge(" (input)", " ", found)
1653 END IF
1654 END IF
1655 END IF
1656 END DO
1657
1658 DO iover = 1, SIZE(used)
1659 IF (.NOT. used(iover)) THEN
1660 CALL cp_abort(__location__, "SCCS| CORE_FILL_KIND: unknown kind label "// &
1661 trim(sccs_control%core_fill_kind_label(iover)))
1662 END IF
1663 END DO
1664 IF (PRESENT(output_unit)) THEN
1665 IF (output_unit > 0) sccs_control%core_fill_printed = .true.
1666 END IF
1667 DEALLOCATE (used)
1668
1669 END SUBROUTINE sccs_core_fill_coefficients
1670
1671END MODULE qs_sccs
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_integrate_function
Definition pw_methods.F:116
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
qs_collocate_density::calculate_rho_core
Definition qs_collocate_density.F:138
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:131
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::low_print_level
integer, parameter, public low_print_level
Definition cp_output_handling.F:77
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1156
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:107
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout, file_access)
...
Definition cp_output_handling.F:935
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
...
Definition cp_realspace_grid_cube.F:48
cp_realspace_grid_init
Definition cp_realspace_grid_init.F:15
cp_realspace_grid_init::init_input_type
subroutine, public init_input_type(input_settings, nsmax, rs_grid_section, ilevel, higher_grid_layout)
parses an input section to assign the proper values to the input type
Definition cp_realspace_grid_init.F:48
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1311
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::sccs_derivative_cd5
integer, parameter, public sccs_derivative_cd5
Definition input_constants.F:1277
input_constants::sccs_saa_andreussi
integer, parameter, public sccs_saa_andreussi
Definition input_constants.F:1272
input_constants::sccs_fattebert_gygi
integer, parameter, public sccs_fattebert_gygi
Definition input_constants.F:1272
input_constants::sccs_derivative_cd7
integer, parameter, public sccs_derivative_cd7
Definition input_constants.F:1277
input_constants::sccs_derivative_fft
integer, parameter, public sccs_derivative_fft
Definition input_constants.F:1277
input_constants::sccs_derivative_cd3
integer, parameter, public sccs_derivative_cd3
Definition input_constants.F:1277
input_constants::sccs_andreussi
integer, parameter, public sccs_andreussi
Definition input_constants.F:1272
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_get_ivals
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
Definition input_section_types.F:1001
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition input_section_types.F:1017
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::oorootpi
real(kind=dp), parameter, public oorootpi
Definition mathconstants.F:115
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition mathconstants.F:113
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_methods::pw_derive
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
Definition pw_methods.F:10149
pw_poisson_methods
Definition pw_poisson_methods.F:13
pw_poisson_methods::pw_func_u_convolution
subroutine, public pw_func_u_convolution(poisson_env, func, convolution, u)
computes the convolution f * u using reciprocal-space multiplication
Definition pw_poisson_methods.F:4306
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_poisson_types::pw_poisson_mt
integer, parameter, public pw_poisson_mt
Definition pw_poisson_types.F:67
pw_poisson_types::pw_poisson_analytic
integer, parameter, public pw_poisson_analytic
Definition pw_poisson_types.F:67
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:483
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_sccs
Self-consistent continuum solvation (SCCS) model implementation.
Definition qs_sccs.F:29
qs_sccs::sccs
subroutine, public sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs, h_stress)
Self-consistent continuum solvation (SCCS) model implementation.
Definition qs_sccs.F:134
qs_sccs::print_sccs_results
subroutine, public print_sccs_results(energy, sccs_control, output_unit)
Print SCCS results.
Definition qs_sccs.F:1528
qs_sccs::sccs_core_fill_coefficients
subroutine, public sccs_core_fill_coefficients(qs_env, calpha, ccore, output_unit)
Exponents and coefficients of the SCCS core fill Gaussians per atomic kind, such that the fill densit...
Definition qs_sccs.F:1569
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
realspace_grid_types
Definition realspace_grid_types.F:18
realspace_grid_types::rs_grid_create
subroutine, public rs_grid_create(rs, desc)
...
Definition realspace_grid_types.F:478
realspace_grid_types::rs_grid_create_descriptor
subroutine, public rs_grid_create_descriptor(desc, pw_grid, input_settings, border_points)
Determine the setup of real space grids - this is divided up into the creation of a descriptor and th...
Definition realspace_grid_types.F:202
realspace_grid_types::rs_grid_release_descriptor
subroutine, public rs_grid_release_descriptor(rs_desc)
releases the given rs grid descriptor (see doc/ReferenceCounting.html)
Definition realspace_grid_types.F:2040
realspace_grid_types::rs_grid_release
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
Definition realspace_grid_types.F:2020
rs_methods
numerical operations on real-space grid
Definition rs_methods.F:14
rs_methods::derive_fdm_cd7
subroutine, public derive_fdm_cd7(f, df, rs_grid)
6th order finite difference derivative of a function on realspace grid
Definition rs_methods.F:198
rs_methods::derive_fdm_cd5
subroutine, public derive_fdm_cd5(f, df, rs_grid)
4th order finite difference derivative of a function on realspace grid
Definition rs_methods.F:129
rs_methods::derive_fdm_cd3
subroutine, public derive_fdm_cd3(f, df, rs_grid)
2nd order finite difference derivative of a function on realspace grid
Definition rs_methods.F:60
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_control_types::sccs_control_type
Definition cp_control_types.F:583
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:123
pw_pool_types::pw_pool_p_type
to create arrays of pools
Definition pw_pool_types.F:135
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:120
realspace_grid_types::realspace_grid_desc_type
Definition realspace_grid_types.F:90
realspace_grid_types::realspace_grid_input_type
Definition realspace_grid_types.F:80
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136