Line data    Source code

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