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             : !> \par History
      10             : !>      JGH [04042007] code refactoring
      11             : ! **************************************************************************************************
      12             : MODULE virial_methods
      13             : 
      14             :    USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
      15             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      16             :                                               get_atomic_kind
      17             :    USE cell_types,                      ONLY: cell_type
      18             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      19             :                                               cp_subsys_type
      20             :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      21             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      22             :    USE kinds,                           ONLY: default_string_length,&
      23             :                                               dp
      24             :    USE mathlib,                         ONLY: det_3x3,&
      25             :                                               jacobi
      26             :    USE message_passing,                 ONLY: mp_comm_type,&
      27             :                                               mp_para_env_type
      28             :    USE particle_list_types,             ONLY: particle_list_type
      29             :    USE particle_types,                  ONLY: particle_type
      30             :    USE virial_types,                    ONLY: virial_type
      31             : #include "./base/base_uses.f90"
      32             : 
      33             :    IMPLICIT NONE
      34             : 
      35             :    PRIVATE
      36             :    PUBLIC:: one_third_sum_diag, virial_evaluate, virial_pair_force, virial_update, &
      37             :             write_stress_tensor, write_stress_tensor_components
      38             : 
      39             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'virial_methods'
      40             : 
      41             : CONTAINS
      42             : ! **************************************************************************************************
      43             : !> \brief Updates the virial given the virial and subsys
      44             : !> \param virial ...
      45             : !> \param subsys ...
      46             : !> \param para_env ...
      47             : !> \par History
      48             : !>      none
      49             : !> \author Teodoro Laino [tlaino] - 03.2008 - Zurich University
      50             : ! **************************************************************************************************
      51        5856 :    SUBROUTINE virial_update(virial, subsys, para_env)
      52             : 
      53             :       TYPE(virial_type), INTENT(INOUT)                   :: virial
      54             :       TYPE(cp_subsys_type), POINTER                      :: subsys
      55             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      56             : 
      57             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
      58             :       TYPE(distribution_1d_type), POINTER                :: local_particles
      59             :       TYPE(particle_list_type), POINTER                  :: particles
      60             : 
      61             :       CALL cp_subsys_get(subsys, local_particles=local_particles, atomic_kinds=atomic_kinds, &
      62        5856 :                          particles=particles)
      63             : 
      64             :       CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, &
      65        5856 :                            virial, para_env)
      66             : 
      67        5856 :    END SUBROUTINE virial_update
      68             : 
      69             : ! **************************************************************************************************
      70             : !> \brief Computes the kinetic part of the pressure tensor and updates
      71             : !>      the full VIRIAL (PV)
      72             : !> \param atomic_kind_set ...
      73             : !> \param particle_set ...
      74             : !> \param local_particles ...
      75             : !> \param virial ...
      76             : !> \param igroup ...
      77             : !> \par History
      78             : !>      none
      79             : !> \author CJM
      80             : ! **************************************************************************************************
      81       61798 :    SUBROUTINE virial_evaluate(atomic_kind_set, particle_set, local_particles, &
      82             :                               virial, igroup)
      83             : 
      84             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      85             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      86             :       TYPE(distribution_1d_type), POINTER                :: local_particles
      87             :       TYPE(virial_type), INTENT(INOUT)                   :: virial
      88             : 
      89             :       CLASS(mp_comm_type), INTENT(IN)                     :: igroup
      90             : 
      91             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = "virial_evaluate"
      92             : 
      93             :       INTEGER                                            :: handle, i, iparticle, iparticle_kind, &
      94             :                                                             iparticle_local, j, nparticle_kind, &
      95             :                                                             nparticle_local
      96             :       REAL(KIND=dp)                                      :: mass
      97             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
      98             : 
      99       61798 :       IF (virial%pv_availability) THEN
     100       19310 :          CALL timeset(routineN, handle)
     101       19310 :          NULLIFY (atomic_kind)
     102       19310 :          nparticle_kind = SIZE(atomic_kind_set)
     103      251030 :          virial%pv_kinetic = 0.0_dp
     104       77240 :          DO i = 1, 3
     105      193100 :             DO j = 1, i
     106      331596 :                DO iparticle_kind = 1, nparticle_kind
     107      215736 :                   atomic_kind => atomic_kind_set(iparticle_kind)
     108      215736 :                   CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     109      215736 :                   nparticle_local = local_particles%n_el(iparticle_kind)
     110     6434382 :                   DO iparticle_local = 1, nparticle_local
     111     6102786 :                      iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     112             :                      virial%pv_kinetic(i, j) = virial%pv_kinetic(i, j) + &
     113     6318522 :                                                mass*particle_set(iparticle)%v(i)*particle_set(iparticle)%v(j)
     114             :                   END DO
     115             :                END DO
     116      173790 :                virial%pv_kinetic(j, i) = virial%pv_kinetic(i, j)
     117             :             END DO
     118             :          END DO
     119             : 
     120       19310 :          CALL igroup%sum(virial%pv_kinetic)
     121             : 
     122             :          ! total virial
     123      251030 :          virial%pv_total = virial%pv_virial + virial%pv_kinetic + virial%pv_constraint
     124             : 
     125       19310 :          CALL timestop(handle)
     126             :       END IF
     127             : 
     128       61798 :    END SUBROUTINE virial_evaluate
     129             : 
     130             : ! **************************************************************************************************
     131             : !> \brief Computes the contribution to the stress tensor from two-body
     132             : !>      pair-wise forces
     133             : !> \param pv_virial ...
     134             : !> \param f0 ...
     135             : !> \param force ...
     136             : !> \param rab ...
     137             : !> \par History
     138             : !>      none
     139             : !> \author JGH
     140             : ! **************************************************************************************************
     141    55011321 :    PURE SUBROUTINE virial_pair_force(pv_virial, f0, force, rab)
     142             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT)      :: pv_virial
     143             :       REAL(KIND=dp), INTENT(IN)                          :: f0
     144             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: force, rab
     145             : 
     146             :       INTEGER                                            :: i, j
     147             : 
     148   220045284 :       DO i = 1, 3
     149   715147173 :          DO j = 1, 3
     150   660135852 :             pv_virial(i, j) = pv_virial(i, j) + f0*force(i)*rab(j)
     151             :          END DO
     152             :       END DO
     153             : 
     154    55011321 :    END SUBROUTINE virial_pair_force
     155             : 
     156             : ! **************************************************************************************************
     157             : !> \brief ...
     158             : !> \param virial ...
     159             : !> \param iw ...
     160             : !> \param cell ...
     161             : !> \param unit_string ...
     162             : !> \par History
     163             : !>      - Revised virial components (14.10.2020, MK)
     164             : !> \author JGH
     165             : ! **************************************************************************************************
     166         112 :    SUBROUTINE write_stress_tensor_components(virial, iw, cell, unit_string)
     167             : 
     168             :       TYPE(virial_type), INTENT(IN)                      :: virial
     169             :       INTEGER, INTENT(IN)                                :: iw
     170             :       TYPE(cell_type), POINTER                           :: cell
     171             :       CHARACTER(LEN=default_string_length), INTENT(IN)   :: unit_string
     172             : 
     173             :       CHARACTER(LEN=*), PARAMETER :: fmt = "(T2,A,T41,2(1X,ES19.11))"
     174             : 
     175             :       REAL(KIND=dp)                                      :: fconv
     176             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv
     177             : 
     178         112 :       IF (iw > 0) THEN
     179         112 :          CPASSERT(ASSOCIATED(cell))
     180             :          ! Conversion factor
     181         112 :          fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, TRIM(ADJUSTL(unit_string)))
     182             :          WRITE (UNIT=iw, FMT="(/,T2,A)") &
     183         112 :             "STRESS| Stress tensor components (GPW/GAPW) ["//TRIM(ADJUSTL(unit_string))//"]"
     184             :          WRITE (UNIT=iw, FMT="(T2,A,T52,A,T70,A)") &
     185         112 :             "STRESS|", "1/3 Trace", "Determinant"
     186        1456 :          pv = fconv*virial%pv_overlap
     187             :          WRITE (UNIT=iw, FMT=fmt) &
     188         112 :             "STRESS| Overlap", one_third_sum_diag(pv), det_3x3(pv)
     189        1456 :          pv = fconv*virial%pv_ekinetic
     190             :          WRITE (UNIT=iw, FMT=fmt) &
     191         112 :             "STRESS| Kinetic energy", one_third_sum_diag(pv), det_3x3(pv)
     192        1456 :          pv = fconv*virial%pv_ppl
     193             :          WRITE (UNIT=iw, FMT=fmt) &
     194         112 :             "STRESS| Local pseudopotential/core", one_third_sum_diag(pv), det_3x3(pv)
     195        1456 :          pv = fconv*virial%pv_ppnl
     196             :          WRITE (UNIT=iw, FMT=fmt) &
     197         112 :             "STRESS| Nonlocal pseudopotential", one_third_sum_diag(pv), det_3x3(pv)
     198        1456 :          pv = fconv*virial%pv_ecore_overlap
     199             :          WRITE (UNIT=iw, FMT=fmt) &
     200         112 :             "STRESS| Core charge overlap", one_third_sum_diag(pv), det_3x3(pv)
     201        1456 :          pv = fconv*virial%pv_ehartree
     202             :          WRITE (UNIT=iw, FMT=fmt) &
     203         112 :             "STRESS| Hartree", one_third_sum_diag(pv), det_3x3(pv)
     204        1456 :          pv = fconv*virial%pv_exc
     205             :          WRITE (UNIT=iw, FMT=fmt) &
     206         112 :             "STRESS| Exchange-correlation", one_third_sum_diag(pv), det_3x3(pv)
     207        1456 :          pv = fconv*virial%pv_exx
     208             :          WRITE (UNIT=iw, FMT=fmt) &
     209         112 :             "STRESS| Exact exchange (EXX)", one_third_sum_diag(pv), det_3x3(pv)
     210        1456 :          pv = fconv*virial%pv_vdw
     211             :          WRITE (UNIT=iw, FMT=fmt) &
     212         112 :             "STRESS| vdW correction", one_third_sum_diag(pv), det_3x3(pv)
     213        1456 :          pv = fconv*virial%pv_mp2
     214             :          WRITE (UNIT=iw, FMT=fmt) &
     215         112 :             "STRESS| Moller-Plesset (MP2)", one_third_sum_diag(pv), det_3x3(pv)
     216        1456 :          pv = fconv*virial%pv_nlcc
     217             :          WRITE (UNIT=iw, FMT=fmt) &
     218         112 :             "STRESS| Nonlinear core correction (NLCC)", one_third_sum_diag(pv), det_3x3(pv)
     219        1456 :          pv = fconv*virial%pv_gapw
     220             :          WRITE (UNIT=iw, FMT=fmt) &
     221         112 :             "STRESS| Local atomic parts (GAPW)", one_third_sum_diag(pv), det_3x3(pv)
     222        1456 :          pv = fconv*virial%pv_lrigpw
     223             :          WRITE (UNIT=iw, FMT=fmt) &
     224         112 :             "STRESS| Resolution-of-the-identity (LRI)", one_third_sum_diag(pv), det_3x3(pv)
     225             :          pv = fconv*(virial%pv_overlap + virial%pv_ekinetic + virial%pv_ppl + virial%pv_ppnl + &
     226             :                      virial%pv_ecore_overlap + virial%pv_ehartree + virial%pv_exc + &
     227             :                      virial%pv_exx + virial%pv_vdw + virial%pv_mp2 + virial%pv_nlcc + &
     228        1456 :                      virial%pv_gapw + virial%pv_lrigpw)
     229             :          WRITE (UNIT=iw, FMT=fmt) &
     230         112 :             "STRESS| Sum of components", one_third_sum_diag(pv), det_3x3(pv)
     231        1456 :          pv = fconv*virial%pv_virial
     232             :          WRITE (UNIT=iw, FMT=fmt) &
     233         112 :             "STRESS| Total", one_third_sum_diag(pv), det_3x3(pv)
     234             :       END IF
     235             : 
     236         112 :    END SUBROUTINE write_stress_tensor_components
     237             : 
     238             : ! **************************************************************************************************
     239             : !> \brief ...
     240             : !> \param a ...
     241             : !> \return ...
     242             : ! **************************************************************************************************
     243        1680 :    PURE FUNCTION one_third_sum_diag(a) RESULT(p)
     244             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: a
     245             :       REAL(KIND=dp)                                      :: p
     246             : 
     247        1680 :       p = (a(1, 1) + a(2, 2) + a(3, 3))/3.0_dp
     248        1680 :    END FUNCTION one_third_sum_diag
     249             : 
     250             : ! **************************************************************************************************
     251             : !> \brief Print stress tensor to output file
     252             : !>
     253             : !> \param pv_virial ...
     254             : !> \param iw ...
     255             : !> \param cell ...
     256             : !> \param unit_string ...
     257             : !> \param numerical ...
     258             : !> \author MK (26.08.2010)
     259             : ! **************************************************************************************************
     260        6074 :    SUBROUTINE write_stress_tensor(pv_virial, iw, cell, unit_string, numerical)
     261             : 
     262             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: pv_virial
     263             :       INTEGER, INTENT(IN)                                :: iw
     264             :       TYPE(cell_type), POINTER                           :: cell
     265             :       CHARACTER(LEN=default_string_length), INTENT(IN)   :: unit_string
     266             :       LOGICAL, INTENT(IN)                                :: numerical
     267             : 
     268             :       REAL(KIND=dp)                                      :: fconv
     269             :       REAL(KIND=dp), DIMENSION(3)                        :: eigval
     270             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: eigvec, stress_tensor
     271             : 
     272        6074 :       IF (iw > 0) THEN
     273        6074 :          CPASSERT(ASSOCIATED(cell))
     274             :          ! Conversion factor
     275        6074 :          fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, TRIM(ADJUSTL(unit_string)))
     276       78962 :          stress_tensor(:, :) = fconv*pv_virial(:, :)
     277        6074 :          IF (numerical) THEN
     278             :             WRITE (UNIT=iw, FMT="(/,T2,A)") &
     279          13 :                "STRESS| Numerical stress tensor ["//TRIM(ADJUSTL(unit_string))//"]"
     280             :          ELSE
     281             :             WRITE (UNIT=iw, FMT="(/,T2,A)") &
     282        6061 :                "STRESS| Analytical stress tensor ["//TRIM(ADJUSTL(unit_string))//"]"
     283             :          END IF
     284             :          WRITE (UNIT=iw, FMT="(T2,A,T14,3(19X,A1))") &
     285        6074 :             "STRESS|", "x", "y", "z"
     286             :          WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
     287        6074 :             "STRESS|      x", stress_tensor(1, 1:3)
     288             :          WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
     289        6074 :             "STRESS|      y", stress_tensor(2, 1:3)
     290             :          WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
     291        6074 :             "STRESS|      z", stress_tensor(3, 1:3)
     292             :          WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.11)") &
     293        6074 :             "STRESS| 1/3 Trace", (stress_tensor(1, 1) + &
     294             :                                   stress_tensor(2, 2) + &
     295       12148 :                                   stress_tensor(3, 3))/3.0_dp
     296             :          WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.11)") &
     297        6074 :             "STRESS| Determinant", det_3x3(stress_tensor(1:3, 1), &
     298             :                                            stress_tensor(1:3, 2), &
     299       12148 :                                            stress_tensor(1:3, 3))
     300        6074 :          eigval(:) = 0.0_dp
     301        6074 :          eigvec(:, :) = 0.0_dp
     302        6074 :          CALL jacobi(stress_tensor, eigval, eigvec)
     303        6074 :          IF (numerical) THEN
     304             :             WRITE (UNIT=iw, FMT="(/,T2,A)") &
     305             :                "STRESS| Eigenvectors and eigenvalues of the numerical stress tensor ["// &
     306          13 :                TRIM(ADJUSTL(unit_string))//"]"
     307             :          ELSE
     308             :             WRITE (UNIT=iw, FMT="(/,T2,A)") &
     309             :                "STRESS| Eigenvectors and eigenvalues of the analytical stress tensor ["// &
     310        6061 :                TRIM(ADJUSTL(unit_string))//"]"
     311             :          END IF
     312             :          WRITE (UNIT=iw, FMT="(T2,A,T14,3(1X,I19))") &
     313        6074 :             "STRESS|", 1, 2, 3
     314             :          WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
     315        6074 :             "STRESS| Eigenvalues", eigval(1:3)
     316             :          WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,F19.12))") &
     317        6074 :             "STRESS|      x", eigvec(1, 1:3)
     318             :          WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,F19.12))") &
     319        6074 :             "STRESS|      y", eigvec(2, 1:3)
     320             :          WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,F19.12))") &
     321        6074 :             "STRESS|      z", eigvec(3, 1:3)
     322             :       END IF
     323             : 
     324        6074 :    END SUBROUTINE write_stress_tensor
     325             : 
     326             : END MODULE virial_methods