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