(git:71c3ab0)
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-2026 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) THEN
98 CALL force_vsite_int(molecule, particle_set)
99 END IF
100 END DO
101 END DO mol
102 ! Intermolecular Virtual Site Constraints
103 IF (do_ext_constraint) THEN
104 IF (gci%nvsite /= 0) THEN
105 CALL force_vsite_ext(gci, particle_set)
106 END IF
107 END IF
108
109 END SUBROUTINE vsite_force_control
110
111! **************************************************************************************************
112!> \brief Intramolecular virtual site
113!> \param molecule ...
114!> \param pos ...
115!> \par History
116!> 12.2008 Marcel Baer
117! **************************************************************************************************
118 SUBROUTINE shake_vsite_int(molecule, pos)
119 TYPE(molecule_type), POINTER :: molecule
120 REAL(kind=dp), INTENT(INOUT) :: pos(:, :)
121
122 INTEGER :: first_atom, nvsite
123 TYPE(molecule_kind_type), POINTER :: molecule_kind
124 TYPE(vsite_constraint_type), POINTER :: vsite_list(:)
125
126 molecule_kind => molecule%molecule_kind
127 CALL get_molecule_kind(molecule_kind, nvsite=nvsite, vsite_list=vsite_list)
128 CALL get_molecule(molecule, first_atom=first_atom)
129 ! Real Shake
130 CALL shake_vsite_low(vsite_list, nvsite, first_atom, pos)
131
132 END SUBROUTINE shake_vsite_int
133
134! **************************************************************************************************
135!> \brief Intramolecular virtual site
136!> \param gci ...
137!> \param pos ...
138!> \par History
139!> 12.2008 Marcel Baer
140! **************************************************************************************************
141 SUBROUTINE shake_vsite_ext(gci, pos)
142
143 TYPE(global_constraint_type), POINTER :: gci
144 REAL(kind=dp), INTENT(INOUT) :: pos(:, :)
145
146 INTEGER :: first_atom, nvsite
147 TYPE(vsite_constraint_type), POINTER :: vsite_list(:)
148
149 first_atom = 1
150 nvsite = gci%nvsite
151 vsite_list => gci%vsite_list
152 ! Real Shake
153 CALL shake_vsite_low(vsite_list, nvsite, first_atom, pos)
154
155 END SUBROUTINE shake_vsite_ext
156
157! **************************************************************************************************
158!> \brief ...
159!> \param vsite_list ...
160!> \param nvsite ...
161!> \param first_atom ...
162!> \param pos ...
163!> \par History
164!> 12.2008 Marcel Bear
165! **************************************************************************************************
166 SUBROUTINE shake_vsite_low(vsite_list, nvsite, first_atom, pos)
167 TYPE(vsite_constraint_type) :: vsite_list(:)
168 INTEGER, INTENT(IN) :: nvsite, first_atom
169 REAL(kind=dp), INTENT(INOUT) :: pos(:, :)
170
171 INTEGER :: iconst, index_a, index_b, index_c, &
172 index_d
173 REAL(kind=dp), DIMENSION(3) :: r1, r2
174
175 DO iconst = 1, nvsite
176 IF (vsite_list(iconst)%restraint%active) cycle
177 index_a = vsite_list(iconst)%a + first_atom - 1
178 index_b = vsite_list(iconst)%b + first_atom - 1
179 index_c = vsite_list(iconst)%c + first_atom - 1
180 index_d = vsite_list(iconst)%d + first_atom - 1
181
182 r1(:) = pos(:, index_b) - pos(:, index_c)
183 r2(:) = pos(:, index_d) - pos(:, index_c)
184 pos(:, index_a) = pos(:, index_c) + vsite_list(iconst)%wbc*r1(:) + &
185 vsite_list(iconst)%wdc*r2(:)
186 END DO
187 END SUBROUTINE shake_vsite_low
188
189! **************************************************************************************************
190!> \brief Intramolecular virtual site
191!> \param molecule ...
192!> \param particle_set ...
193!> \par History
194!> 12.2008 Marcel Bear
195! **************************************************************************************************
196 SUBROUTINE force_vsite_int(molecule, particle_set)
197 TYPE(molecule_type), POINTER :: molecule
198 TYPE(particle_type), POINTER :: particle_set(:)
199
200 INTEGER :: first_atom, iconst, index_a, index_b, &
201 index_c, index_d, nvsite
202 REAL(kind=dp) :: wb, wc, wd
203 TYPE(molecule_kind_type), POINTER :: molecule_kind
204 TYPE(vsite_constraint_type), POINTER :: vsite_list(:)
205
206 molecule_kind => molecule%molecule_kind
207 CALL get_molecule_kind(molecule_kind, nvsite=nvsite, vsite_list=vsite_list)
208 CALL get_molecule(molecule, first_atom=first_atom)
209
210 DO iconst = 1, nvsite
211 IF (vsite_list(iconst)%restraint%active) cycle
212 index_a = vsite_list(iconst)%a + first_atom - 1
213 index_b = vsite_list(iconst)%b + first_atom - 1
214 index_c = vsite_list(iconst)%c + first_atom - 1
215 index_d = vsite_list(iconst)%d + first_atom - 1
216
217 wb = vsite_list(iconst)%wbc
218 wd = vsite_list(iconst)%wdc
219 wc = 1.0_dp - vsite_list(iconst)%wbc - vsite_list(iconst)%wdc
220
221 particle_set(index_b)%f(:) = particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
222 particle_set(index_c)%f(:) = particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
223 particle_set(index_d)%f(:) = particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
224 particle_set(index_a)%f(:) = 0.0_dp
225 END DO
226
227 END SUBROUTINE force_vsite_int
228
229! **************************************************************************************************
230!> \brief Intramolecular virtual site
231!> \param gci ...
232!> \param particle_set ...
233!> \par History
234!> 12.2008 Marcel Bear
235! **************************************************************************************************
236 SUBROUTINE force_vsite_ext(gci, particle_set)
237 TYPE(global_constraint_type), POINTER :: gci
238 TYPE(particle_type), POINTER :: particle_set(:)
239
240 INTEGER :: first_atom, iconst, index_a, index_b, &
241 index_c, index_d, nvsite
242 REAL(kind=dp) :: wb, wc, wd
243 TYPE(vsite_constraint_type), POINTER :: vsite_list(:)
244
245 first_atom = 1
246 nvsite = gci%nvsite
247 vsite_list => gci%vsite_list
248 ! Real Shake
249
250 DO iconst = 1, nvsite
251 IF (vsite_list(iconst)%restraint%active) cycle
252 index_a = vsite_list(iconst)%a + first_atom - 1
253 index_b = vsite_list(iconst)%b + first_atom - 1
254 index_c = vsite_list(iconst)%c + first_atom - 1
255 index_d = vsite_list(iconst)%d + first_atom - 1
256
257 wb = vsite_list(iconst)%wbc
258 wd = vsite_list(iconst)%wdc
259 wc = 1.0_dp - vsite_list(iconst)%wbc - vsite_list(iconst)%wdc
260
261 particle_set(index_b)%f(:) = particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
262 particle_set(index_c)%f(:) = particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
263 particle_set(index_d)%f(:) = particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
264 particle_set(index_a)%f(:) = 0.0_dp
265 END DO
266 END SUBROUTINE force_vsite_ext
267
268END 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:142
constraint_vsite::shake_vsite_int
subroutine, public shake_vsite_int(molecule, pos)
Intramolecular virtual site.
Definition constraint_vsite.F:119
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, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
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:344
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:513
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:332
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:151
molecule_kind_types::vsite_constraint_type
Definition molecule_kind_types.F:134
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:134
molecule_types::molecule_type
Definition molecule_types.F:159
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