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