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