1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief initialize fist environment
8!> \author CJM
9! **************************************************************************************************
10MODULE fist_environment
11   USE atomic_kind_types,               ONLY: atomic_kind_type,&
12                                              get_atomic_kind_set
13   USE bibliography,                    ONLY: Devynck2012,&
14                                              Dick1958,&
15                                              Mitchell1993,&
16                                              cite_reference
17   USE cell_methods,                    ONLY: read_cell,&
18                                              write_cell
19   USE cell_types,                      ONLY: cell_release,&
20                                              cell_type,&
21                                              get_cell
22   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
23                                              cp_logger_type
24   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
25                                              cp_print_key_unit_nr
26   USE cp_para_types,                   ONLY: cp_para_env_type
27   USE cp_subsys_methods,               ONLY: cp_subsys_create
28   USE cp_subsys_types,                 ONLY: cp_subsys_release,&
29                                              cp_subsys_set,&
30                                              cp_subsys_type
31   USE cp_symmetry,                     ONLY: write_symmetry
32   USE distribution_1d_types,           ONLY: distribution_1d_release,&
33                                              distribution_1d_type
34   USE distribution_methods,            ONLY: distribute_molecules_1d
35   USE ewald_environment_types,         ONLY: ewald_env_create,&
36                                              ewald_env_get,&
37                                              ewald_env_release,&
38                                              ewald_env_set,&
39                                              ewald_environment_type,&
40                                              read_ewald_section
41   USE ewald_pw_methods,                ONLY: ewald_pw_grid_update
42   USE ewald_pw_types,                  ONLY: ewald_pw_create,&
43                                              ewald_pw_release,&
44                                              ewald_pw_type
45   USE exclusion_types,                 ONLY: exclusion_type
46   USE fist_efield_types,               ONLY: fist_efield_type,&
47                                              read_efield_section
48   USE fist_energy_types,               ONLY: allocate_fist_energy,&
49                                              fist_energy_type
50   USE fist_environment_types,          ONLY: fist_env_get,&
51                                              fist_env_set,&
52                                              fist_environment_type
53   USE fist_nonbond_env_types,          ONLY: fist_nonbond_env_release,&
54                                              fist_nonbond_env_type
55   USE force_fields,                    ONLY: force_field_control
56   USE header,                          ONLY: fist_header
57   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
58                                              section_vals_type
59   USE kinds,                           ONLY: dp
60   USE molecule_kind_types,             ONLY: molecule_kind_type,&
61                                              write_molecule_kind_set
62   USE molecule_types,                  ONLY: molecule_type
63   USE multipole_types,                 ONLY: create_multipole_type,&
64                                              multipole_type,&
65                                              release_multipole_type
66   USE particle_list_types,             ONLY: particle_list_create,&
67                                              particle_list_release,&
68                                              particle_list_type
69   USE particle_methods,                ONLY: write_fist_particle_coordinates,&
70                                              write_particle_distances,&
71                                              write_structure_data
72   USE particle_types,                  ONLY: particle_type
73   USE qmmm_types_low,                  ONLY: qmmm_env_mm_type
74#include "./base/base_uses.f90"
75
76   IMPLICIT NONE
77
78   PRIVATE
79
80   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_environment'
81   PUBLIC :: fist_init
82
83CONTAINS
84! **************************************************************************************************
85!> \brief reads the input and database file for fist
86!> \param fist_env ...
87!> \param root_section ...
88!> \param para_env ...
89!> \param force_env_section ...
90!> \param subsys_section ...
91!> \param use_motion_section ...
92!> \param prev_subsys ...
93!> \par Used By
94!>      fist_main
95! **************************************************************************************************
96   SUBROUTINE fist_init(fist_env, root_section, para_env, force_env_section, &
97                        subsys_section, use_motion_section, prev_subsys)
98
99      TYPE(fist_environment_type), POINTER               :: fist_env
100      TYPE(section_vals_type), POINTER                   :: root_section
101      TYPE(cp_para_env_type), POINTER                    :: para_env
102      TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
103      LOGICAL, INTENT(IN)                                :: use_motion_section
104      TYPE(cp_subsys_type), OPTIONAL, POINTER            :: prev_subsys
105
106      CHARACTER(len=*), PARAMETER :: routineN = 'fist_init', routineP = moduleN//':'//routineN
107
108      INTEGER                                            :: handle, iw
109      LOGICAL                                            :: qmmm, shell_adiabatic, shell_present
110      REAL(KIND=dp), DIMENSION(3)                        :: abc
111      TYPE(cell_type), POINTER                           :: cell, cell_ref
112      TYPE(cp_logger_type), POINTER                      :: logger
113      TYPE(cp_subsys_type), POINTER                      :: subsys
114      TYPE(ewald_environment_type), POINTER              :: ewald_env
115      TYPE(exclusion_type), DIMENSION(:), POINTER        :: exclusions
116      TYPE(fist_efield_type), POINTER                    :: efield
117      TYPE(fist_nonbond_env_type), POINTER               :: fist_nonbond_env
118      TYPE(particle_list_type), POINTER                  :: core_particles, shell_particles
119      TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, shell_particle_set
120      TYPE(qmmm_env_mm_type), POINTER                    :: qmmm_env
121      TYPE(section_vals_type), POINTER                   :: cell_section, ewald_section, mm_section, &
122                                                            poisson_section
123
124      CALL timeset(routineN, handle)
125      logger => cp_get_default_logger()
126
127      NULLIFY (subsys, cell, cell_ref)
128      NULLIFY (ewald_env, fist_nonbond_env, qmmm_env, cell_section, &
129               poisson_section, shell_particle_set, shell_particles, &
130               core_particle_set, core_particles, exclusions)
131      IF (.NOT. ASSOCIATED(subsys_section)) THEN
132         subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
133      END IF
134      mm_section => section_vals_get_subs_vals(force_env_section, "MM")
135      cell_section => section_vals_get_subs_vals(subsys_section, "CELL")
136      poisson_section => section_vals_get_subs_vals(mm_section, "POISSON")
137      ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
138
139      CALL fist_env_set(fist_env, input=force_env_section)
140
141      iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%PROGRAM_BANNER", &
142                                extension=".mmLog")
143      CALL fist_header(iw)
144      CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%PROGRAM_BANNER")
145
146      CALL read_cell(cell, cell_ref, cell_section=cell_section, para_env=para_env)
147      CALL get_cell(cell, abc=abc)
148
149      ! Print the cell parameters
150      CALL write_cell(cell, subsys_section, cell_ref)
151
152      ! Create the ewald environment
153      CALL ewald_env_create(ewald_env, para_env)
154
155      ! Read the input section and set the ewald_env
156      CALL read_ewald_section(ewald_env, ewald_section)
157      CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
158
159      ! Read the efield section
160      NULLIFY (efield)
161      CALL read_efield_section(mm_section, efield)
162      CALL fist_env_set(fist_env, efield=efield)
163
164      ! Topology
165      CALL fist_env_get(fist_env, qmmm=qmmm, qmmm_env=qmmm_env)
166      CALL cp_subsys_create(subsys, para_env=para_env, root_section=root_section, &
167                            force_env_section=force_env_section, subsys_section=subsys_section, &
168                            qmmm=qmmm, qmmm_env=qmmm_env, exclusions=exclusions, &
169                            use_motion_section=use_motion_section)
170      CALL fist_env_set(fist_env, subsys=subsys, exclusions=exclusions)
171
172      CALL force_field_control(subsys%atomic_kinds%els, subsys%particles%els, &
173                               subsys%molecule_kinds%els, subsys%molecules%els, &
174                               ewald_env, fist_nonbond_env, root_section, para_env, qmmm=qmmm, &
175                               qmmm_env=qmmm_env, subsys_section=subsys_section, &
176                               mm_section=mm_section, shell_particle_set=shell_particle_set, &
177                               core_particle_set=core_particle_set, cell=cell)
178
179      NULLIFY (shell_particles, core_particles)
180      IF (ASSOCIATED(shell_particle_set)) THEN
181         CALL cite_reference(Devynck2012)
182         CALL cite_reference(Mitchell1993)
183         CALL cite_reference(Dick1958)
184         CALL particle_list_create(shell_particles, els_ptr=shell_particle_set)
185      END IF
186      IF (ASSOCIATED(core_particle_set)) THEN
187         CALL particle_list_create(core_particles, els_ptr=core_particle_set)
188      END IF
189      CALL get_atomic_kind_set(atomic_kind_set=subsys%atomic_kinds%els, &
190                               shell_present=shell_present, shell_adiabatic=shell_adiabatic)
191      CALL fist_env_set(fist_env, shell_model=shell_present, &
192                        shell_model_ad=shell_adiabatic)
193      CALL cp_subsys_set(subsys, shell_particles=shell_particles, &
194                         core_particles=core_particles)
195      CALL particle_list_release(shell_particles)
196      CALL particle_list_release(core_particles)
197
198      CALL fist_init_subsys(fist_env, subsys, cell, cell_ref, fist_nonbond_env, ewald_env, &
199                            force_env_section, subsys_section, prev_subsys)
200
201      CALL cell_release(cell)
202      CALL cell_release(cell_ref)
203      CALL ewald_env_release(ewald_env)
204      CALL fist_nonbond_env_release(fist_nonbond_env)
205      CALL cp_subsys_release(subsys)
206
207      CALL timestop(handle)
208
209   END SUBROUTINE fist_init
210
211! **************************************************************************************************
212!> \brief   Read the input and the database files for the setup of the
213!>          FIST environment.
214!> \param fist_env ...
215!> \param subsys ...
216!> \param cell ...
217!> \param cell_ref ...
218!> \param fist_nonbond_env ...
219!> \param ewald_env ...
220!> \param force_env_section ...
221!> \param subsys_section ...
222!> \param prev_subsys ...
223!> \date    22.05.2000
224!> \author  MK
225!> \version 1.0
226! **************************************************************************************************
227   SUBROUTINE fist_init_subsys(fist_env, subsys, cell, cell_ref, fist_nonbond_env, &
228                               ewald_env, force_env_section, subsys_section, &
229                               prev_subsys)
230
231      TYPE(fist_environment_type), POINTER               :: fist_env
232      TYPE(cp_subsys_type), POINTER                      :: subsys
233      TYPE(cell_type), POINTER                           :: cell, cell_ref
234      TYPE(fist_nonbond_env_type), POINTER               :: fist_nonbond_env
235      TYPE(ewald_environment_type), POINTER              :: ewald_env
236      TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
237      TYPE(cp_subsys_type), OPTIONAL, POINTER            :: prev_subsys
238
239      CHARACTER(len=*), PARAMETER :: routineN = 'fist_init_subsys', &
240         routineP = moduleN//':'//routineN
241
242      INTEGER                                            :: handle, max_multipole
243      LOGICAL                                            :: do_multipoles
244      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
245      TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles, &
246                                                            prev_local_molecules
247      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
248      TYPE(fist_energy_type), POINTER                    :: thermo
249      TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set, prev_molecule_kind_set
250      TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
251      TYPE(multipole_type), POINTER                      :: multipoles
252      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
253      TYPE(section_vals_type), POINTER                   :: grid_print_section
254
255      CALL timeset(routineN, handle)
256      NULLIFY (thermo, ewald_pw, local_molecules, local_particles, multipoles)
257      particle_set => subsys%particles%els
258      atomic_kind_set => subsys%atomic_kinds%els
259      molecule_set => subsys%molecules%els
260      molecule_kind_set => subsys%molecule_kinds%els
261
262      IF (PRESENT(prev_subsys)) THEN
263         prev_molecule_kind_set => prev_subsys%molecule_kinds%els
264         prev_local_molecules => prev_subsys%local_molecules
265      ELSE
266         NULLIFY (prev_molecule_kind_set)
267         NULLIFY (prev_local_molecules)
268      ENDIF
269
270      ! Create the fist_energy_type
271      CALL allocate_fist_energy(thermo)
272
273      ! Print the molecule kind set
274      CALL write_molecule_kind_set(molecule_kind_set, subsys_section)
275
276      ! Print the atomic coordinates
277      CALL write_fist_particle_coordinates(particle_set, subsys_section, &
278                                           fist_nonbond_env%charges)
279      CALL write_particle_distances(particle_set, cell, subsys_section)
280      CALL write_structure_data(particle_set, cell=cell, input_section=subsys_section)
281
282      ! Print symmetry information
283      CALL write_symmetry(particle_set, cell, subsys_section)
284
285      ! Distribute molecules and atoms using the new data structures ***
286      CALL distribute_molecules_1d(atomic_kind_set=atomic_kind_set, &
287                                   particle_set=particle_set, &
288                                   local_particles=local_particles, &
289                                   molecule_kind_set=molecule_kind_set, &
290                                   molecule_set=molecule_set, &
291                                   local_molecules=local_molecules, &
292                                   prev_molecule_kind_set=prev_molecule_kind_set, &
293                                   prev_local_molecules=prev_local_molecules, &
294                                   force_env_section=force_env_section)
295
296      ! Create ewald grids
297      grid_print_section => section_vals_get_subs_vals(force_env_section, &
298                                                       "PRINT%GRID_INFORMATION")
299      CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, grid_print_section)
300
301      ! Initialize ewald grids
302      CALL ewald_pw_grid_update(ewald_pw, ewald_env, cell%hmat)
303
304      ! Possibly Initialize the multipole environment
305      CALL ewald_env_get(ewald_env, do_multipoles=do_multipoles, &
306                         max_multipole=max_multipole)
307      IF (do_multipoles) &
308         CALL create_multipole_type(multipoles, particle_set, subsys_section, max_multipole)
309      CALL cp_subsys_set(subsys, multipoles=multipoles, cell=cell)
310
311      ! Set the fist_env
312      CALL fist_env_set(fist_env=fist_env, &
313                        cell_ref=cell_ref, &
314                        local_molecules=local_molecules, &
315                        local_particles=local_particles, &
316                        ewald_env=ewald_env, ewald_pw=ewald_pw, &
317                        fist_nonbond_env=fist_nonbond_env, &
318                        thermo=thermo)
319
320      CALL distribution_1d_release(local_particles)
321      CALL distribution_1d_release(local_molecules)
322      CALL ewald_pw_release(ewald_pw)
323      CALL release_multipole_type(multipoles)
324      CALL timestop(handle)
325
326   END SUBROUTINE fist_init_subsys
327END MODULE fist_environment
328