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