1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ***************************************************************************** 7!> \brief ... 8!> \author ... 9! ***************************************************************************** 10MODULE qs_basis_gradient 11 12 USE atomic_kind_types, ONLY: atomic_kind_type,& 13 get_atomic_kind,& 14 get_atomic_kind_set 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set 17 USE cp_fm_types, ONLY: cp_fm_type 18 USE cp_log_handling, ONLY: cp_get_default_logger,& 19 cp_logger_type 20 USE cp_para_types, ONLY: cp_para_env_type 21 USE dbcsr_api, ONLY: dbcsr_copy,& 22 dbcsr_p_type,& 23 dbcsr_set,& 24 dbcsr_type 25 USE input_section_types, ONLY: section_vals_type 26 USE kinds, ONLY: dp 27 USE particle_types, ONLY: particle_type 28 USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix 29 USE qs_energy_matrix_w, ONLY: qs_energies_compute_w 30 USE qs_environment_types, ONLY: get_qs_env,& 31 qs_environment_type 32 USE qs_force_types, ONLY: allocate_qs_force,& 33 deallocate_qs_force,& 34 qs_force_type,& 35 replicate_qs_force,& 36 zero_qs_force 37 USE qs_kind_types, ONLY: get_qs_kind,& 38 qs_kind_type 39 USE qs_ks_methods, ONLY: qs_ks_allocate_basics,& 40 qs_ks_update_qs_env 41 USE qs_ks_types, ONLY: qs_ks_env_type,& 42 set_ks_env 43 USE qs_mixing_utils, ONLY: mixing_allocate 44 USE qs_mo_methods, ONLY: calculate_density_matrix,& 45 make_basis_sm 46 USE qs_mo_types, ONLY: get_mo_set,& 47 mo_set_p_type 48 USE qs_neighbor_lists, ONLY: build_qs_neighbor_lists 49 USE qs_rho_methods, ONLY: qs_rho_update_rho 50 USE qs_rho_types, ONLY: qs_rho_get,& 51 qs_rho_set,& 52 qs_rho_type 53 USE qs_scf_types, ONLY: qs_scf_env_type 54 USE qs_subsys_types, ONLY: qs_subsys_set,& 55 qs_subsys_type 56 USE qs_update_s_mstruct, ONLY: qs_env_update_s_mstruct 57#include "./base/base_uses.f90" 58 59 IMPLICIT NONE 60 61 PRIVATE 62 63! *** Global parameters *** 64 65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_basis_gradient' 66 67! *** Public subroutines *** 68 69 PUBLIC :: qs_basis_center_gradient, qs_update_basis_center_pos, & 70 return_basis_center_gradient_norm 71 72CONTAINS 73 74! ***************************************************************************** 75! for development of floating basis functions 76! ***************************************************************************** 77!> \brief ... 78!> \param qs_env ... 79! ************************************************************************************************** 80 SUBROUTINE qs_basis_center_gradient(qs_env) 81 82 TYPE(qs_environment_type), POINTER :: qs_env 83 84 CHARACTER(len=*), PARAMETER :: routineN = 'qs_basis_center_gradient', & 85 routineP = moduleN//':'//routineN 86 87 INTEGER :: handle, i, iatom, ikind, img, ispin, & 88 natom, nimg, nkind, nspin 89 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind 90 LOGICAL :: floating, has_unit_metric 91 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient 92 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 93 TYPE(cp_logger_type), POINTER :: logger 94 TYPE(cp_para_env_type), POINTER :: para_env 95 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_w_kp 96 TYPE(dbcsr_type), POINTER :: matrix 97 TYPE(dft_control_type), POINTER :: dft_control 98 TYPE(qs_force_type), DIMENSION(:), POINTER :: basis_force, force, qs_force 99 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 100 TYPE(qs_ks_env_type), POINTER :: ks_env 101 TYPE(qs_scf_env_type), POINTER :: scf_env 102 TYPE(qs_subsys_type), POINTER :: subsys 103 104 CALL timeset(routineN, handle) 105 NULLIFY (logger) 106 logger => cp_get_default_logger() 107 108 ! get the gradient array 109 CALL get_qs_env(qs_env, scf_env=scf_env, natom=natom) 110 IF (ASSOCIATED(scf_env%floating_basis%gradient)) THEN 111 gradient => scf_env%floating_basis%gradient 112 CPASSERT(SIZE(gradient) == 3*natom) 113 ELSE 114 ALLOCATE (gradient(3, natom)) 115 scf_env%floating_basis%gradient => gradient 116 END IF 117 gradient = 0.0_dp 118 119 ! init the force environment 120 CALL get_qs_env(qs_env, force=force, subsys=subsys, atomic_kind_set=atomic_kind_set) 121 IF (ASSOCIATED(force)) THEN 122 qs_force => force 123 ELSE 124 NULLIFY (qs_force) 125 END IF 126 ! Allocate the force data structure 127 nkind = SIZE(atomic_kind_set) 128 ALLOCATE (natom_of_kind(nkind)) 129 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind) 130 CALL allocate_qs_force(basis_force, natom_of_kind) 131 DEALLOCATE (natom_of_kind) 132 CALL qs_subsys_set(subsys, force=basis_force) 133 CALL zero_qs_force(basis_force) 134 135 ! get atom mapping 136 ALLOCATE (atom_of_kind(natom), kind_of(natom)) 137 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of) 138 139 ! allocate energy weighted density matrices, if needed 140 CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric) 141 IF (.NOT. has_unit_metric) THEN 142 CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp) 143 IF (.NOT. ASSOCIATED(matrix_w_kp)) THEN 144 NULLIFY (matrix_w_kp) 145 CALL get_qs_env(qs_env, ks_env=ks_env, matrix_s_kp=matrix_s, dft_control=dft_control) 146 nspin = dft_control%nspins 147 nimg = dft_control%nimages 148 matrix => matrix_s(1, 1)%matrix 149 CALL dbcsr_allocate_matrix_set(matrix_w_kp, nspin, nimg) 150 DO ispin = 1, nspin 151 DO img = 1, nimg 152 ALLOCATE (matrix_w_kp(ispin, img)%matrix) 153 CALL dbcsr_copy(matrix_w_kp(ispin, img)%matrix, matrix, name="W MATRIX") 154 CALL dbcsr_set(matrix_w_kp(ispin, img)%matrix, 0.0_dp) 155 END DO 156 END DO 157 CALL set_ks_env(ks_env, matrix_w_kp=matrix_w_kp) 158 ENDIF 159 ENDIF 160 ! time to compute the w matrix 161 CALL qs_energies_compute_w(qs_env, .TRUE.) 162 163 ! core hamiltonian forces 164 CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.TRUE.) 165 ! Compute grid-based forces 166 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.TRUE.) 167 168 ! replicate forces 169 CALL replicate_qs_force(basis_force, para_env) 170 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set) 171 DO iatom = 1, natom 172 ikind = kind_of(iatom) 173 i = atom_of_kind(iatom) 174 CALL get_qs_kind(qs_kind_set(ikind), floating=floating) 175 IF (floating) gradient(1:3, iatom) = -basis_force(ikind)%total(1:3, i) 176 END DO 177 ! clean up force environment and reinitialize qs_force 178 IF (ASSOCIATED(basis_force)) CALL deallocate_qs_force(basis_force) 179 CALL qs_subsys_set(subsys, force=qs_force) 180 CALL get_qs_env(qs_env, ks_env=ks_env) 181 CALL set_ks_env(ks_env, forces_up_to_date=.FALSE.) 182 183 DEALLOCATE (atom_of_kind, kind_of) 184 185 CALL timestop(handle) 186 187 END SUBROUTINE qs_basis_center_gradient 188 189! ***************************************************************************** 190!> \brief ... returns the norm of the gradient vector, taking only floating 191!> components into account 192!> \param qs_env ... 193!> \return ... 194! ************************************************************************************************** 195 FUNCTION return_basis_center_gradient_norm(qs_env) RESULT(norm) 196 197 TYPE(qs_environment_type), POINTER :: qs_env 198 REAL(KIND=dp) :: norm 199 200 INTEGER :: iatom, ikind, natom, nfloat 201 LOGICAL :: floating 202 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient 203 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 204 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 205 TYPE(qs_scf_env_type), POINTER :: scf_env 206 207 norm = 0.0_dp 208 CALL get_qs_env(qs_env, scf_env=scf_env, particle_set=particle_set, qs_kind_set=qs_kind_set) 209 gradient => scf_env%floating_basis%gradient 210 natom = SIZE(particle_set) 211 nfloat = 0 212 DO iatom = 1, natom 213 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind) 214 CALL get_qs_kind(qs_kind_set(ikind), floating=floating) 215 IF (floating) THEN 216 nfloat = nfloat + 1 217 norm = norm + SUM(ABS(gradient(1:3, iatom))) 218 END IF 219 END DO 220 IF (nfloat > 0) THEN 221 norm = norm/(3.0_dp*REAL(nfloat, KIND=dp)) 222 END IF 223 224 END FUNCTION return_basis_center_gradient_norm 225 226! ***************************************************************************** 227!> \brief move atoms with kind float according to gradient 228!> \param qs_env ... 229! ************************************************************************************************** 230 SUBROUTINE qs_update_basis_center_pos(qs_env) 231 232 TYPE(qs_environment_type), POINTER :: qs_env 233 234 CHARACTER(len=*), PARAMETER :: routineN = 'qs_update_basis_center_pos', & 235 routineP = moduleN//':'//routineN 236 237 INTEGER :: handle, iatom, ikind, natom 238 LOGICAL :: floating 239 REAL(KIND=dp) :: alpha 240 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient 241 TYPE(cp_logger_type), POINTER :: logger 242 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 243 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 244 TYPE(qs_scf_env_type), POINTER :: scf_env 245 246 CALL timeset(routineN, handle) 247 NULLIFY (logger) 248 logger => cp_get_default_logger() 249 250 ! update positions 251 CALL get_qs_env(qs_env, scf_env=scf_env, particle_set=particle_set, qs_kind_set=qs_kind_set) 252 gradient => scf_env%floating_basis%gradient 253 natom = SIZE(particle_set) 254 alpha = 0.50_dp 255 DO iatom = 1, natom 256 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind) 257 CALL get_qs_kind(qs_kind_set(ikind), floating=floating) 258 IF (floating) THEN 259 particle_set(iatom)%r(1:3) = particle_set(iatom)%r(1:3) + alpha*gradient(1:3, iatom) 260 END IF 261 END DO 262 263 CALL qs_basis_reinit_energy(qs_env) 264 265 CALL timestop(handle) 266 267 END SUBROUTINE qs_update_basis_center_pos 268 269! ***************************************************************************** 270!> \brief rebuilds the structures after a floating basis update 271!> \param qs_env ... 272!> \par History 273!> 05.2016 created [JGH] 274!> \author JGH 275! ************************************************************************************************** 276 SUBROUTINE qs_basis_reinit_energy(qs_env) 277 TYPE(qs_environment_type), POINTER :: qs_env 278 279 CHARACTER(len=*), PARAMETER :: routineN = 'qs_basis_reinit_energy', & 280 routineP = moduleN//':'//routineN 281 282 INTEGER :: handle, ispin, nmo 283 TYPE(cp_fm_type), POINTER :: mo_coeff 284 TYPE(cp_para_env_type), POINTER :: para_env 285 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_kp 286 TYPE(dft_control_type), POINTER :: dft_control 287 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 288 TYPE(qs_rho_type), POINTER :: rho 289 TYPE(qs_scf_env_type), POINTER :: scf_env 290 TYPE(section_vals_type), POINTER :: input 291 292 CALL timeset(routineN, handle) 293 294 ! rebuild neighbor lists 295 CALL get_qs_env(qs_env, input=input, para_env=para_env) 296 CALL build_qs_neighbor_lists(qs_env, para_env, molecular=.FALSE., & 297 force_env_section=input) 298 ! update core hamiltonian 299 CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.FALSE.) 300 ! update structures 301 CALL qs_env_update_s_mstruct(qs_env) 302 ! KS matrices 303 CALL qs_ks_allocate_basics(qs_env) 304 ! reinit SCF task matrices 305 NULLIFY (scf_env) 306 CALL get_qs_env(qs_env, scf_env=scf_env, dft_control=dft_control) 307 IF (scf_env%mixing_method > 0) THEN 308 CALL mixing_allocate(qs_env, scf_env%mixing_method, scf_env%p_mix_new, & 309 scf_env%p_delta, dft_control%nspins, & 310 scf_env%mixing_store) 311 ELSE 312 NULLIFY (scf_env%p_mix_new) 313 END IF 314 CALL get_qs_env(qs_env, mos=mos, rho=rho, matrix_s_kp=matrix_s_kp) 315 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 316 DO ispin = 1, SIZE(mos) 317 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, nmo=nmo) 318 ! reorthogonalize MOs 319 CALL make_basis_sm(mo_coeff, nmo, matrix_s_kp(1, 1)%matrix) 320 ! update density matrix 321 CALL calculate_density_matrix(mos(ispin)%mo_set, rho_ao_kp(ispin, 1)%matrix) 322 END DO 323 CALL qs_rho_set(rho, rho_r_valid=.FALSE., drho_r_valid=.FALSE., rho_g_valid=.FALSE., & 324 drho_g_valid=.FALSE., tau_r_valid=.FALSE., tau_g_valid=.FALSE.) 325 CALL qs_rho_update_rho(rho, qs_env) 326 327 CALL timestop(handle) 328 329 END SUBROUTINE qs_basis_reinit_energy 330 331END MODULE qs_basis_gradient 332