1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Split and build its own idependent core_core SE interaction module 8!> \author Teodoro Laino [tlaino] - 05.2009 9!> \par History 10!> Teodoro Laino (05.2009) [tlaino] - create 11! ************************************************************************************************** 12MODULE se_core_core 13 USE atomic_kind_types, ONLY: atomic_kind_type,& 14 get_atomic_kind_set 15 USE atprop_types, ONLY: atprop_array_init,& 16 atprop_type 17 USE cell_types, ONLY: cell_type 18 USE cp_control_types, ONLY: dft_control_type,& 19 semi_empirical_control_type 20 USE cp_para_types, ONLY: cp_para_env_type 21 USE ewald_environment_types, ONLY: ewald_env_get,& 22 ewald_environment_type 23 USE ewald_pw_types, ONLY: ewald_pw_get,& 24 ewald_pw_type 25 USE input_constants, ONLY: & 26 do_method_am1, do_method_mndo, do_method_mndod, do_method_pdg, do_method_pm3, & 27 do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1 28 USE kinds, ONLY: dp 29 USE message_passing, ONLY: mp_sum 30 USE particle_types, ONLY: particle_type 31 USE qs_energy_types, ONLY: qs_energy_type 32 USE qs_environment_types, ONLY: get_qs_env,& 33 qs_environment_type 34 USE qs_force_types, ONLY: qs_force_type 35 USE qs_kind_types, ONLY: get_qs_kind,& 36 qs_kind_type 37 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 38 neighbor_list_iterate,& 39 neighbor_list_iterator_create,& 40 neighbor_list_iterator_p_type,& 41 neighbor_list_iterator_release,& 42 neighbor_list_set_p_type 43 USE semi_empirical_int_arrays, ONLY: rij_threshold 44 USE semi_empirical_integrals, ONLY: corecore,& 45 dcorecore 46 USE semi_empirical_types, ONLY: get_se_param,& 47 se_int_control_type,& 48 se_taper_type,& 49 semi_empirical_p_type,& 50 semi_empirical_type,& 51 setup_se_int_control_type 52 USE semi_empirical_utils, ONLY: finalize_se_taper,& 53 get_se_type,& 54 initialize_se_taper 55 USE virial_methods, ONLY: virial_pair_force 56 USE virial_types, ONLY: virial_type 57#include "./base/base_uses.f90" 58 59 IMPLICIT NONE 60 PRIVATE 61 62 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_core_core' 63 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 64 65 PUBLIC :: se_core_core_interaction 66 67CONTAINS 68 69! ************************************************************************************************** 70!> \brief Evaluates the core-core interactions for NDDO methods 71!> \param qs_env ... 72!> \param para_env ... 73!> \param calculate_forces ... 74!> \date 04.2008 [tlaino] 75!> \author Teodoro Laino [tlaino] - University of Zurich 76! ************************************************************************************************** 77 SUBROUTINE se_core_core_interaction(qs_env, para_env, calculate_forces) 78 TYPE(qs_environment_type), POINTER :: qs_env 79 TYPE(cp_para_env_type), POINTER :: para_env 80 LOGICAL, INTENT(in) :: calculate_forces 81 82 CHARACTER(len=*), PARAMETER :: routineN = 'se_core_core_interaction', & 83 routineP = moduleN//':'//routineN 84 85 INTEGER :: atom_a, atom_b, handle, iab, iatom, & 86 ikind, itype, jatom, jkind, natom, & 87 nkind 88 INTEGER, DIMENSION(:), POINTER :: atom_of_kind 89 LOGICAL :: anag, atener, defined, use_virial 90 LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined 91 REAL(KIND=dp) :: delta, dr1, dr3inv(3), enuc, enucij, & 92 enuclear, r2inv, r3inv, rinv 93 REAL(KIND=dp), DIMENSION(3) :: force_ab, rij 94 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 95 TYPE(atprop_type), POINTER :: atprop 96 TYPE(cell_type), POINTER :: cell 97 TYPE(dft_control_type), POINTER :: dft_control 98 TYPE(ewald_environment_type), POINTER :: ewald_env 99 TYPE(ewald_pw_type), POINTER :: ewald_pw 100 TYPE(neighbor_list_iterator_p_type), & 101 DIMENSION(:), POINTER :: nl_iterator 102 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 103 POINTER :: sab_se 104 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 105 TYPE(qs_energy_type), POINTER :: energy 106 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 107 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 108 TYPE(se_int_control_type) :: se_int_control 109 TYPE(se_taper_type), POINTER :: se_taper 110 TYPE(semi_empirical_control_type), POINTER :: se_control 111 TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_param 112 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b 113 TYPE(virial_type), POINTER :: virial 114 115 enuclear = 0.0_dp 116 NULLIFY (dft_control, cell, force, particle_set, se_control, se_taper, atomic_kind_set, & 117 virial, atprop) 118 119 CALL timeset(routineN, handle) 120 CPASSERT(ASSOCIATED(qs_env)) 121 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper, & 122 virial=virial, atprop=atprop, energy=energy) 123 124 CALL initialize_se_taper(se_taper, coulomb=.TRUE.) 125 ! Parameters 126 se_control => dft_control%qs_control%se_control 127 anag = se_control%analytical_gradients 128 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 129 CALL setup_se_int_control_type(se_int_control, do_ewald_r3=se_control%do_ewald_r3, & 130 do_ewald_gks=se_control%do_ewald_gks, integral_screening=se_control%integral_screening, & 131 shortrange=(se_control%do_ewald .OR. se_control%do_ewald_gks), & 132 max_multipole=se_control%max_multipole, pc_coulomb_int=.FALSE.) 133 134 ! atomic energy decomposition 135 atener = atprop%energy 136 IF (atener) THEN 137 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set) 138 natom = SIZE(particle_set) 139 CALL atprop_array_init(atprop%atecc, natom) 140 END IF 141 142 ! Retrieve some information if GKS ewald scheme is used 143 IF (se_control%do_ewald_gks) THEN 144 CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw) 145 CALL ewald_env_get(ewald_env, alpha=se_int_control%ewald_gks%alpha) 146 CALL ewald_pw_get(ewald_pw, pw_big_pool=se_int_control%ewald_gks%pw_pool, & 147 dg=se_int_control%ewald_gks%dg) 148 ! Virial not implemented 149 CPASSERT(.NOT. use_virial) 150 END IF 151 152 CALL get_qs_env(qs_env=qs_env, sab_se=sab_se, atomic_kind_set=atomic_kind_set, & 153 qs_kind_set=qs_kind_set) 154 155 nkind = SIZE(atomic_kind_set) 156 ! Possibly compute forces 157 IF (calculate_forces) THEN 158 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, force=force) 159 natom = SIZE(particle_set) 160 ALLOCATE (atom_of_kind(natom)) 161 delta = se_control%delta 162 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind) 163 END IF 164 165 itype = get_se_type(dft_control%qs_control%method_id) 166 167 ALLOCATE (se_kind_param(nkind), se_defined(nkind)) 168 DO ikind = 1, nkind 169 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a) 170 se_kind_param(ikind)%se_param => se_kind_a 171 CALL get_se_param(se_kind_a, defined=defined) 172 se_defined(ikind) = defined 173 END DO 174 CALL neighbor_list_iterator_create(nl_iterator, sab_se) 175 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 176 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij) 177 IF (.NOT. se_defined(ikind)) CYCLE 178 IF (.NOT. se_defined(jkind)) CYCLE 179 se_kind_a => se_kind_param(ikind)%se_param 180 se_kind_b => se_kind_param(jkind)%se_param 181 iab = ikind + nkind*(jkind - 1) 182 dr1 = DOT_PRODUCT(rij, rij) 183 enucij = 0._dp 184 IF (dr1 > rij_threshold) THEN 185 SELECT CASE (dft_control%qs_control%method_id) 186 CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, & 187 do_method_rm1, do_method_mndod, do_method_pnnl) 188 189 ! Core-Core energy term 190 CALL corecore(se_kind_a, se_kind_b, rij, enuc=enuc, itype=itype, anag=anag, & 191 se_int_control=se_int_control, se_taper=se_taper) 192 enucij = enucij + enuc 193 ! Residual integral (1/R^3) correction 194 IF (se_int_control%do_ewald_r3) THEN 195 r2inv = 1.0_dp/dr1 196 rinv = SQRT(r2inv) 197 r3inv = rinv**3 198 ! Core-Core term 199 enucij = enucij + se_kind_a%expns3_int(jkind)%expns3%core_core*r3inv 200 END IF 201 202 ! Core-Core Derivatives 203 IF (calculate_forces) THEN 204 atom_a = atom_of_kind(iatom) 205 atom_b = atom_of_kind(jatom) 206 207 CALL dcorecore(se_kind_a, se_kind_b, rij, denuc=force_ab, itype=itype, delta=delta, & 208 anag=anag, se_int_control=se_int_control, se_taper=se_taper) 209 210 ! Residual integral (1/R^3) correction 211 IF (se_int_control%do_ewald_r3) THEN 212 dr3inv = -3.0_dp*rij*r3inv*r2inv 213 ! Derivatives of core-core terms 214 force_ab = force_ab + se_kind_a%expns3_int(jkind)%expns3%core_core*dr3inv 215 END IF 216 IF (use_virial) THEN 217 CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij) 218 END IF 219 220 ! Sum up force components 221 force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1) 222 force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1) 223 224 force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2) 225 force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2) 226 227 force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3) 228 force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3) 229 END IF 230 CASE DEFAULT 231 CPABORT("") 232 END SELECT 233 ELSE 234 IF (se_int_control%do_ewald_gks) THEN 235 ! Core-Core energy term (self term in periodic systems) 236 CALL corecore(se_kind_a, se_kind_b, rij, enuc=enuc, itype=itype, anag=anag, & 237 se_int_control=se_int_control, se_taper=se_taper) 238 enucij = enucij + 0.5_dp*enuc 239 END IF 240 END IF 241 IF (atener) THEN 242 atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*enucij 243 atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*enucij 244 END IF 245 enuclear = enuclear + enucij 246 END DO 247 CALL neighbor_list_iterator_release(nl_iterator) 248 249 DEALLOCATE (se_kind_param, se_defined) 250 251 IF (calculate_forces) THEN 252 DEALLOCATE (atom_of_kind) 253 END IF 254 255 CALL mp_sum(enuclear, para_env%group) 256 energy%core_overlap = enuclear 257 energy%core_overlap0 = enuclear 258 259 CALL finalize_se_taper(se_taper) 260 CALL timestop(handle) 261 END SUBROUTINE se_core_core_interaction 262 263END MODULE se_core_core 264 265