LCOV - code coverage report
Current view: top level - src/motion - md_conserved_quantities.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 94.1 % 219 206
Test Date: 2025-07-25 12:55:17 Functions: 88.9 % 9 8

            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 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        42043 :    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        42043 :       NULLIFY (itimes, para_env, simpar)
      97              : 
      98        42043 :       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        42043 :                       para_env=para_env)
     106              : 
     107        42043 :       CALL get_part_ke(md_env, md_ener, tkind, tshell, para_env)
     108              : 
     109        42043 :       IF (md_ener%nfree /= 0) THEN
     110        41791 :          md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(simpar%nfree, KIND=dp)
     111        41791 :          md_ener%temp_part = md_ener%temp_part*kelvin
     112              :       END IF
     113              : 
     114        42043 :       nfree_qm = calc_nfree_qm(md_env, md_ener)
     115        42043 :       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        42043 :       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        42043 :       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        42043 :       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        31143 :          CALL get_econs_nve(md_env, md_ener, para_env)
     150              :       CASE (nvt_ensemble)
     151         7936 :          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        42043 :          CALL get_econs_nvt_adiabatic(md_env, md_ener, para_env)
     168              :       END SELECT
     169              : 
     170        42043 :       IF (init) THEN
     171              :          ! If the value was not read from input let's set it at the begin of the MD
     172         1731 :          IF (constant == 0.0_dp) THEN
     173         1543 :             constant = md_ener%constant
     174         1543 :             CALL set_md_env(md_env=md_env, constant=constant)
     175              :          END IF
     176              :       ELSE
     177        40312 :          CALL get_md_env(md_env=md_env, constant=constant)
     178        40312 :          md_ener%delta_cons = (md_ener%constant - constant)/REAL(natom, KIND=dp)*kelvin
     179              :       END IF
     180              : 
     181        42043 :    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        84399 :    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        84399 :       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        84399 :       NULLIFY (qmmm_env, qmmmx_env, subsys, particles, force_env, force_env_section)
     204        84399 :       nfree_qm = 0
     205              : 
     206        84399 :       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        84399 :                          force_env_section=force_env_section)
     212              : 
     213        84399 :       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        84399 :       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        84399 :       CPASSERT(.NOT. (ASSOCIATED(qmmm_env) .AND. ASSOCIATED(qmmmx_env)))
     236        84399 :    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        31143 :    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        31143 :       NULLIFY (force_env, thermostat_coeff, thermostat_shell)
     256              : 
     257              :       CALL get_md_env(md_env, force_env=force_env, thermostat_coeff=thermostat_coeff, &
     258        31143 :                       thermostat_shell=thermostat_shell)
     259        31143 :       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        31143 :                                    md_ener%thermostat_shell_kin, para_env)
     263        31143 :       md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
     264              : 
     265        31143 :    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         7936 :    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         7936 :       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         7936 :                       thermostat_coeff=thermostat_coeff, thermostat_shell=thermostat_shell)
     319              :       CALL get_thermostat_energies(thermostat_part, md_ener%thermostat_part_pot, &
     320         7936 :                                    md_ener%thermostat_part_kin, para_env)
     321              :       md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%ekin_shell + &
     322         7936 :                          md_ener%thermostat_part_kin + md_ener%thermostat_part_pot
     323              : 
     324              :       CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
     325         7936 :                                    md_ener%thermostat_shell_kin, para_env)
     326         7936 :       md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
     327              : 
     328         7936 :    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        42043 :    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        42043 :       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        42043 :       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        42043 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
     468        42043 :                                                             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        42043 :       NULLIFY (force_env, qmmm_env, qmmmx_env, subsys, force_env_section)
     475        42043 :       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        42043 :                          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        42043 :                          core_particles=core_particles)
     487              : 
     488        42043 :       nparticle_kind = atomic_kinds%n_els
     489        42043 :       atomic_kind_set => atomic_kinds%els
     490              : 
     491        42043 :       ekin_s = zero
     492        42043 :       ekin_c = zero
     493        42043 :       ekin_com = zero
     494        42043 :       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        42043 :       particle_set => particles%els
     502        42043 :       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       146194 :          DO iparticle_kind = 1, nparticle_kind
     576       106571 :             atomic_kind => atomic_kind_set(iparticle_kind)
     577       106571 :             CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     578       106571 :             nparticle_local = local_particles%n_el(iparticle_kind)
     579      3264622 :             DO iparticle_local = 1, nparticle_local
     580      3118428 :                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      3118428 :                            + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
     586              : 
     587              :                !vcom
     588      3118428 :                md_ener%vcom(1) = md_ener%vcom(1) + particle_set(iparticle)%v(1)*mass
     589      3118428 :                md_ener%vcom(2) = md_ener%vcom(2) + particle_set(iparticle)%v(2)*mass
     590      3118428 :                md_ener%vcom(3) = md_ener%vcom(3) + particle_set(iparticle)%v(3)*mass
     591      3118428 :                md_ener%total_mass = md_ener%total_mass + mass
     592              : 
     593      3118428 :                md_ener%ekin = md_ener%ekin + ekin_com
     594      3224999 :                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        39623 :          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        42043 :       CALL group%sum(md_ener%ekin)
     608        42043 :       CALL group%sum(md_ener%vcom)
     609        42043 :       CALL group%sum(md_ener%total_mass)
     610       168172 :       md_ener%vcom = md_ener%vcom/md_ener%total_mass
     611              :       !
     612              :       ! Compute the QM/MM kinetic energy
     613              : 
     614        42043 :       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        42043 :       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        42043 :       IF (ASSOCIATED(qmmm_env) .AND. ASSOCIATED(qmmmx_env)) &
     641            0 :          CPABORT("get_part_ke: qmmm bug")
     642        42043 :    END SUBROUTINE get_part_ke
     643              : 
     644              : ! **************************************************************************************************
     645              : 
     646              : END MODULE md_conserved_quantities
        

Generated by: LCOV version 2.0-1