(git:ed6f26b)
Loading...
Searching...
No Matches
constraint_vsite.F
Go to the documentation of this file.
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 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! **************************************************************************************************
14MODULE 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
42CONTAINS
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
266END 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, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
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
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
molecule_kind_list_types::molecule_kind_list_type
represent a list of objects
Definition molecule_kind_list_types.F:45
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_kind_types::vsite_constraint_type
Definition molecule_kind_types.F:142
molecule_list_types::molecule_list_type
represent a list of objects
Definition molecule_list_types.F:46
molecule_types::global_constraint_type
Definition molecule_types.F:93
molecule_types::molecule_type
Definition molecule_types.F:110
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35