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