Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2023 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Contains types used for a Dimer Method calculations
10 : !> \par History
11 : !> Luca Bellucci 11.2017 added kdimer and beta
12 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
13 : ! **************************************************************************************************
14 : MODULE dimer_types
15 :
16 : USE cell_types, ONLY: use_perd_x,&
17 : use_perd_xy,&
18 : use_perd_xyz,&
19 : use_perd_xz,&
20 : use_perd_y,&
21 : use_perd_yz,&
22 : use_perd_z
23 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
24 : USE cp_subsys_types, ONLY: cp_subsys_get,&
25 : cp_subsys_type
26 : USE global_types, ONLY: global_environment_type
27 : USE input_constants, ONLY: do_first_rotation_step
28 : USE input_section_types, ONLY: section_vals_get,&
29 : section_vals_get_subs_vals,&
30 : section_vals_type,&
31 : section_vals_val_get
32 : USE kinds, ONLY: dp
33 : USE molecule_kind_list_types, ONLY: molecule_kind_list_type
34 : USE molecule_kind_types, ONLY: fixd_constraint_type,&
35 : get_molecule_kind,&
36 : molecule_kind_type
37 : #include "../base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 : PRIVATE
41 :
42 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dimer_types'
44 :
45 : PUBLIC :: dimer_env_type, &
46 : dimer_env_create, &
47 : dimer_env_retain, &
48 : dimer_env_release, &
49 : dimer_fixed_atom_control
50 :
51 : ! **************************************************************************************************
52 : !> \brief Type containing all informations abour the rotation of the Dimer
53 : !> \par History
54 : !> none
55 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
56 : ! **************************************************************************************************
57 : TYPE dimer_rotational_type
58 : ! Rotational parameters
59 : INTEGER :: rotation_step
60 : LOGICAL :: interpolate_gradient
61 : REAL(KIND=dp) :: angle_tol, angle1, angle2, dCdp, curvature
62 : REAL(KIND=dp), POINTER, DIMENSION(:) :: g0, g1, g1p
63 : END TYPE dimer_rotational_type
64 :
65 : ! **************************************************************************************************
66 : !> \brief Type containing all informations abour the translation of the Dimer
67 : !> \par History
68 : !> none
69 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
70 : ! **************************************************************************************************
71 : TYPE dimer_translational_type
72 : ! Translational parameters
73 : REAL(KIND=dp), POINTER, DIMENSION(:) :: tls_vec
74 : END TYPE dimer_translational_type
75 :
76 : ! **************************************************************************************************
77 : !> \brief Conjugate Directions type
78 : !> \par History
79 : !> none
80 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
81 : ! **************************************************************************************************
82 : TYPE dimer_cg_rot_type
83 : REAL(KIND=dp) :: norm_theta, norm_theta_old, norm_h
84 : REAL(KIND=dp), POINTER, DIMENSION(:) :: nvec_old
85 : END TYPE dimer_cg_rot_type
86 :
87 : ! **************************************************************************************************
88 : !> \brief Defines the environment for a Dimer Method calculation
89 : !> \par History
90 : !> none
91 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
92 : ! **************************************************************************************************
93 : TYPE dimer_env_type
94 : INTEGER :: ref_count
95 : REAL(KIND=dp) :: dr
96 : REAL(KIND=dp), POINTER, DIMENSION(:) :: nvec
97 : REAL(KIND=dp) :: beta
98 : TYPE(dimer_rotational_type) :: rot
99 : TYPE(dimer_translational_type) :: tsl
100 : TYPE(dimer_cg_rot_type) :: cg_rot
101 : LOGICAL :: kdimer
102 : END TYPE dimer_env_type
103 :
104 : CONTAINS
105 :
106 : ! **************************************************************************************************
107 : !> \brief ...
108 : !> \param dimer_env ...
109 : !> \param subsys ...
110 : !> \param globenv ...
111 : !> \param dimer_section ...
112 : !> \par History
113 : !> Luca Bellucci 11.2017 added K-DIMER and BETA
114 : !> 2016/03/03 [LTong] changed input natom to subsys
115 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
116 : ! **************************************************************************************************
117 22 : SUBROUTINE dimer_env_create(dimer_env, subsys, globenv, dimer_section)
118 : TYPE(dimer_env_type), POINTER :: dimer_env
119 : TYPE(cp_subsys_type), POINTER :: subsys
120 : TYPE(global_environment_type), POINTER :: globenv
121 : TYPE(section_vals_type), POINTER :: dimer_section
122 :
123 : INTEGER :: i, isize, j, k, n_rep_val, natom, unit_nr
124 : LOGICAL :: explicit
125 : REAL(KIND=dp) :: norm, xval(3)
126 22 : REAL(KIND=dp), DIMENSION(:), POINTER :: array
127 : TYPE(section_vals_type), POINTER :: nvec_section
128 :
129 44 : unit_nr = cp_logger_get_default_io_unit()
130 22 : CPASSERT(.NOT. ASSOCIATED(dimer_env))
131 22 : ALLOCATE (dimer_env)
132 22 : dimer_env%ref_count = 1
133 : ! Setup NVEC
134 22 : NULLIFY (dimer_env%nvec, dimer_env%rot%g0, dimer_env%rot%g1, dimer_env%rot%g1p, &
135 22 : dimer_env%tsl%tls_vec)
136 : ! get natom
137 22 : CALL cp_subsys_get(subsys=subsys, natom=natom)
138 : ! Allocate the working arrays
139 66 : ALLOCATE (dimer_env%nvec(natom*3))
140 66 : ALLOCATE (dimer_env%rot%g0(natom*3))
141 66 : ALLOCATE (dimer_env%rot%g1(natom*3))
142 66 : ALLOCATE (dimer_env%rot%g1p(natom*3))
143 : ! Check if the dimer vector is available in the input or not..
144 22 : nvec_section => section_vals_get_subs_vals(dimer_section, "DIMER_VECTOR")
145 22 : CALL section_vals_get(nvec_section, explicit=explicit)
146 22 : IF (explicit) THEN
147 4 : IF (unit_nr > 0) WRITE (unit_nr, *) "Reading Dimer Vector from file!"
148 4 : NULLIFY (array)
149 4 : CALL section_vals_val_get(nvec_section, "_DEFAULT_KEYWORD_", n_rep_val=n_rep_val)
150 4 : isize = 0
151 10 : DO i = 1, n_rep_val
152 6 : CALL section_vals_val_get(nvec_section, "_DEFAULT_KEYWORD_", r_vals=array, i_rep_val=i)
153 34 : DO j = 1, SIZE(array)
154 24 : isize = isize + 1
155 30 : dimer_env%nvec(isize) = array(j)
156 : END DO
157 : END DO
158 4 : CPASSERT(isize == SIZE(dimer_env%nvec))
159 : ELSE
160 18 : CALL globenv%gaussian_rng_stream%fill(dimer_env%nvec)
161 : END IF
162 : ! Check for translation in the dimer vector and remove them
163 22 : IF (natom > 1) THEN
164 20 : xval = 0.0_dp
165 90 : DO j = 1, natom
166 300 : DO k = 1, 3
167 210 : i = (j - 1)*3 + k
168 280 : xval(k) = xval(k) + dimer_env%nvec(i)
169 : END DO
170 : END DO
171 : ! Subtract net translations
172 80 : xval = xval/REAL(natom*3, KIND=dp)
173 90 : DO j = 1, natom
174 300 : DO k = 1, 3
175 210 : i = (j - 1)*3 + k
176 280 : dimer_env%nvec(i) = dimer_env%nvec(i) - xval(k)
177 : END DO
178 : END DO
179 : END IF
180 : ! set nvec components to zero for the corresponding constraints
181 22 : CALL dimer_fixed_atom_control(dimer_env%nvec, subsys)
182 238 : norm = SQRT(SUM(dimer_env%nvec**2))
183 22 : IF (norm <= EPSILON(0.0_dp)) &
184 0 : CPABORT("The norm of the dimer vector is 0! Calculation cannot proceed further.")
185 238 : dimer_env%nvec = dimer_env%nvec/norm
186 22 : dimer_env%rot%rotation_step = do_first_rotation_step
187 22 : CALL section_vals_val_get(dimer_section, "DR", r_val=dimer_env%dr)
188 : CALL section_vals_val_get(dimer_section, "INTERPOLATE_GRADIENT", &
189 22 : l_val=dimer_env%rot%interpolate_gradient)
190 : CALL section_vals_val_get(dimer_section, "ANGLE_TOLERANCE", &
191 22 : r_val=dimer_env%rot%angle_tol)
192 : CALL section_vals_val_get(dimer_section, "K-DIMER", &
193 22 : l_val=dimer_env%kdimer)
194 : CALL section_vals_val_get(dimer_section, "BETA", &
195 22 : r_val=dimer_env%beta)
196 : ! initialise values
197 22 : dimer_env%cg_rot%norm_h = 1.0_dp
198 22 : dimer_env%cg_rot%norm_theta = 0.0_dp
199 22 : dimer_env%cg_rot%norm_theta_old = 0.0_dp
200 238 : dimer_env%rot%g0 = 0.0_dp
201 238 : dimer_env%rot%g1 = 0.0_dp
202 238 : dimer_env%rot%g1p = 0.0_dp
203 66 : ALLOCATE (dimer_env%cg_rot%nvec_old(natom*3))
204 66 : END SUBROUTINE dimer_env_create
205 :
206 : ! **************************************************************************************************
207 : !> \brief ...
208 : !> \param dimer_env ...
209 : !> \par History
210 : !> none
211 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
212 : ! **************************************************************************************************
213 22 : SUBROUTINE dimer_env_retain(dimer_env)
214 : TYPE(dimer_env_type), POINTER :: dimer_env
215 :
216 22 : CPASSERT(ASSOCIATED(dimer_env))
217 22 : CPASSERT(dimer_env%ref_count > 0)
218 22 : dimer_env%ref_count = dimer_env%ref_count + 1
219 22 : END SUBROUTINE dimer_env_retain
220 :
221 : ! **************************************************************************************************
222 : !> \brief ...
223 : !> \param dimer_env ...
224 : !> \par History
225 : !> none
226 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
227 : ! **************************************************************************************************
228 1863 : SUBROUTINE dimer_env_release(dimer_env)
229 : TYPE(dimer_env_type), POINTER :: dimer_env
230 :
231 1863 : IF (ASSOCIATED(dimer_env)) THEN
232 44 : CPASSERT(dimer_env%ref_count > 0)
233 44 : dimer_env%ref_count = dimer_env%ref_count - 1
234 44 : IF (dimer_env%ref_count == 0) THEN
235 22 : IF (ASSOCIATED(dimer_env%nvec)) THEN
236 22 : DEALLOCATE (dimer_env%nvec)
237 : END IF
238 22 : IF (ASSOCIATED(dimer_env%rot%g0)) THEN
239 22 : DEALLOCATE (dimer_env%rot%g0)
240 : END IF
241 22 : IF (ASSOCIATED(dimer_env%rot%g1)) THEN
242 22 : DEALLOCATE (dimer_env%rot%g1)
243 : END IF
244 22 : IF (ASSOCIATED(dimer_env%rot%g1p)) THEN
245 22 : DEALLOCATE (dimer_env%rot%g1p)
246 : END IF
247 22 : IF (ASSOCIATED(dimer_env%cg_rot%nvec_old)) THEN
248 22 : DEALLOCATE (dimer_env%cg_rot%nvec_old)
249 : END IF
250 : ! No need to deallocate tls_vec (just a pointer to aother local array)
251 22 : NULLIFY (dimer_env%tsl%tls_vec)
252 22 : DEALLOCATE (dimer_env)
253 : END IF
254 : END IF
255 1863 : END SUBROUTINE dimer_env_release
256 :
257 : ! **************************************************************************************************
258 : !> \brief Set parts of a given array vec to zero according to fixed atom constraints.
259 : !> When atoms are (partially) fixed then the relevant components of
260 : !> nvec should be set to zero. Furthermore, the relevant components
261 : !> of the gradient in CG should also be set to zero.
262 : !> \param vec : vector to be modified
263 : !> \param subsys : subsys type object used by CP2k
264 : !> \par History
265 : !> 2016/03/03 [LTong] created
266 : !> \author Lianheng Tong [LTong]
267 : ! **************************************************************************************************
268 410 : SUBROUTINE dimer_fixed_atom_control(vec, subsys)
269 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vec
270 : TYPE(cp_subsys_type), POINTER :: subsys
271 :
272 : INTEGER :: ii, ikind, ind, iparticle, nfixed_atoms, &
273 : nkinds
274 410 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
275 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
276 410 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
277 : TYPE(molecule_kind_type), POINTER :: molecule_kind
278 :
279 410 : NULLIFY (molecule_kinds, molecule_kind, fixd_list)
280 :
281 : ! need to get constraint information from molecule information
282 : CALL cp_subsys_get(subsys=subsys, &
283 410 : molecule_kinds=molecule_kinds)
284 410 : molecule_kind_set => molecule_kinds%els
285 :
286 : ! get total number of fixed atoms
287 : ! nkinds is the kinds of molecules, not atoms
288 410 : nkinds = molecule_kinds%n_els
289 1026 : DO ikind = 1, nkinds
290 616 : molecule_kind => molecule_kind_set(ikind)
291 : CALL get_molecule_kind(molecule_kind, &
292 : nfixd=nfixed_atoms, &
293 616 : fixd_list=fixd_list)
294 1026 : IF (ASSOCIATED(fixd_list)) THEN
295 224 : DO ii = 1, nfixed_atoms
296 224 : IF (.NOT. fixd_list(ii)%restraint%active) THEN
297 56 : iparticle = fixd_list(ii)%fixd
298 56 : ind = (iparticle - 1)*3
299 : ! apply constraint to nvec
300 56 : SELECT CASE (fixd_list(ii)%itype)
301 : CASE (use_perd_x)
302 0 : vec(ind + 1) = 0.0_dp
303 : CASE (use_perd_y)
304 0 : vec(ind + 2) = 0.0_dp
305 : CASE (use_perd_z)
306 0 : vec(ind + 3) = 0.0_dp
307 : CASE (use_perd_xy)
308 0 : vec(ind + 1) = 0.0_dp
309 0 : vec(ind + 2) = 0.0_dp
310 : CASE (use_perd_xz)
311 0 : vec(ind + 1) = 0.0_dp
312 0 : vec(ind + 3) = 0.0_dp
313 : CASE (use_perd_yz)
314 0 : vec(ind + 2) = 0.0_dp
315 0 : vec(ind + 3) = 0.0_dp
316 : CASE (use_perd_xyz)
317 56 : vec(ind + 1) = 0.0_dp
318 56 : vec(ind + 2) = 0.0_dp
319 56 : vec(ind + 3) = 0.0_dp
320 : END SELECT
321 : END IF ! .NOT.fixd_list(ii)%restraint%active
322 : END DO ! ii
323 : END IF ! ASSOCIATED(fixd_list)
324 : END DO ! ikind
325 410 : END SUBROUTINE dimer_fixed_atom_control
326 :
327 0 : END MODULE dimer_types
|
