1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of the derivative of the QMMM Hamiltonian integral 8!> matrix <a|\sum_i q_i|b> for semi-empirical methods 9!> \author Teodoro Laino - 04.2007 [tlaino] 10! ************************************************************************************************** 11MODULE qmmm_se_forces 12 USE atomic_kind_types, ONLY: atomic_kind_type,& 13 get_atomic_kind 14 USE cell_types, ONLY: cell_type,& 15 pbc 16 USE cp_control_types, ONLY: dft_control_type 17 USE cp_para_types, ONLY: cp_para_env_type 18 USE dbcsr_api, ONLY: dbcsr_get_block_p,& 19 dbcsr_p_type 20 USE input_constants, ONLY: & 21 do_method_am1, do_method_mndo, do_method_mndod, do_method_pchg, do_method_pdg, & 22 do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1 23 USE kinds, ONLY: dp 24 USE message_passing, ONLY: mp_sum 25 USE multipole_types, ONLY: do_multipole_none 26 USE particle_types, ONLY: particle_type 27 USE qmmm_types_low, ONLY: qmmm_env_qm_type,& 28 qmmm_pot_p_type,& 29 qmmm_pot_type 30 USE qmmm_util, ONLY: spherical_cutoff_factor 31 USE qs_environment_types, ONLY: get_qs_env,& 32 qs_environment_type 33 USE qs_kind_types, ONLY: get_qs_kind,& 34 qs_kind_type 35 USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type 36 USE qs_rho_types, ONLY: qs_rho_get,& 37 qs_rho_type 38 USE semi_empirical_int_arrays, ONLY: se_orbital_pointer 39 USE semi_empirical_integrals, ONLY: dcorecore,& 40 drotnuc 41 USE semi_empirical_types, ONLY: get_se_param,& 42 se_int_control_type,& 43 se_taper_type,& 44 semi_empirical_create,& 45 semi_empirical_release,& 46 semi_empirical_type,& 47 setup_se_int_control_type 48 USE semi_empirical_utils, ONLY: get_se_type,& 49 se_param_set_default 50#include "./base/base_uses.f90" 51 52 IMPLICIT NONE 53 54 PRIVATE 55 56 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_se_forces' 57 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 58 PUBLIC :: deriv_se_qmmm_matrix 59 60CONTAINS 61 62! ************************************************************************************************** 63!> \brief Constructs the derivative w.r.t. 1-el semi-empirical hamiltonian 64!> QMMM terms 65!> \param qs_env ... 66!> \param qmmm_env ... 67!> \param particles_mm ... 68!> \param mm_cell ... 69!> \param para_env ... 70!> \param calc_force ... 71!> \param Forces ... 72!> \param Forces_added_charges ... 73!> \author Teodoro Laino 04.2007 [created] 74! ************************************************************************************************** 75 SUBROUTINE deriv_se_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, & 76 calc_force, Forces, Forces_added_charges) 77 78 TYPE(qs_environment_type), POINTER :: qs_env 79 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env 80 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm 81 TYPE(cell_type), POINTER :: mm_cell 82 TYPE(cp_para_env_type), POINTER :: para_env 83 LOGICAL, INTENT(in), OPTIONAL :: calc_force 84 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces, Forces_added_charges 85 86 CHARACTER(len=*), PARAMETER :: routineN = 'deriv_se_qmmm_matrix', & 87 routineP = moduleN//':'//routineN 88 89 INTEGER :: handle, i, iatom, ikind, iqm, ispin, & 90 itype, natom, natorb_a, nkind, & 91 number_qm_atoms 92 INTEGER, DIMENSION(:), POINTER :: list 93 LOGICAL :: anag, defined, found 94 REAL(KIND=dp) :: delta, enuclear 95 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Forces_QM, p_block_a 96 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 97 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p 98 TYPE(dft_control_type), POINTER :: dft_control 99 TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm 100 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 101 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc 102 TYPE(qs_rho_type), POINTER :: rho 103 TYPE(se_int_control_type) :: se_int_control 104 TYPE(se_taper_type), POINTER :: se_taper 105 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm 106 107 CALL timeset(routineN, handle) 108 IF (calc_force) THEN 109 NULLIFY (rho, atomic_kind_set, qs_kind_set, se_taper) 110 NULLIFY (se_kind_a, se_kind_mm, particles_qm) 111 CALL get_qs_env(qs_env=qs_env, & 112 rho=rho, & 113 se_taper=se_taper, & 114 atomic_kind_set=atomic_kind_set, & 115 qs_kind_set=qs_kind_set, & 116 ks_qmmm_env=ks_qmmm_env_loc, & 117 dft_control=dft_control, & 118 particle_set=particles_qm, & 119 natom=number_qm_atoms) 120 SELECT CASE (dft_control%qs_control%method_id) 121 CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, & 122 do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl) 123 ! Go on with the calculation.. 124 CASE DEFAULT 125 ! Otherwise stop.. 126 CPABORT("Method not available") 127 END SELECT 128 anag = dft_control%qs_control%se_control%analytical_gradients 129 delta = dft_control%qs_control%se_control%delta 130 ! Setup SE integral control type 131 CALL setup_se_int_control_type( & 132 se_int_control, shortrange=.FALSE., do_ewald_r3=.FALSE., & 133 do_ewald_gks=.FALSE., integral_screening=dft_control%qs_control%se_control%integral_screening, & 134 max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.) 135 136 ! Create a fake semi-empirical type to handle the classical atom 137 ALLOCATE (Forces_QM(3, number_qm_atoms)) 138 CALL semi_empirical_create(se_kind_mm) 139 CALL se_param_set_default(se_kind_mm, 0, do_method_pchg) 140 itype = get_se_type(se_kind_mm%typ) 141 nkind = SIZE(atomic_kind_set) 142 enuclear = 0.0_dp 143 Forces_QM = 0.0_dp 144 CALL qs_rho_get(rho, rho_ao=matrix_p) 145 146 DO ispin = 1, dft_control%nspins 147 iqm = 0 148 Kinds: DO ikind = 1, nkind 149 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=list) 150 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a) 151 CALL get_se_param(se_kind_a, & 152 defined=defined, & 153 natorb=natorb_a) 154 IF (.NOT. defined .OR. natorb_a < 1) CYCLE 155 Atoms: DO i = 1, SIZE(list) 156 iqm = iqm + 1 157 iatom = list(i) 158 ! Give back block 159 NULLIFY (p_block_a) 160 CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, & 161 row=iatom, col=iatom, BLOCK=p_block_a, found=found) 162 163 IF (ASSOCIATED(p_block_a)) THEN 164 ! Expand derivative of geometrical factors 165 CALL deriv_se_qmmm_matrix_low(p_block_a, & 166 se_kind_a, & 167 se_kind_mm, & 168 qmmm_env%Potentials, & 169 particles_mm, & 170 qmmm_env%mm_atom_chrg, & 171 qmmm_env%mm_atom_index, & 172 mm_cell, & 173 iatom, & 174 itype, & 175 Forces, & 176 Forces_QM(:, iqm), & 177 se_taper, & 178 se_int_control, & 179 anag, & 180 delta, & 181 qmmm_env%spherical_cutoff, & 182 particles_qm) 183 ! Possibly added charges 184 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN 185 CALL deriv_se_qmmm_matrix_low(p_block_a, & 186 se_kind_a, & 187 se_kind_mm, & 188 qmmm_env%added_charges%potentials, & 189 qmmm_env%added_charges%added_particles, & 190 qmmm_env%added_charges%mm_atom_chrg, & 191 qmmm_env%added_charges%mm_atom_index, & 192 mm_cell, & 193 iatom, & 194 itype, & 195 Forces_added_charges, & 196 Forces_QM(:, iqm), & 197 se_taper, & 198 se_int_control, & 199 anag, & 200 delta, & 201 qmmm_env%spherical_cutoff, & 202 particles_qm) 203 END IF 204 END IF 205 END DO Atoms 206 END DO Kinds 207 END DO 208 CPASSERT(iqm == number_qm_atoms) 209 ! Transfer QM gradients to the QM particles.. 210 CALL mp_sum(Forces_QM, para_env%group) 211 iqm = 0 212 DO ikind = 1, nkind 213 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list) 214 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a) 215 CALL get_se_param(se_kind_a, & 216 defined=defined, & 217 natorb=natorb_a) 218 IF (.NOT. defined .OR. natorb_a < 1) CYCLE 219 DO i = 1, SIZE(list) 220 iqm = iqm + 1 221 iatom = qmmm_env%qm_atom_index(list(i)) 222 particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm) 223 END DO 224 END DO 225 ! MM forces will be handled directly from the QMMM module in the same way 226 ! as for GPW/GAPW methods 227 DEALLOCATE (Forces_QM) 228 CALL semi_empirical_release(se_kind_mm) 229 230 END IF 231 CALL timestop(handle) 232 END SUBROUTINE deriv_se_qmmm_matrix 233 234! ************************************************************************************************** 235!> \brief Low Level : Computes derivatives of the 1-el semi-empirical QMMM 236!> hamiltonian block w.r.t. MM and QM coordinates 237!> \param p_block_a ... 238!> \param se_kind_a ... 239!> \param se_kind_mm ... 240!> \param potentials ... 241!> \param particles_mm ... 242!> \param mm_charges ... 243!> \param mm_atom_index ... 244!> \param mm_cell ... 245!> \param IndQM ... 246!> \param itype ... 247!> \param forces ... 248!> \param forces_qm ... 249!> \param se_taper ... 250!> \param se_int_control ... 251!> \param anag ... 252!> \param delta ... 253!> \param qmmm_spherical_cutoff ... 254!> \param particles_qm ... 255!> \author Teodoro Laino 04.2007 [created] 256! ************************************************************************************************** 257 SUBROUTINE deriv_se_qmmm_matrix_low(p_block_a, se_kind_a, se_kind_mm, & 258 potentials, particles_mm, mm_charges, mm_atom_index, & 259 mm_cell, IndQM, itype, forces, forces_qm, se_taper, & 260 se_int_control, anag, delta, qmmm_spherical_cutoff, particles_qm) 261 262 REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block_a 263 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm 264 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials 265 TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm 266 REAL(KIND=dp), DIMENSION(:), POINTER :: mm_charges 267 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 268 TYPE(cell_type), POINTER :: mm_cell 269 INTEGER, INTENT(IN) :: IndQM, itype 270 REAL(KIND=dp), DIMENSION(:, :), POINTER :: forces 271 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: forces_qm 272 TYPE(se_taper_type), POINTER :: se_taper 273 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 274 LOGICAL, INTENT(IN) :: anag 275 REAL(KIND=dp), INTENT(IN) :: delta, qmmm_spherical_cutoff(2) 276 TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm 277 278 CHARACTER(len=*), PARAMETER :: routineN = 'deriv_se_qmmm_matrix_low', & 279 routineP = moduleN//':'//routineN 280 281 INTEGER :: handle, i1, i1L, i2, Imm, Imp, IndMM, & 282 Ipot, j1, j1L 283 REAL(KIND=dp) :: rt1, rt2, rt3, sph_chrg_factor 284 REAL(KIND=dp), DIMENSION(3) :: denuc, force_ab, r_pbc, rij 285 REAL(KIND=dp), DIMENSION(3, 45) :: de1b 286 TYPE(qmmm_pot_type), POINTER :: Pot 287 288 CALL timeset(routineN, handle) 289 ! Loop Over MM atoms - parallelization over MM atoms... 290 ! Loop over Pot stores atoms with the same charge 291 MainLoopPot: DO Ipot = 1, SIZE(Potentials) 292 Pot => Potentials(Ipot)%Pot 293 ! Loop over atoms belonging to this type 294 LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index) 295 Imm = Pot%mm_atom_index(Imp) 296 IndMM = mm_atom_index(Imm) 297 r_pbc = pbc(particles_mm(IndMM)%r - particles_qm(IndQM)%r, mm_cell) 298 rt1 = r_pbc(1) 299 rt2 = r_pbc(2) 300 rt3 = r_pbc(3) 301 rij = (/rt1, rt2, rt3/) 302 se_kind_mm%zeff = mm_charges(Imm) 303 ! Computes the screening factor for the spherical cutoff 304 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN 305 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor) 306 se_kind_mm%zeff = se_kind_mm%zeff*sph_chrg_factor 307 END IF 308 IF (ABS(se_kind_mm%zeff) <= EPSILON(0.0_dp)) CYCLE 309 ! Integrals derivatives involving QM - MM atoms 310 CALL drotnuc(se_kind_a, se_kind_mm, rij, itype=itype, de1b=de1b, & 311 se_int_control=se_int_control, anag=anag, delta=delta, & 312 se_taper=se_taper) 313 CALL dcorecore(se_kind_a, se_kind_mm, rij, itype=itype, denuc=denuc, & 314 se_int_control=se_int_control, anag=anag, delta=delta, & 315 se_taper=se_taper) 316 ! Nucler - Nuclear term 317 force_ab(1:3) = -denuc(1:3) 318 ! Force contribution from the QMMM Hamiltonian 319 i2 = 0 320 DO i1L = 1, se_kind_a%natorb 321 i1 = se_orbital_pointer(i1L) 322 DO j1L = 1, i1L - 1 323 j1 = se_orbital_pointer(j1L) 324 i2 = i2 + 1 325 force_ab = force_ab - 2.0_dp*de1b(:, i2)*p_block_a(i1, j1) 326 END DO 327 j1 = se_orbital_pointer(j1L) 328 i2 = i2 + 1 329 force_ab = force_ab - de1b(:, i2)*p_block_a(i1, j1) 330 END DO 331 ! The array of QM forces are really the forces 332 forces_qm(:) = forces_qm(:) - force_ab 333 ! The one of MM atoms are instead gradients 334 forces(:, Imm) = forces(:, Imm) - force_ab 335 END DO LoopMM 336 END DO MainLoopPot 337 CALL timestop(handle) 338 END SUBROUTINE deriv_se_qmmm_matrix_low 339 340END MODULE qmmm_se_forces 341