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