1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Routines to handle the virtual site constraint/restraint
8!> \par History
9!>      Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
10!>                                       (patch by Marcel Baer)
11! **************************************************************************************************
12MODULE constraint_vsite
13   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
14                                              cp_subsys_type
15   USE distribution_1d_types,           ONLY: distribution_1d_type
16   USE force_env_types,                 ONLY: force_env_get,&
17                                              force_env_type
18   USE kinds,                           ONLY: dp
19   USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
20   USE molecule_kind_types,             ONLY: get_molecule_kind,&
21                                              molecule_kind_type,&
22                                              vsite_constraint_type
23   USE molecule_list_types,             ONLY: molecule_list_type
24   USE molecule_types,                  ONLY: get_molecule,&
25                                              global_constraint_type,&
26                                              molecule_type
27   USE particle_list_types,             ONLY: particle_list_type
28   USE particle_types,                  ONLY: particle_type
29#include "./base/base_uses.f90"
30
31   IMPLICIT NONE
32
33   PRIVATE
34   PUBLIC :: shake_vsite_int, &
35             shake_vsite_ext, &
36             vsite_force_control
37
38   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_vsite'
39
40CONTAINS
41
42! **************************************************************************************************
43!> \brief control force distribution for virtual sites
44!> \param force_env ...
45!> \date 12.2008
46!> \par History
47!>      - none
48!> \author Marcel Baer
49! **************************************************************************************************
50   SUBROUTINE vsite_force_control(force_env)
51      TYPE(force_env_type), POINTER                      :: force_env
52
53      INTEGER                                            :: i, ikind, imol, nconstraint, nkind, &
54                                                            nmol_per_kind, nvsitecon
55      LOGICAL                                            :: do_ext_constraint
56      TYPE(cp_subsys_type), POINTER                      :: subsys
57      TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
58      TYPE(global_constraint_type), POINTER              :: gci
59      TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
60      TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
61      TYPE(molecule_kind_type), POINTER                  :: molecule_kind
62      TYPE(molecule_list_type), POINTER                  :: molecules
63      TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
64      TYPE(molecule_type), POINTER                       :: molecule
65      TYPE(particle_list_type), POINTER                  :: particles
66      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
67
68      NULLIFY (gci, subsys, local_molecules, local_particles, &
69               molecule_kinds)
70
71      CALL force_env_get(force_env=force_env, subsys=subsys)
72
73      CALL cp_subsys_get(subsys=subsys, local_particles=local_particles, &
74                         particles=particles, local_molecules=local_molecules, &
75                         molecule_kinds=molecule_kinds, gci=gci, molecules=molecules)
76
77      molecule_kind_set => molecule_kinds%els
78      molecule_set => molecules%els
79      particle_set => particles%els
80      nkind = SIZE(molecule_kind_set)
81      ! Intermolecular Virtual Site Constraints
82      do_ext_constraint = .FALSE.
83      IF (ASSOCIATED(gci)) THEN
84         do_ext_constraint = (gci%ntot /= 0)
85      END IF
86      ! Intramolecular Virtual Site Constraints
87      MOL: DO ikind = 1, nkind
88         nmol_per_kind = local_molecules%n_el(ikind)
89         DO imol = 1, nmol_per_kind
90            i = local_molecules%list(ikind)%array(imol)
91            molecule => molecule_set(i)
92            molecule_kind => molecule%molecule_kind
93            CALL get_molecule_kind(molecule_kind, nconstraint=nconstraint, nvsite=nvsitecon)
94            IF (nconstraint == 0) CYCLE
95            IF (nvsitecon /= 0) &
96               CALL force_vsite_int(molecule, particle_set)
97         END DO
98      END DO MOL
99      ! Intermolecular Virtual Site Constraints
100      IF (do_ext_constraint) THEN
101         IF (gci%nvsite /= 0) &
102            CALL force_vsite_ext(gci, particle_set)
103      END IF
104
105   END SUBROUTINE vsite_force_control
106
107! **************************************************************************************************
108!> \brief Intramolecular virtual site
109!> \param molecule ...
110!> \param pos ...
111!> \par History
112!>      12.2008 Marcel Baer
113! **************************************************************************************************
114   SUBROUTINE shake_vsite_int(molecule, pos)
115      TYPE(molecule_type), POINTER                       :: molecule
116      REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :)
117
118      INTEGER                                            :: first_atom, nvsite
119      TYPE(molecule_kind_type), POINTER                  :: molecule_kind
120      TYPE(vsite_constraint_type), POINTER               :: vsite_list(:)
121
122      molecule_kind => molecule%molecule_kind
123      CALL get_molecule_kind(molecule_kind, nvsite=nvsite, vsite_list=vsite_list)
124      CALL get_molecule(molecule, first_atom=first_atom)
125      ! Real Shake
126      CALL shake_vsite_low(vsite_list, nvsite, first_atom, pos)
127
128   END SUBROUTINE shake_vsite_int
129
130! **************************************************************************************************
131!> \brief Intramolecular virtual site
132!> \param gci ...
133!> \param pos ...
134!> \par History
135!>      12.2008 Marcel Baer
136! **************************************************************************************************
137   SUBROUTINE shake_vsite_ext(gci, pos)
138
139      TYPE(global_constraint_type), POINTER              :: gci
140      REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :)
141
142      INTEGER                                            :: first_atom, nvsite
143      TYPE(vsite_constraint_type), POINTER               :: vsite_list(:)
144
145      first_atom = 1
146      nvsite = gci%nvsite
147      vsite_list => gci%vsite_list
148      ! Real Shake
149      CALL shake_vsite_low(vsite_list, nvsite, first_atom, pos)
150
151   END SUBROUTINE shake_vsite_ext
152
153! **************************************************************************************************
154!> \brief ...
155!> \param vsite_list ...
156!> \param nvsite ...
157!> \param first_atom ...
158!> \param pos ...
159!> \par History
160!>      12.2008 Marcel Bear
161! **************************************************************************************************
162   SUBROUTINE shake_vsite_low(vsite_list, nvsite, first_atom, pos)
163      TYPE(vsite_constraint_type)                        :: vsite_list(:)
164      INTEGER, INTENT(IN)                                :: nvsite, first_atom
165      REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :)
166
167      INTEGER                                            :: iconst, index_a, index_b, index_c, &
168                                                            index_d
169      REAL(KIND=dp), DIMENSION(3)                        :: r1, r2
170
171      DO iconst = 1, nvsite
172         IF (vsite_list(iconst)%restraint%active) CYCLE
173         index_a = vsite_list(iconst)%a + first_atom - 1
174         index_b = vsite_list(iconst)%b + first_atom - 1
175         index_c = vsite_list(iconst)%c + first_atom - 1
176         index_d = vsite_list(iconst)%d + first_atom - 1
177
178         r1(:) = pos(:, index_b) - pos(:, index_c)
179         r2(:) = pos(:, index_d) - pos(:, index_c)
180         pos(:, index_a) = pos(:, index_c) + vsite_list(iconst)%wbc*r1(:) + &
181                           vsite_list(iconst)%wdc*r2(:)
182      END DO
183   END SUBROUTINE shake_vsite_low
184
185! **************************************************************************************************
186!> \brief Intramolecular virtual site
187!> \param molecule ...
188!> \param particle_set ...
189!> \par History
190!>      12.2008 Marcel Bear
191! **************************************************************************************************
192   SUBROUTINE force_vsite_int(molecule, particle_set)
193      TYPE(molecule_type), POINTER                       :: molecule
194      TYPE(particle_type), POINTER                       :: particle_set(:)
195
196      INTEGER                                            :: first_atom, iconst, index_a, index_b, &
197                                                            index_c, index_d, nvsite
198      REAL(KIND=dp)                                      :: wb, wc, wd
199      TYPE(molecule_kind_type), POINTER                  :: molecule_kind
200      TYPE(vsite_constraint_type), POINTER               :: vsite_list(:)
201
202      molecule_kind => molecule%molecule_kind
203      CALL get_molecule_kind(molecule_kind, nvsite=nvsite, vsite_list=vsite_list)
204      CALL get_molecule(molecule, first_atom=first_atom)
205
206      DO iconst = 1, nvsite
207         IF (vsite_list(iconst)%restraint%active) CYCLE
208         index_a = vsite_list(iconst)%a + first_atom - 1
209         index_b = vsite_list(iconst)%b + first_atom - 1
210         index_c = vsite_list(iconst)%c + first_atom - 1
211         index_d = vsite_list(iconst)%d + first_atom - 1
212
213         wb = vsite_list(iconst)%wbc
214         wd = vsite_list(iconst)%wdc
215         wc = 1.0_dp - vsite_list(iconst)%wbc - vsite_list(iconst)%wdc
216
217         particle_set(index_b)%f(:) = particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
218         particle_set(index_c)%f(:) = particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
219         particle_set(index_d)%f(:) = particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
220         particle_set(index_a)%f(:) = 0.0_dp
221      END DO
222
223   END SUBROUTINE force_vsite_int
224
225! **************************************************************************************************
226!> \brief Intramolecular virtual site
227!> \param gci ...
228!> \param particle_set ...
229!> \par History
230!>      12.2008 Marcel Bear
231! **************************************************************************************************
232   SUBROUTINE force_vsite_ext(gci, particle_set)
233      TYPE(global_constraint_type), POINTER              :: gci
234      TYPE(particle_type), POINTER                       :: particle_set(:)
235
236      INTEGER                                            :: first_atom, iconst, index_a, index_b, &
237                                                            index_c, index_d, nvsite
238      REAL(KIND=dp)                                      :: wb, wc, wd
239      TYPE(vsite_constraint_type), POINTER               :: vsite_list(:)
240
241      first_atom = 1
242      nvsite = gci%nvsite
243      vsite_list => gci%vsite_list
244      ! Real Shake
245
246      DO iconst = 1, nvsite
247         IF (vsite_list(iconst)%restraint%active) CYCLE
248         index_a = vsite_list(iconst)%a + first_atom - 1
249         index_b = vsite_list(iconst)%b + first_atom - 1
250         index_c = vsite_list(iconst)%c + first_atom - 1
251         index_d = vsite_list(iconst)%d + first_atom - 1
252
253         wb = vsite_list(iconst)%wbc
254         wd = vsite_list(iconst)%wdc
255         wc = 1.0_dp - vsite_list(iconst)%wbc - vsite_list(iconst)%wdc
256
257         particle_set(index_b)%f(:) = particle_set(index_b)%f(:) + wb*particle_set(index_a)%f(:)
258         particle_set(index_c)%f(:) = particle_set(index_c)%f(:) + wc*particle_set(index_a)%f(:)
259         particle_set(index_d)%f(:) = particle_set(index_d)%f(:) + wd*particle_set(index_a)%f(:)
260         particle_set(index_a)%f(:) = 0.0_dp
261      END DO
262   END SUBROUTINE force_vsite_ext
263
264END MODULE constraint_vsite
265