1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Util mixed_environment
8!> \author Teodoro Laino [tlaino] - 02.2011
9! **************************************************************************************************
10MODULE mixed_environment_utils
11
12   USE cp_result_methods,               ONLY: cp_results_erase,&
13                                              get_results,&
14                                              put_results,&
15                                              test_for_result
16   USE cp_result_types,                 ONLY: cp_result_p_type,&
17                                              cp_result_type
18   USE input_section_types,             ONLY: section_vals_get,&
19                                              section_vals_get_subs_vals,&
20                                              section_vals_type,&
21                                              section_vals_val_get
22   USE kinds,                           ONLY: default_string_length,&
23                                              dp
24   USE mixed_energy_types,              ONLY: mixed_force_type
25   USE particle_list_types,             ONLY: particle_list_type
26   USE virial_types,                    ONLY: virial_p_type,&
27                                              virial_type,&
28                                              zero_virial
29#include "./base/base_uses.f90"
30
31   IMPLICIT NONE
32
33   PRIVATE
34
35   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mixed_environment_utils'
36
37   PUBLIC :: mixed_map_forces, &
38             get_subsys_map_index
39
40CONTAINS
41
42! **************************************************************************************************
43!> \brief Maps forces between the different force_eval sections/environments
44!> \param particles_mix ...
45!> \param virial_mix ...
46!> \param results_mix ...
47!> \param global_forces ...
48!> \param virials ...
49!> \param results ...
50!> \param factor ...
51!> \param iforce_eval ...
52!> \param nforce_eval ...
53!> \param map_index ...
54!> \param mapping_section ...
55!> \param overwrite ...
56!> \author Teodoro Laino - University of Zurich [tlaino] - 05.2007
57! **************************************************************************************************
58   SUBROUTINE mixed_map_forces(particles_mix, virial_mix, results_mix, global_forces, &
59                               virials, results, factor, iforce_eval, nforce_eval, map_index, &
60                               mapping_section, overwrite)
61
62      TYPE(particle_list_type), POINTER                  :: particles_mix
63      TYPE(virial_type), POINTER                         :: virial_mix
64      TYPE(cp_result_type), POINTER                      :: results_mix
65      TYPE(mixed_force_type), DIMENSION(:), POINTER      :: global_forces
66      TYPE(virial_p_type), DIMENSION(:), POINTER         :: virials
67      TYPE(cp_result_p_type), DIMENSION(:), POINTER      :: results
68      REAL(KIND=dp), INTENT(IN)                          :: factor
69      INTEGER, INTENT(IN)                                :: iforce_eval, nforce_eval
70      INTEGER, DIMENSION(:), POINTER                     :: map_index
71      TYPE(section_vals_type), POINTER                   :: mapping_section
72      LOGICAL, INTENT(IN)                                :: overwrite
73
74      CHARACTER(len=*), PARAMETER :: routineN = 'mixed_map_forces', &
75         routineP = moduleN//':'//routineN
76
77      CHARACTER(LEN=default_string_length)               :: description
78      INTEGER                                            :: iparticle, jparticle, natom, nres
79      LOGICAL                                            :: dip_exists
80      REAL(KIND=dp), DIMENSION(3)                        :: dip_mix, dip_tmp
81
82! Get Mapping index array
83
84      natom = SIZE(global_forces(iforce_eval)%forces, 2)
85      CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, map_index)
86      DO iparticle = 1, natom
87         jparticle = map_index(iparticle)
88         IF (overwrite) THEN
89            particles_mix%els(jparticle)%f(:) = factor*global_forces(iforce_eval)%forces(:, iparticle)
90         ELSE
91            particles_mix%els(jparticle)%f(:) = particles_mix%els(jparticle)%f(:) + &
92                                                factor*global_forces(iforce_eval)%forces(:, iparticle)
93         END IF
94      END DO
95      ! Mixing Virial
96      IF (virial_mix%pv_availability) THEN
97         IF (overwrite) CALL zero_virial(virial_mix, reset=.FALSE.)
98         virial_mix%pv_total = virial_mix%pv_total + factor*virials(iforce_eval)%virial%pv_total
99         virial_mix%pv_kinetic = virial_mix%pv_kinetic + factor*virials(iforce_eval)%virial%pv_kinetic
100         virial_mix%pv_virial = virial_mix%pv_virial + factor*virials(iforce_eval)%virial%pv_virial
101         virial_mix%pv_xc = virial_mix%pv_xc + factor*virials(iforce_eval)%virial%pv_xc
102         virial_mix%pv_fock_4c = virial_mix%pv_fock_4c + factor*virials(iforce_eval)%virial%pv_fock_4c
103         virial_mix%pv_constraint = virial_mix%pv_constraint + factor*virials(iforce_eval)%virial%pv_constraint
104      END IF
105      ! Deallocate map_index array
106      IF (ASSOCIATED(map_index)) THEN
107         DEALLOCATE (map_index)
108      END IF
109
110      ! Collect Requested Results info
111      description = '[DIPOLE]'
112      IF (overwrite) CALL cp_results_erase(results_mix)
113
114      dip_exists = test_for_result(results=results(iforce_eval)%results, description=description)
115      IF (dip_exists) THEN
116         CALL get_results(results=results_mix, description=description, n_rep=nres)
117         CPASSERT(nres <= 1)
118         dip_mix = 0.0_dp
119         IF (nres == 1) CALL get_results(results=results_mix, description=description, values=dip_mix)
120         CALL get_results(results=results(iforce_eval)%results, description=description, n_rep=nres)
121         CALL get_results(results=results(iforce_eval)%results, description=description, &
122                          values=dip_tmp, nval=nres)
123         dip_mix = dip_mix + factor*dip_tmp
124         CALL cp_results_erase(results=results_mix, description=description)
125         CALL put_results(results=results_mix, description=description, values=dip_mix)
126      END IF
127
128   END SUBROUTINE mixed_map_forces
129
130! **************************************************************************************************
131!> \brief performs mapping of the subsystems of different force_eval
132!> \param mapping_section ...
133!> \param natom ...
134!> \param iforce_eval ...
135!> \param nforce_eval ...
136!> \param map_index ...
137!> \param force_eval_embed ...
138!> \author Teodoro Laino - University of Zurich [tlaino] - 05.2007
139! **************************************************************************************************
140   SUBROUTINE get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, map_index, &
141                                   force_eval_embed)
142
143      TYPE(section_vals_type), POINTER                   :: mapping_section
144      INTEGER, INTENT(IN)                                :: natom, iforce_eval, nforce_eval
145      INTEGER, DIMENSION(:), POINTER                     :: map_index
146      LOGICAL, OPTIONAL                                  :: force_eval_embed
147
148      CHARACTER(len=*), PARAMETER :: routineN = 'get_subsys_map_index', &
149         routineP = moduleN//':'//routineN
150
151      INTEGER                                            :: i, iatom, ival, j, jval, k, n_rep, &
152                                                            n_rep_loc, n_rep_map, n_rep_sys, tmp
153      INTEGER, DIMENSION(:), POINTER                     :: index_glo, index_loc, list
154      LOGICAL                                            :: check, explicit
155      TYPE(section_vals_type), POINTER                   :: fragments_loc, fragments_sys, &
156                                                            map_force_ev, map_full_sys
157
158      CPASSERT(.NOT. ASSOCIATED(map_index))
159      ALLOCATE (map_index(natom))
160      CALL section_vals_get(mapping_section, explicit=explicit)
161      IF (.NOT. explicit) THEN
162         ! Standard Mapping.. subsys are assumed to have the same structure
163         DO i = 1, natom
164            map_index(i) = i
165         END DO
166      ELSE
167         ! Mapping systems with different structures
168         IF (.NOT. PRESENT(force_eval_embed)) THEN
169            map_full_sys => section_vals_get_subs_vals(mapping_section, "FORCE_EVAL_MIXED")
170         ELSE
171            map_full_sys => section_vals_get_subs_vals(mapping_section, "FORCE_EVAL_EMBED")
172         ENDIF
173         map_force_ev => section_vals_get_subs_vals(mapping_section, "FORCE_EVAL")
174         CALL section_vals_get(map_full_sys, explicit=explicit)
175         CPASSERT(explicit)
176         CALL section_vals_get(map_force_ev, explicit=explicit, n_repetition=n_rep)
177         CPASSERT(explicit)
178         CPASSERT(n_rep == nforce_eval)
179         DO i = 1, n_rep
180            CALL section_vals_val_get(map_force_ev, "_SECTION_PARAMETERS_", i_rep_section=i, i_val=ival)
181            IF (ival == iforce_eval) EXIT
182         END DO
183         CPASSERT(i <= nforce_eval)
184         fragments_sys => section_vals_get_subs_vals(map_full_sys, "FRAGMENT")
185         fragments_loc => section_vals_get_subs_vals(map_force_ev, "FRAGMENT", i_rep_section=i)
186         !Perform few check on the structure of the input mapping section. as provided by the user
187         CALL section_vals_get(fragments_loc, n_repetition=n_rep_loc)
188         CALL section_vals_get(fragments_sys, explicit=explicit, n_repetition=n_rep_sys)
189         CPASSERT(explicit)
190         CPASSERT(n_rep_sys >= n_rep_loc)
191         IF (n_rep_loc == 0) THEN
192            NULLIFY (list)
193            ! We expect an easier syntax in this case..
194            CALL section_vals_val_get(map_force_ev, "DEFINE_FRAGMENTS", i_rep_section=i, n_rep_val=n_rep_map)
195            check = (n_rep_map /= 0)
196            CPASSERT(check)
197            CALL section_vals_val_get(map_force_ev, "DEFINE_FRAGMENTS", i_rep_section=i, i_vals=list)
198            CPASSERT(SIZE(list) > 0)
199            iatom = 0
200            DO i = 1, SIZE(list)
201               jval = list(i)
202               DO j = 1, n_rep_sys
203                  CALL section_vals_val_get(fragments_sys, "_SECTION_PARAMETERS_", i_rep_section=j, i_val=tmp)
204                  IF (tmp == jval) EXIT
205               END DO
206               CALL section_vals_val_get(fragments_sys, "_DEFAULT_KEYWORD_", i_rep_section=j, i_vals=index_glo)
207               DO k = 0, index_glo(2) - index_glo(1)
208                  iatom = iatom + 1
209                  CPASSERT(iatom <= natom)
210                  map_index(iatom) = index_glo(1) + k
211               END DO
212            END DO
213            check = (iatom == natom)
214            CPASSERT(check)
215         ELSE
216            ! General syntax..
217            !Loop over the fragment of the force_eval
218            DO i = 1, n_rep_loc
219               CALL section_vals_val_get(fragments_loc, "_SECTION_PARAMETERS_", i_rep_section=i, i_val=ival)
220               CALL section_vals_val_get(fragments_loc, "MAP", i_rep_section=i, i_val=jval)
221               ! Index corresponding to the mixed_force_eval fragment
222               DO j = 1, n_rep_sys
223                  CALL section_vals_val_get(fragments_sys, "_SECTION_PARAMETERS_", i_rep_section=j, i_val=tmp)
224                  IF (tmp == jval) EXIT
225               END DO
226               CPASSERT(j <= n_rep_sys)
227               CALL section_vals_val_get(fragments_loc, "_DEFAULT_KEYWORD_", i_rep_section=i, i_vals=index_loc)
228               CALL section_vals_val_get(fragments_sys, "_DEFAULT_KEYWORD_", i_rep_section=j, i_vals=index_glo)
229               check = ((index_loc(2) - index_loc(1)) == (index_glo(2) - index_glo(1)))
230               CPASSERT(check)
231               ! Now let's build the real mapping
232               DO k = 0, index_loc(2) - index_loc(1)
233                  map_index(index_loc(1) + k) = index_glo(1) + k
234               END DO
235            END DO
236         END IF
237      END IF
238
239   END SUBROUTINE get_subsys_map_index
240
241END MODULE mixed_environment_utils
242