(git:e5fdd81)
constraint_vsite.F
Go to the documentation of this file.
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  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  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  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  NULLIFY (gci, subsys, local_molecules, local_particles, &
71  molecule_kinds)
72 
73  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  molecule_kinds=molecule_kinds, gci=gci, molecules=molecules)
78 
79  molecule_kind_set => molecule_kinds%els
80  molecule_set => molecules%els
81  particle_set => particles%els
82  nkind = SIZE(molecule_kind_set)
83  ! Intermolecular Virtual Site Constraints
84  do_ext_constraint = .false.
85  IF (ASSOCIATED(gci)) THEN
86  do_ext_constraint = (gci%ntot /= 0)
87  END IF
88  ! Intramolecular Virtual Site Constraints
89  mol: DO ikind = 1, nkind
90  nmol_per_kind = local_molecules%n_el(ikind)
91  DO imol = 1, nmol_per_kind
92  i = local_molecules%list(ikind)%array(imol)
93  molecule => molecule_set(i)
94  molecule_kind => molecule%molecule_kind
95  CALL get_molecule_kind(molecule_kind, nconstraint=nconstraint, nvsite=nvsitecon)
96  IF (nconstraint == 0) cycle
97  IF (nvsitecon /= 0) &
98  CALL force_vsite_int(molecule, particle_set)
99  END DO
100  END DO mol
101  ! Intermolecular Virtual Site Constraints
102  IF (do_ext_constraint) THEN
103  IF (gci%nvsite /= 0) &
104  CALL force_vsite_ext(gci, particle_set)
105  END IF
106 
107  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  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  TYPE(vsite_constraint_type), POINTER :: vsite_list(:)
123 
124  molecule_kind => molecule%molecule_kind
125  CALL get_molecule_kind(molecule_kind, nvsite=nvsite, vsite_list=vsite_list)
126  CALL get_molecule(molecule, first_atom=first_atom)
127  ! Real Shake
128  CALL shake_vsite_low(vsite_list, nvsite, first_atom, pos)
129 
130  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  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  first_atom = 1
148  nvsite = gci%nvsite
149  vsite_list => gci%vsite_list
150  ! Real Shake
151  CALL shake_vsite_low(vsite_list, nvsite, first_atom, pos)
152 
153  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  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  DO iconst = 1, nvsite
174  IF (vsite_list(iconst)%restraint%active) cycle
175  index_a = vsite_list(iconst)%a + first_atom - 1
176  index_b = vsite_list(iconst)%b + first_atom - 1
177  index_c = vsite_list(iconst)%c + first_atom - 1
178  index_d = vsite_list(iconst)%d + first_atom - 1
179 
180  r1(:) = pos(:, index_b) - pos(:, index_c)
181  r2(:) = pos(:, index_d) - pos(:, index_c)
182  pos(:, index_a) = pos(:, index_c) + vsite_list(iconst)%wbc*r1(:) + &
183  vsite_list(iconst)%wdc*r2(:)
184  END DO
185  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  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  TYPE(vsite_constraint_type), POINTER :: vsite_list(:)
203 
204  molecule_kind => molecule%molecule_kind
205  CALL get_molecule_kind(molecule_kind, nvsite=nvsite, vsite_list=vsite_list)
206  CALL get_molecule(molecule, first_atom=first_atom)
207 
208  DO iconst = 1, nvsite
209  IF (vsite_list(iconst)%restraint%active) cycle
210  index_a = vsite_list(iconst)%a + first_atom - 1
211  index_b = vsite_list(iconst)%b + first_atom - 1
212  index_c = vsite_list(iconst)%c + first_atom - 1
213  index_d = vsite_list(iconst)%d + first_atom - 1
214 
215  wb = vsite_list(iconst)%wbc
216  wd = vsite_list(iconst)%wdc
217  wc = 1.0_dp - vsite_list(iconst)%wbc - vsite_list(iconst)%wdc
218 
219  particle_set(index_b)%f(:) = particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
220  particle_set(index_c)%f(:) = particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
221  particle_set(index_d)%f(:) = particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
222  particle_set(index_a)%f(:) = 0.0_dp
223  END DO
224 
225  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  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  TYPE(vsite_constraint_type), POINTER :: vsite_list(:)
242 
243  first_atom = 1
244  nvsite = gci%nvsite
245  vsite_list => gci%vsite_list
246  ! Real Shake
247 
248  DO iconst = 1, nvsite
249  IF (vsite_list(iconst)%restraint%active) cycle
250  index_a = vsite_list(iconst)%a + first_atom - 1
251  index_b = vsite_list(iconst)%b + first_atom - 1
252  index_c = vsite_list(iconst)%c + first_atom - 1
253  index_d = vsite_list(iconst)%d + first_atom - 1
254 
255  wb = vsite_list(iconst)%wbc
256  wd = vsite_list(iconst)%wdc
257  wc = 1.0_dp - vsite_list(iconst)%wbc - vsite_list(iconst)%wdc
258 
259  particle_set(index_b)%f(:) = particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
260  particle_set(index_c)%f(:) = particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
261  particle_set(index_d)%f(:) = particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
262  particle_set(index_a)%f(:) = 0.0_dp
263  END DO
264  END SUBROUTINE force_vsite_ext
265 
266 END MODULE constraint_vsite
constraint_vsite
Routines to handle the virtual site constraint/restraint.
Definition: constraint_vsite.F:14
constraint_vsite::shake_vsite_ext
subroutine, public shake_vsite_ext(gci, pos)
Intramolecular virtual site.
Definition: constraint_vsite.F:140
constraint_vsite::shake_vsite_int
subroutine, public shake_vsite_int(molecule, pos)
Intramolecular virtual site.
Definition: constraint_vsite.F:117
constraint_vsite::vsite_force_control
subroutine, public vsite_force_control(force_env)
control force distribution for virtual sites
Definition: constraint_vsite.F:53
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition: cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition: cp_subsys_types.F:345
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition: distribution_1d_types.F:24
force_env_types
Interface for the force calculations.
Definition: force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
Definition: force_env_types.F:329
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
molecule_kind_list_types
represent a simple array based list of the given type
Definition: molecule_kind_list_types.F:15
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition: molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition: molecule_kind_types.F:510
molecule_list_types
represent a simple array based list of the given type
Definition: molecule_list_types.F:15
molecule_types
Define the data structure for the molecule information.
Definition: molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition: molecule_types.F:277
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19