LCOV - code coverage report
Current view: top level - src - qs_basis_gradient.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 0.0 % 126 0
Test Date: 2025-07-25 12:55:17 Functions: 0.0 % 4 0

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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_api,                    ONLY: dbcsr_copy,&
      19              :                                               dbcsr_p_type,&
      20              :                                               dbcsr_set,&
      21              :                                               dbcsr_type
      22              :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set
      23              :    USE cp_fm_types,                     ONLY: cp_fm_type
      24              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      25              :                                               cp_logger_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_environment_types,            ONLY: get_qs_env,&
      33              :                                               qs_environment_type
      34              :    USE qs_force_types,                  ONLY: allocate_qs_force,&
      35              :                                               deallocate_qs_force,&
      36              :                                               qs_force_type,&
      37              :                                               replicate_qs_force,&
      38              :                                               zero_qs_force
      39              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      40              :                                               qs_kind_type
      41              :    USE qs_ks_methods,                   ONLY: qs_ks_allocate_basics,&
      42              :                                               qs_ks_update_qs_env
      43              :    USE qs_ks_types,                     ONLY: get_ks_env,&
      44              :                                               qs_ks_env_type,&
      45              :                                               set_ks_env
      46              :    USE qs_matrix_w,                     ONLY: compute_matrix_w
      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            0 :       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 compute_matrix_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
        

Generated by: LCOV version 2.0-1