1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief
8!> \author
9! **************************************************************************************************
10MODULE shell_opt
11
12   USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
13   USE atomic_kind_types,               ONLY: atomic_kind_type,&
14                                              get_atomic_kind
15   USE cg_optimizer,                    ONLY: geoopt_cg
16   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
17                                              cp_logger_type
18   USE cp_output_handling,              ONLY: cp_add_iter_level,&
19                                              cp_rm_iter_level
20   USE cp_para_types,                   ONLY: cp_para_env_type
21   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
22                                              cp_subsys_type
23   USE distribution_1d_types,           ONLY: distribution_1d_type
24   USE force_env_types,                 ONLY: force_env_get,&
25                                              force_env_type
26   USE global_types,                    ONLY: global_environment_type
27   USE gopt_f_types,                    ONLY: gopt_f_create,&
28                                              gopt_f_release,&
29                                              gopt_f_type
30   USE gopt_param_types,                ONLY: gopt_param_read,&
31                                              gopt_param_release,&
32                                              gopt_param_type
33   USE input_constants,                 ONLY: default_shellcore_method_id
34   USE input_section_types,             ONLY: section_vals_get,&
35                                              section_vals_get_subs_vals,&
36                                              section_vals_type
37   USE integrator_utils,                ONLY: tmp_variables_type
38   USE kinds,                           ONLY: dp
39   USE message_passing,                 ONLY: mp_sum
40   USE particle_types,                  ONLY: particle_type
41   USE shell_potential_types,           ONLY: shell_kind_type
42#include "../base/base_uses.f90"
43
44   IMPLICIT NONE
45   PRIVATE
46   PUBLIC :: optimize_shell_core
47   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'shell_opt'
48
49CONTAINS
50
51! **************************************************************************************************
52!> \brief Optimize shell-core positions along an MD run
53!> \param force_env ...
54!> \param particle_set ...
55!> \param shell_particle_set ...
56!> \param core_particle_set ...
57!> \param globenv ...
58!> \param tmp ...
59!> \param check ...
60!> \author
61! **************************************************************************************************
62
63   SUBROUTINE optimize_shell_core(force_env, particle_set, shell_particle_set, core_particle_set, globenv, tmp, check)
64      TYPE(force_env_type), POINTER                      :: force_env
65      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set, shell_particle_set, &
66                                                            core_particle_set
67      TYPE(global_environment_type), POINTER             :: globenv
68      TYPE(tmp_variables_type), OPTIONAL, POINTER        :: tmp
69      LOGICAL, INTENT(IN), OPTIONAL                      :: check
70
71      CHARACTER(len=*), PARAMETER :: routineN = 'optimize_shell_core', &
72         routineP = moduleN//':'//routineN
73
74      INTEGER                                            :: handle, i, iat, nshell
75      LOGICAL                                            :: do_update, explicit, my_check, optimize
76      REAL(dp), DIMENSION(:), POINTER                    :: dvec_sc, dvec_sc_0
77      TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
78      TYPE(cp_logger_type), POINTER                      :: logger
79      TYPE(cp_para_env_type), POINTER                    :: para_env
80      TYPE(cp_subsys_type), POINTER                      :: subsys
81      TYPE(distribution_1d_type), POINTER                :: local_particles
82      TYPE(gopt_f_type), POINTER                         :: gopt_env
83      TYPE(gopt_param_type), POINTER                     :: gopt_param
84      TYPE(section_vals_type), POINTER                   :: force_env_section, geo_section, &
85                                                            root_section
86
87      NULLIFY (logger)
88      logger => cp_get_default_logger()
89
90      CPASSERT(ASSOCIATED(force_env))
91      CPASSERT(ASSOCIATED(globenv))
92
93      NULLIFY (gopt_param, force_env_section, gopt_env, dvec_sc, dvec_sc_0, root_section, geo_section)
94      root_section => force_env%root_section
95      force_env_section => force_env%force_env_section
96      geo_section => section_vals_get_subs_vals(root_section, "MOTION%SHELL_OPT")
97
98      CALL section_vals_get(geo_section, explicit=explicit)
99      IF (.NOT. explicit) RETURN
100
101      CALL timeset(routineN, handle)
102
103      optimize = .FALSE.
104      my_check = .FALSE.
105      IF (PRESENT(check)) my_check = check
106      IF (my_check) THEN
107         NULLIFY (subsys, para_env, atomic_kinds, local_particles)
108         CALL force_env_get(force_env=force_env, subsys=subsys, para_env=para_env)
109         CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles)
110         CALL check_shell_core_distance(atomic_kinds, local_particles, particle_set, shell_particle_set, &
111                                        core_particle_set, para_env, optimize)
112
113         IF (.NOT. optimize) THEN
114            CALL timestop(handle)
115            RETURN
116         END IF
117      END IF
118
119      nshell = SIZE(shell_particle_set)
120      ALLOCATE (dvec_sc(3*nshell))
121      ALLOCATE (dvec_sc_0(3*nshell))
122      DO i = 1, nshell
123         dvec_sc(1 + 3*(i - 1)) = core_particle_set(i)%r(1) - shell_particle_set(i)%r(1)
124         dvec_sc(2 + 3*(i - 1)) = core_particle_set(i)%r(2) - shell_particle_set(i)%r(2)
125         dvec_sc(3 + 3*(i - 1)) = core_particle_set(i)%r(3) - shell_particle_set(i)%r(3)
126      END DO
127      dvec_sc_0 = dvec_sc
128
129      CALL gopt_param_read(gopt_param, geo_section, type_id=default_shellcore_method_id)
130      CALL gopt_f_create(gopt_env, gopt_param, force_env=force_env, globenv=globenv, &
131                         geo_opt_section=geo_section)
132
133      CALL cp_add_iter_level(logger%iter_info, "SHELL_OPT")
134      gopt_env%eval_opt_geo = .FALSE.
135      CALL geoopt_cg(force_env, gopt_param, globenv, &
136                     geo_section, gopt_env, dvec_sc, do_update=do_update)
137      IF (.NOT. do_update) THEN
138         DO i = 1, nshell
139            shell_particle_set(i)%r(1) = -dvec_sc_0(1 + 3*(i - 1)) + core_particle_set(i)%r(1)
140            shell_particle_set(i)%r(2) = -dvec_sc_0(2 + 3*(i - 1)) + core_particle_set(i)%r(2)
141            shell_particle_set(i)%r(3) = -dvec_sc_0(3 + 3*(i - 1)) + core_particle_set(i)%r(3)
142         END DO
143      END IF
144      CALL cp_rm_iter_level(logger%iter_info, "SHELL_OPT")
145
146      CALL gopt_f_release(gopt_env)
147      CALL gopt_param_release(gopt_param)
148      DEALLOCATE (dvec_sc)
149      DEALLOCATE (dvec_sc_0)
150
151      IF (PRESENT(tmp)) THEN
152         DO i = 1, nshell
153            iat = shell_particle_set(i)%atom_index
154            tmp%shell_vel(1:3, i) = tmp%vel(1:3, iat)
155            tmp%core_vel(1:3, i) = tmp%vel(1:3, iat)
156         END DO
157      ELSE
158         DO i = 1, nshell
159            iat = shell_particle_set(i)%atom_index
160            shell_particle_set(i)%v(1:3) = particle_set(iat)%v(1:3)
161            core_particle_set(i)%v(1:3) = particle_set(iat)%v(1:3)
162         END DO
163      END IF
164
165      CALL timestop(handle)
166
167   END SUBROUTINE optimize_shell_core
168
169! **************************************************************************************************
170!> \brief Check shell_core_distance
171!> \param atomic_kinds ...
172!> \param local_particles ...
173!> \param particle_set ...
174!> \param shell_particle_set ...
175!> \param core_particle_set ...
176!> \param para_env ...
177!> \param optimize ...
178!> \par History
179!>      none
180!> \author MI (October 2008)
181!>     I soliti ignoti
182! **************************************************************************************************
183   SUBROUTINE check_shell_core_distance(atomic_kinds, local_particles, particle_set, &
184                                        shell_particle_set, core_particle_set, para_env, optimize)
185
186      TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
187      TYPE(distribution_1d_type), POINTER                :: local_particles
188      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set, shell_particle_set, &
189                                                            core_particle_set
190      TYPE(cp_para_env_type), POINTER                    :: para_env
191      LOGICAL, INTENT(INOUT)                             :: optimize
192
193      CHARACTER(LEN=*), PARAMETER :: routineN = 'check_shell_core_distance', &
194         routineP = moduleN//':'//routineN
195
196      INTEGER                                            :: ikind, iparticle, iparticle_local, &
197                                                            itest, nkind, nparticle_local, &
198                                                            shell_index
199      LOGICAL                                            :: is_shell
200      REAL(dp)                                           :: dsc, rc(3), rs(3)
201      TYPE(atomic_kind_type), POINTER                    :: atomic_kind
202      TYPE(shell_kind_type), POINTER                     :: shell
203
204      nkind = atomic_kinds%n_els
205      itest = 0
206      DO ikind = 1, nkind
207         NULLIFY (atomic_kind)
208         atomic_kind => atomic_kinds%els(ikind)
209         CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_shell, shell=shell)
210         IF (is_shell) THEN
211            IF (shell%max_dist > 0.0_dp) THEN
212               nparticle_local = local_particles%n_el(ikind)
213               DO iparticle_local = 1, nparticle_local
214                  iparticle = local_particles%list(ikind)%array(iparticle_local)
215                  shell_index = particle_set(iparticle)%shell_index
216
217                  rc(:) = core_particle_set(shell_index)%r(:)
218                  rs(:) = shell_particle_set(shell_index)%r(:)
219                  dsc = SQRT((rc(1) - rs(1))**2 + (rc(2) - rs(2))**2 + (rc(3) - rs(3))**2)
220                  IF (dsc > shell%max_dist) THEN
221                     itest = 1
222                  END IF
223               END DO
224            END IF
225         END IF
226      END DO
227
228      CALL mp_sum(itest, para_env%group)
229      IF (itest > 0) optimize = .TRUE.
230
231   END SUBROUTINE check_shell_core_distance
232END MODULE shell_opt
233