1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Initialize a small environment for a particular calculation
8!> \par History
9!>      5.2004 created [fawzi]
10!>      9.2007 cleaned [tlaino] - University of Zurich
11!> \author Teodoro Laino
12! **************************************************************************************************
13MODULE cp_subsys_methods
14   USE atomic_kind_list_types,          ONLY: atomic_kind_list_create,&
15                                              atomic_kind_list_release,&
16                                              atomic_kind_list_type
17   USE atomic_kind_types,               ONLY: atomic_kind_type
18   USE atprop_types,                    ONLY: atprop_create
19   USE cell_types,                      ONLY: cell_retain,&
20                                              cell_type
21   USE colvar_methods,                  ONLY: colvar_read
22   USE cp_para_env,                     ONLY: cp_para_env_retain
23   USE cp_para_types,                   ONLY: cp_para_env_type
24   USE cp_result_types,                 ONLY: cp_result_create
25   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
26                                              cp_subsys_set,&
27                                              cp_subsys_type
28   USE exclusion_types,                 ONLY: exclusion_type
29   USE input_constants,                 ONLY: do_conn_off,&
30                                              do_stress_analytical,&
31                                              do_stress_diagonal_anal,&
32                                              do_stress_diagonal_numer,&
33                                              do_stress_none,&
34                                              do_stress_numerical
35   USE input_section_types,             ONLY: section_vals_get,&
36                                              section_vals_get_subs_vals,&
37                                              section_vals_type,&
38                                              section_vals_val_get
39   USE kinds,                           ONLY: default_string_length,&
40                                              dp
41   USE molecule_kind_list_types,        ONLY: molecule_kind_list_create,&
42                                              molecule_kind_list_release,&
43                                              molecule_kind_list_type
44   USE molecule_kind_types,             ONLY: molecule_kind_type
45   USE molecule_list_types,             ONLY: molecule_list_create,&
46                                              molecule_list_release,&
47                                              molecule_list_type
48   USE molecule_types,                  ONLY: molecule_type
49   USE particle_list_types,             ONLY: particle_list_create,&
50                                              particle_list_release,&
51                                              particle_list_type
52   USE particle_types,                  ONLY: particle_type
53   USE qmmm_types_low,                  ONLY: qmmm_env_mm_type
54   USE string_table,                    ONLY: id2str,&
55                                              s2s,&
56                                              str2id
57   USE topology,                        ONLY: connectivity_control,&
58                                              topology_control
59   USE topology_connectivity_util,      ONLY: topology_connectivity_pack
60   USE topology_coordinate_util,        ONLY: topology_coordinate_pack
61   USE topology_types,                  ONLY: deallocate_topology,&
62                                              init_topology,&
63                                              topology_parameters_type
64   USE topology_util,                   ONLY: check_subsys_element
65   USE virial_types,                    ONLY: virial_create,&
66                                              virial_set
67#include "./base/base_uses.f90"
68
69   IMPLICIT NONE
70   PRIVATE
71
72   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
73   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_subsys_methods'
74
75   PUBLIC :: create_small_subsys, cp_subsys_create
76
77CONTAINS
78
79! **************************************************************************************************
80!> \brief Creates allocates and fills subsys from given input.
81!> \param subsys ...
82!> \param para_env ...
83!> \param root_section ...
84!> \param force_env_section ...
85!> \param subsys_section ...
86!> \param use_motion_section ...
87!> \param qmmm ...
88!> \param qmmm_env ...
89!> \param exclusions ...
90!> \param elkind ...
91!> \author Ole Schuett
92! **************************************************************************************************
93   SUBROUTINE cp_subsys_create(subsys, para_env, &
94                               root_section, force_env_section, subsys_section, &
95                               use_motion_section, qmmm, qmmm_env, exclusions, elkind)
96      TYPE(cp_subsys_type), POINTER                      :: subsys
97      TYPE(cp_para_env_type), POINTER                    :: para_env
98      TYPE(section_vals_type), POINTER                   :: root_section
99      TYPE(section_vals_type), OPTIONAL, POINTER         :: force_env_section, subsys_section
100      LOGICAL, INTENT(IN), OPTIONAL                      :: use_motion_section
101      LOGICAL, OPTIONAL                                  :: qmmm
102      TYPE(qmmm_env_mm_type), OPTIONAL, POINTER          :: qmmm_env
103      TYPE(exclusion_type), DIMENSION(:), OPTIONAL, &
104         POINTER                                         :: exclusions
105      LOGICAL, INTENT(IN), OPTIONAL                      :: elkind
106
107      CHARACTER(len=*), PARAMETER :: routineN = 'cp_subsys_create', &
108         routineP = moduleN//':'//routineN
109
110      INTEGER                                            :: stress_tensor
111      INTEGER, DIMENSION(:), POINTER                     :: seed_vals
112      LOGICAL                                            :: atomic_energy, atomic_stress, &
113                                                            my_use_motion_section, &
114                                                            pv_availability, pv_diagonal, &
115                                                            pv_numerical
116      TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
117      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
118      TYPE(molecule_kind_list_type), POINTER             :: mol_kinds
119      TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
120      TYPE(molecule_list_type), POINTER                  :: mols
121      TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
122      TYPE(particle_list_type), POINTER                  :: particles
123      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
124      TYPE(section_vals_type), POINTER                   :: colvar_section, my_force_env_section, &
125                                                            my_subsys_section
126
127      CPASSERT(.NOT. ASSOCIATED(subsys))
128      ALLOCATE (subsys)
129
130      CALL cp_para_env_retain(para_env)
131      subsys%para_env => para_env
132
133      my_use_motion_section = .FALSE.
134      IF (PRESENT(use_motion_section)) &
135         my_use_motion_section = use_motion_section
136
137      my_force_env_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
138      IF (PRESENT(force_env_section)) &
139         my_force_env_section => force_env_section
140
141      my_subsys_section => section_vals_get_subs_vals(my_force_env_section, "SUBSYS")
142      IF (PRESENT(subsys_section)) &
143         my_subsys_section => subsys_section
144
145      CALL section_vals_val_get(my_subsys_section, "SEED", i_vals=seed_vals)
146      IF (SIZE(seed_vals) == 1) THEN
147         subsys%seed(:, :) = REAL(seed_vals(1), KIND=dp)
148      ELSE IF (SIZE(seed_vals) == 6) THEN
149         subsys%seed(1:3, 1:2) = RESHAPE(REAL(seed_vals(:), KIND=dp), (/3, 2/))
150      ELSE
151         CPABORT("Supply exactly 1 or 6 arguments for SEED in &SUBSYS only!")
152      END IF
153
154      colvar_section => section_vals_get_subs_vals(my_subsys_section, "COLVAR")
155
156      CALL cp_subsys_read_colvar(subsys, colvar_section)
157
158      !   *** Read the particle coordinates and allocate the atomic kind, ***
159      !   *** the molecule kind, and the molecule data structures         ***
160      CALL topology_control(atomic_kind_set, particle_set, molecule_kind_set, molecule_set, &
161                            subsys%colvar_p, subsys%gci, root_section, para_env, &
162                            force_env_section=my_force_env_section, &
163                            subsys_section=my_subsys_section, use_motion_section=my_use_motion_section, &
164                            qmmm=qmmm, qmmm_env=qmmm_env, exclusions=exclusions, elkind=elkind)
165
166      CALL particle_list_create(particles, els_ptr=particle_set)
167      CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
168      CALL molecule_list_create(mols, els_ptr=molecule_set)
169      CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
170
171      CALL cp_subsys_set(subsys, particles=particles, atomic_kinds=atomic_kinds, &
172                         molecules=mols, molecule_kinds=mol_kinds)
173
174      CALL particle_list_release(particles)
175      CALL atomic_kind_list_release(atomic_kinds)
176      CALL molecule_list_release(mols)
177      CALL molecule_kind_list_release(mol_kinds)
178
179      ! Should we compute the virial?
180      CALL section_vals_val_get(my_force_env_section, "STRESS_TENSOR", i_val=stress_tensor)
181      SELECT CASE (stress_tensor)
182      CASE (do_stress_none)
183         pv_availability = .FALSE.
184         pv_numerical = .FALSE.
185         pv_diagonal = .FALSE.
186      CASE (do_stress_analytical)
187         pv_availability = .TRUE.
188         pv_numerical = .FALSE.
189         pv_diagonal = .FALSE.
190      CASE (do_stress_numerical)
191         pv_availability = .TRUE.
192         pv_numerical = .TRUE.
193         pv_diagonal = .FALSE.
194      CASE (do_stress_diagonal_anal)
195         pv_availability = .TRUE.
196         pv_numerical = .FALSE.
197         pv_diagonal = .TRUE.
198      CASE (do_stress_diagonal_numer)
199         pv_availability = .TRUE.
200         pv_numerical = .TRUE.
201         pv_diagonal = .TRUE.
202      END SELECT
203
204      CALL virial_create(subsys%virial)
205      CALL virial_set(virial=subsys%virial, &
206                      pv_availability=pv_availability, &
207                      pv_numer=pv_numerical, &
208                      pv_diagonal=pv_diagonal)
209
210      ! Should we compute atomic properties?
211      CALL atprop_create(subsys%atprop)
212      CALL section_vals_val_get(my_force_env_section, "PROPERTIES%ATOMIC%ENERGY", l_val=atomic_energy)
213      subsys%atprop%energy = atomic_energy
214      CALL section_vals_val_get(my_force_env_section, "PROPERTIES%ATOMIC%PRESSURE", l_val=atomic_stress)
215      IF (atomic_stress) THEN
216         CPASSERT(pv_availability)
217         CPASSERT(.NOT. pv_numerical)
218      END IF
219      subsys%atprop%stress = atomic_stress
220
221      CALL cp_result_create(subsys%results)
222   END SUBROUTINE cp_subsys_create
223
224! **************************************************************************************************
225!> \brief reads the colvar section of the colvar
226!> \param subsys ...
227!> \param colvar_section ...
228!> \par History
229!>      2006.01 Joost VandeVondele
230! **************************************************************************************************
231   SUBROUTINE cp_subsys_read_colvar(subsys, colvar_section)
232      TYPE(cp_subsys_type), POINTER                      :: subsys
233      TYPE(section_vals_type), POINTER                   :: colvar_section
234
235      CHARACTER(len=*), PARAMETER :: routineN = 'cp_subsys_read_colvar', &
236         routineP = moduleN//':'//routineN
237
238      INTEGER                                            :: ig, ncol
239
240      CALL section_vals_get(colvar_section, n_repetition=ncol)
241      ALLOCATE (subsys%colvar_p(ncol))
242      DO ig = 1, ncol
243         NULLIFY (subsys%colvar_p(ig)%colvar)
244         CALL colvar_read(subsys%colvar_p(ig)%colvar, ig, colvar_section, subsys%para_env)
245      ENDDO
246   END SUBROUTINE cp_subsys_read_colvar
247
248! **************************************************************************************************
249!> \brief updates the molecule information of the given subsys
250!> \param small_subsys the subsys to create
251!> \param big_subsys the superset of small_subsys
252!> \param small_cell ...
253!> \param small_para_env the parallel environment for the new (small)
254!>        subsys
255!> \param sub_atom_index indexes of the atoms that should be in small_subsys
256!> \param sub_atom_kind_name ...
257!> \param para_env ...
258!> \param force_env_section ...
259!> \param subsys_section ...
260!> \param ignore_outside_box ...
261!> \par History
262!>      05.2004 created [fawzi]
263!> \author Fawzi Mohamed, Teodoro Laino
264!> \note
265!>      not really ready to be used with different para_envs for the small
266!>      and big part
267! **************************************************************************************************
268   SUBROUTINE create_small_subsys(small_subsys, big_subsys, small_cell, &
269                                  small_para_env, sub_atom_index, sub_atom_kind_name, &
270                                  para_env, force_env_section, subsys_section, ignore_outside_box)
271
272      TYPE(cp_subsys_type), POINTER                      :: small_subsys, big_subsys
273      TYPE(cell_type), POINTER                           :: small_cell
274      TYPE(cp_para_env_type), POINTER                    :: small_para_env
275      INTEGER, DIMENSION(:), INTENT(in)                  :: sub_atom_index
276      CHARACTER(len=default_string_length), &
277         DIMENSION(:), INTENT(in)                        :: sub_atom_kind_name
278      TYPE(cp_para_env_type), POINTER                    :: para_env
279      TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
280      LOGICAL, INTENT(in), OPTIONAL                      :: ignore_outside_box
281
282      CHARACTER(len=*), PARAMETER :: routineN = 'create_small_subsys', &
283         routineP = moduleN//':'//routineN
284
285      CHARACTER(len=default_string_length)               :: my_element, strtmp1
286      INTEGER                                            :: iat, id_, nat
287      TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
288      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
289      TYPE(molecule_kind_list_type), POINTER             :: mol_kinds
290      TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
291      TYPE(molecule_list_type), POINTER                  :: mols
292      TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
293      TYPE(particle_list_type), POINTER                  :: particles
294      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
295      TYPE(topology_parameters_type)                     :: topology
296
297      NULLIFY (mol_kinds, mols, particles, atomic_kinds, atomic_kind_set, particle_set, &
298               molecule_kind_set, molecule_set, particles, atomic_kinds)
299
300      CPASSERT(.NOT. ASSOCIATED(small_subsys))
301      CPASSERT(ASSOCIATED(big_subsys))
302      IF (big_subsys%para_env%group /= small_para_env%group) &
303         CPABORT("big_subsys%para_env%group==small_para_env%group")
304
305      !-----------------------------------------------------------------------------
306      !-----------------------------------------------------------------------------
307      ! 1. Initialize the topology structure type
308      !-----------------------------------------------------------------------------
309      CALL init_topology(topology)
310
311      !-----------------------------------------------------------------------------
312      !-----------------------------------------------------------------------------
313      ! 2. Get the cell info
314      !-----------------------------------------------------------------------------
315      topology%cell => small_cell
316      CALL cell_retain(small_cell)
317
318      !-----------------------------------------------------------------------------
319      !-----------------------------------------------------------------------------
320      ! 3. Initialize atom coords from the bigger system
321      !-----------------------------------------------------------------------------
322      nat = SIZE(sub_atom_index)
323      topology%natoms = nat
324      CPASSERT(.NOT. ASSOCIATED(topology%atom_info%r))
325      CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_atmname))
326      CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_molname))
327      CPASSERT(.NOT. ASSOCIATED(topology%atom_info%id_resname))
328      CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_mass))
329      CPASSERT(.NOT. ASSOCIATED(topology%atom_info%atm_charge))
330      ALLOCATE (topology%atom_info%r(3, nat), topology%atom_info%id_atmname(nat), &
331                topology%atom_info%id_molname(nat), topology%atom_info%id_resname(nat), &
332                topology%atom_info%id_element(nat), topology%atom_info%atm_mass(nat), &
333                topology%atom_info%atm_charge(nat))
334
335      CALL cp_subsys_get(big_subsys, particles=particles)
336      DO iat = 1, nat
337         topology%atom_info%r(:, iat) = particles%els(sub_atom_index(iat))%r
338         topology%atom_info%id_atmname(iat) = str2id(s2s(sub_atom_kind_name(iat)))
339         topology%atom_info%id_molname(iat) = topology%atom_info%id_atmname(iat)
340         topology%atom_info%id_resname(iat) = topology%atom_info%id_atmname(iat)
341         !
342         ! Defining element
343         !
344         id_ = INDEX(id2str(topology%atom_info%id_atmname(iat)), "_") - 1
345         IF (id_ == -1) id_ = LEN_TRIM(id2str(topology%atom_info%id_atmname(iat)))
346         strtmp1 = id2str(topology%atom_info%id_atmname(iat))
347         strtmp1 = strtmp1(1:id_)
348         CALL check_subsys_element(strtmp1, strtmp1, my_element, &
349                                   subsys_section, use_mm_map_first=.FALSE.)
350         topology%atom_info%id_element(iat) = str2id(s2s(my_element))
351         topology%atom_info%atm_mass(iat) = 0._dp
352         topology%atom_info%atm_charge(iat) = 0._dp
353      END DO
354      topology%conn_type = do_conn_off
355
356      !-----------------------------------------------------------------------------
357      !-----------------------------------------------------------------------------
358      ! 4. Read in or generate the molecular connectivity
359      !-----------------------------------------------------------------------------
360      CALL connectivity_control(topology, para_env, subsys_section=subsys_section, &
361                                force_env_section=force_env_section)
362
363      !-----------------------------------------------------------------------------
364      !-----------------------------------------------------------------------------
365      ! 5. Pack everything into the molecular types
366      !-----------------------------------------------------------------------------
367      CALL topology_connectivity_pack(molecule_kind_set, molecule_set, &
368                                      topology, subsys_section=subsys_section)
369
370      !-----------------------------------------------------------------------------
371      !-----------------------------------------------------------------------------
372      ! 6. Pack everything into the atomic types
373      !-----------------------------------------------------------------------------
374      CALL topology_coordinate_pack(particle_set, atomic_kind_set, &
375                                    molecule_kind_set, molecule_set, topology, subsys_section=subsys_section, &
376                                    force_env_section=force_env_section, ignore_outside_box=ignore_outside_box)
377
378      !-----------------------------------------------------------------------------
379      !-----------------------------------------------------------------------------
380      ! 7. Cleanup the topology structure type
381      !-----------------------------------------------------------------------------
382      CALL deallocate_topology(topology)
383
384      !-----------------------------------------------------------------------------
385      !-----------------------------------------------------------------------------
386      ! 8. Allocate new subsys
387      !-----------------------------------------------------------------------------
388      ALLOCATE (small_subsys)
389      CALL cp_para_env_retain(para_env)
390      small_subsys%para_env => para_env
391      CALL particle_list_create(particles, els_ptr=particle_set)
392      CALL atomic_kind_list_create(atomic_kinds, els_ptr=atomic_kind_set)
393      CALL molecule_list_create(mols, els_ptr=molecule_set)
394      CALL molecule_kind_list_create(mol_kinds, els_ptr=molecule_kind_set)
395      CALL cp_subsys_set(small_subsys, particles=particles, atomic_kinds=atomic_kinds, &
396                         molecules=mols, molecule_kinds=mol_kinds)
397      CALL particle_list_release(particles)
398      CALL atomic_kind_list_release(atomic_kinds)
399      CALL molecule_list_release(mols)
400      CALL molecule_kind_list_release(mol_kinds)
401
402      CALL virial_create(small_subsys%virial)
403      CALL atprop_create(small_subsys%atprop)
404      CALL cp_result_create(small_subsys%results)
405   END SUBROUTINE create_small_subsys
406
407END MODULE cp_subsys_methods
408