1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 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 CPABORT("No stress tensor is implemented for SCCS models yet") 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="(/,T2,A,/,/,T2,A)") & 341 "The dielectric function is written in cube file format to the file:", TRIM(filename) 342 END IF 343 CALL cp_pw_to_cube(eps_elec, cube_unit, TRIM(filename), particles=particles, & 344 stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), & 345 mpi_io=mpi_io) 346 CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io) 347 END IF 348 349 ! Calculate the (quantum) volume and surface of the solute cavity 350 cavity_surface = 0.0_dp 351 cavity_volume = 0.0_dp 352 353 IF (should_output .AND. (ABS(eps0 - 1.0_dp) > tol)) THEN 354 355 ! Initialise the switching function theta 356 theta => work_r3d(4)%pw 357 CALL pw_zero(theta) 358 359 ! Calculate the (quantum) volume of the solute cavity 360 f = 1.0_dp/(eps0 - 1.0_dp) 361!$OMP PARALLEL DO DEFAULT(NONE) & 362!$OMP PRIVATE(i,j,k) & 363!$OMP SHARED(eps0,eps_elec,f,lb,theta,ub) 364 DO k = lb(3), ub(3) 365 DO j = lb(2), ub(2) 366 DO i = lb(1), ub(1) 367 theta%cr3d(i, j, k) = f*(eps0 - eps_elec%cr3d(i, j, k)) 368 END DO 369 END DO 370 END DO 371!$OMP END PARALLEL DO 372 cavity_volume = accurate_sum(theta%cr3d)*dvol 373 CALL mp_sum(cavity_volume, para_env%group) 374 375 ! Calculate the derivative of the electronic density in r-space 376 ! TODO: Could be retrieved from the QS environment 377 DO i = 1, 3 378 NULLIFY (drho_elec(i)%pw) 379 CALL pw_pool_create_pw(auxbas_pw_pool, & 380 drho_elec(i)%pw, & 381 use_data=REALDATA3D, & 382 in_space=REALSPACE) 383 END DO 384 CALL derive(rho_elec, drho_elec, sccs_derivative_fft, pw_env, input, para_env) 385 386 ! Calculate the norm of the gradient of the electronic density in r-space 387 norm_drho_elec => work_r3d(5)%pw 388!$OMP PARALLEL DO DEFAULT(NONE) & 389!$OMP PRIVATE(i,j,k) & 390!$OMP SHARED(drho_elec,lb,norm_drho_elec,ub) 391 DO k = lb(3), ub(3) 392 DO j = lb(2), ub(2) 393 DO i = lb(1), ub(1) 394 norm_drho_elec%cr3d(i, j, k) = SQRT(drho_elec(1)%pw%cr3d(i, j, k)* & 395 drho_elec(1)%pw%cr3d(i, j, k) + & 396 drho_elec(2)%pw%cr3d(i, j, k)* & 397 drho_elec(2)%pw%cr3d(i, j, k) + & 398 drho_elec(3)%pw%cr3d(i, j, k)* & 399 drho_elec(3)%pw%cr3d(i, j, k)) 400 END DO 401 END DO 402 END DO 403!$OMP END PARALLEL DO 404 405 ! Optional output of the norm of the density gradient in cube file format 406 filename = "DENSITY_GRADIENT" 407 cube_path = TRIM(print_path)//"%"//TRIM(filename) 408 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), & 409 cp_p_file)) THEN 410 append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND") 411 my_pos_cube = "REWIND" 412 IF (append_cube) my_pos_cube = "APPEND" 413 mpi_io = .TRUE. 414 cube_unit = cp_print_key_unit_nr(logger, input, TRIM(cube_path), & 415 extension=".cube", middle_name=TRIM(filename), & 416 file_position=my_pos_cube, log_filename=.FALSE., & 417 mpi_io=mpi_io, fout=mpi_filename) 418 IF (output_unit > 0) THEN 419 IF (.NOT. mpi_io) THEN 420 INQUIRE (UNIT=cube_unit, NAME=filename) 421 ELSE 422 filename = mpi_filename 423 END IF 424 WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & 425 "The norm of the density gradient is written in cube file format to the file:", TRIM(filename) 426 END IF 427 CALL cp_pw_to_cube(norm_drho_elec, cube_unit, TRIM(filename), particles=particles, & 428 stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), & 429 mpi_io=mpi_io) 430 CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io) 431 END IF 432 433 ! Calculate the (quantum) surface of the solute cavity 434 SELECT CASE (sccs_control%method_id) 435 CASE (sccs_andreussi) 436 CALL surface_andreussi(rho_elec, norm_drho_elec, theta, eps0, & 437 sccs_control%rho_max, sccs_control%rho_min, & 438 sccs_control%delta_rho) 439 CASE (sccs_fattebert_gygi) 440 CALL surface_fattebert_gygi(rho_elec, norm_drho_elec, theta, eps0, & 441 sccs_control%beta, sccs_control%rho_zero, & 442 sccs_control%delta_rho) 443 CASE DEFAULT 444 CPABORT("Invalid method specified for SCCS model") 445 END SELECT 446 cavity_surface = accurate_sum(theta%cr3d)*dvol 447 CALL mp_sum(cavity_surface, para_env%group) 448 449 IF (sccs_control%gamma_solvent > 0.0_dp) THEN 450 451 CALL cp_warn(__LOCATION__, & 452 "Experimental feature requested: The contribution from the "// & 453 "cavitation energy is added as an extra potential term to the "// & 454 "Kohn-Sham potential") 455 456 ! Calculate the second derivatives of the electronic density in r-space 457 DO di = 1, 3 458 DO dj = 1, 3 459 NULLIFY (d2rho_elec(dj, di)%pw) 460 CALL pw_pool_create_pw(auxbas_pw_pool, & 461 d2rho_elec(dj, di)%pw, & 462 use_data=REALDATA3D, & 463 in_space=REALSPACE) 464 END DO 465 CALL derive(drho_elec(di)%pw, d2rho_elec(:, di), sccs_derivative_fft, pw_env, & 466 input, para_env) 467 END DO 468 469 ! Calculate the contribution of the cavitation energy to the Kohn-Sham potential 470!$OMP PARALLEL DO DEFAULT(NONE) & 471!$OMP PRIVATE(di,dj,i,j,k,norm_drho2) & 472!$OMP SHARED(d2rho_elec,drho_elec,lb,norm_drho_elec,sccs_control) & 473!$OMP SHARED(theta,ub,v_sccs) 474 DO k = lb(3), ub(3) 475 DO j = lb(2), ub(2) 476 DO i = lb(1), ub(1) 477 norm_drho2 = norm_drho_elec%cr3d(i, j, k)*norm_drho_elec%cr3d(i, j, k) 478 DO di = 1, 3 479 DO dj = 1, 3 480 v_sccs%cr3d(i, j, k) = sccs_control%gamma_solvent*theta%cr3d(i, j, k)* & 481 (drho_elec(di)%pw%cr3d(i, j, k)* & 482 drho_elec(dj)%pw%cr3d(i, j, k)* & 483 d2rho_elec(di, dj)%pw%cr3d(i, j, k)/norm_drho2 - & 484 d2rho_elec(di, di)%pw%cr3d(i, j, k))/norm_drho2 485 END DO 486 END DO 487 END DO 488 END DO 489 END DO 490!$OMP END PARALLEL DO 491 CALL pw_scale(v_sccs, dvol) 492 DO di = 1, 3 493 DO dj = 1, 3 494 CALL pw_pool_give_back_pw(auxbas_pw_pool, d2rho_elec(dj, di)%pw) 495 END DO 496 END DO 497 END IF 498 499 ! Release storage 500 NULLIFY (theta) 501 DO i = 1, 3 502 CALL pw_pool_give_back_pw(auxbas_pw_pool, drho_elec(i)%pw) 503 END DO 504 505 END IF 506 507 ! Retrieve the total charge density (core + elec) of the solute in r-space 508 rho_solute => work_r3d(4)%pw 509 CALL pw_zero(rho_solute) 510 CALL pw_transfer(rho_tot_gspace, rho_solute) 511 tot_rho_solute = accurate_sum(rho_solute%cr3d)*dvol 512 CALL mp_sum(tot_rho_solute, para_env%group) 513 514 ! Reassign work storage to rho_tot_zero, because rho_elec is no longer needed 515 rho_tot_zero => rho_elec 516 NULLIFY (rho_elec) 517 518 ! Build the initial (rho_iter = 0) total charge density (solute plus polarisation) in r-space 519 ! eps_elec <- ln(eps_elec) 520!$OMP PARALLEL DO DEFAULT(NONE) & 521!$OMP PRIVATE(i,j,k) & 522!$OMP SHARED(eps_elec,lb,message,output_unit,para_env,ub) & 523!$OMP SHARED(rho_solute,rho_tot_zero) 524 DO k = lb(3), ub(3) 525 DO j = lb(2), ub(2) 526 DO i = lb(1), ub(1) 527 IF (eps_elec%cr3d(i, j, k) < 1.0E-12_dp) THEN 528 WRITE (UNIT=message, FMT="(A,ES12.3,A,3(I0,A))") & 529 "SCCS| Invalid dielectric function value ", eps_elec%cr3d(i, j, k), & 530 " encountered at grid point (", i, ",", j, ",", k, ")" 531 CPABORT(message) 532 END IF 533 rho_tot_zero%cr3d(i, j, k) = rho_solute%cr3d(i, j, k)/eps_elec%cr3d(i, j, k) 534 eps_elec%cr3d(i, j, k) = LOG(eps_elec%cr3d(i, j, k)) 535 END DO 536 END DO 537 END DO 538!$OMP END PARALLEL DO 539 540 ! Reassign pointers due to new content 541 ln_eps_elec => eps_elec 542 NULLIFY (eps_elec) 543 544 ! Build the derivative of ln_eps_elec 545 DO i = 1, 3 546 NULLIFY (dln_eps_elec(i)%pw) 547 CALL pw_pool_create_pw(auxbas_pw_pool, & 548 dln_eps_elec(i)%pw, & 549 use_data=REALDATA3D, & 550 in_space=REALSPACE) 551 CALL pw_zero(dln_eps_elec(i)%pw) 552 END DO 553 CALL derive(ln_eps_elec, dln_eps_elec, sccs_control%derivative_method, pw_env, input, para_env) 554 555 ! Print header for the SCCS cycle 556 IF (should_output .AND. (output_unit > 0)) THEN 557 IF (print_level > low_print_level) THEN 558 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 559 "SCCS| Total electronic charge density ", -tot_rho_elec, & 560 "SCCS| Total charge density (solute) ", -tot_rho_solute 561 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.3)") & 562 "SCCS| Surface of the solute cavity [bohr^2]", cavity_surface, & 563 "SCCS| [angstrom^2]", & 564 cp_unit_from_cp2k(cavity_surface, "angstrom^2"), & 565 "SCCS| Volume of the solute cavity [bohr^3]", cavity_volume, & 566 "SCCS| [angstrom^3]", & 567 cp_unit_from_cp2k(cavity_volume, "angstrom^3"), & 568 "SCCS| Volume of the cell [bohr^3]", cell_volume, & 569 "SCCS| [angstrom^3]", & 570 cp_unit_from_cp2k(cell_volume, "angstrom^3") 571 WRITE (UNIT=output_unit, FMT="(T3,A)") & 572 "SCCS|", & 573 "SCCS| Step Average residual Maximum residual" 574 END IF 575 END IF 576 577 ! Get storage for the derivative of the total potential (field) in r-space 578 DO i = 1, 3 579 NULLIFY (dphi_tot(i)%pw) 580 CALL pw_pool_create_pw(auxbas_pw_pool, & 581 dphi_tot(i)%pw, & 582 use_data=REALDATA3D, & 583 in_space=REALSPACE) 584 END DO 585 586 ! Reassign work storage to rho_tot, because ln_eps_elec is no longer needed 587 rho_tot => ln_eps_elec 588 NULLIFY (ln_eps_elec) 589 590 ! Initialise the total electronic density in r-space rho_tot with rho_tot_zero + rho_iter_zero 591 CALL pw_copy(rho_tot_zero, rho_tot) 592 CALL pw_axpy(rho_iter_old, rho_tot) 593 594 ! Main SCCS iteration loop 595 iter = 0 596 597 iter_loop: DO 598 599 ! Increment iteration counter 600 iter = iter + 1 601 602 ! Check if the requested maximum number of SCCS iterations is reached 603 IF (iter > sccs_control%max_iter) THEN 604 IF (output_unit > 0) THEN 605 WRITE (UNIT=output_unit, FMT="(T3,A,/,T3,A,I0,A)") & 606 "SCCS| Maximum number of SCCS iterations reached", & 607 "SCCS| Iteration cycle did not converge in ", sccs_control%max_iter, " steps" 608 END IF 609 EXIT iter_loop 610 END IF 611 612 ! Calculate derivative of the current total potential in r-space 613 CALL pw_poisson_solve(poisson_env=poisson_env, & 614 density=rho_tot, & 615 dvhartree=dphi_tot) 616 617 ! Update total charge density (solute plus polarisation) in r-space 618 ! based on the iterated polarisation charge density 619 f = 0.5_dp/twopi 620 rho_delta_avg = 0.0_dp 621 rho_delta_max = 0.0_dp 622!$OMP PARALLEL DO DEFAULT(NONE) & 623!$OMP PRIVATE(i,j,k,rho_delta,rho_iter_new) & 624!$OMP SHARED(dln_eps_elec,dphi_tot,f,lb,rho_iter_old,ub) & 625!$OMP SHARED(rho_tot,rho_tot_zero,sccs_control) & 626!$OMP REDUCTION(+:rho_delta_avg) & 627!$OMP REDUCTION(MAX:rho_delta_max) 628 DO k = lb(3), ub(3) 629 DO j = lb(2), ub(2) 630 DO i = lb(1), ub(1) 631 rho_iter_new = (dln_eps_elec(1)%pw%cr3d(i, j, k)*dphi_tot(1)%pw%cr3d(i, j, k) + & 632 dln_eps_elec(2)%pw%cr3d(i, j, k)*dphi_tot(2)%pw%cr3d(i, j, k) + & 633 dln_eps_elec(3)%pw%cr3d(i, j, k)*dphi_tot(3)%pw%cr3d(i, j, k))*f 634 rho_iter_new = rho_iter_old%cr3d(i, j, k) + & 635 sccs_control%mixing*(rho_iter_new - rho_iter_old%cr3d(i, j, k)) 636 rho_delta = ABS(rho_iter_new - rho_iter_old%cr3d(i, j, k)) 637 rho_delta_max = MAX(rho_delta, rho_delta_max) 638 rho_delta_avg = rho_delta_avg + rho_delta 639 rho_tot%cr3d(i, j, k) = rho_tot_zero%cr3d(i, j, k) + rho_iter_new 640 rho_iter_old%cr3d(i, j, k) = rho_iter_new 641 END DO 642 END DO 643 END DO 644!$OMP END PARALLEL DO 645 646 CALL mp_sum(rho_delta_avg, para_env%group) 647 rho_delta_avg = rho_delta_avg/REAL(ngpts, KIND=dp) 648 CALL mp_max(rho_delta_max, para_env%group) 649 650 IF (should_output .AND. (output_unit > 0)) THEN 651 IF (print_level > low_print_level) THEN 652 IF ((ABS(rho_delta_avg) < 1.0E-8_dp) .OR. & 653 (ABS(rho_delta_avg) >= 1.0E5_dp)) THEN 654 WRITE (UNIT=output_unit, FMT="(T3,A,I6,4X,ES16.4,4X,ES16.4)") & 655 "SCCS| ", iter, rho_delta_avg, rho_delta_max 656 ELSE 657 WRITE (UNIT=output_unit, FMT="(T3,A,I6,4X,F16.8,4X,F16.8)") & 658 "SCCS| ", iter, rho_delta_avg, rho_delta_max 659 END IF 660 END IF 661 END IF 662 663 ! Check if the SCCS iteration cycle is converged to the requested tolerance 664 IF (rho_delta_max <= sccs_control%eps_sccs) THEN 665 IF (should_output .AND. (output_unit > 0)) THEN 666 WRITE (UNIT=output_unit, FMT="(T3,A,I0,A)") & 667 "SCCS| Iteration cycle converged in ", iter, " steps" 668 END IF 669 EXIT iter_loop 670 END IF 671 672 END DO iter_loop 673 674 ! Calculate the total Hartree energy, potential, and its derivatives of 675 ! the solute and the implicit solvent 676 CALL pw_transfer(rho_tot, rho_tot_gspace) 677 IF (calculate_stress_tensor) THEN 678 CALL pw_poisson_solve(poisson_env=poisson_env, & 679 density=rho_tot_gspace, & 680 ehartree=energy%sccs_hartree, & 681 vhartree=v_hartree_gspace, & 682 dvhartree=dphi_tot, & 683 h_stress=h_stress) 684 ELSE 685 CALL pw_poisson_solve(poisson_env=poisson_env, & 686 density=rho_tot_gspace, & 687 ehartree=energy%sccs_hartree, & 688 vhartree=v_hartree_gspace, & 689 dvhartree=dphi_tot) 690 END IF 691 phi_pol => work_r3d(5)%pw 692 CALL pw_transfer(v_hartree_gspace, phi_pol) 693 694 ! Calculate the Hartree energy and potential of the solute only 695 phi_solute => rho_tot_zero 696 CALL pw_zero(phi_solute) 697 NULLIFY (rho_tot_zero) 698 CALL pw_poisson_solve(poisson_env=poisson_env, & 699 density=rho_solute, & 700 ehartree=energy%hartree, & 701 vhartree=phi_solute) 702 703 ! Calculate the polarisation potential 704 ! phi_pol = phi_tot - phi_solute 705 CALL pw_axpy(phi_solute, phi_pol, alpha=-1.0_dp) 706 707 ! Calculate the polarisation charge 708 ! rho_pol = rho_tot - rho_solute 709 CALL pw_axpy(rho_solute, rho_tot, alpha=-1.0_dp) 710 polarisation_charge = accurate_sum(rho_tot%cr3d)*dvol 711 CALL mp_sum(polarisation_charge, para_env%group) 712 713 filename = "POLARISATION_POTENTIAL" 714 cube_path = TRIM(print_path)//"%"//TRIM(filename) 715 IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), & 716 cp_p_file)) THEN 717 append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND") 718 my_pos_cube = "REWIND" 719 IF (append_cube) my_pos_cube = "APPEND" 720 mpi_io = .TRUE. 721 cube_unit = cp_print_key_unit_nr(logger, input, TRIM(cube_path), & 722 extension=".cube", middle_name=TRIM(filename), & 723 file_position=my_pos_cube, log_filename=.FALSE., & 724 mpi_io=mpi_io, fout=mpi_filename) 725 IF (output_unit > 0) THEN 726 IF (.NOT. mpi_io) THEN 727 INQUIRE (UNIT=cube_unit, NAME=filename) 728 ELSE 729 filename = mpi_filename 730 END IF 731 WRITE (UNIT=output_unit, FMT="(/,T2,A,/,/,T2,A)") & 732 "The SCCS polarisation potential is written in cube file format to the file:", TRIM(filename) 733 END IF 734 CALL cp_pw_to_cube(phi_pol, cube_unit, TRIM(filename), particles=particles, & 735 stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), mpi_io=mpi_io) 736 CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io) 737 END IF 738 739 ! Calculate SCCS polarisation energy 740 energy%sccs_pol = 0.5_dp*pw_integral_ab(rho_solute, phi_pol) 741 742 ! Calculate the Makov-Payne energy correction for charged systems with PBC 743 ! Madelung energy of a simple cubic lattice of point charges immersed in neutralising jellium 744 energy%sccs_mpc = -0.5_dp*2.8373_dp*tot_rho_solute**2/MINVAL(abc(1:3)) 745 746 ! Calculate additional solvation terms 747 energy%sccs_cav = sccs_control%gamma_solvent*cavity_surface 748 energy%sccs_dis = sccs_control%beta_solvent*cavity_volume 749 energy%sccs_rep = sccs_control%alpha_solvent*cavity_surface 750 751 IF (should_output .AND. (output_unit > 0)) THEN 752 WRITE (UNIT=output_unit, FMT="(T3,A)") & 753 "SCCS|" 754 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.12)") & 755 "SCCS| Polarisation charge", polarisation_charge 756 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.12)") & 757 "SCCS| Hartree energy of the solute only [Hartree]", energy%hartree, & 758 "SCCS| Hartree energy of solute and solvent [Hartree]", energy%sccs_hartree 759 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.12,/,T3,A,T61,F20.3)") & 760 "SCCS| Polarisation energy [Hartree]", energy%sccs_pol, & 761 "SCCS| [kcal/mol]", & 762 cp_unit_from_cp2k(energy%sccs_pol, "kcalmol"), & 763 "SCCS| Makov-Payne correction [Hartree]", energy%sccs_mpc, & 764 "SCCS| [kcal/mol]", & 765 cp_unit_from_cp2k(energy%sccs_mpc, "kcalmol"), & 766 "SCCS| Cavitation energy [Hartree]", energy%sccs_cav, & 767 "SCCS| [kcal/mol]", & 768 cp_unit_from_cp2k(energy%sccs_cav, "kcalmol"), & 769 "SCCS| Dispersion free energy [Hartree]", energy%sccs_dis, & 770 "SCCS| [kcal/mol]", & 771 cp_unit_from_cp2k(energy%sccs_dis, "kcalmol"), & 772 "SCCS| Repulsion free energy [Hartree]", energy%sccs_rep, & 773 "SCCS| [kcal/mol]", & 774 cp_unit_from_cp2k(energy%sccs_rep, "kcalmol") 775 END IF 776 777 ! Calculate SCCS contribution to the Kohn-Sham potential 778 f = -0.25_dp*dvol/twopi 779!$OMP PARALLEL DO DEFAULT(NONE) & 780!$OMP PRIVATE(dphi2,i,j,k) & 781!$OMP SHARED(calculate_stress_tensor,f,deps_elec,dphi_tot) & 782!$OMP SHARED(dvol,lb,ub,v_sccs) 783 DO k = lb(3), ub(3) 784 DO j = lb(2), ub(2) 785 DO i = lb(1), ub(1) 786 dphi2 = dphi_tot(1)%pw%cr3d(i, j, k)*dphi_tot(1)%pw%cr3d(i, j, k) + & 787 dphi_tot(2)%pw%cr3d(i, j, k)*dphi_tot(2)%pw%cr3d(i, j, k) + & 788 dphi_tot(3)%pw%cr3d(i, j, k)*dphi_tot(3)%pw%cr3d(i, j, k) 789 v_sccs%cr3d(i, j, k) = v_sccs%cr3d(i, j, k) + f*deps_elec%cr3d(i, j, k)*dphi2 790 END DO 791 END DO 792 END DO 793!$OMP END PARALLEL DO 794 795 ! Release work storage 796 NULLIFY (phi_solute) 797 NULLIFY (rho_tot) 798 NULLIFY (deps_elec) 799 NULLIFY (rho_solute) 800 DO i = 1, SIZE(work_r3d) 801 CALL pw_pool_give_back_pw(auxbas_pw_pool, work_r3d(i)%pw) 802 END DO 803 DO i = 1, 3 804 CALL pw_pool_give_back_pw(auxbas_pw_pool, dln_eps_elec(i)%pw) 805 CALL pw_pool_give_back_pw(auxbas_pw_pool, dphi_tot(i)%pw) 806 END DO 807 808 ! Release the SCCS printout environment 809 CALL cp_print_key_finished_output(output_unit, logger, input, TRIM(print_path), & 810 ignore_should_output=should_output) 811 812 CALL timestop(handle) 813 814 END SUBROUTINE sccs 815 816! ************************************************************************************************** 817!> \brief Calculate the smoothed dielectric function of Andreussi et al. 818!> \param rho_elec ... 819!> \param eps_elec ... 820!> \param deps_elec ... 821!> \param epsilon_solvent ... 822!> \param rho_max ... 823!> \param rho_min ... 824!> \par History: 825!> - Creation (16.10.2013,MK) 826!> - Finite difference of isosurfaces implemented (21.12.2013,MK) 827!> \author Matthias Krack (MK) 828!> \version 1.1 829! ************************************************************************************************** 830 SUBROUTINE andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, rho_max, & 831 rho_min) 832 833 TYPE(pw_type), POINTER :: rho_elec, eps_elec, deps_elec 834 REAL(KIND=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min 835 836 CHARACTER(LEN=*), PARAMETER :: routineN = 'andreussi', routineP = moduleN//':'//routineN 837 REAL(KIND=dp), PARAMETER :: tol = 1.0E-12_dp 838 839 INTEGER :: handle, i, j, k 840 INTEGER, DIMENSION(3) :: lb, ub 841 REAL(KIND=dp) :: diff, dq, dt, f, ln_rho_max, ln_rho_min, & 842 q, rho, t, x, y 843 844 CALL timeset(routineN, handle) 845 846 f = LOG(epsilon_solvent)/twopi 847 diff = rho_max - rho_min 848 IF (diff > tol) THEN 849 ln_rho_max = LOG(rho_max) 850 ln_rho_min = LOG(rho_min) 851 q = twopi/(ln_rho_max - ln_rho_min) 852 dq = -f*q 853 END IF 854 855 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3) 856 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3) 857 858 ! Calculate the dielectric function and its derivative 859!$OMP PARALLEL DO DEFAULT(NONE) & 860!$OMP PRIVATE(dt,i,j,k,rho,t,x,y) & 861!$OMP SHARED(deps_elec,diff,dq,eps_elec,epsilon_solvent,f,lb,ub) & 862!$OMP SHARED(ln_rho_max,rho_elec,q,rho_max,rho_min) 863 DO k = lb(3), ub(3) 864 DO j = lb(2), ub(2) 865 DO i = lb(1), ub(1) 866 rho = rho_elec%cr3d(i, j, k) 867 IF (rho < rho_min) THEN 868 eps_elec%cr3d(i, j, k) = epsilon_solvent 869 deps_elec%cr3d(i, j, k) = 0.0_dp 870 ELSE IF (rho <= rho_max) THEN 871 IF (diff > tol) THEN 872 x = LOG(rho) 873 y = q*(ln_rho_max - x) 874 t = f*(y - SIN(y)) 875 eps_elec%cr3d(i, j, k) = EXP(t) 876 dt = dq*(1.0_dp - COS(y)) 877 deps_elec%cr3d(i, j, k) = eps_elec%cr3d(i, j, k)*dt/rho 878 ELSE 879 eps_elec%cr3d(i, j, k) = 1.0_dp 880 deps_elec%cr3d(i, j, k) = 0.0_dp 881 END IF 882 ELSE 883 eps_elec%cr3d(i, j, k) = 1.0_dp 884 deps_elec%cr3d(i, j, k) = 0.0_dp 885 END IF 886 END DO 887 END DO 888 END DO 889!$OMP END PARALLEL DO 890 891 CALL timestop(handle) 892 893 END SUBROUTINE andreussi 894 895! ************************************************************************************************** 896!> \brief Calculate the smoothed dielectric function of Fattebert and Gygi 897!> \param rho_elec ... 898!> \param eps_elec ... 899!> \param deps_elec ... 900!> \param epsilon_solvent ... 901!> \param beta ... 902!> \param rho_zero ... 903!> \par History: 904!> - Creation (15.10.2013,MK) 905!> \author Matthias Krack (MK) 906!> \version 1.0 907! ************************************************************************************************** 908 SUBROUTINE fattebert_gygi(rho_elec, eps_elec, deps_elec, epsilon_solvent, beta, & 909 rho_zero) 910 911 TYPE(pw_type), POINTER :: rho_elec, eps_elec, deps_elec 912 REAL(KIND=dp), INTENT(IN) :: epsilon_solvent, beta, rho_zero 913 914 CHARACTER(LEN=*), PARAMETER :: routineN = 'fattebert_gygi', routineP = moduleN//':'//routineN 915 REAL(KIND=dp), PARAMETER :: tol = 1.0E-12_dp 916 917 INTEGER :: handle, i, j, k 918 INTEGER, DIMENSION(3) :: lb, ub 919 REAL(KIND=dp) :: df, f, p, q, rho, s, t, twobeta 920 921 CALL timeset(routineN, handle) 922 923 df = (1.0_dp - epsilon_solvent)/rho_zero 924 f = 0.5_dp*(epsilon_solvent - 1.0_dp) 925 q = 1.0_dp/rho_zero 926 twobeta = 2.0_dp*beta 927 928 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3) 929 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3) 930 931 ! Calculate the smoothed dielectric function and its derivative 932!$OMP PARALLEL DO DEFAULT(NONE) & 933!$OMP PRIVATE(i,j,k,p,rho,s,t) & 934!$OMP SHARED(df,deps_elec,eps_elec,epsilon_solvent,f,lb,ub) & 935!$OMP SHARED(q,rho_elec,twobeta) 936 DO k = lb(3), ub(3) 937 DO j = lb(2), ub(2) 938 DO i = lb(1), ub(1) 939 rho = rho_elec%cr3d(i, j, k) 940 IF (rho < tol) THEN 941 eps_elec%cr3d(i, j, k) = epsilon_solvent 942 deps_elec%cr3d(i, j, k) = 0.0_dp 943 ELSE 944 s = rho*q 945 p = s**twobeta 946 t = 1.0_dp/(1.0_dp + p) 947 eps_elec%cr3d(i, j, k) = 1.0_dp + f*(1.0_dp + (1.0_dp - p)*t) 948 deps_elec%cr3d(i, j, k) = df*twobeta*t*t*p/s 949 END IF 950 END DO 951 END DO 952 END DO 953!$OMP END PARALLEL DO 954 955 CALL timestop(handle) 956 957 END SUBROUTINE fattebert_gygi 958 959! ************************************************************************************************** 960!> \brief Build the numerical derivative of a function on realspace grid 961!> \param f ... 962!> \param df ... 963!> \param method ... 964!> \param pw_env ... 965!> \param input ... 966!> \param para_env ... 967!> \par History: 968!> - Creation (15.11.2013,MK) 969!> \author Matthias Krack (MK) 970!> \version 1.0 971! ************************************************************************************************** 972 SUBROUTINE derive(f, df, method, pw_env, input, para_env) 973 974 TYPE(pw_type), POINTER :: f 975 TYPE(pw_p_type), DIMENSION(3), INTENT(INOUT) :: df 976 INTEGER, INTENT(IN) :: method 977 TYPE(pw_env_type), POINTER :: pw_env 978 TYPE(section_vals_type), POINTER :: input 979 TYPE(cp_para_env_type), POINTER :: para_env 980 981 CHARACTER(LEN=*), PARAMETER :: routineN = 'derive', routineP = moduleN//':'//routineN 982 983 INTEGER :: border_points, handle, i 984 INTEGER, DIMENSION(3) :: lb, n, ub 985 TYPE(pw_p_type), DIMENSION(2) :: work_g1d 986 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 987 TYPE(realspace_grid_desc_type), POINTER :: rs_desc 988 TYPE(realspace_grid_input_type) :: input_settings 989 TYPE(realspace_grid_type), POINTER :: rs_grid 990 TYPE(section_vals_type), POINTER :: rs_grid_section 991 992 CALL timeset(routineN, handle) 993 994 CPASSERT(ASSOCIATED(f)) 995 CPASSERT(ASSOCIATED(pw_env)) 996 CPASSERT(ASSOCIATED(para_env)) 997 998 ! Perform method specific setup 999 SELECT CASE (method) 1000 CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7) 1001 NULLIFY (rs_desc) 1002 NULLIFY (rs_grid) 1003 rs_grid_section => section_vals_get_subs_vals(input, "DFT%MGRID%RS_GRID") 1004 SELECT CASE (method) 1005 CASE (sccs_derivative_cd3) 1006 border_points = 1 1007 CASE (sccs_derivative_cd5) 1008 border_points = 2 1009 CASE (sccs_derivative_cd7) 1010 border_points = 3 1011 END SELECT 1012 CALL init_input_type(input_settings, 2*border_points + 1, rs_grid_section, & 1013 1, (/-1, -1, -1/)) 1014 CALL rs_grid_create_descriptor(rs_desc, f%pw_grid, input_settings, & 1015 border_points=border_points) 1016 CALL rs_grid_create(rs_grid, rs_desc) 1017!MK CALL rs_grid_print(rs_grid,6) 1018 CASE (sccs_derivative_fft) 1019 lb(1:3) = f%pw_grid%bounds_local(1, 1:3) 1020 ub(1:3) = f%pw_grid%bounds_local(2, 1:3) 1021 NULLIFY (auxbas_pw_pool) 1022 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 1023 ! Get work storage for the 1d grids in g-space (derivative calculation) 1024 DO i = 1, SIZE(work_g1d) 1025 NULLIFY (work_g1d(i)%pw) 1026 CALL pw_pool_create_pw(auxbas_pw_pool, & 1027 work_g1d(i)%pw, & 1028 use_data=COMPLEXDATA1D, & 1029 in_space=RECIPROCALSPACE) 1030 END DO 1031 END SELECT 1032 1033 ! Calculate the derivatives 1034 SELECT CASE (method) 1035 CASE (sccs_derivative_cd3) 1036 CALL derive_fdm_cd3(f, df, rs_grid) 1037 CASE (sccs_derivative_cd5) 1038 CALL derive_fdm_cd5(f, df, rs_grid) 1039 CASE (sccs_derivative_cd7) 1040 CALL derive_fdm_cd7(f, df, rs_grid) 1041 CASE (sccs_derivative_fft) 1042 ! FFT 1043 CALL pw_transfer(f, work_g1d(1)%pw) 1044 DO i = 1, 3 1045 n(:) = 0 1046 n(i) = 1 1047 CALL pw_copy(work_g1d(1)%pw, work_g1d(2)%pw) 1048 CALL pw_derive(work_g1d(2)%pw, n(:)) 1049 CALL pw_transfer(work_g1d(2)%pw, df(i)%pw) 1050 END DO 1051 CASE DEFAULT 1052 CPABORT("Invalid derivative method for SCCS specified") 1053 END SELECT 1054 1055 ! Perform method specific cleanup 1056 SELECT CASE (method) 1057 CASE (sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7) 1058 CALL rs_grid_release(rs_grid) 1059 CALL rs_grid_release_descriptor(rs_desc) 1060 CASE (sccs_derivative_fft) 1061 DO i = 1, SIZE(work_g1d) 1062 CALL pw_pool_give_back_pw(auxbas_pw_pool, work_g1d(i)%pw) 1063 END DO 1064 END SELECT 1065 1066 CALL timestop(handle) 1067 1068 END SUBROUTINE derive 1069 1070! ************************************************************************************************** 1071!> \brief Calculate the finite difference between two isosurfaces of the 1072!> electronic density. The smoothed dielectric function of 1073!> Andreussi et al. is used as switching function eventually 1074!> defining the quantum volume and surface of the cavity. 1075!> \param rho_elec ... 1076!> \param norm_drho_elec ... 1077!> \param dtheta ... 1078!> \param epsilon_solvent ... 1079!> \param rho_max ... 1080!> \param rho_min ... 1081!> \param delta_rho ... 1082!> \par History: 1083!> - Creation (21.12.2013,MK) 1084!> \author Matthias Krack (MK) 1085!> \version 1.0 1086! ************************************************************************************************** 1087 SUBROUTINE surface_andreussi(rho_elec, norm_drho_elec, dtheta, & 1088 epsilon_solvent, rho_max, rho_min, delta_rho) 1089 1090 TYPE(pw_type), POINTER :: rho_elec, norm_drho_elec, dtheta 1091 REAL(KIND=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min, & 1092 delta_rho 1093 1094 CHARACTER(LEN=*), PARAMETER :: routineN = 'surface_andreussi', & 1095 routineP = moduleN//':'//routineN 1096 REAL(KIND=dp), PARAMETER :: tol = 1.0E-12_dp 1097 1098 INTEGER :: handle, i, j, k, l 1099 INTEGER, DIMENSION(3) :: lb, ub 1100 REAL(KIND=dp) :: diff, e, eps_elec, f, ln_rho_max, & 1101 ln_rho_min, q, rho, t, x, y 1102 REAL(KIND=dp), DIMENSION(2) :: theta 1103 1104 CALL timeset(routineN, handle) 1105 1106 e = epsilon_solvent - 1.0_dp 1107 f = LOG(epsilon_solvent)/twopi 1108 diff = rho_max - rho_min 1109 IF (diff > tol) THEN 1110 ln_rho_max = LOG(rho_max) 1111 ln_rho_min = LOG(rho_min) 1112 q = twopi/(ln_rho_max - ln_rho_min) 1113 END IF 1114 1115 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3) 1116 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3) 1117 1118 ! Calculate finite difference between two isosurfaces 1119!$OMP PARALLEL DO DEFAULT(NONE) & 1120!$OMP PRIVATE(eps_elec,i,j,k,l,rho,t,theta,x,y) & 1121!$OMP SHARED(delta_rho,diff,dtheta,e,epsilon_solvent,f,lb) & 1122!$OMP SHARED(ln_rho_max,norm_drho_elec,rho_elec,q,rho_max,rho_min,ub) 1123 DO k = lb(3), ub(3) 1124 DO j = lb(2), ub(2) 1125 DO i = lb(1), ub(1) 1126 DO l = 1, 2 1127 rho = rho_elec%cr3d(i, j, k) + (REAL(l, KIND=dp) - 1.5_dp)*delta_rho 1128 IF (rho < rho_min) THEN 1129 eps_elec = epsilon_solvent 1130 ELSE IF (rho <= rho_max) THEN 1131 IF (diff > tol) THEN 1132 x = LOG(rho) 1133 y = q*(ln_rho_max - x) 1134 t = f*(y - SIN(y)) 1135 eps_elec = EXP(t) 1136 ELSE 1137 eps_elec = 1.0_dp 1138 END IF 1139 ELSE 1140 eps_elec = 1.0_dp 1141 END IF 1142 theta(l) = (epsilon_solvent - eps_elec)/e 1143 END DO 1144 dtheta%cr3d(i, j, k) = (theta(2) - theta(1))*norm_drho_elec%cr3d(i, j, k)/delta_rho 1145 END DO 1146 END DO 1147 END DO 1148!$OMP END PARALLEL DO 1149 1150 CALL timestop(handle) 1151 1152 END SUBROUTINE surface_andreussi 1153 1154! ************************************************************************************************** 1155!> \brief Calculate the finite difference between two isosurfaces of the 1156!> the electronic density. The smoothed dielectric function of 1157!> Fattebert and Gygi is used as switching function eventually 1158!> defining the quantum volume and surface of the cavity. 1159!> \param rho_elec ... 1160!> \param norm_drho_elec ... 1161!> \param dtheta ... 1162!> \param epsilon_solvent ... 1163!> \param beta ... 1164!> \param rho_zero ... 1165!> \param delta_rho ... 1166!> \par History: 1167!> - Creation (21.12.2013,MK) 1168!> \author Matthias Krack (MK) 1169!> \version 1.0 1170! ************************************************************************************************** 1171 SUBROUTINE surface_fattebert_gygi(rho_elec, norm_drho_elec, dtheta, & 1172 epsilon_solvent, beta, rho_zero, delta_rho) 1173 1174 TYPE(pw_type), POINTER :: rho_elec, norm_drho_elec, dtheta 1175 REAL(KIND=dp), INTENT(IN) :: epsilon_solvent, beta, rho_zero, & 1176 delta_rho 1177 1178 CHARACTER(LEN=*), PARAMETER :: routineN = 'surface_fattebert_gygi', & 1179 routineP = moduleN//':'//routineN 1180 REAL(KIND=dp), PARAMETER :: tol = 1.0E-12_dp 1181 1182 INTEGER :: handle, i, j, k, l 1183 INTEGER, DIMENSION(3) :: lb, ub 1184 REAL(KIND=dp) :: e, eps_elec, f, p, q, rho, s, t, twobeta 1185 REAL(KIND=dp), DIMENSION(2) :: theta 1186 1187 CALL timeset(routineN, handle) 1188 1189 e = epsilon_solvent - 1.0_dp 1190 f = 0.5_dp*e 1191 q = 1.0_dp/rho_zero 1192 twobeta = 2.0_dp*beta 1193 1194 lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3) 1195 ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3) 1196 1197 ! Calculate finite difference between two isosurfaces 1198!$OMP PARALLEL DO DEFAULT(NONE) & 1199!$OMP PRIVATE(eps_elec,i,j,k,l,p,rho,s,t,theta) & 1200!$OMP SHARED(delta_rho,dtheta,e,epsilon_solvent,f,lb) & 1201!$OMP SHARED(norm_drho_elec,q,rho_elec,twobeta,ub) 1202 DO k = lb(3), ub(3) 1203 DO j = lb(2), ub(2) 1204 DO i = lb(1), ub(1) 1205 DO l = 1, 2 1206 rho = rho_elec%cr3d(i, j, k) + (REAL(l, KIND=dp) - 1.5_dp)*delta_rho 1207 IF (rho < tol) THEN 1208 eps_elec = epsilon_solvent 1209 ELSE 1210 s = rho*q 1211 p = s**twobeta 1212 t = 1.0_dp/(1.0_dp + p) 1213 eps_elec = 1.0_dp + f*(1.0_dp + (1.0_dp - p)*t) 1214 END IF 1215 theta(l) = (epsilon_solvent - eps_elec)/e 1216 END DO 1217 dtheta%cr3d(i, j, k) = (theta(2) - theta(1))*norm_drho_elec%cr3d(i, j, k)/delta_rho 1218 END DO 1219 END DO 1220 END DO 1221!$OMP END PARALLEL DO 1222 1223 CALL timestop(handle) 1224 1225 END SUBROUTINE surface_fattebert_gygi 1226 1227END MODULE qs_sccs 1228