Line data Source code
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 : ! **************************************************************************************************
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 102231 : 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 102231 : 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 102231 : 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 102231 : NULLIFY (gci, subsys, local_molecules, local_particles, &
71 102231 : molecule_kinds)
72 :
73 102231 : 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 102231 : molecule_kinds=molecule_kinds, gci=gci, molecules=molecules)
78 :
79 102231 : molecule_kind_set => molecule_kinds%els
80 102231 : molecule_set => molecules%els
81 102231 : particle_set => particles%els
82 102231 : nkind = SIZE(molecule_kind_set)
83 : ! Intermolecular Virtual Site Constraints
84 102231 : do_ext_constraint = .FALSE.
85 102231 : IF (ASSOCIATED(gci)) THEN
86 102231 : do_ext_constraint = (gci%ntot /= 0)
87 : END IF
88 : ! Intramolecular Virtual Site Constraints
89 1757216 : MOL: DO ikind = 1, nkind
90 1654985 : nmol_per_kind = local_molecules%n_el(ikind)
91 3832923 : DO imol = 1, nmol_per_kind
92 2075707 : i = local_molecules%list(ikind)%array(imol)
93 2075707 : molecule => molecule_set(i)
94 2075707 : molecule_kind => molecule%molecule_kind
95 2075707 : CALL get_molecule_kind(molecule_kind, nconstraint=nconstraint, nvsite=nvsitecon)
96 2075707 : IF (nconstraint == 0) CYCLE
97 3929798 : IF (nvsitecon /= 0) THEN
98 1059 : CALL force_vsite_int(molecule, particle_set)
99 : END IF
100 : END DO
101 : END DO MOL
102 : ! Intermolecular Virtual Site Constraints
103 102231 : IF (do_ext_constraint) THEN
104 1226 : IF (gci%nvsite /= 0) THEN
105 0 : CALL force_vsite_ext(gci, particle_set)
106 : END IF
107 : END IF
108 :
109 102231 : 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 838 : 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 838 : TYPE(vsite_constraint_type), POINTER :: vsite_list(:)
125 :
126 838 : molecule_kind => molecule%molecule_kind
127 838 : CALL get_molecule_kind(molecule_kind, nvsite=nvsite, vsite_list=vsite_list)
128 838 : CALL get_molecule(molecule, first_atom=first_atom)
129 : ! Real Shake
130 838 : CALL shake_vsite_low(vsite_list, nvsite, first_atom, pos)
131 :
132 838 : 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 0 : 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 0 : first_atom = 1
150 0 : nvsite = gci%nvsite
151 0 : vsite_list => gci%vsite_list
152 : ! Real Shake
153 0 : CALL shake_vsite_low(vsite_list, nvsite, first_atom, pos)
154 :
155 0 : 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 838 : 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 1676 : DO iconst = 1, nvsite
176 838 : IF (vsite_list(iconst)%restraint%active) CYCLE
177 838 : index_a = vsite_list(iconst)%a + first_atom - 1
178 838 : index_b = vsite_list(iconst)%b + first_atom - 1
179 838 : index_c = vsite_list(iconst)%c + first_atom - 1
180 838 : index_d = vsite_list(iconst)%d + first_atom - 1
181 :
182 3352 : r1(:) = pos(:, index_b) - pos(:, index_c)
183 3352 : r2(:) = pos(:, index_d) - pos(:, index_c)
184 : pos(:, index_a) = pos(:, index_c) + vsite_list(iconst)%wbc*r1(:) + &
185 4190 : vsite_list(iconst)%wdc*r2(:)
186 : END DO
187 838 : 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 1059 : 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 1059 : TYPE(vsite_constraint_type), POINTER :: vsite_list(:)
205 :
206 1059 : molecule_kind => molecule%molecule_kind
207 1059 : CALL get_molecule_kind(molecule_kind, nvsite=nvsite, vsite_list=vsite_list)
208 1059 : CALL get_molecule(molecule, first_atom=first_atom)
209 :
210 2118 : DO iconst = 1, nvsite
211 1059 : IF (vsite_list(iconst)%restraint%active) CYCLE
212 1059 : index_a = vsite_list(iconst)%a + first_atom - 1
213 1059 : index_b = vsite_list(iconst)%b + first_atom - 1
214 1059 : index_c = vsite_list(iconst)%c + first_atom - 1
215 1059 : index_d = vsite_list(iconst)%d + first_atom - 1
216 :
217 1059 : wb = vsite_list(iconst)%wbc
218 1059 : wd = vsite_list(iconst)%wdc
219 1059 : wc = 1.0_dp - vsite_list(iconst)%wbc - vsite_list(iconst)%wdc
220 :
221 4236 : particle_set(index_b)%f(:) = particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
222 4236 : particle_set(index_c)%f(:) = particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
223 4236 : particle_set(index_d)%f(:) = particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
224 5295 : particle_set(index_a)%f(:) = 0.0_dp
225 : END DO
226 :
227 1059 : 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 0 : 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 0 : TYPE(vsite_constraint_type), POINTER :: vsite_list(:)
244 :
245 0 : first_atom = 1
246 0 : nvsite = gci%nvsite
247 0 : vsite_list => gci%vsite_list
248 : ! Real Shake
249 :
250 0 : DO iconst = 1, nvsite
251 0 : IF (vsite_list(iconst)%restraint%active) CYCLE
252 0 : index_a = vsite_list(iconst)%a + first_atom - 1
253 0 : index_b = vsite_list(iconst)%b + first_atom - 1
254 0 : index_c = vsite_list(iconst)%c + first_atom - 1
255 0 : index_d = vsite_list(iconst)%d + first_atom - 1
256 :
257 0 : wb = vsite_list(iconst)%wbc
258 0 : wd = vsite_list(iconst)%wdc
259 0 : wc = 1.0_dp - vsite_list(iconst)%wbc - vsite_list(iconst)%wdc
260 :
261 0 : particle_set(index_b)%f(:) = particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
262 0 : particle_set(index_c)%f(:) = particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
263 0 : particle_set(index_d)%f(:) = particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
264 0 : particle_set(index_a)%f(:) = 0.0_dp
265 : END DO
266 0 : END SUBROUTINE force_vsite_ext
267 :
268 : END MODULE constraint_vsite
|
