diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 23e8f0c7a2..07462bc5e2 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -4285,6 +4285,15 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="IN_MEMORY", & + description="Precompute gradients due to constraint during"// & + " initial formation of constraint and store them in memory. Does"// & + " nothing if forces are not calculated.", & + usage="IN_MEMORY", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_cdft_control_section ! ************************************************************************************************** @@ -4356,6 +4365,40 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="PRINT_DENSITY", & + description="Logical to control printing of Hirshfeld densities to .cube file.", & + usage="PRINT_DENSITY TRUE", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ATOMS_MEMORY", & + description="Number of atomic gradients to store in memory.", & + usage="ATOMS_MEMORY", & + n_var=1, type_of_var=integer_t, & + default_i_val=80) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="USE_ATOMIC_CUTOFF", & + description="Logical to control use of ATOMIC_CUTOFF.", & + usage="USE_ATOMIC_CUTOFF TRUE", & + default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="EPS_CUTOFF", & + description="Numerical cutoff for calculation of weight function.", & + usage="EPS_CUTOFF {real} ", default_r_val=1.0e-12_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ATOMIC_CUTOFF", & + description="Numerical cutoff for calculation of Hirshfeld densities.", & + usage="ATOMIC_CUTOFF {real} ", default_r_val=1.0e-12_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_hirshfeld_constraint_section ! ************************************************************************************************** diff --git a/src/qs_cdft_methods.F b/src/qs_cdft_methods.F index 972ec0a6d0..0b158a817a 100644 --- a/src/qs_cdft_methods.F +++ b/src/qs_cdft_methods.F @@ -27,9 +27,7 @@ MODULE qs_cdft_methods USE cp_realspace_grid_cube, ONLY: cp_cube_to_pw USE grid_api, ONLY: GRID_FUNC_AB,& collocate_pgf_product - USE hirshfeld_types, ONLY: get_hirshfeld_info,& - hirshfeld_type,& - set_hirshfeld_info + USE hirshfeld_types, ONLY: hirshfeld_type USE input_constants, ONLY: cdft_alpha_constraint,& cdft_beta_constraint,& cdft_charge_constraint,& @@ -59,6 +57,7 @@ MODULE qs_cdft_methods hirshfeld_constraint_type USE qs_cdft_utils, ONLY: becke_constraint_init,& cdft_constraint_print,& + cdft_print_hirshfeld_density,& hfun_scale,& hirshfeld_constraint_init USE qs_energy_types, ONLY: qs_energy_type @@ -124,7 +123,7 @@ CONTAINS ! Integrate the Becke constraint CALL cdft_constraint_integrate(qs_env) ! Calculate forces - IF (calculate_forces) CALL becke_constraint_force(qs_env) + IF (calculate_forces) CALL cdft_constraint_force(qs_env) END IF CALL timestop(handle) @@ -579,193 +578,571 @@ CONTAINS END SUBROUTINE becke_constraint_low ! ************************************************************************************************** -!> \brief Calculates atomic forces due to a Becke constraint -!> \param qs_env the qs_env where to build the gradients -!> \par History -!> Created 01.2007 [fschiff] -!> Extended functionality 12/15-12/16 [Nico Holmberg] +!> \brief Driver routine for calculating a Hirshfeld constraint +!> \param qs_env ... +!> \param calc_pot ... +!> \param calculate_forces ... ! ************************************************************************************************** - SUBROUTINE becke_constraint_force(qs_env) + SUBROUTINE hirshfeld_constraint(qs_env, calc_pot, calculate_forces) TYPE(qs_environment_type), POINTER :: qs_env + LOGICAL :: calc_pot, calculate_forces - CHARACTER(len=*), PARAMETER :: routineN = 'becke_constraint_force' + CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint' - INTEGER :: handle, i, iatom, igroup, ikind, ind(3), & - ispin, j, k, natom, nvar - INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of - INTEGER, DIMENSION(2, 3) :: bo - REAL(kind=dp) :: dvol, eps_cavity, sign, th - REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: strength - REAL(KIND=dp), DIMENSION(:), POINTER :: cutoffs - TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set - TYPE(becke_constraint_type), POINTER :: becke_control + INTEGER :: handle + TYPE(cdft_control_type), POINTER :: cdft_control + TYPE(dft_control_type), POINTER :: dft_control + + CALL timeset(routineN, handle) + CALL get_qs_env(qs_env, dft_control=dft_control) + cdft_control => dft_control%qs_control%cdft_control + IF (dft_control%qs_control%cdft .AND. cdft_control%type == outer_scf_hirshfeld_constraint) THEN + IF (calc_pot) THEN + ! Initialize the Hirshfeld constraint environment + CALL hirshfeld_constraint_init(qs_env) + ! Calculate the Hirshfeld weight function and possibly the gradients + CALL hirshfeld_constraint_low(qs_env) + END IF + ! Integrate the Hirshfeld constraint + CALL cdft_constraint_integrate(qs_env) + ! Calculate forces + IF (calculate_forces) CALL cdft_constraint_force(qs_env) + END IF + CALL timestop(handle) + + END SUBROUTINE hirshfeld_constraint + + ! ************************************************************************************************** + !> \brief Calculates Hirshfeld constraints + ! ************************************************************************************************** +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param just_gradients ... +! ************************************************************************************************** + SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients) + TYPE(qs_environment_type), POINTER :: qs_env + LOGICAL, OPTIONAL :: just_gradients + + CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint_low', & + routineP = moduleN//':'//routineN + + INTEGER :: atom_a, atoms_memory, atoms_memory_num, handle, i, iatom, iex, igroup, ikind, & + ithread, j, k, natom, npme, nthread, num_atoms, num_species, numexp, subpatch_pattern + INTEGER, ALLOCATABLE, DIMENSION(:) :: num_species_small + INTEGER, DIMENSION(2, 3) :: bo + INTEGER, DIMENSION(3) :: lb_pw, lb_rs, npts, ub_pw, ub_rs + INTEGER, DIMENSION(:), POINTER :: atom_list, cores + LOGICAL :: my_just_gradients + LOGICAL, ALLOCATABLE, DIMENSION(:) :: compute_charge, is_constraint + REAL(kind=dp) :: alpha, coef, eps_rho_rspace, exp_eval, & + prefactor, radius + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coefficients + REAL(kind=dp), DIMENSION(3) :: dr_pw, dr_rs, origin, r2, r_pbc, ra + REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cdft_control_type), POINTER :: cdft_control - TYPE(cdft_group_type), DIMENSION(:), POINTER :: group TYPE(cell_type), POINTER :: cell TYPE(cp_para_env_type), POINTER :: para_env TYPE(dft_control_type), POINTER :: dft_control + TYPE(hirshfeld_constraint_type), POINTER :: hirshfeld_control + TYPE(hirshfeld_type), POINTER :: hirshfeld_env TYPE(particle_type), DIMENSION(:), POINTER :: particle_set - TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r - TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(pw_env_type), POINTER :: pw_env + TYPE(pw_p_type), DIMENSION(:), POINTER :: pw_single_dr, rho_r + TYPE(pw_pool_type), POINTER :: auxbas_pw_pool TYPE(qs_rho_type), POINTER :: rho + TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc + TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: rs_single, rs_single_charge, rs_single_dr + TYPE(realspace_grid_type), POINTER :: rs_rho_all, rs_rho_constr + + NULLIFY (atom_list, atomic_kind_set, dft_control, pw_single_dr, & + hirshfeld_env, particle_set, pw_env, auxbas_pw_pool, para_env, & + auxbas_rs_desc, cdft_control, pab, rs_rho_all, & + hirshfeld_control, cell, rho_r, rho, rs_single, rs_single_dr, & + rs_rho_constr, rs_single_charge) CALL timeset(routineN, handle) - NULLIFY (atomic_kind_set, cell, para_env, dft_control, particle_set, & - rho, rho_r, force, cutoffs, becke_control, group) - CALL get_qs_env(qs_env, & atomic_kind_set=atomic_kind_set, & - natom=natom, & particle_set=particle_set, & + natom=natom, & cell=cell, & rho=rho, & - force=force, & dft_control=dft_control, & - para_env=para_env) + para_env=para_env, & + pw_env=pw_env) CALL qs_rho_get(rho, rho_r=rho_r) - th = 1.0e-8_dp + num_atoms = natom + cdft_control => dft_control%qs_control%cdft_control - becke_control => cdft_control%becke_control - group => cdft_control%group - nvar = SIZE(cdft_control%target) - ALLOCATE (strength(nvar)) - strength(:) = cdft_control%strength(:) - cutoffs => becke_control%cutoffs - eps_cavity = becke_control%eps_cavity - ALLOCATE (atom_of_kind(natom)) - ALLOCATE (kind_of(natom)) + hirshfeld_control => cdft_control%hirshfeld_control + hirshfeld_env => hirshfeld_control%hirshfeld_env - CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & - atom_of_kind=atom_of_kind, & - kind_of=kind_of) - DO igroup = 1, SIZE(group) - ALLOCATE (group(igroup)%integrated(3, natom)) - group(igroup)%integrated = 0.0_dp - END DO - bo = group(1)%weight%pw%pw_grid%bounds_local - dvol = group(1)%weight%pw%pw_grid%dvol + ! Check if only gradient should be calculated, if gradients should be precomputed + my_just_gradients = .FALSE. + IF (PRESENT(just_gradients)) my_just_gradients = just_gradients + IF (my_just_gradients) THEN + cdft_control%in_memory = .TRUE. + cdft_control%atomic_charges = .FALSE. + hirshfeld_control%print_density = .FALSE. + END IF - IF (.NOT. becke_control%in_memory) & - ! Gradients were not precomputed => calculate them now - CALL becke_constraint_low(qs_env, just_gradients=.TRUE.) - ! Integrate gradients - IF (.NOT. ASSOCIATED(becke_control%cavity_mat)) THEN - ! No external control - DO k = bo(1, 1), bo(2, 1) - DO j = bo(1, 2), bo(2, 2) - DO i = bo(1, 3), bo(2, 3) - ! First check if this grid point should be skipped - IF (becke_control%cavity_confine) THEN - IF (becke_control%cavity%pw%cr3d(k, j, i) < eps_cavity) CYCLE + ALLOCATE (coefficients(natom)) + ALLOCATE (is_constraint(natom)) + + subpatch_pattern = 0 + eps_rho_rspace = dft_control%qs_control%eps_rho_rspace + radius = 100.0_dp + + ! For each CDFT group + DO igroup = 1, SIZE(cdft_control%group) + + IF (.NOT. ASSOCIATED(rs_rho_all)) THEN + CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, auxbas_rs_grid=rs_rho_all, & + auxbas_pw_pool=auxbas_pw_pool) + CALL rs_grid_retain(rs_rho_all) + CALL rs_grid_zero(rs_rho_all) + END IF + + ! Get size of rs and pw grids + npts = rho_r(1)%pw%pw_grid%npts + dr_pw = rho_r(1)%pw%pw_grid%dr + dr_rs(1) = rs_rho_all%desc%dh(1, 1) + dr_rs(2) = rs_rho_all%desc%dh(2, 2) + dr_rs(3) = rs_rho_all%desc%dh(3, 3) + origin(:) = (dr_pw(:)*npts(:))*0.5_dp + lb_pw(1:3) = rho_r(1)%pw%pw_grid%bounds_local(1, 1:3) + ub_pw(1:3) = rho_r(1)%pw%pw_grid%bounds_local(2, 1:3) + lb_rs(1) = LBOUND(rs_rho_all%r(:, :, :), 1) + lb_rs(2) = LBOUND(rs_rho_all%r(:, :, :), 2) + lb_rs(3) = LBOUND(rs_rho_all%r(:, :, :), 3) + ub_rs(1) = UBOUND(rs_rho_all%r(:, :, :), 1) + ub_rs(2) = UBOUND(rs_rho_all%r(:, :, :), 2) + ub_rs(3) = UBOUND(rs_rho_all%r(:, :, :), 3) + bo = cdft_control%group(igroup)%weight%pw%pw_grid%bounds_local + + ! Coefficients + coefficients(:) = 0.0_dp + is_constraint = .FALSE. + DO i = 1, SIZE(cdft_control%group(igroup)%atoms) + coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i) + is_constraint(cdft_control%group(igroup)%atoms(i)) = .TRUE. + END DO + + ! rs_rho_constr: Sum of isolated Gaussian densities over constraint atoms in this constraint group + CALL rs_grid_create(rs_rho_constr, auxbas_rs_desc) + CALL rs_grid_zero(rs_rho_constr) + + ! rs_single: Gaussian density over single atoms when required + IF (hirshfeld_control%print_density .AND. igroup == 1) THEN + ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic(cdft_control%natoms)) + ALLOCATE (rs_single(cdft_control%natoms)) + DO i = 1, cdft_control%natoms + NULLIFY (rs_single(i)%rs_grid) + CALL rs_grid_create(rs_single(i)%rs_grid, auxbas_rs_desc) + CALL rs_grid_zero(rs_single(i)%rs_grid) + END DO + END IF + + ! Setup pw + cdft_control%group(igroup)%weight%pw%cr3d = 0.0_dp + + CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%group(igroup)%hw_rho_total_constraint%pw, & + use_data=REALDATA3D, in_space=REALSPACE) + cdft_control%group(igroup)%hw_rho_total_constraint%pw%cr3d = 1.0_dp + + IF (igroup == 1) THEN + CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%group(igroup)%hw_rho_total%pw, & + use_data=REALDATA3D, in_space=REALSPACE) + cdft_control%group(igroup)%hw_rho_total%pw%cr3d = 1.0_dp + + IF (hirshfeld_control%print_density) THEN + DO iatom = 1, cdft_control%natoms + CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%group(igroup)%hw_rho_atomic(iatom)%pw, & + use_data=REALDATA3D, in_space=REALSPACE) + cdft_control%group(igroup)%hw_rho_atomic(iatom)%pw%cr3d = 1.0_dp + END DO + END IF + END IF + + IF (cdft_control%atomic_charges .AND. igroup == 1) THEN + ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge(cdft_control%natoms)) + ALLOCATE (rs_single_charge(cdft_control%natoms)) + ALLOCATE (compute_charge(natom)) + compute_charge = .FALSE. + + DO i = 1, cdft_control%natoms + NULLIFY (rs_single_charge(i)%rs_grid) + CALL rs_grid_create(rs_single_charge(i)%rs_grid, auxbas_rs_desc) + CALL rs_grid_zero(rs_single_charge(i)%rs_grid) + compute_charge(cdft_control%atoms(i)) = .TRUE. + END DO + + DO iatom = 1, cdft_control%natoms + CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%group(igroup)%hw_rho_atomic_charge(iatom)%pw, & + use_data=REALDATA3D, in_space=REALSPACE) + cdft_control%group(igroup)%hw_rho_atomic_charge(iatom)%pw%cr3d = 1.0_dp + END DO + END IF + + ALLOCATE (pab(1, 1)) + nthread = 1 + ithread = 0 + + DO ikind = 1, SIZE(atomic_kind_set) + numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp + IF (numexp <= 0) CYCLE + CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list) + ALLOCATE (cores(num_species)) + + DO iex = 1, numexp + alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex) + coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex) + npme = 0 + cores = 0 + DO iatom = 1, num_species + atom_a = atom_list(iatom) + ra(:) = pbc(particle_set(atom_a)%r, cell) + IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN + IF (MODULO(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN + npme = npme + 1 + cores(npme) = iatom + END IF + ELSE + npme = npme + 1 + cores(npme) = iatom END IF - ind = (/k, j, i/) - DO igroup = 1, SIZE(group) - DO iatom = 1, natom - DO ispin = 1, dft_control%nspins - SELECT CASE (group(igroup)%constraint_type) - CASE (cdft_charge_constraint) - sign = 1.0_dp - CASE (cdft_magnetization_constraint) - IF (ispin == 1) THEN - sign = 1.0_dp - ELSE - sign = -1.0_dp - END IF - CASE (cdft_alpha_constraint) - sign = 1.0_dp - IF (ispin == 2) CYCLE - CASE (cdft_beta_constraint) - sign = 1.0_dp - IF (ispin == 1) CYCLE - CASE DEFAULT - CPABORT("Unknown constraint type.") - END SELECT - group(igroup)%integrated(:, iatom) = group(igroup)%integrated(:, iatom) + sign* & - group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i)* & - rho_r(ispin)%pw%cr3d(k, j, i)*dvol - END DO - END DO - END DO - END DO - END DO - END DO - ELSE - DO k = LBOUND(becke_control%cavity_mat, 1), UBOUND(becke_control%cavity_mat, 1) - DO j = LBOUND(becke_control%cavity_mat, 2), UBOUND(becke_control%cavity_mat, 2) - DO i = LBOUND(becke_control%cavity_mat, 3), UBOUND(becke_control%cavity_mat, 3) - ! First check if this grid point should be skipped - IF (becke_control%cavity_mat(k, j, i) < eps_cavity) CYCLE - DO igroup = 1, SIZE(group) - DO iatom = 1, natom - DO ispin = 1, dft_control%nspins - SELECT CASE (group(igroup)%constraint_type) - CASE (cdft_charge_constraint) - sign = 1.0_dp - CASE (cdft_magnetization_constraint) - IF (ispin == 1) THEN - sign = 1.0_dp - ELSE - sign = -1.0_dp - END IF - CASE (cdft_alpha_constraint) - sign = 1.0_dp - IF (ispin == 2) CYCLE - CASE (cdft_beta_constraint) - sign = 1.0_dp - IF (ispin == 1) CYCLE - CASE DEFAULT - CPABORT("Unknown constraint type.") - END SELECT - group(igroup)%integrated(:, iatom) = group(igroup)%integrated(:, iatom) + sign* & - group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i)* & - rho_r(ispin)%pw%cr3d(k, j, i)*dvol - END DO + END DO + DO j = 1, npme + iatom = cores(j) + atom_a = atom_list(iatom) + pab(1, 1) = coef*hirshfeld_env%charges(atom_a) + ra(:) = pbc(particle_set(atom_a)%r, cell) + + IF (hirshfeld_control%use_atomic_cutoff) THEN + radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, & + ra=ra, rb=ra, rp=ra, & + zetp=alpha, eps=hirshfeld_control%atomic_cutoff, & + pab=pab, o1=0, o2=0, & ! without map_consistent + prefactor=1.0_dp, cutoff=0.0_dp) + ! IF (para_env%ionode) PRINT *, "radius", radius + END IF + + IF (igroup == 1) THEN + CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & + (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, & + rs_rho_all, radius=radius, & + ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & + subpatch_pattern=subpatch_pattern) + END IF + + IF (is_constraint(atom_a)) THEN + CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & + (/0.0_dp, 0.0_dp, 0.0_dp/), coefficients(atom_a), & + pab, 0, 0, rs_rho_constr, & + radius=radius, & + ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & + subpatch_pattern=subpatch_pattern) + END IF + + IF (hirshfeld_control%print_density .AND. igroup == 1) THEN + IF (is_constraint(atom_a)) THEN + DO iatom = 1, cdft_control%natoms + IF (atom_a == cdft_control%atoms(iatom)) EXIT END DO + CPASSERT(iatom <= cdft_control%natoms) + CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & + (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, & + rs_single(iatom)%rs_grid, radius=radius, & + ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & + subpatch_pattern=subpatch_pattern) + END IF + END IF + + IF (cdft_control%atomic_charges .AND. igroup == 1) THEN + IF (compute_charge(atom_a)) THEN + DO iatom = 1, cdft_control%natoms + IF (atom_a == cdft_control%atoms(iatom)) EXIT + END DO + CPASSERT(iatom <= cdft_control%natoms) + CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & + (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, & + rs_single_charge(iatom)%rs_grid, radius=radius, & + ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & + subpatch_pattern=subpatch_pattern) + END IF + END IF + + END DO + END DO + DEALLOCATE (cores) + END DO + DEALLOCATE (pab) + + IF (igroup == 1) THEN + CALL rs_pw_transfer(rs_rho_all, cdft_control%group(igroup)%hw_rho_total%pw, rs2pw) + IF (.NOT. cdft_control%in_memory) CALL rs_grid_release(rs_rho_all) + END IF + + CALL rs_pw_transfer(rs_rho_constr, cdft_control%group(igroup)%hw_rho_total_constraint%pw, rs2pw) + CALL rs_grid_release(rs_rho_constr) + + ! Calculate weight function + CALL hfun_scale(cdft_control%group(igroup)%weight%pw%cr3d, & + cdft_control%group(igroup)%hw_rho_total_constraint%pw%cr3d, & + cdft_control%group(1)%hw_rho_total%pw%cr3d, divide=.TRUE., & + small=hirshfeld_control%eps_cutoff) + + ! Calculate charges + IF (cdft_control%atomic_charges .AND. igroup == 1) THEN + DO i = 1, cdft_control%natoms + CALL rs_pw_transfer(rs_single_charge(i)%rs_grid, cdft_control%group(igroup)%hw_rho_atomic_charge(i)%pw, rs2pw) + CALL hfun_scale(cdft_control%charge(i)%pw%cr3d, & + cdft_control%group(igroup)%hw_rho_atomic_charge(i)%pw%cr3d, & + cdft_control%group(igroup)%hw_rho_total%pw%cr3d, divide=.TRUE., & + small=hirshfeld_control%eps_cutoff) + END DO + END IF + + ! Print atomic densities if requested + IF (hirshfeld_control%print_density .AND. igroup == 1) THEN + DO i = 1, cdft_control%natoms + CALL rs_pw_transfer(rs_single(i)%rs_grid, cdft_control%group(igroup)%hw_rho_atomic(i)%pw, rs2pw) + END DO + CALL cdft_print_hirshfeld_density(qs_env) + END IF + + END DO + + DO igroup = 1, SIZE(cdft_control%group) + + CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%group(igroup)%hw_rho_total_constraint%pw) + + IF (.NOT. cdft_control%in_memory .AND. igroup == 1) THEN + CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%group(1)%hw_rho_total%pw) + END IF + + IF (hirshfeld_control%print_density .AND. igroup == 1) THEN + DO i = 1, cdft_control%natoms + CALL rs_grid_release(rs_single(i)%rs_grid) + CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%group(igroup)%hw_rho_atomic(i)%pw) + END DO + DEALLOCATE (rs_single) + DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic) + END IF + + IF (cdft_control%atomic_charges .AND. igroup == 1) THEN + DO i = 1, cdft_control%natoms + CALL rs_grid_release(rs_single_charge(i)%rs_grid) + CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%group(igroup)%hw_rho_atomic_charge(i)%pw) + END DO + DEALLOCATE (rs_single_charge) + DEALLOCATE (compute_charge) + DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge) + END IF + + END DO + + IF (cdft_control%in_memory) THEN + DO igroup = 1, SIZE(cdft_control%group) + ALLOCATE (cdft_control%group(igroup)%gradients_x(1*natom, lb_pw(1):ub_pw(1), & + lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3))) + cdft_control%group(igroup)%gradients_x(:, :, :, :) = 0.0_dp + END DO + END IF + + IF (cdft_control%in_memory) THEN + DO igroup = 1, SIZE(cdft_control%group) + + ALLOCATE (pab(1, 1)) + nthread = 1 + ithread = 0 + atoms_memory = hirshfeld_control%atoms_memory + + DO ikind = 1, SIZE(atomic_kind_set) + numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp + IF (numexp <= 0) CYCLE + CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list) + + ALLOCATE (pw_single_dr(num_species)) + ALLOCATE (rs_single_dr(num_species)) + + DO i = 1, num_species + CALL pw_pool_create_pw(auxbas_pw_pool, pw_single_dr(i)%pw, use_data=REALDATA3D, & + in_space=REALSPACE) + pw_single_dr(i)%pw%cr3d(:, :, :) = 0.0_dp + END DO + + atoms_memory_num = SIZE((/(j, j=1, num_species, atoms_memory)/)) + + ! Can't store all pw grids, therefore split into groups of size atom_memory + ! Ideally this code should be re-written to be more memory efficient + IF (num_species > atoms_memory) THEN + ALLOCATE (num_species_small(atoms_memory_num + 1)) + num_species_small(1:atoms_memory_num) = (/(j, j=1, num_species, atoms_memory)/) + num_species_small(atoms_memory_num + 1) = num_species + ELSE + ALLOCATE (num_species_small(2)) + num_species_small(:) = (/1, num_species/) + END IF + + DO k = 1, SIZE(num_species_small) - 1 + IF (num_species > atoms_memory) THEN + ALLOCATE (cores(num_species_small(k + 1) - (num_species_small(k) - 1))) + ELSE + ALLOCATE (cores(num_species)) + END IF + + DO i = num_species_small(k), num_species_small(k + 1) + NULLIFY (rs_single_dr(i)%rs_grid) + CALL rs_grid_create(rs_single_dr(i)%rs_grid, auxbas_rs_desc) + CALL rs_grid_zero(rs_single_dr(i)%rs_grid) + END DO + DO iex = 1, numexp + + alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex) + coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex) + prefactor = 2.0_dp*alpha + npme = 0 + cores = 0 + + DO iatom = 1, SIZE(cores) + atom_a = atom_list(iatom + (num_species_small(k) - 1)) + ra(:) = pbc(particle_set(atom_a)%r, cell) + + IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN + IF (MODULO(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN + npme = npme + 1 + cores(npme) = iatom + END IF + ELSE + npme = npme + 1 + cores(npme) = iatom + END IF + END DO + DO j = 1, npme + iatom = cores(j) + atom_a = atom_list(iatom + (num_species_small(k) - 1)) + pab(1, 1) = coef*hirshfeld_env%charges(atom_a) + ra(:) = pbc(particle_set(atom_a)%r, cell) + subpatch_pattern = 0 + + ! Calculate cutoff + IF (hirshfeld_control%use_atomic_cutoff) THEN + radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, & + ra=ra, rb=ra, rp=ra, & + zetp=alpha, eps=hirshfeld_control%atomic_cutoff, & + pab=pab, o1=0, o2=0, & ! without map_consistent + prefactor=1.0_dp, cutoff=0.0_dp) + END IF + + CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & + (/0.0_dp, 0.0_dp, 0.0_dp/), prefactor, & + pab, 0, 0, rs_single_dr(iatom + (num_species_small(k) - 1))%rs_grid, & + radius=radius, & + ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & + subpatch_pattern=subpatch_pattern) + + END DO + END DO + + DO iatom = num_species_small(k), num_species_small(k + 1) + CALL rs_pw_transfer(rs_single_dr(iatom)%rs_grid, pw_single_dr(iatom)%pw, rs2pw) + CALL rs_grid_release(rs_single_dr(iatom)%rs_grid) + END DO + + DEALLOCATE (cores) + END DO + + DO iatom = 1, num_species + atom_a = atom_list(iatom) + cdft_control%group(igroup)%gradients_x(atom_a, :, :, :) = pw_single_dr(iatom)%pw%cr3d(:, :, :) + CALL pw_pool_give_back_pw(auxbas_pw_pool, pw_single_dr(iatom)%pw) + END DO + + DEALLOCATE (rs_single_dr) + DEALLOCATE (num_species_small) + DEALLOCATE (pw_single_dr) + END DO + DEALLOCATE (pab) + END DO + END IF + + IF (cdft_control%in_memory) THEN + DO igroup = 1, SIZE(cdft_control%group) + ALLOCATE (cdft_control%group(igroup)%gradients_y(1*num_atoms, lb_pw(1):ub_pw(1), & + lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3))) + ALLOCATE (cdft_control%group(igroup)%gradients_z(1*num_atoms, lb_pw(1):ub_pw(1), & + lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3))) + cdft_control%group(igroup)%gradients_y(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :) + cdft_control%group(igroup)%gradients_z(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :) + END DO + END IF + + ! Calculate gradient if requested + IF (cdft_control%in_memory) THEN + + DO igroup = 1, SIZE(cdft_control%group) + DO k = lb_pw(3), ub_pw(3) + DO j = lb_pw(2), ub_pw(2) + DO i = lb_pw(1), ub_pw(1) + DO iatom = 1, natom + + ra(:) = particle_set(iatom)%r + + IF (cdft_control%group(igroup)%hw_rho_total%pw%cr3d(i, j, k) > hirshfeld_control%eps_cutoff) THEN + + exp_eval = (coefficients(iatom) - & + cdft_control%group(igroup)%weight%pw%cr3d(i, j, k))/ & + cdft_control%group(igroup)%hw_rho_total%pw%cr3d(i, j, k) + + r2 = (/i*dr_pw(1), j*dr_pw(2), k*dr_pw(3)/) + origin + r_pbc = pbc(ra, r2, cell) + + ! Store gradient d/dR_x w, including term: (r_x - R_x) + cdft_control%group(igroup)%gradients_x(iatom, i, j, k) = & + cdft_control%group(igroup)%gradients_x(iatom, i, j, k)* & + r_pbc(1)*exp_eval + + ! Store gradient d/dR_y w, including term: (r_y - R_y) + cdft_control%group(igroup)%gradients_y(iatom, i, j, k) = & + cdft_control%group(igroup)%gradients_y(iatom, i, j, k)* & + r_pbc(2)*exp_eval + + ! Store gradient d/dR_z w, including term:(r_z - R_z) + cdft_control%group(igroup)%gradients_z(iatom, i, j, k) = & + cdft_control%group(igroup)%gradients_z(iatom, i, j, k)* & + r_pbc(3)*exp_eval + + END IF + END DO END DO END DO END DO + + IF (igroup == 1) THEN + CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%group(1)%hw_rho_total%pw) + END IF + CALL rs_grid_release(rs_rho_all) END DO END IF - IF (.NOT. cdft_control%transfer_pot) THEN - DO igroup = 1, SIZE(group) - DEALLOCATE (group(igroup)%gradients) - END DO - END IF - DO igroup = 1, SIZE(group) - CALL mp_sum(group(igroup)%integrated, para_env%group) - END DO - ! Update force only on master process. Otherwise force due to constraint becomes multiplied - ! by the number of processes when the final force%rho_elec is constructed in qs_force - ! by mp_summing [the final integrated(:,:) is distributed on all processors] - IF (para_env%ionode) THEN - DO igroup = 1, SIZE(group) - DO iatom = 1, natom - ikind = kind_of(iatom) - i = atom_of_kind(iatom) - force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + group(igroup)%integrated(:, iatom)*strength(igroup) - END DO - END DO - END IF - DEALLOCATE (strength) - DEALLOCATE (atom_of_kind) - DEALLOCATE (kind_of) - DO igroup = 1, SIZE(group) - DEALLOCATE (group(igroup)%integrated) - END DO - NULLIFY (group) + + IF (ALLOCATED(coefficients)) DEALLOCATE (coefficients) + IF (ALLOCATED(is_constraint)) DEALLOCATE (is_constraint) CALL timestop(handle) - END SUBROUTINE becke_constraint_force + END SUBROUTINE hirshfeld_constraint_low + ! ************************************************************************************************** + !> \brief Calculates the value of a CDFT constraint by integrating the product of the CDFT + !> weight function and the realspace electron density + ! ************************************************************************************************** ! ************************************************************************************************** -!> \brief Calculates the value of a CDFT constraint by integrating the product of the CDFT -!> weight function and the realspace electron density -!> \param qs_env the qs_env where to build the constraint -!> \par History -!> Created 3.2017 [Nico Holmberg] -!> Generalized to any CDFT weight function 9.2018 [Nico Holmberg] +!> \brief ... +!> \param qs_env ... ! ************************************************************************************************** SUBROUTINE cdft_constraint_integrate(qs_env) TYPE(qs_environment_type), POINTER :: qs_env @@ -951,12 +1328,280 @@ CONTAINS END SUBROUTINE cdft_constraint_integrate + ! ************************************************************************************************** + !> \brief Calculates atomic forces due to a CDFT constraint (Becke or Hirshfeld) + ! ************************************************************************************************** ! ************************************************************************************************** -!> \brief Prepare CDFT fragment constraints. Fragment densities are read from cube files, multiplied -!> by the CDFT weight functions and integrated over the realspace grid. -!> \param qs_env the qs_env where to build the constraint -!> \par History -!> Created 9.2018 [Nico Holmberg] +!> \brief ... +!> \param qs_env ... +! ************************************************************************************************** + SUBROUTINE cdft_constraint_force(qs_env) + TYPE(qs_environment_type), POINTER :: qs_env + + CHARACTER(len=*), PARAMETER :: routineN = 'cdft_constraint_force', & + routineP = moduleN//':'//routineN + + INTEGER :: handle, i, iatom, igroup, ikind, ispin, & + j, k, natom, nvar + INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of + INTEGER, DIMENSION(2, 3) :: bo + INTEGER, DIMENSION(3) :: lb, ub + REAL(kind=dp) :: dvol, eps_cavity, sign + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: strength + REAL(KIND=dp), DIMENSION(:), POINTER :: cutoffs + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(becke_constraint_type), POINTER :: becke_control + TYPE(cdft_control_type), POINTER :: cdft_control + TYPE(cdft_group_type), DIMENSION(:), POINTER :: group + TYPE(cell_type), POINTER :: cell + TYPE(cp_para_env_type), POINTER :: para_env + TYPE(dft_control_type), POINTER :: dft_control + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r + TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(qs_rho_type), POINTER :: rho + + CALL timeset(routineN, handle) + NULLIFY (atomic_kind_set, cell, para_env, dft_control, particle_set, & + rho, rho_r, force, cutoffs, becke_control, group) + + CALL get_qs_env(qs_env, & + atomic_kind_set=atomic_kind_set, & + natom=natom, & + particle_set=particle_set, & + cell=cell, & + rho=rho, & + force=force, & + dft_control=dft_control, & + para_env=para_env) + CALL qs_rho_get(rho, rho_r=rho_r) + + cdft_control => dft_control%qs_control%cdft_control + becke_control => cdft_control%becke_control + group => cdft_control%group + nvar = SIZE(cdft_control%target) + ALLOCATE (strength(nvar)) + strength(:) = cdft_control%strength(:) + cutoffs => cdft_control%becke_control%cutoffs + eps_cavity = cdft_control%becke_control%eps_cavity + ALLOCATE (atom_of_kind(natom)) + ALLOCATE (kind_of(natom)) + + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & + atom_of_kind=atom_of_kind, & + kind_of=kind_of) + DO igroup = 1, SIZE(cdft_control%group) + ALLOCATE (cdft_control%group(igroup)%integrated(3, natom)) + cdft_control%group(igroup)%integrated = 0.0_dp + END DO + + lb(1:3) = rho_r(1)%pw%pw_grid%bounds_local(1, 1:3) + ub(1:3) = rho_r(1)%pw%pw_grid%bounds_local(2, 1:3) + bo = cdft_control%group(1)%weight%pw%pw_grid%bounds_local + dvol = cdft_control%group(1)%weight%pw%pw_grid%dvol + + IF (cdft_control%type == outer_scf_becke_constraint) THEN + IF (.NOT. cdft_control%becke_control%in_memory) THEN + CALL becke_constraint_low(qs_env, just_gradients=.TRUE.) + END IF + + ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN + IF (.NOT. cdft_control%in_memory) THEN + CALL hirshfeld_constraint_low(qs_env, just_gradients=.TRUE.) + END IF + END IF + + ! If no Becke Gaussian confinement + IF (.NOT. ASSOCIATED(becke_control%cavity_mat)) THEN + ! No external control + DO k = bo(1, 1), bo(2, 1) + DO j = bo(1, 2), bo(2, 2) + DO i = bo(1, 3), bo(2, 3) + ! First check if this grid point should be skipped + IF (cdft_control%becke_control%cavity_confine) THEN + IF (cdft_control%becke_control%cavity%pw%cr3d(k, j, i) < eps_cavity) CYCLE + END IF + + DO igroup = 1, SIZE(cdft_control%group) + DO iatom = 1, natom + DO ispin = 1, dft_control%nspins + + SELECT CASE (cdft_control%group(igroup)%constraint_type) + CASE (cdft_charge_constraint) + sign = 1.0_dp + CASE (cdft_magnetization_constraint) + IF (ispin == 1) THEN + sign = 1.0_dp + ELSE + sign = -1.0_dp + END IF + CASE (cdft_alpha_constraint) + sign = 1.0_dp + IF (ispin == 2) CYCLE + CASE (cdft_beta_constraint) + sign = 1.0_dp + IF (ispin == 1) CYCLE + CASE DEFAULT + CPABORT("Unknown constraint type.") + END SELECT + + IF (cdft_control%type == outer_scf_becke_constraint) THEN + + cdft_control%group(igroup)%integrated(:, iatom) = & + cdft_control%group(igroup)%integrated(:, iatom) + sign* & + cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) & + *rho_r(ispin)%pw%cr3d(k, j, i) & + *dvol + + ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN + + cdft_control%group(igroup)%integrated(1, iatom) = & + cdft_control%group(igroup)%integrated(1, iatom) + sign* & + cdft_control%group(igroup)%gradients_x(iatom, k, j, i) & + *rho_r(ispin)%pw%cr3d(k, j, i) & + *dvol + + cdft_control%group(igroup)%integrated(2, iatom) = & + cdft_control%group(igroup)%integrated(2, iatom) + sign* & + cdft_control%group(igroup)%gradients_y(iatom, k, j, i) & + *rho_r(ispin)%pw%cr3d(k, j, i) & + *dvol + + cdft_control%group(igroup)%integrated(3, iatom) = & + cdft_control%group(igroup)%integrated(3, iatom) + sign* & + cdft_control%group(igroup)%gradients_z(iatom, k, j, i) & + *rho_r(ispin)%pw%cr3d(k, j, i) & + *dvol + + END IF + + END DO + END DO + END DO + END DO + END DO + END DO + + ! If Becke Gaussian confinement + ELSE + DO k = LBOUND(cdft_control%becke_control%cavity_mat, 1), UBOUND(cdft_control%becke_control%cavity_mat, 1) + DO j = LBOUND(cdft_control%becke_control%cavity_mat, 2), UBOUND(cdft_control%becke_control%cavity_mat, 2) + DO i = LBOUND(cdft_control%becke_control%cavity_mat, 3), UBOUND(cdft_control%becke_control%cavity_mat, 3) + + ! First check if this grid point should be skipped + IF (cdft_control%becke_control%cavity_mat(k, j, i) < eps_cavity) CYCLE + + DO igroup = 1, SIZE(group) + DO iatom = 1, natom + DO ispin = 1, dft_control%nspins + SELECT CASE (group(igroup)%constraint_type) + CASE (cdft_charge_constraint) + sign = 1.0_dp + CASE (cdft_magnetization_constraint) + IF (ispin == 1) THEN + sign = 1.0_dp + ELSE + sign = -1.0_dp + END IF + CASE (cdft_alpha_constraint) + sign = 1.0_dp + IF (ispin == 2) CYCLE + CASE (cdft_beta_constraint) + sign = 1.0_dp + IF (ispin == 1) CYCLE + CASE DEFAULT + CPABORT("Unknown constraint type.") + END SELECT + + ! Integrate gradient of weight function + IF (cdft_control%type == outer_scf_becke_constraint) THEN + + cdft_control%group(igroup)%integrated(:, iatom) = & + cdft_control%group(igroup)%integrated(:, iatom) + sign* & + cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) & + *rho_r(ispin)%pw%cr3d(k, j, i) & + *dvol + + ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN + + cdft_control%group(igroup)%integrated(1, iatom) = & + cdft_control%group(igroup)%integrated(1, iatom) + sign* & + cdft_control%group(igroup)%gradients_x(iatom, k, j, i) & + *rho_r(ispin)%pw%cr3d(k, j, i) & + *dvol + + cdft_control%group(igroup)%integrated(2, iatom) = & + cdft_control%group(igroup)%integrated(2, iatom) + sign* & + cdft_control%group(igroup)%gradients_y(iatom, k, j, i) & + *rho_r(ispin)%pw%cr3d(k, j, i) & + *dvol + + cdft_control%group(igroup)%integrated(3, iatom) = & + cdft_control%group(igroup)%integrated(3, iatom) + sign* & + cdft_control%group(igroup)%gradients_z(iatom, k, j, i) & + *rho_r(ispin)%pw%cr3d(k, j, i) & + *dvol + + END IF + + END DO + END DO + END DO + END DO + END DO + END DO + END IF + + IF (.NOT. cdft_control%transfer_pot) THEN + IF (cdft_control%type == outer_scf_becke_constraint) THEN + DO igroup = 1, SIZE(group) + DEALLOCATE (cdft_control%group(igroup)%gradients) + END DO + ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN + DO igroup = 1, SIZE(group) + DEALLOCATE (cdft_control%group(igroup)%gradients_x) + DEALLOCATE (cdft_control%group(igroup)%gradients_y) + DEALLOCATE (cdft_control%group(igroup)%gradients_z) + END DO + END IF + END IF + + DO igroup = 1, SIZE(group) + CALL mp_sum(group(igroup)%integrated, para_env%group) + END DO + + ! Update force only on master process. Otherwise force due to constraint becomes multiplied + ! by the number of processes when the final force%rho_elec is constructed in qs_force + ! by mp_summing [the final integrated(:,:) is distributed on all processors] + IF (para_env%ionode) THEN + DO igroup = 1, SIZE(group) + DO iatom = 1, natom + ikind = kind_of(iatom) + i = atom_of_kind(iatom) + force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + group(igroup)%integrated(:, iatom)*strength(igroup) + END DO + END DO + END IF + + DEALLOCATE (strength) + DEALLOCATE (atom_of_kind) + DEALLOCATE (kind_of) + DO igroup = 1, SIZE(group) + DEALLOCATE (group(igroup)%integrated) + END DO + NULLIFY (group) + + CALL timestop(handle) + + END SUBROUTINE cdft_constraint_force + + ! ************************************************************************************************** + !> \brief Prepare CDFT fragment constraints. Fragment densities are read from cube files, multiplied + !> by the CDFT weight functions and integrated over the realspace grid. + ! ************************************************************************************************** +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... ! ************************************************************************************************** SUBROUTINE prepare_fragment_constraint(qs_env) TYPE(qs_environment_type), POINTER :: qs_env @@ -1106,259 +1751,4 @@ CONTAINS END SUBROUTINE prepare_fragment_constraint -! ************************************************************************************************** -!> \brief Driver routine for calculating a Gaussian Hirshfeld constraint -!> \param qs_env the qs_env where to build the constraint -!> \param calc_pot if the potential needs to be recalculated or just integrated -!> \param calculate_forces logical if potential has to be calculated or only_energy -!> \par History -!> Created 09.2018 [Nico Holmberg] -! ************************************************************************************************** - SUBROUTINE hirshfeld_constraint(qs_env, calc_pot, calculate_forces) - TYPE(qs_environment_type), POINTER :: qs_env - LOGICAL :: calc_pot, calculate_forces - - CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint' - - INTEGER :: handle - TYPE(cdft_control_type), POINTER :: cdft_control - TYPE(dft_control_type), POINTER :: dft_control - - CALL timeset(routineN, handle) - CALL get_qs_env(qs_env, dft_control=dft_control) - cdft_control => dft_control%qs_control%cdft_control - IF (dft_control%qs_control%cdft .AND. cdft_control%type == outer_scf_hirshfeld_constraint) THEN - IF (calc_pot) THEN - ! Initialize the Hirshfeld constraint environment - CALL hirshfeld_constraint_init(qs_env) - ! Calculate the Hirshfeld weight function and possibly the gradients - CALL hirshfeld_constraint_low(qs_env) - END IF - ! Integrate the Hirshfeld constraint - CALL cdft_constraint_integrate(qs_env) - ! Calculate forces - IF (calculate_forces) CPABORT("Hirshfeld force NYI.") - !CALL hirshfeld_constraint_force(qs_env) - END IF - CALL timestop(handle) - - END SUBROUTINE hirshfeld_constraint - -! ************************************************************************************************** -!> \brief Calculates Gaussian Hirshfeld constraints -!> \param qs_env the qs_env where to build the constraint -!> \author Nico Holmberg (09.2018) -! ************************************************************************************************** - SUBROUTINE hirshfeld_constraint_low(qs_env) - TYPE(qs_environment_type), POINTER :: qs_env - - CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint_low' - - INTEGER :: atom_a, handle, i, iatom, iex, igroup, & - ikind, ithread, j, natom, npme, & - nthread, numexp, subpatch_pattern - INTEGER, DIMENSION(:), POINTER :: atom_list, cores - LOGICAL, ALLOCATABLE, DIMENSION(:) :: compute_charge, is_constraint - REAL(kind=dp) :: alpha, coef, dvol, eps_rho_rspace, & - radius, radius_constr - REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coefficients - REAL(kind=dp), DIMENSION(3) :: ra - REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab - TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set - TYPE(cdft_control_type), POINTER :: cdft_control - TYPE(cell_type), POINTER :: cell - TYPE(dft_control_type), POINTER :: dft_control - TYPE(hirshfeld_constraint_type), POINTER :: hirshfeld_control - TYPE(hirshfeld_type), POINTER :: hirshfeld_env - TYPE(particle_type), DIMENSION(:), POINTER :: particle_set - TYPE(pw_env_type), POINTER :: pw_env - TYPE(pw_p_type) :: tmp - TYPE(pw_p_type), POINTER :: fnorm - TYPE(pw_pool_type), POINTER :: auxbas_pw_pool - TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc - TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: rs_single - TYPE(realspace_grid_type), POINTER :: rs_rho_all, rs_rho_constr - - NULLIFY (atom_list, cores, pab, atomic_kind_set, cell, dft_control, & - hirshfeld_env, particle_set, pw_env, fnorm, auxbas_pw_pool, & - auxbas_rs_desc, rs_rho_all, rs_rho_constr, cdft_control, & - hirshfeld_control, rs_single) - - CALL timeset(routineN, handle) - CALL get_qs_env(qs_env, & - atomic_kind_set=atomic_kind_set, & - cell=cell, & - particle_set=particle_set, & - natom=natom, & - dft_control=dft_control, & - pw_env=pw_env) - - cdft_control => dft_control%qs_control%cdft_control - hirshfeld_control => cdft_control%hirshfeld_control - hirshfeld_env => hirshfeld_control%hirshfeld_env - - ! Build weight function - ALLOCATE (is_constraint(natom)) - ALLOCATE (coefficients(natom)) - IF (cdft_control%atomic_charges) THEN - ALLOCATE (compute_charge(natom)) - compute_charge = .FALSE. - END IF - - DO igroup = 1, SIZE(cdft_control%group) - dvol = cdft_control%group(igroup)%weight%pw%pw_grid%dvol - cdft_control%group(igroup)%weight%pw%cr3d = 0.0_dp - ! Keep track of those atoms that are active in this constraint group - coefficients(:) = 0.0_dp - is_constraint = .FALSE. - DO i = 1, SIZE(cdft_control%group(igroup)%atoms) - coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i) - is_constraint(cdft_control%group(igroup)%atoms(i)) = .TRUE. - END DO - ! Calculate atomic Hirshfeld weight functions: rs_rho_constr / rs_rho_all - ! rs_rho_all: Superposition of isolated Gaussian densities over all atoms in system - ! rs_rho_constr: Sum of isolated Gaussian densities over constraint atoms in this constraint group - IF (.NOT. ASSOCIATED(rs_rho_all)) THEN - CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, auxbas_rs_grid=rs_rho_all, & - auxbas_pw_pool=auxbas_pw_pool) - CALL rs_grid_retain(rs_rho_all) - CALL rs_grid_zero(rs_rho_all) - END IF - CALL rs_grid_create(rs_rho_constr, auxbas_rs_desc) - CALL rs_grid_zero(rs_rho_constr) - ! Compute Gaussian density over single atoms (rs_single) when atomic charges are requested - IF (cdft_control%atomic_charges .AND. igroup == 1) THEN - ALLOCATE (rs_single(cdft_control%natoms)) - DO i = 1, cdft_control%natoms - NULLIFY (rs_single(i)%rs_grid) - CALL rs_grid_create(rs_single(i)%rs_grid, auxbas_rs_desc) - CALL rs_grid_zero(rs_single(i)%rs_grid) - compute_charge(cdft_control%atoms(i)) = .TRUE. - END DO - END IF - - ! Collocate Gaussians - eps_rho_rspace = dft_control%qs_control%eps_rho_rspace - ALLOCATE (pab(1, 1)) - nthread = 1 - ithread = 0 - - DO ikind = 1, SIZE(atomic_kind_set) - numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp - IF (numexp <= 0) CYCLE - CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list) - ALLOCATE (cores(natom)) - - DO iex = 1, numexp - alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex) - coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex) - npme = 0 - cores = 0 - DO iatom = 1, natom - atom_a = atom_list(iatom) - ra(:) = pbc(particle_set(atom_a)%r, cell) - IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN - ! replicated realspace grid, split the atoms up between procs - IF (MODULO(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN - npme = npme + 1 - cores(npme) = iatom - END IF - ELSE - npme = npme + 1 - cores(npme) = iatom - END IF - END DO - DO j = 1, npme - iatom = cores(j) - atom_a = atom_list(iatom) - pab(1, 1) = coef*hirshfeld_env%charges(atom_a) - ra(:) = pbc(particle_set(atom_a)%r, cell) - subpatch_pattern = 0 - radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, & - ra=ra, rb=ra, rp=ra, & - zetp=alpha, eps=eps_rho_rspace, & - pab=pab, o1=0, o2=0, & ! without map_consistent - prefactor=1.0_dp, cutoff=0.0_dp) - - CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & - (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, & - rs_rho_all, radius=radius, & - ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & - subpatch_pattern=subpatch_pattern) - IF (is_constraint(atom_a)) THEN - radius_constr = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, & - ra=ra, rb=ra, rp=ra, & - zetp=alpha, eps=eps_rho_rspace, & - pab=pab, o1=0, o2=0, & ! without map_consistent - prefactor=coefficients(atom_a), cutoff=0.0_dp) - - CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & - (/0.0_dp, 0.0_dp, 0.0_dp/), coefficients(atom_a), & - pab, 0, 0, rs_rho_constr, & - radius=radius_constr, & - ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & - subpatch_pattern=subpatch_pattern) - END IF - IF (cdft_control%atomic_charges .AND. igroup == 1) THEN - IF (compute_charge(atom_a)) THEN - DO iatom = 1, cdft_control%natoms - IF (atom_a == cdft_control%atoms(iatom)) EXIT - END DO - CPASSERT(iatom <= cdft_control%natoms) - CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & - (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, & - rs_single(iatom)%rs_grid, radius=radius, & - ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & - subpatch_pattern=subpatch_pattern) - END IF - END IF - END DO - END DO - DEALLOCATE (cores) - END DO - DEALLOCATE (pab) - - ! Transfer rs_rho_all to the correct grid and save it - CALL get_hirshfeld_info(hirshfeld_env, fnorm=fnorm) - IF (igroup == 1) THEN - IF (ASSOCIATED(fnorm)) THEN - CALL pw_pool_give_back_pw(auxbas_pw_pool, fnorm%pw) - END IF - ALLOCATE (fnorm) - CALL pw_pool_create_pw(auxbas_pw_pool, fnorm%pw, use_data=REALDATA3D, & - in_space=REALSPACE) - CALL set_hirshfeld_info(hirshfeld_env, fnorm=fnorm) - CALL rs_pw_transfer(rs_rho_all, fnorm%pw, rs2pw) - CALL rs_grid_release(rs_rho_all) - END IF - ! Compute CDFT weight function - CALL pw_pool_create_pw(auxbas_pw_pool, tmp%pw, use_data=REALDATA3D, & - in_space=REALSPACE) - CALL rs_pw_transfer(rs_rho_constr, tmp%pw, rs2pw) - CALL rs_grid_release(rs_rho_constr) - CALL hfun_scale(cdft_control%group(igroup)%weight%pw%cr3d, tmp%pw%cr3d, & - fnorm%pw%cr3d, divide=.TRUE.) - CALL pw_pool_give_back_pw(auxbas_pw_pool, tmp%pw) - ! Compute atomic weight functions if charges are needed - IF (cdft_control%atomic_charges .AND. igroup == 1) THEN - CALL pw_pool_create_pw(auxbas_pw_pool, tmp%pw, use_data=REALDATA3D, & - in_space=REALSPACE) - DO i = 1, cdft_control%natoms - CALL rs_pw_transfer(rs_single(i)%rs_grid, tmp%pw, rs2pw) - CALL rs_grid_release(rs_single(i)%rs_grid) - CALL hfun_scale(cdft_control%charge(i)%pw%cr3d, tmp%pw%cr3d, & - fnorm%pw%cr3d, divide=.TRUE.) - END DO - CALL pw_pool_give_back_pw(auxbas_pw_pool, tmp%pw) - DEALLOCATE (rs_single) - DEALLOCATE (compute_charge) - END IF - END DO - - DEALLOCATE (is_constraint) - DEALLOCATE (coefficients) - CALL timestop(handle) - - END SUBROUTINE hirshfeld_constraint_low - END MODULE qs_cdft_methods diff --git a/src/qs_cdft_types.F b/src/qs_cdft_types.F index c1df944f43..d61a734f0c 100644 --- a/src/qs_cdft_types.F +++ b/src/qs_cdft_types.F @@ -18,6 +18,8 @@ MODULE qs_cdft_types USE hirshfeld_types, ONLY: hirshfeld_type,& release_hirshfeld_type USE input_constants, ONLY: becke_cutoff_global,& + outer_scf_becke_constraint,& + outer_scf_hirshfeld_constraint,& outer_scf_none,& radius_single,& shape_function_gaussian @@ -98,12 +100,17 @@ MODULE qs_cdft_types !> \param radius Gaussian radius parameter !> \param shape_function the constraint type: atomic density or single Gaussian !> \param use_bohr determines whether to use angstrom or bohr units for the radii of Gaussians +!> \param use_atomic_cutoff Logical to control use of ATOMIC_CUTOFF +!> \param atomic_cutoff Numerical cutoff for calculation of Hirshfeld densities +!> \param atoms_memory Number of atomic gradients to store in memory +!> \param eps_cutoff Numerical cutoff for calculation of weight function +!> \param print_density Logical to control printing of Hirshfeld densities to .cube file !> \param hirshfeld_env auxiliary type storing information about the Gaussians ! ************************************************************************************************** TYPE hirshfeld_constraint_type - INTEGER :: gaussian_shape, shape_function - LOGICAL :: use_bohr - REAL(KIND=dp) :: radius + INTEGER :: gaussian_shape, shape_function, atoms_memory + LOGICAL :: use_bohr, print_density, use_atomic_cutoff + REAL(KIND=dp) :: radius, eps_cutoff, atomic_cutoff REAL(KIND=dp), DIMENSION(:), POINTER :: radii TYPE(hirshfeld_type), POINTER :: hirshfeld_env END TYPE hirshfeld_constraint_type @@ -154,6 +161,7 @@ MODULE qs_cdft_types !> in a mixed_env that holds multiple CDFT states !> \param first_iteration a flag to mark the first iteration for printing of additional data !> \param print_weight logical which determines if CDFT weight functions should be saved to a file +!> \param in_memory logical which determines if the gradients of the Becke potential should be !> \param is_constraint list of logicals which determines if an atom is included in a constraint group !> \param strength Lagrangian multipliers of the constraints !> \param target target values of the constraints @@ -194,8 +202,19 @@ MODULE qs_cdft_types ! Atomic gradients of the weight function at every grid point REAL(KIND=dp), POINTER, & DIMENSION(:, :, :, :) :: gradients + REAL(KIND=dp), POINTER, & + DIMENSION(:, :, :, :) :: gradients_x + REAL(KIND=dp), POINTER, & + DIMENSION(:, :, :, :) :: gradients_y + REAL(KIND=dp), POINTER, & + DIMENSION(:, :, :, :) :: gradients_z ! The weight function of this constraint group TYPE(pw_p_type) :: weight + TYPE(pw_p_type), DIMENSION(:), POINTER :: hw_rho_atomic + TYPE(pw_p_type), DIMENSION(:), POINTER :: hw_rho_atomic_dr + TYPE(pw_p_type), DIMENSION(:), POINTER :: hw_rho_atomic_charge + TYPE(pw_p_type) :: hw_rho_total_constraint + TYPE(pw_p_type) :: hw_rho_total END TYPE cdft_group_type TYPE cdft_control_type @@ -214,7 +233,7 @@ MODULE qs_cdft_types atomic_charges, fragment_density, & fragments_integrated, flip_fragment(2), & transfer_pot, external_control, & - first_iteration, print_weight + first_iteration, print_weight, in_memory LOGICAL, POINTER, DIMENSION(:) :: is_constraint REAL(KIND=dp), DIMENSION(:), POINTER :: strength, TARGET, value REAL(KIND=dp), POINTER, & @@ -371,6 +390,7 @@ CONTAINS cdft_control%should_purge = .FALSE. cdft_control%purge_history = .FALSE. cdft_control%calculate_metric = .FALSE. + cdft_control%in_memory = .FALSE. cdft_control%purge_freq = 0 cdft_control%nbad_conv = 0 cdft_control%purge_offset = 0 @@ -451,8 +471,17 @@ CONTAINS DEALLOCATE (cdft_control%group(i)%coeff) IF (ALLOCATED(cdft_control%group(i)%d_sum_const_dR)) & DEALLOCATE (cdft_control%group(i)%d_sum_const_dR) - IF (ASSOCIATED(cdft_control%group(i)%gradients)) & - DEALLOCATE (cdft_control%group(i)%gradients) + IF (cdft_control%type == outer_scf_becke_constraint) THEN + IF (ASSOCIATED(cdft_control%group(i)%gradients)) & + DEALLOCATE (cdft_control%group(i)%gradients) + ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN + IF (ASSOCIATED(cdft_control%group(i)%gradients_x)) & + DEALLOCATE (cdft_control%group(i)%gradients_x) + IF (ASSOCIATED(cdft_control%group(i)%gradients_y)) & + DEALLOCATE (cdft_control%group(i)%gradients_y) + IF (ASSOCIATED(cdft_control%group(i)%gradients_z)) & + DEALLOCATE (cdft_control%group(i)%gradients_z) + END IF IF (ASSOCIATED(cdft_control%group(i)%integrated)) & DEALLOCATE (cdft_control%group(i)%integrated) END DO @@ -518,8 +547,13 @@ CONTAINS ALLOCATE (hirshfeld_control) hirshfeld_control%use_bohr = .FALSE. + hirshfeld_control%print_density = .FALSE. + hirshfeld_control%use_atomic_cutoff = .TRUE. hirshfeld_control%radius = 3.0_dp + hirshfeld_control%eps_cutoff = 1.0e-12_dp + hirshfeld_control%atomic_cutoff = 1.0e-12_dp hirshfeld_control%shape_function = shape_function_gaussian + hirshfeld_control%atoms_memory = 80 hirshfeld_control%gaussian_shape = radius_single NULLIFY (hirshfeld_control%hirshfeld_env) NULLIFY (hirshfeld_control%radii) diff --git a/src/qs_cdft_utils.F b/src/qs_cdft_utils.F index 5fe00b22a0..00b493762e 100644 --- a/src/qs_cdft_utils.F +++ b/src/qs_cdft_utils.F @@ -86,7 +86,8 @@ MODULE qs_cdft_utils ! *** Public subroutines *** PUBLIC :: becke_constraint_init, read_becke_section, read_cdft_control_section - PUBLIC :: hfun_scale, hirshfeld_constraint_init, cdft_constraint_print + PUBLIC :: hfun_scale, hirshfeld_constraint_init, cdft_constraint_print, & + cdft_print_hirshfeld_density, cdft_print_weight_function CONTAINS @@ -564,7 +565,6 @@ CONTAINS ALLOCATE (cdft_control%group(k)%atoms(natoms)) ALLOCATE (cdft_control%group(k)%coeff(natoms)) NULLIFY (cdft_control%group(k)%weight%pw) - NULLIFY (cdft_control%group(k)%gradients) NULLIFY (cdft_control%group(k)%integrated) tot_natoms = tot_natoms + natoms ! Now parse @@ -726,16 +726,19 @@ CONTAINS TYPE(qs_control_type), INTENT(INOUT) :: qs_control TYPE(section_vals_type), POINTER :: cdft_control_section + INTEGER :: k, nvar LOGICAL :: exists TYPE(cdft_control_type), POINTER :: cdft_control - TYPE(section_vals_type), POINTER :: becke_constraint_section, & + TYPE(section_vals_type), POINTER :: becke_constraint_section, group_section, & hirshfeld_constraint_section, & outer_scf_section, print_section NULLIFY (outer_scf_section, hirshfeld_constraint_section, becke_constraint_section, & - print_section) + print_section, group_section) cdft_control => qs_control%cdft_control CPASSERT(ASSOCIATED(cdft_control)) + group_section => section_vals_get_subs_vals(cdft_control_section, "ATOM_GROUP") + CALL section_vals_get(group_section, n_repetition=nvar, explicit=exists) CALL section_vals_val_get(cdft_control_section, "TYPE_OF_CONSTRAINT", & i_val=qs_control%cdft_control%type) @@ -771,11 +774,19 @@ CONTAINS becke_constraint_section => section_vals_get_subs_vals(cdft_control_section, "BECKE_CONSTRAINT") CALL section_vals_get(becke_constraint_section, explicit=exists) IF (.NOT. exists) CPABORT("BECKE_CONSTRAINT section is missing.") + DO k = 1, nvar + NULLIFY (cdft_control%group(k)%gradients) + END DO CALL read_becke_section(cdft_control, becke_constraint_section) CASE (outer_scf_hirshfeld_constraint) hirshfeld_constraint_section => section_vals_get_subs_vals(cdft_control_section, "HIRSHFELD_CONSTRAINT") CALL section_vals_get(hirshfeld_constraint_section, explicit=exists) IF (.NOT. exists) CPABORT("HIRSHFELD_CONSTRAINT section is missing.") + DO k = 1, nvar + NULLIFY (cdft_control%group(k)%gradients_x) + NULLIFY (cdft_control%group(k)%gradients_y) + NULLIFY (cdft_control%group(k)%gradients_z) + END DO CALL read_hirshfeld_constraint_section(cdft_control, hirshfeld_constraint_section) CASE DEFAULT CPABORT("Unknown constraint type.") @@ -796,7 +807,6 @@ CONTAINS !> \brief reads the input parameters needed for Hirshfeld constraint !> \param cdft_control the cdft_control which holds the Hirshfeld constraint !> \param hirshfeld_section the input section for a Hirshfeld constraint -!> \author Nico Holmberg [09.2018] ! ************************************************************************************************** SUBROUTINE read_hirshfeld_constraint_section(cdft_control, hirshfeld_section) TYPE(cdft_control_type), INTENT(INOUT) :: cdft_control @@ -814,6 +824,11 @@ CONTAINS CALL section_vals_val_get(hirshfeld_section, "GAUSSIAN_SHAPE", i_val=hirshfeld_control%gaussian_shape) CALL section_vals_val_get(hirshfeld_section, "GAUSSIAN_RADIUS", r_val=hirshfeld_control%radius) CALL section_vals_val_get(hirshfeld_section, "USE_BOHR", l_val=hirshfeld_control%use_bohr) + CALL section_vals_val_get(hirshfeld_section, "USE_ATOMIC_CUTOFF", l_val=hirshfeld_control%use_atomic_cutoff) + CALL section_vals_val_get(hirshfeld_section, "PRINT_DENSITY", l_val=hirshfeld_control%print_density) + CALL section_vals_val_get(hirshfeld_section, "EPS_CUTOFF", r_val=hirshfeld_control%eps_cutoff) + CALL section_vals_val_get(hirshfeld_section, "ATOMIC_CUTOFF", r_val=hirshfeld_control%atomic_cutoff) + IF (.NOT. hirshfeld_control%use_bohr) THEN hirshfeld_control%radius = cp_unit_from_cp2k(hirshfeld_control%radius, "angstrom") END IF @@ -837,19 +852,21 @@ CONTAINS END SUBROUTINE read_hirshfeld_constraint_section -! ************************************************************************************************** + ! ************************************************************************************************** !> \brief Calculate fout = fun1/fun2 or fout = fun1*fun2 !> \param fout the output 3D potential !> \param fun1 the first input 3D potential !> \param fun2 the second input 3D potential !> \param divide logical that decides whether to divide or multiply the input potentials +!> \param small customisable parameter to determine lower bound of division ! ************************************************************************************************** - SUBROUTINE hfun_scale(fout, fun1, fun2, divide) + SUBROUTINE hfun_scale(fout, fun1, fun2, divide, small) REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: fout REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: fun1, fun2 LOGICAL, INTENT(IN) :: divide + REAL(KIND=dp), INTENT(IN) :: small - REAL(KIND=dp), PARAMETER :: small = 1.0e-12_dp + CHARACTER(len=*), PARAMETER :: routineN = 'hfun_scale', routineP = moduleN//':'//routineN INTEGER :: i1, i2, i3, n1, n2, n3 @@ -1250,9 +1267,7 @@ CONTAINS ! ************************************************************************************************** !> \brief Prints CDFT weight functions to cube files -!> \param qs_env the qs_env where to build the constraint -!> \par History -!> Created 9.2018 [Nico Holmberg] +!> \param qs_env ... ! ************************************************************************************************** SUBROUTINE cdft_print_weight_function(qs_env) TYPE(qs_environment_type), POINTER :: qs_env @@ -1300,4 +1315,82 @@ CONTAINS END SUBROUTINE cdft_print_weight_function + ! ************************************************************************************************** +!> \brief Prints Hirshfeld weight function and promolecule density +!> \param qs_env ... +! ************************************************************************************************** + SUBROUTINE cdft_print_hirshfeld_density(qs_env) + TYPE(qs_environment_type), POINTER :: qs_env + + CHARACTER(len=*), PARAMETER :: routineN = 'cdft_print_hirshfeld_density', & + routineP = moduleN//':'//routineN + + CHARACTER(LEN=default_path_length) :: middle_name + INTEGER :: iatom, igroup, unit_nr + LOGICAL :: mpi_io + TYPE(cdft_control_type), POINTER :: cdft_control + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_para_env_type), POINTER :: para_env + TYPE(dft_control_type), POINTER :: dft_control + TYPE(particle_list_type), POINTER :: particles + TYPE(pw_env_type), POINTER :: pw_env + TYPE(qs_subsys_type), POINTER :: subsys + TYPE(section_vals_type), POINTER :: cdft_constraint_section + + NULLIFY (cdft_constraint_section, logger, particles, dft_control, & + para_env, subsys, cdft_control, pw_env) + logger => cp_get_default_logger() + + CALL get_qs_env(qs_env, subsys=subsys, para_env=para_env, dft_control=dft_control, pw_env=pw_env) + CALL qs_subsys_get(subsys, particles=particles) + cdft_control => dft_control%qs_control%cdft_control + cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT") + + mpi_io = .TRUE. + + DO igroup = 1, SIZE(cdft_control%group) + + middle_name = "hw_rho_total"//TRIM(ADJUSTL(cp_to_string(igroup))) + unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io, & + file_position="REWIND", middle_name=middle_name, extension=".cube") + + CALL cp_pw_to_cube(cdft_control%group(igroup)%hw_rho_total%pw, unit_nr, "CDFT Weight Function", mpi_io=mpi_io, & + particles=particles, stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE")) + + CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io) + + END DO + + DO igroup = 1, SIZE(cdft_control%group) + + middle_name = "hw_rho_total_constraint_"//TRIM(ADJUSTL(cp_to_string(igroup))) + unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io, & + file_position="REWIND", middle_name=middle_name, extension=".cube") + + CALL cp_pw_to_cube(cdft_control%group(igroup)%hw_rho_total_constraint%pw, unit_nr, & + "CDFT Weight Function", mpi_io=mpi_io, particles=particles, & + stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE")) + + CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io) + + END DO + + DO igroup = 1, SIZE(cdft_control%group) + DO iatom = 1, (cdft_control%natoms) + + middle_name = "hw_rho_atomic_"//TRIM(ADJUSTL(cp_to_string(iatom))) + unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io, & + file_position="REWIND", middle_name=middle_name, extension=".cube") + + CALL cp_pw_to_cube(cdft_control%group(igroup)%hw_rho_atomic(iatom)%pw, unit_nr, & + "CDFT Weight Function", mpi_io=mpi_io, particles=particles, & + stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE")) + + CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io) + + END DO + END DO + + END SUBROUTINE cdft_print_hirshfeld_density + END MODULE qs_cdft_utils diff --git a/tests/QS/regtest-cdft-hirshfeld-2/TEST_FILES b/tests/QS/regtest-cdft-hirshfeld-2/TEST_FILES index 6813493d96..531ba506cc 100644 --- a/tests/QS/regtest-cdft-hirshfeld-2/TEST_FILES +++ b/tests/QS/regtest-cdft-hirshfeld-2/TEST_FILES @@ -3,16 +3,15 @@ # see regtest/TEST_FILES # HeH-noconstraint.inp 1 2e-13 -3.01067446615063 -HeH-cdft-state-1.inp 71 5e-08 0.161058337228 -HeH-cdft-state-2.inp 71 1e-09 1.849710970513 +HeH-cdft-state-1.inp 71 1e-7 0.161058337228 +HeH-cdft-state-2.inp 71 1e-7 1.849710970513 # These tests compute the electronic coupling and related quantities # Coupling via rotation -HeH-mixed-cdft-1.inp 73 5e-10 204.357138287654 +HeH-mixed-cdft-1.inp 73 1e-7 204.357138287654 # Coupling via Lowdin orthogonalization -HeH-mixed-cdft-2.inp 76 5e-10 154.149819062673 +HeH-mixed-cdft-2.inp 76 1e-7 154.149819062673 # CDFT-CI -HeH-mixed-cdft-3.inp 77 3e-12 -2.80477609770094 +HeH-mixed-cdft-3.inp 77 1e-7 -2.80477609770094 # Rerun some tests using parallel mode (NGROUPS 2) # Electronic coupling and related -HeH-mixed-cdft-8.inp 73 5e-10 204.357138287651 -#EOF +HeH-mixed-cdft-8.inp 73 1e-7 204.357138287651 \ No newline at end of file diff --git a/tests/QS/regtest-cdft-hirshfeld-3/HeH-cdft-1.inp b/tests/QS/regtest-cdft-hirshfeld-3/HeH-cdft-1.inp new file mode 100644 index 0000000000..b20c465203 --- /dev/null +++ b/tests/QS/regtest-cdft-hirshfeld-3/HeH-cdft-1.inp @@ -0,0 +1,32 @@ +@SET RESTART_WFN TRUE +@SET WFN_FILE HeH-noconstraint-1_0.wfn +@SET PROJECT_NAME HeH-cdft-1 +@SET WRITE_WFN 0 +@SET CHARGE 1 +@SET WRITE_CUBE FALSE +@SET CENTER_SYS ON +@SET XYZFILE HeH.xyz + +@SET HIRSHFELD_ACTIVE TRUE +@SET HIRSHFELD_FRAGMENT FALSE +@SET PRINT_WEIGHT_FUNCTION FALSE +@SET MAX_SCF 0 + +@SET HIRSHFELD_TARGET 1.0 +@SET HIRSHFELD_STR 0.186894372937 +@SET HIRSHFELD_SHAPEFN DENSITY + +&GLOBAL + PROJECT ${PROJECT_NAME} + RUN_TYPE ENERGY_FORCE + PRINT_LEVEL MEDIUM + PREFERRED_DIAG_LIBRARY SL +&END GLOBAL + +&FORCE_EVAL + METHOD QS + @SET HIRSHFELD_STR ${HIRSHFELD_STR} + @SET HIRSHFELD_TARGET ${HIRSHFELD_TARGET} + @include dft-common-params.inc + @include subsys.inc +&END FORCE_EVAL diff --git a/tests/QS/regtest-cdft-hirshfeld-3/HeH-noconstraint.inp b/tests/QS/regtest-cdft-hirshfeld-3/HeH-noconstraint.inp new file mode 100644 index 0000000000..9311b49c0d --- /dev/null +++ b/tests/QS/regtest-cdft-hirshfeld-3/HeH-noconstraint.inp @@ -0,0 +1,23 @@ +@SET RESTART_WFN FALSE +@SET WFN_FILE HeH-noconstraint-1_0.wfn +@SET PROJECT_NAME HeH-noconstraint +@SET WRITE_WFN 1 +@SET CHARGE 1 +@SET WRITE_CUBE FALSE +@SET CENTER_SYS ON +@SET XYZFILE HeH.xyz + +@SET HIRSHFELD_ACTIVE FALSE + +&GLOBAL + PROJECT ${PROJECT_NAME} + RUN_TYPE ENERGY_FORCE + PRINT_LEVEL MEDIUM + PREFERRED_DIAG_LIBRARY SL +&END GLOBAL + +&FORCE_EVAL + METHOD QS + @include dft-common-params.inc + @include subsys.inc +&END FORCE_EVAL diff --git a/tests/QS/regtest-cdft-hirshfeld-3/HeH.xyz b/tests/QS/regtest-cdft-hirshfeld-3/HeH.xyz new file mode 100644 index 0000000000..493b1a9a40 --- /dev/null +++ b/tests/QS/regtest-cdft-hirshfeld-3/HeH.xyz @@ -0,0 +1,4 @@ +2 + +He 3.000000 3.000000 2.600000 +H 3.000000 3.000000 3.400000 diff --git a/tests/QS/regtest-cdft-hirshfeld-3/TEST_FILES b/tests/QS/regtest-cdft-hirshfeld-3/TEST_FILES new file mode 100644 index 0000000000..87d1edcd80 --- /dev/null +++ b/tests/QS/regtest-cdft-hirshfeld-3/TEST_FILES @@ -0,0 +1,7 @@ +# runs are executed in the same order as in this file +# the second field tells which test should be run in order to compare with the last available output +# see regtest/TEST_FILES +# +HeH-noconstraint.inp 72 1e-7 0.1450972684448 +HeH-cdft-1.inp 72 1e-7 0.1552195046628 +#EOF diff --git a/tests/QS/regtest-cdft-hirshfeld-3/dft-common-params.inc b/tests/QS/regtest-cdft-hirshfeld-3/dft-common-params.inc new file mode 100644 index 0000000000..f660b284f6 --- /dev/null +++ b/tests/QS/regtest-cdft-hirshfeld-3/dft-common-params.inc @@ -0,0 +1,82 @@ +&DFT + @IF ( ${HIRSHFELD_ACTIVE} == TRUE ) + @include hirshfeld_qs.inc + @ENDIF + @IF ( ${HIRSHFELD_ACTIVE} == FALSE ) + &QS + METHOD GPW + EPS_DEFAULT 1.0E-12 + EXTRAPOLATION ASPC + EXTRAPOLATION_ORDER 3 + &END QS + @ENDIF + BASIS_SET_FILE_NAME BASIS_MOLOPT + POTENTIAL_FILE_NAME POTENTIAL + @IF ( ${RESTART_WFN} == TRUE ) + WFN_RESTART_FILE_NAME ${WFN_FILE} + @ENDIF + LSD + CHARGE ${CHARGE} + &MGRID + CUTOFF 100 + NGRIDS 5 + &END MGRID + &SCF + @IF ( ${RESTART_WFN} == TRUE ) + SCF_GUESS RESTART + @ENDIF + @IF ( ${RESTART_WFN} == FALSE ) + SCF_GUESS ATOMIC + @ENDIF + EPS_SCF 1.0E-5 + CHOLESKY INVERSE_DBCSR + MAX_SCF 20 + &OT ON + ! DIIS convergence might not be constant with different mpiranks => use CG + MINIMIZER CG + PRECONDITIONER FULL_ALL + ALGORITHM IRAC + &END OT + &OUTER_SCF ON + EPS_SCF 1.0E-5 + MAX_SCF 2 + &END + &PRINT + &RESTART + FILENAME ./${PROJECT_NAME} + BACKUP_COPIES 0 + COMMON_ITERATION_LEVELS 1 + &EACH + JUST_ENERGY ${WRITE_WFN} + QS_SCF 0 + &END EACH + &END RESTART + &RESTART_HISTORY OFF + &END RESTART_HISTORY + &END PRINT + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &XC_GRID + XC_DERIV SPLINE2 + XC_SMOOTH_RHO NONE + &END XC_GRID + &END XC + &PRINT + &MULLIKEN OFF + &END MULLIKEN + &HIRSHFELD OFF + &END HIRSHFELD + @IF ( ${WRITE_CUBE} == TRUE ) + &E_DENSITY_CUBE ON + STRIDE 1 1 1 + &END E_DENSITY_CUBE + @ENDIF + &END PRINT +&END DFT +&PRINT + &FORCES ON + NDIGITS 13 + &END +&END diff --git a/tests/QS/regtest-cdft-hirshfeld-3/hirshfeld_qs.inc b/tests/QS/regtest-cdft-hirshfeld-3/hirshfeld_qs.inc new file mode 100644 index 0000000000..42fb4be0c9 --- /dev/null +++ b/tests/QS/regtest-cdft-hirshfeld-3/hirshfeld_qs.inc @@ -0,0 +1,46 @@ +&QS + METHOD GPW + EPS_DEFAULT 1.0E-12 + EXTRAPOLATION ASPC + EXTRAPOLATION_ORDER 3 + &CDFT + TYPE_OF_CONSTRAINT HIRSHFELD + STRENGTH ${HIRSHFELD_STR} + TARGET ${HIRSHFELD_TARGET} + + &ATOM_GROUP + ATOMS 1 2 + COEFF 1 -1 + &END ATOM_GROUP + + &OUTER_SCF ON + TYPE CDFT_CONSTRAINT + MAX_SCF 50 + + EPS_SCF 1.0E-3 + OPTIMIZER NEWTON_LS + + &CDFT_OPT ON + CONTINUE_LS TRUE + &END CDFT_OPT + &END + + &HIRSHFELD_CONSTRAINT + SHAPE_FUNCTION DENSITY + &END HIRSHFELD_CONSTRAINT + + &PROGRAM_RUN_INFO ON + ! Print weight function to cube file + @IF ( ${PRINT_WEIGHT_FUNCTION} == TRUE ) + &WEIGHT_FUNCTION + &END WEIGHT_FUNCTION + @ENDIF + &EACH + QS_SCF 1 + &END EACH + COMMON_ITERATION_LEVELS 2 + ADD_LAST NUMERIC + FILENAME ./${PROJECT_NAME} + &END PROGRAM_RUN_INFO + &END CDFT +&END QS diff --git a/tests/QS/regtest-cdft-hirshfeld-3/subsys.inc b/tests/QS/regtest-cdft-hirshfeld-3/subsys.inc new file mode 100644 index 0000000000..877394e3ec --- /dev/null +++ b/tests/QS/regtest-cdft-hirshfeld-3/subsys.inc @@ -0,0 +1,20 @@ +&SUBSYS + &CELL + PERIODIC XYZ + ABC 6. 6. 6. + &END CELL + &TOPOLOGY + COORD_FILE_FORMAT XYZ + COORD_FILE_NAME ${XYZFILE} + &CENTER_COORDINATES ${CENTER_SYS} + &END CENTER_COORDINATES + &END TOPOLOGY + &KIND He + BASIS_SET SZV-MOLOPT-SR-GTH + POTENTIAL GTH-PBE-q2 + &END KIND + &KIND H + BASIS_SET SZV-MOLOPT-SR-GTH + POTENTIAL GTH-PBE-q1 + &END KIND +&END SUBSYS diff --git a/tests/QS/regtest-cdft-hirshfeld/TEST_FILES b/tests/QS/regtest-cdft-hirshfeld/TEST_FILES index 433351895b..e599d9ea53 100644 --- a/tests/QS/regtest-cdft-hirshfeld/TEST_FILES +++ b/tests/QS/regtest-cdft-hirshfeld/TEST_FILES @@ -6,8 +6,7 @@ HeH-noconstraint.inp 1 2e-13 He+-noconstraint.inp 1 2e-13 -2.05007934458372 H-noconstraint.inp 1 4e-13 -0.45734465780293 # These tests use different constraint formalisms so their value differs (see outputted charges) -HeH-cdft-1.inp 71 3e-11 1.135281643193 -HeH-cdft-2.inp 71 2e-11 1.681172869379 -HeH-cdft-3.inp 71 8e-08 1.197176189824 -HeH-cdft-4.inp 71 2e-11 1.567706385507 -#EOF +HeH-cdft-1.inp 71 1e-7 1.135281643193 +HeH-cdft-2.inp 71 1e-7 1.681172869379 +HeH-cdft-3.inp 71 1e-7 1.197176189824 +HeH-cdft-4.inp 71 1e-7 1.567706385507 \ No newline at end of file diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 3d4bfc4a6c..f22e1c03e7 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -44,6 +44,7 @@ QS/regtest-negf QS/regtest-negf-fft fftw3 QS/regtest-ri-rpa-axk libint QS/regtest-ri-rpa-rse libint +QS/regtest-cdft-hirshfeld-3 QS/regtest-cdft-hirshfeld-2 parallel mpiranks>1 QS/regtest-cdft-hirshfeld SIRIUS/regtest-1 sirius