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