LCOV - code coverage report
Current view: top level - src/motion - md_conserved_quantities.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 206 219 94.1 %
Date: 2024-04-20 06:29:22 Functions: 8 9 88.9 %

          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 computes the conserved quantities for a given md ensemble
      10             : !>      and also kinetic energies, thermo/barostat stuff
      11             : !> \author gtb, 05.02.2003
      12             : ! **************************************************************************************************
      13             : MODULE md_conserved_quantities
      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 barostat_utils,                  ONLY: get_baro_energies
      18             :    USE cell_types,                      ONLY: cell_type
      19             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      20             :                                               cp_subsys_type
      21             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      22             :    USE extended_system_types,           ONLY: npt_info_type
      23             :    USE force_env_types,                 ONLY: force_env_get,&
      24             :                                               force_env_type
      25             :    USE input_constants,                 ONLY: &
      26             :         isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
      27             :         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
      28             :         npt_ia_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, reftraj_ensemble
      29             :    USE input_section_types,             ONLY: section_vals_type,&
      30             :                                               section_vals_val_get
      31             :    USE kinds,                           ONLY: dp
      32             :    USE mathconstants,                   ONLY: zero
      33             :    USE md_ener_types,                   ONLY: md_ener_type,&
      34             :                                               zero_md_ener
      35             :    USE md_environment_types,            ONLY: get_md_env,&
      36             :                                               md_environment_type,&
      37             :                                               set_md_env
      38             :    USE message_passing,                 ONLY: mp_comm_type,&
      39             :                                               mp_para_env_type
      40             :    USE particle_list_types,             ONLY: particle_list_type
      41             :    USE particle_types,                  ONLY: particle_type
      42             :    USE physcon,                         ONLY: kelvin
      43             :    USE qmmm_types,                      ONLY: qmmm_env_type
      44             :    USE qmmm_types_low,                  ONLY: force_mixing_label_QM_dynamics
      45             :    USE qmmmx_types,                     ONLY: qmmmx_env_type
      46             :    USE shell_potential_types,           ONLY: shell_kind_type
      47             :    USE simpar_types,                    ONLY: simpar_type
      48             :    USE thermostat_types,                ONLY: thermostat_type
      49             :    USE thermostat_utils,                ONLY: get_thermostat_energies
      50             : #include "../base/base_uses.f90"
      51             : 
      52             :    IMPLICIT NONE
      53             : 
      54             :    PRIVATE
      55             : 
      56             :    PUBLIC :: compute_conserved_quantity, calc_nfree_qm
      57             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_conserved_quantities'
      58             : 
      59             : CONTAINS
      60             : 
      61             : ! **************************************************************************************************
      62             : !> \brief calculates conserved quantity.
      63             : !> \param md_env ...
      64             : !> \param md_ener ...
      65             : !> \param tkind ...
      66             : !> \param tshell ...
      67             : !> \param natom ...
      68             : !> \par Input Arguments
      69             : !>     md_env is the md_environment
      70             : !>     epot is the total potential energy
      71             : !> \par Output Arguments
      72             : !>     cons is the conserved quantity
      73             : !> \par Output Optional Arguments
      74             : !>     cons_rel : relative cons. quantity (to the first md step)
      75             : !>     ekin : kinetic energy of particles
      76             : !>     temp : temperature
      77             : !>     temp_qm : temperature of the QM system in a QM/MM calculation
      78             : !> \par History
      79             : !>      none
      80             : !> \author gloria
      81             : ! **************************************************************************************************
      82       44539 :    SUBROUTINE compute_conserved_quantity(md_env, md_ener, tkind, tshell, &
      83             :                                          natom)
      84             :       TYPE(md_environment_type), POINTER                 :: md_env
      85             :       TYPE(md_ener_type), POINTER                        :: md_ener
      86             :       LOGICAL, INTENT(IN)                                :: tkind, tshell
      87             :       INTEGER, INTENT(IN)                                :: natom
      88             : 
      89             :       INTEGER                                            :: ikind, nfree_qm, nkind
      90             :       INTEGER, POINTER                                   :: itimes
      91             :       LOGICAL                                            :: init
      92             :       REAL(KIND=dp), POINTER                             :: constant
      93             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      94             :       TYPE(simpar_type), POINTER                         :: simpar
      95             : 
      96       44539 :       NULLIFY (itimes, para_env, simpar)
      97             : 
      98       44539 :       CALL zero_md_ener(md_ener, tkind, tshell)
      99             : 
     100             :       CALL get_md_env(md_env=md_env, &
     101             :                       constant=constant, &
     102             :                       itimes=itimes, &
     103             :                       init=init, &
     104             :                       simpar=simpar, &
     105       44539 :                       para_env=para_env)
     106             : 
     107       44539 :       CALL get_part_ke(md_env, md_ener, tkind, tshell, para_env)
     108             : 
     109       44539 :       IF (md_ener%nfree /= 0) THEN
     110       44285 :          md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(simpar%nfree, KIND=dp)
     111       44285 :          md_ener%temp_part = md_ener%temp_part*kelvin
     112             :       END IF
     113             : 
     114       44539 :       nfree_qm = calc_nfree_qm(md_env, md_ener)
     115       44539 :       IF (nfree_qm > 0) THEN
     116        1390 :          md_ener%temp_qm = 2.0_dp*md_ener%ekin_qm/REAL(nfree_qm, KIND=dp)
     117        1390 :          md_ener%temp_qm = md_ener%temp_qm*kelvin
     118             :       END IF
     119             : 
     120       44539 :       IF (md_ener%nfree_shell > 0) THEN
     121        2420 :          md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/REAL(md_ener%nfree_shell, KIND=dp)
     122        2420 :          md_ener%temp_shell = md_ener%temp_shell*kelvin
     123             :       END IF
     124             : 
     125       44539 :       IF (tkind) THEN
     126         902 :          nkind = SIZE(md_ener%temp_kind)
     127        2714 :          DO ikind = 1, nkind
     128             :             md_ener%temp_kind(ikind) = 2.0_dp* &
     129        1812 :                                        md_ener%ekin_kind(ikind)/REAL(md_ener%nfree_kind(ikind), KIND=dp)
     130        2714 :             md_ener%temp_kind(ikind) = md_ener%temp_kind(ikind)*kelvin
     131             :          END DO
     132         902 :          IF (tshell) THEN
     133        1134 :             DO ikind = 1, nkind
     134             :                md_ener%temp_shell_kind(ikind) = 2.0_dp* &
     135         756 :                                                 md_ener%ekin_shell_kind(ikind)/REAL(md_ener%nfree_shell_kind(ikind), KIND=dp)
     136        1134 :                md_ener%temp_shell_kind(ikind) = md_ener%temp_shell_kind(ikind)*kelvin
     137             :             END DO
     138             :          END IF
     139             :       END IF
     140             : 
     141       44539 :       SELECT CASE (simpar%ensemble)
     142             :       CASE DEFAULT
     143           0 :          CPABORT('Unknown ensemble')
     144             :       CASE (isokin_ensemble)
     145           8 :          md_ener%constant = md_ener%ekin
     146             :       CASE (reftraj_ensemble) ! no constant of motion available
     147           0 :          md_ener%constant = md_ener%epot
     148             :       CASE (nve_ensemble)
     149       33855 :          CALL get_econs_nve(md_env, md_ener, para_env)
     150             :       CASE (nvt_ensemble)
     151        7720 :          CALL get_econs_nvt(md_env, md_ener, para_env)
     152             :       CASE (npt_i_ensemble, npt_f_ensemble, npt_ia_ensemble)
     153        2332 :          CALL get_econs_npt(md_env, md_ener, para_env)
     154        2332 :          md_ener%temp_baro = md_ener%temp_baro*kelvin
     155             :       CASE (nph_uniaxial_ensemble)
     156          44 :          CALL get_econs_nph_uniaxial(md_env, md_ener)
     157          44 :          md_ener%temp_baro = md_ener%temp_baro*kelvin
     158             :       CASE (nph_uniaxial_damped_ensemble)
     159          22 :          CALL get_econs_nph_uniaxial(md_env, md_ener)
     160          22 :          md_ener%temp_baro = md_ener%temp_baro*kelvin
     161             :       CASE (langevin_ensemble)
     162         264 :          md_ener%constant = md_ener%ekin + md_ener%epot
     163             :       CASE (npe_f_ensemble, npe_i_ensemble)
     164         294 :          CALL get_econs_npe(md_env, md_ener, para_env)
     165         294 :          md_ener%temp_baro = md_ener%temp_baro*kelvin
     166             :       CASE (nvt_adiabatic_ensemble)
     167       44539 :          CALL get_econs_nvt_adiabatic(md_env, md_ener, para_env)
     168             :       END SELECT
     169             : 
     170       44539 :       IF (init) THEN
     171             :          ! If the value was not read from input let's set it at the begin of the MD
     172        1758 :          IF (constant == 0.0_dp) THEN
     173        1570 :             constant = md_ener%constant
     174        1570 :             CALL set_md_env(md_env=md_env, constant=constant)
     175             :          END IF
     176             :       ELSE
     177       42781 :          CALL get_md_env(md_env=md_env, constant=constant)
     178       42781 :          md_ener%delta_cons = (md_ener%constant - constant)/REAL(natom, KIND=dp)*kelvin
     179             :       END IF
     180             : 
     181       44539 :    END SUBROUTINE compute_conserved_quantity
     182             : 
     183             : ! **************************************************************************************************
     184             : !> \brief Calculates the number of QM degress of freedom
     185             : !> \param md_env ...
     186             : !> \param md_ener ...
     187             : !> \return ...
     188             : ! **************************************************************************************************
     189       89426 :    FUNCTION calc_nfree_qm(md_env, md_ener) RESULT(nfree_qm)
     190             :       TYPE(md_environment_type), POINTER                 :: md_env
     191             :       TYPE(md_ener_type), POINTER                        :: md_ener
     192             :       INTEGER                                            :: nfree_qm
     193             : 
     194             :       INTEGER                                            :: ip
     195       89426 :       INTEGER, POINTER                                   :: cur_indices(:), cur_labels(:)
     196             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     197             :       TYPE(force_env_type), POINTER                      :: force_env
     198             :       TYPE(particle_list_type), POINTER                  :: particles
     199             :       TYPE(qmmm_env_type), POINTER                       :: qmmm_env
     200             :       TYPE(qmmmx_env_type), POINTER                      :: qmmmx_env
     201             :       TYPE(section_vals_type), POINTER                   :: force_env_section
     202             : 
     203       89426 :       NULLIFY (qmmm_env, qmmmx_env, subsys, particles, force_env, force_env_section)
     204       89426 :       nfree_qm = 0
     205             : 
     206       89426 :       CALL get_md_env(md_env, force_env=force_env)
     207             :       CALL force_env_get(force_env, &
     208             :                          subsys=subsys, &
     209             :                          qmmm_env=qmmm_env, &
     210             :                          qmmmx_env=qmmmx_env, &
     211       89426 :                          force_env_section=force_env_section)
     212             : 
     213       89426 :       IF (ASSOCIATED(qmmm_env)) THEN ! conventional QM/MM
     214        2756 :          CALL cp_subsys_get(subsys, particles=particles)
     215             :          ! The degrees of freedom for the quantum part of the system
     216             :          ! are set to 3*Number of QM atoms and to simpar%nfree in case all the MM
     217             :          ! system is treated at QM level (not really QM/MM, just for consistency).
     218             :          ! The degree of freedom will not be correct if 1-3 atoms are treated only
     219             :          ! MM. In this case we should take care of rotations
     220        2756 :          nfree_qm = 3*SIZE(qmmm_env%qm%qm_atom_index)
     221        2756 :          IF (nfree_qm == 3*(particles%n_els)) nfree_qm = md_ener%nfree
     222             :       END IF
     223             : 
     224       89426 :       IF (ASSOCIATED(qmmmx_env)) THEN ! doing force mixing
     225          64 :          CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%INDICES", i_vals=cur_indices)
     226          64 :          CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%LABELS", i_vals=cur_labels)
     227          64 :          nfree_qm = 0
     228        3520 :          DO ip = 1, SIZE(cur_indices)
     229        3520 :             IF (cur_labels(ip) >= force_mixing_label_QM_dynamics) THEN ! this is a QM atom
     230        1002 :                nfree_qm = nfree_qm + 3
     231             :             END IF
     232             :          END DO
     233             :       END IF
     234             : 
     235       89426 :       CPASSERT(.NOT. (ASSOCIATED(qmmm_env) .AND. ASSOCIATED(qmmmx_env)))
     236       89426 :    END FUNCTION calc_nfree_qm
     237             : 
     238             : ! **************************************************************************************************
     239             : !> \brief calculates conserved quantity for nvt ensemble
     240             : !> \param md_env ...
     241             : !> \param md_ener ...
     242             : !> \param para_env ...
     243             : !> \par History
     244             : !>      none
     245             : !> \author gloria
     246             : ! **************************************************************************************************
     247       33855 :    SUBROUTINE get_econs_nve(md_env, md_ener, para_env)
     248             :       TYPE(md_environment_type), POINTER                 :: md_env
     249             :       TYPE(md_ener_type), INTENT(inout)                  :: md_ener
     250             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     251             : 
     252             :       TYPE(force_env_type), POINTER                      :: force_env
     253             :       TYPE(thermostat_type), POINTER                     :: thermostat_coeff, thermostat_shell
     254             : 
     255       33855 :       NULLIFY (force_env, thermostat_coeff, thermostat_shell)
     256             : 
     257             :       CALL get_md_env(md_env, force_env=force_env, thermostat_coeff=thermostat_coeff, &
     258       33855 :                       thermostat_shell=thermostat_shell)
     259       33855 :       md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%ekin_shell
     260             : 
     261             :       CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
     262       33855 :                                    md_ener%thermostat_shell_kin, para_env)
     263       33855 :       md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
     264             : 
     265       33855 :    END SUBROUTINE get_econs_nve
     266             : 
     267             : ! **************************************************************************************************
     268             : !> \brief calculates conserved quantity for nvt ensemble
     269             : !> \param md_env ...
     270             : !> \param md_ener ...
     271             : !> \param para_env ...
     272             : !> \par History
     273             : !>      none
     274             : !> \author gloria
     275             : ! **************************************************************************************************
     276           0 :    SUBROUTINE get_econs_nvt_adiabatic(md_env, md_ener, para_env)
     277             :       TYPE(md_environment_type), POINTER                 :: md_env
     278             :       TYPE(md_ener_type), INTENT(inout)                  :: md_ener
     279             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     280             : 
     281             :       TYPE(force_env_type), POINTER                      :: force_env
     282             :       TYPE(thermostat_type), POINTER                     :: thermostat_fast, thermostat_slow
     283             : 
     284           0 :       NULLIFY (force_env, thermostat_fast, thermostat_slow)
     285             :       CALL get_md_env(md_env, force_env=force_env, thermostat_fast=thermostat_fast, &
     286           0 :                       thermostat_slow=thermostat_slow)
     287             :       CALL get_thermostat_energies(thermostat_fast, md_ener%thermostat_fast_pot, &
     288           0 :                                    md_ener%thermostat_fast_kin, para_env)
     289             :       md_ener%constant = md_ener%ekin + md_ener%epot + &
     290           0 :                          md_ener%thermostat_fast_kin + md_ener%thermostat_fast_pot
     291             :       CALL get_thermostat_energies(thermostat_slow, md_ener%thermostat_slow_pot, &
     292           0 :                                    md_ener%thermostat_slow_kin, para_env)
     293             :       md_ener%constant = md_ener%constant + &
     294           0 :                          md_ener%thermostat_slow_kin + md_ener%thermostat_slow_pot
     295             : 
     296           0 :    END SUBROUTINE get_econs_nvt_adiabatic
     297             : 
     298             : ! **************************************************************************************************
     299             : !> \brief calculates conserved quantity for nvt ensemble
     300             : !> \param md_env ...
     301             : !> \param md_ener ...
     302             : !> \param para_env ...
     303             : !> \par History
     304             : !>      none
     305             : !> \author gloria
     306             : ! **************************************************************************************************
     307        7720 :    SUBROUTINE get_econs_nvt(md_env, md_ener, para_env)
     308             :       TYPE(md_environment_type), POINTER                 :: md_env
     309             :       TYPE(md_ener_type), INTENT(inout)                  :: md_ener
     310             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     311             : 
     312             :       TYPE(force_env_type), POINTER                      :: force_env
     313             :       TYPE(thermostat_type), POINTER                     :: thermostat_coeff, thermostat_part, &
     314             :                                                             thermostat_shell
     315             : 
     316        7720 :       NULLIFY (force_env, thermostat_part, thermostat_coeff, thermostat_shell)
     317             :       CALL get_md_env(md_env, force_env=force_env, thermostat_part=thermostat_part, &
     318        7720 :                       thermostat_coeff=thermostat_coeff, thermostat_shell=thermostat_shell)
     319             :       CALL get_thermostat_energies(thermostat_part, md_ener%thermostat_part_pot, &
     320        7720 :                                    md_ener%thermostat_part_kin, para_env)
     321             :       md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%ekin_shell + &
     322        7720 :                          md_ener%thermostat_part_kin + md_ener%thermostat_part_pot
     323             : 
     324             :       CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
     325        7720 :                                    md_ener%thermostat_shell_kin, para_env)
     326        7720 :       md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
     327             : 
     328        7720 :    END SUBROUTINE get_econs_nvt
     329             : 
     330             : ! **************************************************************************************************
     331             : !> \brief calculates conserved quantity for npe ensemble
     332             : !> \param md_env ...
     333             : !> \param md_ener ...
     334             : !> \param para_env ...
     335             : !> \par History
     336             : !>      none
     337             : !> \author  marcella (02-2008)
     338             : ! **************************************************************************************************
     339         294 :    SUBROUTINE get_econs_npe(md_env, md_ener, para_env)
     340             :       TYPE(md_environment_type), POINTER                 :: md_env
     341             :       TYPE(md_ener_type), INTENT(inout)                  :: md_ener
     342             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     343             : 
     344             :       INTEGER                                            :: nfree
     345             :       TYPE(cell_type), POINTER                           :: box
     346         294 :       TYPE(npt_info_type), POINTER                       :: npt(:, :)
     347             :       TYPE(simpar_type), POINTER                         :: simpar
     348             :       TYPE(thermostat_type), POINTER                     :: thermostat_baro, thermostat_shell
     349             : 
     350         294 :       NULLIFY (thermostat_baro, thermostat_shell, npt)
     351             :       CALL get_md_env(md_env, thermostat_baro=thermostat_baro, &
     352             :                       simpar=simpar, npt=npt, cell=box, &
     353         294 :                       thermostat_shell=thermostat_shell)
     354             :       CALL get_baro_energies(box, simpar, npt, md_ener%baro_kin, &
     355         294 :                              md_ener%baro_pot)
     356         294 :       nfree = SIZE(npt, 1)*SIZE(npt, 2)
     357         294 :       md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/nfree
     358             : 
     359             :       md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%ekin_shell &
     360         294 :                          + md_ener%baro_kin + md_ener%baro_pot
     361             : 
     362             :       CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
     363         294 :                                    md_ener%thermostat_shell_kin, para_env)
     364             :       md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + &
     365         294 :                          md_ener%thermostat_shell_pot
     366             : 
     367         294 :    END SUBROUTINE get_econs_npe
     368             : 
     369             : ! **************************************************************************************************
     370             : !> \brief calculates conserved quantity for npt ensemble
     371             : !> \param md_env ...
     372             : !> \param md_ener ...
     373             : !> \param para_env ...
     374             : !> \par History
     375             : !>      none
     376             : !> \author gloria
     377             : ! **************************************************************************************************
     378        2332 :    SUBROUTINE get_econs_npt(md_env, md_ener, para_env)
     379             :       TYPE(md_environment_type), POINTER                 :: md_env
     380             :       TYPE(md_ener_type), INTENT(inout)                  :: md_ener
     381             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     382             : 
     383             :       INTEGER                                            :: nfree
     384             :       TYPE(cell_type), POINTER                           :: box
     385        2332 :       TYPE(npt_info_type), POINTER                       :: npt(:, :)
     386             :       TYPE(simpar_type), POINTER                         :: simpar
     387             :       TYPE(thermostat_type), POINTER                     :: thermostat_baro, thermostat_part, &
     388             :                                                             thermostat_shell
     389             : 
     390        2332 :       NULLIFY (thermostat_baro, thermostat_part, thermostat_shell, npt, simpar, box)
     391             :       CALL get_md_env(md_env, thermostat_part=thermostat_part, thermostat_baro=thermostat_baro, &
     392        2332 :                       simpar=simpar, npt=npt, cell=box, thermostat_shell=thermostat_shell)
     393             :       CALL get_thermostat_energies(thermostat_part, md_ener%thermostat_part_pot, &
     394        2332 :                                    md_ener%thermostat_part_kin, para_env)
     395             :       CALL get_thermostat_energies(thermostat_baro, md_ener%thermostat_baro_pot, &
     396        2332 :                                    md_ener%thermostat_baro_kin, para_env)
     397        2332 :       CALL get_baro_energies(box, simpar, npt, md_ener%baro_kin, md_ener%baro_pot)
     398        2332 :       nfree = SIZE(npt, 1)*SIZE(npt, 2)
     399        2332 :       md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/nfree
     400             : 
     401             :       md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%ekin_shell &
     402             :                          + md_ener%thermostat_part_kin + md_ener%thermostat_part_pot &
     403             :                          + md_ener%thermostat_baro_kin + md_ener%thermostat_baro_pot &
     404        2332 :                          + md_ener%baro_kin + md_ener%baro_pot
     405             : 
     406             :       CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
     407        2332 :                                    md_ener%thermostat_shell_kin, para_env)
     408        2332 :       md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
     409             : 
     410        2332 :    END SUBROUTINE get_econs_npt
     411             : 
     412             : ! **************************************************************************************************
     413             : !> \brief calculates conserved quantity for nph_uniaxial
     414             : !> \param md_env ...
     415             : !> \param md_ener ...
     416             : !> \par History
     417             : !>      none
     418             : !> \author cjm
     419             : ! **************************************************************************************************
     420          66 :    SUBROUTINE get_econs_nph_uniaxial(md_env, md_ener)
     421             :       TYPE(md_environment_type), POINTER                 :: md_env
     422             :       TYPE(md_ener_type), INTENT(inout)                  :: md_ener
     423             : 
     424             :       TYPE(cell_type), POINTER                           :: box
     425          66 :       TYPE(npt_info_type), POINTER                       :: npt(:, :)
     426             :       TYPE(simpar_type), POINTER                         :: simpar
     427             : 
     428          66 :       CALL get_md_env(md_env, simpar=simpar, npt=npt, cell=box)
     429             : 
     430          66 :       CALL get_baro_energies(box, simpar, npt, md_ener%baro_kin, md_ener%baro_pot)
     431          66 :       md_ener%temp_baro = 2.0_dp*md_ener%baro_kin
     432          66 :       md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%baro_kin + md_ener%baro_pot
     433          66 :    END SUBROUTINE get_econs_nph_uniaxial
     434             : 
     435             : ! **************************************************************************************************
     436             : !> \brief Calculates kinetic energy of particles
     437             : !> \param md_env ...
     438             : !> \param md_ener ...
     439             : !> \param tkind ...
     440             : !> \param tshell ...
     441             : !> \param group ...
     442             : !> \par History
     443             : !>      none
     444             : !> \author CJM
     445             : ! **************************************************************************************************
     446       44539 :    SUBROUTINE get_part_ke(md_env, md_ener, tkind, tshell, group)
     447             :       TYPE(md_environment_type), POINTER                 :: md_env
     448             :       TYPE(md_ener_type), POINTER                        :: md_ener
     449             :       LOGICAL, INTENT(IN)                                :: tkind, tshell
     450             : 
     451             :       CLASS(mp_comm_type), INTENT(IN)                     :: group
     452             : 
     453             :       INTEGER                                            :: i, iparticle, iparticle_kind, &
     454             :                                                             iparticle_local, nparticle_kind, &
     455             :                                                             nparticle_local, shell_index
     456       44539 :       INTEGER, POINTER                                   :: cur_indices(:), cur_labels(:)
     457             :       LOGICAL                                            :: is_shell
     458             :       REAL(KIND=dp)                                      :: ekin_c, ekin_com, ekin_s, mass
     459             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     460       44539 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     461             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     462             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     463             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     464             :       TYPE(force_env_type), POINTER                      :: force_env
     465             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
     466             :                                                             shell_particles
     467       44539 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
     468       44539 :                                                             shell_particle_set
     469             :       TYPE(qmmm_env_type), POINTER                       :: qmmm_env
     470             :       TYPE(qmmmx_env_type), POINTER                      :: qmmmx_env
     471             :       TYPE(section_vals_type), POINTER                   :: force_env_section
     472             :       TYPE(shell_kind_type), POINTER                     :: shell
     473             : 
     474       44539 :       NULLIFY (force_env, qmmm_env, qmmmx_env, subsys, force_env_section)
     475       44539 :       CALL get_md_env(md_env, force_env=force_env)
     476             :       CALL force_env_get(force_env, &
     477             :                          subsys=subsys, &
     478             :                          qmmm_env=qmmm_env, &
     479             :                          qmmmx_env=qmmmx_env, &
     480       44539 :                          force_env_section=force_env_section)
     481             : 
     482             :       CALL cp_subsys_get(subsys=subsys, &
     483             :                          atomic_kinds=atomic_kinds, &
     484             :                          local_particles=local_particles, &
     485             :                          particles=particles, shell_particles=shell_particles, &
     486       44539 :                          core_particles=core_particles)
     487             : 
     488       44539 :       nparticle_kind = atomic_kinds%n_els
     489       44539 :       atomic_kind_set => atomic_kinds%els
     490             : 
     491       44539 :       ekin_s = zero
     492       44539 :       ekin_c = zero
     493       44539 :       ekin_com = zero
     494       44539 :       IF (tkind) THEN
     495        2714 :          md_ener%nfree_kind = 0
     496         902 :          IF (tshell) THEN
     497        1134 :             md_ener%nfree_shell_kind = 0
     498             :          END IF
     499             :       END IF
     500             : 
     501       44539 :       particle_set => particles%els
     502       44539 :       IF (tshell) THEN
     503        2420 :          shell_particle_set => shell_particles%els
     504        2420 :          core_particle_set => core_particles%els
     505        7238 :          DO iparticle_kind = 1, nparticle_kind
     506        4818 :             atomic_kind => atomic_kind_set(iparticle_kind)
     507             :             CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, &
     508        4818 :                                  shell_active=is_shell, shell=shell)
     509        4818 :             nparticle_local = local_particles%n_el(iparticle_kind)
     510        7238 :             IF (is_shell) THEN
     511      124659 :                DO iparticle_local = 1, nparticle_local
     512      119883 :                   iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     513      119883 :                   shell_index = particle_set(iparticle)%shell_index
     514             :                   !ekin
     515             :                   ekin_com = 0.5_dp*mass* &
     516             :                              (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
     517             :                               + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
     518      119883 :                               + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
     519             :                   !vcom
     520      119883 :                   md_ener%vcom(1) = md_ener%vcom(1) + particle_set(iparticle)%v(1)*mass
     521      119883 :                   md_ener%vcom(2) = md_ener%vcom(2) + particle_set(iparticle)%v(2)*mass
     522      119883 :                   md_ener%vcom(3) = md_ener%vcom(3) + particle_set(iparticle)%v(3)*mass
     523      119883 :                   md_ener%total_mass = md_ener%total_mass + mass
     524             : 
     525      119883 :                   md_ener%ekin = md_ener%ekin + ekin_com
     526             :                   ekin_c = 0.5_dp*shell%mass_core* &
     527             :                            (core_particle_set(shell_index)%v(1)*core_particle_set(shell_index)%v(1) &
     528             :                             + core_particle_set(shell_index)%v(2)*core_particle_set(shell_index)%v(2) &
     529      119883 :                             + core_particle_set(shell_index)%v(3)*core_particle_set(shell_index)%v(3))
     530             :                   ekin_s = 0.5_dp*shell%mass_shell* &
     531             :                            (shell_particle_set(shell_index)%v(1)*shell_particle_set(shell_index)%v(1) &
     532             :                             + shell_particle_set(shell_index)%v(2)*shell_particle_set(shell_index)%v(2) &
     533      119883 :                             + shell_particle_set(shell_index)%v(3)*shell_particle_set(shell_index)%v(3))
     534      119883 :                   md_ener%ekin_shell = md_ener%ekin_shell + ekin_c + ekin_s - ekin_com
     535             : 
     536      124659 :                   IF (tkind) THEN
     537       18144 :                      md_ener%ekin_kind(iparticle_kind) = md_ener%ekin_kind(iparticle_kind) + ekin_com
     538       18144 :                      md_ener%nfree_kind(iparticle_kind) = md_ener%nfree_kind(iparticle_kind) + 3
     539             :                      md_ener%ekin_shell_kind(iparticle_kind) = md_ener%ekin_shell_kind(iparticle_kind) + &
     540       18144 :                                                                ekin_c + ekin_s - ekin_com
     541       18144 :                      md_ener%nfree_shell_kind(iparticle_kind) = md_ener%nfree_shell_kind(iparticle_kind) + 3
     542             :                   END IF
     543             : 
     544             :                END DO ! iparticle_local
     545             :             ELSE
     546         714 :                DO iparticle_local = 1, nparticle_local
     547         672 :                   iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     548             :                   ekin_com = 0.5_dp*mass* &
     549             :                              (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
     550             :                               + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
     551         672 :                               + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
     552             :                   !vcom
     553         672 :                   md_ener%vcom(1) = md_ener%vcom(1) + particle_set(iparticle)%v(1)*mass
     554         672 :                   md_ener%vcom(2) = md_ener%vcom(2) + particle_set(iparticle)%v(2)*mass
     555         672 :                   md_ener%vcom(3) = md_ener%vcom(3) + particle_set(iparticle)%v(3)*mass
     556         672 :                   md_ener%total_mass = md_ener%total_mass + mass
     557             : 
     558         672 :                   md_ener%ekin = md_ener%ekin + ekin_com
     559         714 :                   IF (tkind) THEN
     560           0 :                      md_ener%ekin_kind(iparticle_kind) = md_ener%ekin_kind(iparticle_kind) + ekin_com
     561           0 :                      md_ener%nfree_kind(iparticle_kind) = md_ener%nfree_kind(iparticle_kind) + 3
     562             :                   END IF
     563             :                END DO ! iparticle_local
     564             :             END IF
     565             :          END DO ! iparticle_kind
     566        2420 :          IF (tkind) THEN
     567        1890 :             CALL group%sum(md_ener%ekin_kind)
     568        1890 :             CALL group%sum(md_ener%nfree_kind)
     569        1890 :             CALL group%sum(md_ener%ekin_shell_kind)
     570        1890 :             CALL group%sum(md_ener%nfree_shell_kind)
     571             :          END IF
     572             :          ! sum all contributions to energy over calculated parts on all processors
     573        2420 :          CALL group%sum(md_ener%ekin_shell)
     574             :       ELSE
     575      150746 :          DO iparticle_kind = 1, nparticle_kind
     576      108627 :             atomic_kind => atomic_kind_set(iparticle_kind)
     577      108627 :             CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     578      108627 :             nparticle_local = local_particles%n_el(iparticle_kind)
     579     3278996 :             DO iparticle_local = 1, nparticle_local
     580     3128250 :                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     581             :                ! ekin
     582             :                ekin_com = 0.5_dp*mass* &
     583             :                           (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
     584             :                            + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
     585     3128250 :                            + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
     586             : 
     587             :                !vcom
     588     3128250 :                md_ener%vcom(1) = md_ener%vcom(1) + particle_set(iparticle)%v(1)*mass
     589     3128250 :                md_ener%vcom(2) = md_ener%vcom(2) + particle_set(iparticle)%v(2)*mass
     590     3128250 :                md_ener%vcom(3) = md_ener%vcom(3) + particle_set(iparticle)%v(3)*mass
     591     3128250 :                md_ener%total_mass = md_ener%total_mass + mass
     592             : 
     593     3128250 :                md_ener%ekin = md_ener%ekin + ekin_com
     594     3236877 :                IF (tkind) THEN
     595      134462 :                   md_ener%ekin_kind(iparticle_kind) = md_ener%ekin_kind(iparticle_kind) + ekin_com
     596      134462 :                   md_ener%nfree_kind(iparticle_kind) = md_ener%nfree_kind(iparticle_kind) + 3
     597             :                END IF
     598             :             END DO
     599             :          END DO ! iparticle_kind
     600       42119 :          IF (tkind) THEN
     601        2636 :             CALL group%sum(md_ener%ekin_kind)
     602        2636 :             CALL group%sum(md_ener%nfree_kind)
     603             :          END IF
     604             :       END IF
     605             : 
     606             :       ! sum all contributions to energy over calculated parts on all processors
     607       44539 :       CALL group%sum(md_ener%ekin)
     608       44539 :       CALL group%sum(md_ener%vcom)
     609       44539 :       CALL group%sum(md_ener%total_mass)
     610      178156 :       md_ener%vcom = md_ener%vcom/md_ener%total_mass
     611             :       !
     612             :       ! Compute the QM/MM kinetic energy
     613             : 
     614       44539 :       IF (ASSOCIATED(qmmm_env)) THEN ! conventional QM/MM
     615        7952 :          DO i = 1, SIZE(qmmm_env%qm%qm_atom_index)
     616        6574 :             iparticle = qmmm_env%qm%qm_atom_index(i)
     617        6574 :             mass = particle_set(iparticle)%atomic_kind%mass
     618             :             md_ener%ekin_qm = md_ener%ekin_qm + 0.5_dp*mass* &
     619             :                               (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
     620             :                                + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
     621        7952 :                                + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
     622             :          END DO
     623             :       END IF
     624             : 
     625       44539 :       IF (ASSOCIATED(qmmmx_env)) THEN ! doing force mixing
     626          12 :          CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%INDICES", i_vals=cur_indices)
     627          12 :          CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%LABELS", i_vals=cur_labels)
     628        1530 :          DO i = 1, SIZE(cur_indices)
     629        1530 :             IF (cur_labels(i) >= force_mixing_label_QM_dynamics) THEN ! this is a QM atom
     630         318 :                iparticle = cur_indices(i)
     631         318 :                mass = particle_set(iparticle)%atomic_kind%mass
     632             :                md_ener%ekin_qm = md_ener%ekin_qm + 0.5_dp*mass* &
     633             :                                  (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
     634             :                                   + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
     635         318 :                                   + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
     636             :             END IF
     637             :          END DO
     638             :       END IF
     639             : 
     640       44539 :       IF (ASSOCIATED(qmmm_env) .AND. ASSOCIATED(qmmmx_env)) &
     641           0 :          CPABORT("get_part_ke: qmmm bug")
     642       44539 :    END SUBROUTINE get_part_ke
     643             : 
     644             : ! **************************************************************************************************
     645             : 
     646             : END MODULE md_conserved_quantities

Generated by: LCOV version 1.15