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