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 Routines to handle the virtual site constraint/restraint
      10             : !> \par History
      11             : !>      Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
      12             : !>                                       (patch by Marcel Baer)
      13             : ! **************************************************************************************************
      14             : MODULE constraint_vsite
      15             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      16             :                                               cp_subsys_type
      17             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      18             :    USE force_env_types,                 ONLY: force_env_get,&
      19             :                                               force_env_type
      20             :    USE kinds,                           ONLY: dp
      21             :    USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
      22             :    USE molecule_kind_types,             ONLY: get_molecule_kind,&
      23             :                                               molecule_kind_type,&
      24             :                                               vsite_constraint_type
      25             :    USE molecule_list_types,             ONLY: molecule_list_type
      26             :    USE molecule_types,                  ONLY: get_molecule,&
      27             :                                               global_constraint_type,&
      28             :                                               molecule_type
      29             :    USE particle_list_types,             ONLY: particle_list_type
      30             :    USE particle_types,                  ONLY: particle_type
      31             : #include "./base/base_uses.f90"
      32             : 
      33             :    IMPLICIT NONE
      34             : 
      35             :    PRIVATE
      36             :    PUBLIC :: shake_vsite_int, &
      37             :              shake_vsite_ext, &
      38             :              vsite_force_control
      39             : 
      40             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_vsite'
      41             : 
      42             : CONTAINS
      43             : 
      44             : ! **************************************************************************************************
      45             : !> \brief control force distribution for virtual sites
      46             : !> \param force_env ...
      47             : !> \date 12.2008
      48             : !> \par History
      49             : !>      - none
      50             : !> \author Marcel Baer
      51             : ! **************************************************************************************************
      52       98012 :    SUBROUTINE vsite_force_control(force_env)
      53             :       TYPE(force_env_type), POINTER                      :: force_env
      54             : 
      55             :       INTEGER                                            :: i, ikind, imol, nconstraint, nkind, &
      56             :                                                             nmol_per_kind, nvsitecon
      57             :       LOGICAL                                            :: do_ext_constraint
      58             :       TYPE(cp_subsys_type), POINTER                      :: subsys
      59             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
      60             :       TYPE(global_constraint_type), POINTER              :: gci
      61             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
      62       98012 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
      63             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
      64             :       TYPE(molecule_list_type), POINTER                  :: molecules
      65       98012 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
      66             :       TYPE(molecule_type), POINTER                       :: molecule
      67             :       TYPE(particle_list_type), POINTER                  :: particles
      68             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      69             : 
      70       98012 :       NULLIFY (gci, subsys, local_molecules, local_particles, &
      71       98012 :                molecule_kinds)
      72             : 
      73       98012 :       CALL force_env_get(force_env=force_env, subsys=subsys)
      74             : 
      75             :       CALL cp_subsys_get(subsys=subsys, local_particles=local_particles, &
      76             :                          particles=particles, local_molecules=local_molecules, &
      77       98012 :                          molecule_kinds=molecule_kinds, gci=gci, molecules=molecules)
      78             : 
      79       98012 :       molecule_kind_set => molecule_kinds%els
      80       98012 :       molecule_set => molecules%els
      81       98012 :       particle_set => particles%els
      82       98012 :       nkind = SIZE(molecule_kind_set)
      83             :       ! Intermolecular Virtual Site Constraints
      84       98012 :       do_ext_constraint = .FALSE.
      85       98012 :       IF (ASSOCIATED(gci)) THEN
      86       98012 :          do_ext_constraint = (gci%ntot /= 0)
      87             :       END IF
      88             :       ! Intramolecular Virtual Site Constraints
      89     1830105 :       MOL: DO ikind = 1, nkind
      90     1732093 :          nmol_per_kind = local_molecules%n_el(ikind)
      91     3947905 :          DO imol = 1, nmol_per_kind
      92     2117800 :             i = local_molecules%list(ikind)%array(imol)
      93     2117800 :             molecule => molecule_set(i)
      94     2117800 :             molecule_kind => molecule%molecule_kind
      95     2117800 :             CALL get_molecule_kind(molecule_kind, nconstraint=nconstraint, nvsite=nvsitecon)
      96     2117800 :             IF (nconstraint == 0) CYCLE
      97      199106 :             IF (nvsitecon /= 0) &
      98     3850952 :                CALL force_vsite_int(molecule, particle_set)
      99             :          END DO
     100             :       END DO MOL
     101             :       ! Intermolecular Virtual Site Constraints
     102       98012 :       IF (do_ext_constraint) THEN
     103        1226 :          IF (gci%nvsite /= 0) &
     104           0 :             CALL force_vsite_ext(gci, particle_set)
     105             :       END IF
     106             : 
     107       98012 :    END SUBROUTINE vsite_force_control
     108             : 
     109             : ! **************************************************************************************************
     110             : !> \brief Intramolecular virtual site
     111             : !> \param molecule ...
     112             : !> \param pos ...
     113             : !> \par History
     114             : !>      12.2008 Marcel Baer
     115             : ! **************************************************************************************************
     116         838 :    SUBROUTINE shake_vsite_int(molecule, pos)
     117             :       TYPE(molecule_type), POINTER                       :: molecule
     118             :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :)
     119             : 
     120             :       INTEGER                                            :: first_atom, nvsite
     121             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     122         838 :       TYPE(vsite_constraint_type), POINTER               :: vsite_list(:)
     123             : 
     124         838 :       molecule_kind => molecule%molecule_kind
     125         838 :       CALL get_molecule_kind(molecule_kind, nvsite=nvsite, vsite_list=vsite_list)
     126         838 :       CALL get_molecule(molecule, first_atom=first_atom)
     127             :       ! Real Shake
     128         838 :       CALL shake_vsite_low(vsite_list, nvsite, first_atom, pos)
     129             : 
     130         838 :    END SUBROUTINE shake_vsite_int
     131             : 
     132             : ! **************************************************************************************************
     133             : !> \brief Intramolecular virtual site
     134             : !> \param gci ...
     135             : !> \param pos ...
     136             : !> \par History
     137             : !>      12.2008 Marcel Baer
     138             : ! **************************************************************************************************
     139           0 :    SUBROUTINE shake_vsite_ext(gci, pos)
     140             : 
     141             :       TYPE(global_constraint_type), POINTER              :: gci
     142             :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :)
     143             : 
     144             :       INTEGER                                            :: first_atom, nvsite
     145             :       TYPE(vsite_constraint_type), POINTER               :: vsite_list(:)
     146             : 
     147           0 :       first_atom = 1
     148           0 :       nvsite = gci%nvsite
     149           0 :       vsite_list => gci%vsite_list
     150             :       ! Real Shake
     151           0 :       CALL shake_vsite_low(vsite_list, nvsite, first_atom, pos)
     152             : 
     153           0 :    END SUBROUTINE shake_vsite_ext
     154             : 
     155             : ! **************************************************************************************************
     156             : !> \brief ...
     157             : !> \param vsite_list ...
     158             : !> \param nvsite ...
     159             : !> \param first_atom ...
     160             : !> \param pos ...
     161             : !> \par History
     162             : !>      12.2008 Marcel Bear
     163             : ! **************************************************************************************************
     164         838 :    SUBROUTINE shake_vsite_low(vsite_list, nvsite, first_atom, pos)
     165             :       TYPE(vsite_constraint_type)                        :: vsite_list(:)
     166             :       INTEGER, INTENT(IN)                                :: nvsite, first_atom
     167             :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :)
     168             : 
     169             :       INTEGER                                            :: iconst, index_a, index_b, index_c, &
     170             :                                                             index_d
     171             :       REAL(KIND=dp), DIMENSION(3)                        :: r1, r2
     172             : 
     173        1676 :       DO iconst = 1, nvsite
     174         838 :          IF (vsite_list(iconst)%restraint%active) CYCLE
     175         838 :          index_a = vsite_list(iconst)%a + first_atom - 1
     176         838 :          index_b = vsite_list(iconst)%b + first_atom - 1
     177         838 :          index_c = vsite_list(iconst)%c + first_atom - 1
     178         838 :          index_d = vsite_list(iconst)%d + first_atom - 1
     179             : 
     180        3352 :          r1(:) = pos(:, index_b) - pos(:, index_c)
     181        3352 :          r2(:) = pos(:, index_d) - pos(:, index_c)
     182             :          pos(:, index_a) = pos(:, index_c) + vsite_list(iconst)%wbc*r1(:) + &
     183        4190 :                            vsite_list(iconst)%wdc*r2(:)
     184             :       END DO
     185         838 :    END SUBROUTINE shake_vsite_low
     186             : 
     187             : ! **************************************************************************************************
     188             : !> \brief Intramolecular virtual site
     189             : !> \param molecule ...
     190             : !> \param particle_set ...
     191             : !> \par History
     192             : !>      12.2008 Marcel Bear
     193             : ! **************************************************************************************************
     194        1059 :    SUBROUTINE force_vsite_int(molecule, particle_set)
     195             :       TYPE(molecule_type), POINTER                       :: molecule
     196             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     197             : 
     198             :       INTEGER                                            :: first_atom, iconst, index_a, index_b, &
     199             :                                                             index_c, index_d, nvsite
     200             :       REAL(KIND=dp)                                      :: wb, wc, wd
     201             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     202        1059 :       TYPE(vsite_constraint_type), POINTER               :: vsite_list(:)
     203             : 
     204        1059 :       molecule_kind => molecule%molecule_kind
     205        1059 :       CALL get_molecule_kind(molecule_kind, nvsite=nvsite, vsite_list=vsite_list)
     206        1059 :       CALL get_molecule(molecule, first_atom=first_atom)
     207             : 
     208        2118 :       DO iconst = 1, nvsite
     209        1059 :          IF (vsite_list(iconst)%restraint%active) CYCLE
     210        1059 :          index_a = vsite_list(iconst)%a + first_atom - 1
     211        1059 :          index_b = vsite_list(iconst)%b + first_atom - 1
     212        1059 :          index_c = vsite_list(iconst)%c + first_atom - 1
     213        1059 :          index_d = vsite_list(iconst)%d + first_atom - 1
     214             : 
     215        1059 :          wb = vsite_list(iconst)%wbc
     216        1059 :          wd = vsite_list(iconst)%wdc
     217        1059 :          wc = 1.0_dp - vsite_list(iconst)%wbc - vsite_list(iconst)%wdc
     218             : 
     219        4236 :          particle_set(index_b)%f(:) = particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
     220        4236 :          particle_set(index_c)%f(:) = particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
     221        4236 :          particle_set(index_d)%f(:) = particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
     222        5295 :          particle_set(index_a)%f(:) = 0.0_dp
     223             :       END DO
     224             : 
     225        1059 :    END SUBROUTINE force_vsite_int
     226             : 
     227             : ! **************************************************************************************************
     228             : !> \brief Intramolecular virtual site
     229             : !> \param gci ...
     230             : !> \param particle_set ...
     231             : !> \par History
     232             : !>      12.2008 Marcel Bear
     233             : ! **************************************************************************************************
     234           0 :    SUBROUTINE force_vsite_ext(gci, particle_set)
     235             :       TYPE(global_constraint_type), POINTER              :: gci
     236             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     237             : 
     238             :       INTEGER                                            :: first_atom, iconst, index_a, index_b, &
     239             :                                                             index_c, index_d, nvsite
     240             :       REAL(KIND=dp)                                      :: wb, wc, wd
     241           0 :       TYPE(vsite_constraint_type), POINTER               :: vsite_list(:)
     242             : 
     243           0 :       first_atom = 1
     244           0 :       nvsite = gci%nvsite
     245           0 :       vsite_list => gci%vsite_list
     246             :       ! Real Shake
     247             : 
     248           0 :       DO iconst = 1, nvsite
     249           0 :          IF (vsite_list(iconst)%restraint%active) CYCLE
     250           0 :          index_a = vsite_list(iconst)%a + first_atom - 1
     251           0 :          index_b = vsite_list(iconst)%b + first_atom - 1
     252           0 :          index_c = vsite_list(iconst)%c + first_atom - 1
     253           0 :          index_d = vsite_list(iconst)%d + first_atom - 1
     254             : 
     255           0 :          wb = vsite_list(iconst)%wbc
     256           0 :          wd = vsite_list(iconst)%wdc
     257           0 :          wc = 1.0_dp - vsite_list(iconst)%wbc - vsite_list(iconst)%wdc
     258             : 
     259           0 :          particle_set(index_b)%f(:) = particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
     260           0 :          particle_set(index_c)%f(:) = particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
     261           0 :          particle_set(index_d)%f(:) = particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
     262           0 :          particle_set(index_a)%f(:) = 0.0_dp
     263             :       END DO
     264           0 :    END SUBROUTINE force_vsite_ext
     265             : 
     266             : END MODULE constraint_vsite