1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \par History
8!>      09.2004 created [tlaino]
9!> \author Teodoro Laino
10! **************************************************************************************************
11MODULE qmmm_util
12   USE cell_types,                      ONLY: cell_type
13   USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
14   USE cp_subsys_types,                 ONLY: cp_subsys_type
15   USE fist_environment_types,          ONLY: fist_env_get
16   USE force_env_types,                 ONLY: force_env_type,&
17                                              use_qmmm,&
18                                              use_qmmmx
19   USE input_constants,                 ONLY: do_qmmm_wall_none,&
20                                              do_qmmm_wall_quadratic,&
21                                              do_qmmm_wall_reflective
22   USE input_section_types,             ONLY: section_vals_get,&
23                                              section_vals_get_subs_vals,&
24                                              section_vals_type,&
25                                              section_vals_val_get
26   USE kinds,                           ONLY: dp
27   USE mathconstants,                   ONLY: pi
28   USE particle_methods,                ONLY: write_fist_particle_coordinates,&
29                                              write_qs_particle_coordinates
30   USE particle_types,                  ONLY: particle_type
31   USE qmmm_types,                      ONLY: qmmm_env_type
32   USE qs_energy_types,                 ONLY: qs_energy_type
33   USE qs_environment_types,            ONLY: get_qs_env
34   USE qs_kind_types,                   ONLY: qs_kind_type
35#include "./base/base_uses.f90"
36
37   IMPLICIT NONE
38   PRIVATE
39
40   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
41   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_util'
42   PUBLIC :: apply_qmmm_walls_reflective, &
43             apply_qmmm_walls, &
44             apply_qmmm_translate, &
45             apply_qmmm_wrap, &
46             apply_qmmm_unwrap, &
47             spherical_cutoff_factor
48
49CONTAINS
50
51! **************************************************************************************************
52!> \brief Apply QM quadratic walls in order to avoid QM atoms escaping from
53!>      the QM Box
54!> \param qmmm_env ...
55!> \par History
56!>      02.2008 created
57!> \author Benjamin G Levine
58! **************************************************************************************************
59   SUBROUTINE apply_qmmm_walls(qmmm_env)
60      TYPE(qmmm_env_type), POINTER                       :: qmmm_env
61
62      CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_walls', &
63         routineP = moduleN//':'//routineN
64
65      INTEGER                                            :: iwall_type
66      LOGICAL                                            :: do_qmmm_force_mixing, explicit
67      TYPE(section_vals_type), POINTER                   :: qmmmx_section, walls_section
68
69      walls_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "QMMM%WALLS")
70      qmmmx_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "QMMM%FORCE_MIXING")
71      CALL section_vals_get(qmmmx_section, explicit=do_qmmm_force_mixing)
72      CALL section_vals_get(walls_section, explicit=explicit)
73      IF (explicit) THEN
74         CALL section_vals_val_get(walls_section, "TYPE", i_val=iwall_type)
75         SELECT CASE (iwall_type)
76         CASE (do_qmmm_wall_quadratic)
77            IF (do_qmmm_force_mixing) THEN
78               CALL cp_warn(__LOCATION__, &
79                            "Quadratic walls for QM/MM are not implemented (or useful), when "// &
80                            "force mixing is active.  Skipping!")
81            ELSE
82               CALL apply_qmmm_walls_quadratic(qmmm_env, walls_section)
83            ENDIF
84         CASE (do_qmmm_wall_reflective)
85            ! Do nothing.. reflective walls are applied directly in the integrator
86         END SELECT
87      ENDIF
88
89   END SUBROUTINE apply_qmmm_walls
90
91! **************************************************************************************************
92!> \brief Apply reflective QM walls in order to avoid QM atoms escaping from
93!>      the QM Box
94!> \param force_env ...
95!> \par History
96!>      08.2007 created [tlaino] - Zurich University
97!> \author Teodoro Laino
98! **************************************************************************************************
99   SUBROUTINE apply_qmmm_walls_reflective(force_env)
100      TYPE(force_env_type), POINTER                      :: force_env
101
102      CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_walls_reflective', &
103         routineP = moduleN//':'//routineN
104
105      INTEGER                                            :: ip, iwall_type, qm_index
106      INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index
107      LOGICAL                                            :: explicit, is_x(2), is_y(2), is_z(2)
108      REAL(KIND=dp), DIMENSION(3)                        :: coord, qm_cell_diag, skin
109      REAL(KIND=dp), DIMENSION(:), POINTER               :: list
110      TYPE(cell_type), POINTER                           :: mm_cell, qm_cell
111      TYPE(cp_subsys_type), POINTER                      :: subsys_mm, subsys_qm
112      TYPE(particle_type), DIMENSION(:), POINTER         :: particles_mm
113      TYPE(section_vals_type), POINTER                   :: walls_section
114
115      NULLIFY (subsys_mm, subsys_qm, qm_atom_index, particles_mm, qm_cell, mm_cell, &
116               walls_section)
117
118      IF (force_env%in_use /= use_qmmm .AND. force_env%in_use /= use_qmmmx) RETURN
119
120      walls_section => section_vals_get_subs_vals(force_env%root_section, "FORCE_EVAL%QMMM%WALLS")
121      CALL section_vals_get(walls_section, explicit=explicit)
122      IF (explicit) THEN
123         NULLIFY (list)
124         CALL section_vals_val_get(walls_section, "WALL_SKIN", r_vals=list)
125         CALL section_vals_val_get(walls_section, "TYPE", i_val=iwall_type)
126         skin(:) = list(:)
127      ELSE
128         ![NB]
129         iwall_type = do_qmmm_wall_reflective
130         skin(:) = 0.0_dp
131      END IF
132
133      IF (force_env%in_use == use_qmmmx) THEN
134         IF (iwall_type /= do_qmmm_wall_none) &
135            CALL cp_warn(__LOCATION__, &
136                         "Reflective walls for QM/MM are not implemented (or useful) when "// &
137                         "force mixing is active.  Skipping!")
138         RETURN
139      END IF
140
141      ! from here on we can be sure that it's conventional QM/MM
142      CPASSERT(ASSOCIATED(force_env%qmmm_env))
143      CPASSERT(force_env%qmmm_env%ref_count > 0)
144
145      CALL fist_env_get(force_env%qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
146      CALL get_qs_env(force_env%qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm)
147      qm_atom_index => force_env%qmmm_env%qm%qm_atom_index
148      CPASSERT(ASSOCIATED(qm_atom_index))
149
150      qm_cell_diag = (/qm_cell%hmat(1, 1), &
151                       qm_cell%hmat(2, 2), &
152                       qm_cell%hmat(3, 3)/)
153      particles_mm => subsys_mm%particles%els
154      DO ip = 1, SIZE(qm_atom_index)
155         qm_index = qm_atom_index(ip)
156         coord = particles_mm(qm_index)%r
157         IF (ANY(coord < skin) .OR. ANY(coord > (qm_cell_diag - skin))) THEN
158            IF (explicit) THEN
159               IF (iwall_type == do_qmmm_wall_reflective) THEN
160                  ! Apply Walls
161                  is_x(1) = (coord(1) < skin(1))
162                  is_x(2) = (coord(1) > (qm_cell_diag(1) - skin(1)))
163                  is_y(1) = (coord(2) < skin(2))
164                  is_y(2) = (coord(2) > (qm_cell_diag(2) - skin(2)))
165                  is_z(1) = (coord(3) < skin(3))
166                  is_z(2) = (coord(3) > (qm_cell_diag(3) - skin(3)))
167                  IF (ANY(is_x)) THEN
168                     ! X coordinate
169                     IF (is_x(1)) THEN
170                        particles_mm(qm_index)%v(1) = ABS(particles_mm(qm_index)%v(1))
171                     ELSE IF (is_x(2)) THEN
172                        particles_mm(qm_index)%v(1) = -ABS(particles_mm(qm_index)%v(1))
173                     END IF
174                  END IF
175                  IF (ANY(is_y)) THEN
176                     ! Y coordinate
177                     IF (is_y(1)) THEN
178                        particles_mm(qm_index)%v(2) = ABS(particles_mm(qm_index)%v(2))
179                     ELSE IF (is_y(2)) THEN
180                        particles_mm(qm_index)%v(2) = -ABS(particles_mm(qm_index)%v(2))
181                     END IF
182                  END IF
183                  IF (ANY(is_z)) THEN
184                     ! Z coordinate
185                     IF (is_z(1)) THEN
186                        particles_mm(qm_index)%v(3) = ABS(particles_mm(qm_index)%v(3))
187                     ELSE IF (is_z(2)) THEN
188                        particles_mm(qm_index)%v(3) = -ABS(particles_mm(qm_index)%v(3))
189                     END IF
190                  END IF
191               ENDIF
192            ELSE
193               ! Otherwise print a warning and continue crossing cp2k's finger..
194               CALL cp_warn(__LOCATION__, &
195                            "One or few QM atoms are within the SKIN of the quantum box. Check your run "// &
196                            "and you may possibly consider: the activation of the QMMM WALLS "// &
197                            "around the QM box, switching ON the centering of the QM box or increase "// &
198                            "the size of the QM cell. CP2K CONTINUE but results could be meaningless. ")
199            END IF
200         END IF
201      END DO
202
203   END SUBROUTINE apply_qmmm_walls_reflective
204
205! **************************************************************************************************
206!> \brief Apply QM quadratic walls in order to avoid QM atoms escaping from
207!>      the QM Box
208!> \param qmmm_env ...
209!> \param walls_section ...
210!> \par History
211!>      02.2008 created
212!> \author Benjamin G Levine
213! **************************************************************************************************
214   SUBROUTINE apply_qmmm_walls_quadratic(qmmm_env, walls_section)
215      TYPE(qmmm_env_type), POINTER                       :: qmmm_env
216      TYPE(section_vals_type), POINTER                   :: walls_section
217
218      CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_walls_quadratic', &
219         routineP = moduleN//':'//routineN
220
221      INTEGER                                            :: ip, qm_index
222      INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index
223      LOGICAL                                            :: is_x(2), is_y(2), is_z(2)
224      REAL(KIND=dp)                                      :: k, wallenergy, wallforce
225      REAL(KIND=dp), DIMENSION(3)                        :: coord, qm_cell_diag, skin
226      REAL(KIND=dp), DIMENSION(:), POINTER               :: list
227      TYPE(cell_type), POINTER                           :: mm_cell, qm_cell
228      TYPE(cp_subsys_type), POINTER                      :: subsys_mm, subsys_qm
229      TYPE(particle_type), DIMENSION(:), POINTER         :: particles_mm
230      TYPE(qs_energy_type), POINTER                      :: energy
231
232      NULLIFY (list)
233      CALL section_vals_val_get(walls_section, "WALL_SKIN", r_vals=list)
234      CALL section_vals_val_get(walls_section, "K", r_val=k)
235      CPASSERT(ASSOCIATED(qmmm_env))
236      CPASSERT(qmmm_env%ref_count > 0)
237
238      CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
239      CALL get_qs_env(qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm)
240
241      qm_atom_index => qmmm_env%qm%qm_atom_index
242      CPASSERT(ASSOCIATED(qm_atom_index))
243
244      skin(:) = list(:)
245
246      qm_cell_diag = (/qm_cell%hmat(1, 1), &
247                       qm_cell%hmat(2, 2), &
248                       qm_cell%hmat(3, 3)/)
249      particles_mm => subsys_mm%particles%els
250      wallenergy = 0.0_dp
251      DO ip = 1, SIZE(qm_atom_index)
252         qm_index = qm_atom_index(ip)
253         coord = particles_mm(qm_index)%r
254         IF (ANY(coord < skin) .OR. ANY(coord > (qm_cell_diag - skin))) THEN
255            is_x(1) = (coord(1) < skin(1))
256            is_x(2) = (coord(1) > (qm_cell_diag(1) - skin(1)))
257            is_y(1) = (coord(2) < skin(2))
258            is_y(2) = (coord(2) > (qm_cell_diag(2) - skin(2)))
259            is_z(1) = (coord(3) < skin(3))
260            is_z(2) = (coord(3) > (qm_cell_diag(3) - skin(3)))
261            IF (is_x(1)) THEN
262               wallforce = 2.0_dp*k*(skin(1) - coord(1))
263               particles_mm(qm_index)%f(1) = particles_mm(qm_index)%f(1) + &
264                                             wallforce
265               wallenergy = wallenergy + wallforce*(skin(1) - coord(1))*0.5_dp
266            ENDIF
267            IF (is_x(2)) THEN
268               wallforce = 2.0_dp*k*((qm_cell_diag(1) - skin(1)) - coord(1))
269               particles_mm(qm_index)%f(1) = particles_mm(qm_index)%f(1) + &
270                                             wallforce
271               wallenergy = wallenergy + wallforce*((qm_cell_diag(1) - skin(1)) - &
272                                                    coord(1))*0.5_dp
273            ENDIF
274            IF (is_y(1)) THEN
275               wallforce = 2.0_dp*k*(skin(2) - coord(2))
276               particles_mm(qm_index)%f(2) = particles_mm(qm_index)%f(2) + &
277                                             wallforce
278               wallenergy = wallenergy + wallforce*(skin(2) - coord(2))*0.5_dp
279            ENDIF
280            IF (is_y(2)) THEN
281               wallforce = 2.0_dp*k*((qm_cell_diag(2) - skin(2)) - coord(2))
282               particles_mm(qm_index)%f(2) = particles_mm(qm_index)%f(2) + &
283                                             wallforce
284               wallenergy = wallenergy + wallforce*((qm_cell_diag(2) - skin(2)) - &
285                                                    coord(2))*0.5_dp
286            ENDIF
287            IF (is_z(1)) THEN
288               wallforce = 2.0_dp*k*(skin(3) - coord(3))
289               particles_mm(qm_index)%f(3) = particles_mm(qm_index)%f(3) + &
290                                             wallforce
291               wallenergy = wallenergy + wallforce*(skin(3) - coord(3))*0.5_dp
292            ENDIF
293            IF (is_z(2)) THEN
294               wallforce = 2.0_dp*k*((qm_cell_diag(3) - skin(3)) - coord(3))
295               particles_mm(qm_index)%f(3) = particles_mm(qm_index)%f(3) + &
296                                             wallforce
297               wallenergy = wallenergy + wallforce*((qm_cell_diag(3) - skin(3)) - &
298                                                    coord(3))*0.5_dp
299            ENDIF
300         ENDIF
301      ENDDO
302
303      CALL get_qs_env(qs_env=qmmm_env%qs_env, energy=energy)
304      energy%total = energy%total + wallenergy
305
306   END SUBROUTINE apply_qmmm_walls_quadratic
307
308! **************************************************************************************************
309!> \brief wrap positions (with mm periodicity)
310!> \param subsys_mm ...
311!> \param mm_cell ...
312!> \param subsys_qm ...
313!> \param qm_atom_index ...
314!> \param saved_pos ...
315! **************************************************************************************************
316   SUBROUTINE apply_qmmm_wrap(subsys_mm, mm_cell, subsys_qm, qm_atom_index, saved_pos)
317      TYPE(cp_subsys_type), POINTER                      :: subsys_mm
318      TYPE(cell_type), POINTER                           :: mm_cell
319      TYPE(cp_subsys_type), OPTIONAL, POINTER            :: subsys_qm
320      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: qm_atom_index
321      REAL(dp), ALLOCATABLE                              :: saved_pos(:, :)
322
323      CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_wrap', &
324         routineP = moduleN//':'//routineN
325
326      INTEGER                                            :: i_dim, ip
327      REAL(dp)                                           :: r_lat(3)
328
329      ALLOCATE (saved_pos(3, subsys_mm%particles%n_els))
330      DO ip = 1, subsys_mm%particles%n_els
331         saved_pos(1:3, ip) = subsys_mm%particles%els(ip)%r(1:3)
332         r_lat = MATMUL(mm_cell%h_inv, subsys_mm%particles%els(ip)%r)
333         DO i_dim = 1, 3
334            IF (mm_cell%perd(i_dim) /= 1) THEN
335               r_lat(i_dim) = 0.0_dp
336            END IF
337         END DO
338         subsys_mm%particles%els(ip)%r = subsys_mm%particles%els(ip)%r - MATMUL(mm_cell%hmat, FLOOR(r_lat))
339      END DO
340
341      IF (PRESENT(subsys_qm) .AND. PRESENT(qm_atom_index)) THEN
342         DO ip = 1, SIZE(qm_atom_index)
343            subsys_qm%particles%els(ip)%r = subsys_mm%particles%els(qm_atom_index(ip))%r
344         END DO
345      END IF
346   END SUBROUTINE apply_qmmm_wrap
347
348! **************************************************************************************************
349!> \brief ...
350!> \param subsys_mm ...
351!> \param subsys_qm ...
352!> \param qm_atom_index ...
353!> \param saved_pos ...
354! **************************************************************************************************
355   SUBROUTINE apply_qmmm_unwrap(subsys_mm, subsys_qm, qm_atom_index, saved_pos)
356      TYPE(cp_subsys_type), POINTER                      :: subsys_mm
357      TYPE(cp_subsys_type), OPTIONAL, POINTER            :: subsys_qm
358      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: qm_atom_index
359      REAL(dp), ALLOCATABLE                              :: saved_pos(:, :)
360
361      CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_unwrap', &
362         routineP = moduleN//':'//routineN
363
364      INTEGER                                            :: ip
365
366      DO ip = 1, subsys_mm%particles%n_els
367         subsys_mm%particles%els(ip)%r(1:3) = saved_pos(1:3, ip)
368      END DO
369
370      IF (PRESENT(subsys_qm) .AND. PRESENT(qm_atom_index)) THEN
371         DO ip = 1, SIZE(qm_atom_index)
372            subsys_qm%particles%els(ip)%r = subsys_mm%particles%els(qm_atom_index(ip))%r
373         END DO
374      END IF
375
376      DEALLOCATE (saved_pos)
377   END SUBROUTINE apply_qmmm_unwrap
378
379! **************************************************************************************************
380!> \brief Apply translation to the full system in order to center the QM
381!>      system into the QM box
382!> \param qmmm_env ...
383!> \par History
384!>      08.2007 created [tlaino] - Zurich University
385!> \author Teodoro Laino
386! **************************************************************************************************
387   SUBROUTINE apply_qmmm_translate(qmmm_env)
388      TYPE(qmmm_env_type), POINTER                       :: qmmm_env
389
390      CHARACTER(len=*), PARAMETER :: routineN = 'apply_qmmm_translate', &
391         routineP = moduleN//':'//routineN
392
393      INTEGER                                            :: bigger_ip, i_dim, ip, max_ip, min_ip, &
394                                                            smaller_ip, tmp_ip, unit_nr
395      INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index
396      LOGICAL, ALLOCATABLE                               :: avoid(:)
397      REAL(DP) :: bigger_lat_dv, center_p(3), lat_dv, lat_dv3(3), lat_min(3), lat_p(3), &
398         max_coord_lat(3), min_coord_lat(3), smaller_lat_dv
399      REAL(DP), POINTER                                  :: charges(:)
400      REAL(KIND=dp), DIMENSION(3)                        :: max_coord, min_coord, transl_v
401      TYPE(cell_type), POINTER                           :: mm_cell, qm_cell
402      TYPE(cp_subsys_type), POINTER                      :: subsys_mm, subsys_qm
403      TYPE(particle_type), DIMENSION(:), POINTER         :: particles_mm, particles_qm
404      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
405      TYPE(section_vals_type), POINTER                   :: subsys_section
406
407      NULLIFY (subsys_mm, subsys_qm, qm_atom_index, particles_mm, particles_qm, &
408               subsys_section, qm_cell, mm_cell, qs_kind_set)
409
410      CPASSERT(ASSOCIATED(qmmm_env))
411      CPASSERT(qmmm_env%ref_count > 0)
412
413      CALL fist_env_get(qmmm_env%fist_env, cell=mm_cell, subsys=subsys_mm)
414      CALL get_qs_env(qmmm_env%qs_env, cell=qm_cell, cp_subsys=subsys_qm)
415      qm_atom_index => qmmm_env%qm%qm_atom_index
416      CPASSERT(ASSOCIATED(qm_atom_index))
417
418      particles_qm => subsys_qm%particles%els
419      particles_mm => subsys_mm%particles%els
420      IF (.NOT. qmmm_env%qm%center_qm_subsys0) qmmm_env%qm%do_translate = .FALSE.
421      IF (qmmm_env%qm%do_translate) THEN
422         IF (.NOT. qmmm_env%qm%center_qm_subsys_pbc_aware) THEN
423            ! naive coordinate based min-max
424            min_coord = HUGE(0.0_dp)
425            max_coord = -HUGE(0.0_dp)
426            DO ip = 1, SIZE(qm_atom_index)
427               min_coord = MIN(min_coord, particles_mm(qm_atom_index(ip))%r)
428               max_coord = MAX(max_coord, particles_mm(qm_atom_index(ip))%r)
429            END DO
430         ELSE
431            !! periodic based min max (uses complex number based mean)
432            center_p = qmmm_pbc_aware_mean(particles_mm, mm_cell, qm_atom_index)
433            ALLOCATE (avoid(SIZE(qm_atom_index)))
434            DO i_dim = 1, 3
435               IF (mm_cell%perd(i_dim) /= 1) THEN
436                  ! find absolute min and max positions (along i_dim direction) in lattice coordinates
437                  min_coord_lat(i_dim) = HUGE(0.0_dp)
438                  max_coord_lat(i_dim) = -HUGE(0.0_dp)
439                  DO ip = 1, SIZE(qm_atom_index)
440                     lat_p = MATMUL(mm_cell%h_inv, particles_mm(qm_atom_index(ip))%r)
441                     min_coord_lat(i_dim) = MIN(lat_p(i_dim), min_coord_lat(i_dim))
442                     max_coord_lat(i_dim) = MAX(lat_p(i_dim), max_coord_lat(i_dim))
443                  END DO
444               ELSE
445                  ! find min_ip closest to (pbc-aware) mean pos
446                  avoid = .FALSE.
447                  min_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, avoid, center_p, i_dim, 0)
448                  avoid(min_ip) = .TRUE.
449                  ! find max_ip closest to min_ip
450                  max_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, avoid, &
451                                             particles_mm(qm_atom_index(min_ip))%r, i_dim, 0, lat_dv)
452                  avoid(max_ip) = .TRUE.
453                  ! switch min and max if necessary
454                  IF (lat_dv < 0.0) THEN
455                     tmp_ip = min_ip
456                     min_ip = max_ip
457                     max_ip = tmp_ip
458                  ENDIF
459                  ! loop over all other atoms
460                  DO WHILE (.NOT. ALL(avoid))
461                     ! find smaller below min, bigger after max
462                     smaller_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, &
463                                                    avoid, particles_mm(qm_atom_index(min_ip))%r, i_dim, -1, smaller_lat_dv)
464                     bigger_ip = qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, &
465                                                   avoid, particles_mm(qm_atom_index(max_ip))%r, i_dim, 1, bigger_lat_dv)
466                     ! move min or max, not both
467                     IF (ABS(smaller_lat_dv) < ABS(bigger_lat_dv)) THEN
468                        min_ip = smaller_ip
469                        avoid(min_ip) = .TRUE.
470                     ELSE
471                        max_ip = bigger_ip
472                        avoid(max_ip) = .TRUE.
473                     END IF
474                  END DO
475                  ! find min and max coordinates in lattice positions (i_dim ! only)
476                  lat_dv3 = qmmm_lat_dv(mm_cell, particles_mm(qm_atom_index(min_ip))%r, particles_mm(qm_atom_index(max_ip))%r)
477                  IF (lat_dv3(i_dim) < 0.0_dp) lat_dv3(i_dim) = lat_dv3(i_dim) + 1.0_dp
478                  lat_min = MATMUL(mm_cell%h_inv, particles_mm(qm_atom_index(min_ip))%r)
479                  min_coord_lat(i_dim) = lat_min(i_dim)
480                  max_coord_lat(i_dim) = lat_min(i_dim) + lat_dv3(i_dim)
481               END IF ! periodic
482            END DO ! i_dim
483            ! min and max coordinates from lattice positions to Cartesian
484            min_coord = MATMUL(mm_cell%hmat, min_coord_lat)
485            max_coord = MATMUL(mm_cell%hmat, max_coord_lat)
486            DEALLOCATE (avoid)
487         END IF ! pbc aware center
488         transl_v = (max_coord + min_coord)/2.0_dp
489
490         !
491         ! The first time we always translate all the system in order
492         ! to centre the QM system in the box.
493         !
494         transl_v(:) = transl_v(:) - SUM(qm_cell%hmat, 2)/2.0_dp
495
496         IF (ANY(qmmm_env%qm%utrasl /= 1.0_dp)) THEN
497            transl_v = REAL(FLOOR(transl_v/qmmm_env%qm%utrasl), KIND=dp)* &
498                       qmmm_env%qm%utrasl
499         END IF
500         qmmm_env%qm%transl_v = qmmm_env%qm%transl_v + transl_v
501         particles_mm => subsys_mm%particles%els
502         DO ip = 1, subsys_mm%particles%n_els
503            particles_mm(ip)%r = particles_mm(ip)%r - transl_v
504         END DO
505         IF (qmmm_env%qm%added_shells%num_mm_atoms .GT. 0) THEN
506            DO ip = 1, qmmm_env%qm%added_shells%num_mm_atoms
507               qmmm_env%qm%added_shells%added_particles(ip)%r = qmmm_env%qm%added_shells%added_particles(ip)%r - transl_v
508               qmmm_env%qm%added_shells%added_cores(ip)%r = qmmm_env%qm%added_shells%added_cores(ip)%r - transl_v
509            END DO
510         END IF
511         unit_nr = cp_logger_get_default_io_unit()
512         IF (unit_nr > 0) WRITE (unit=unit_nr, fmt='(/1X,A)') &
513            " Translating the system in order to center the QM fragment in the QM box."
514         IF (.NOT. qmmm_env%qm%center_qm_subsys) qmmm_env%qm%do_translate = .FALSE.
515      END IF
516      particles_mm => subsys_mm%particles%els
517      DO ip = 1, SIZE(qm_atom_index)
518         particles_qm(ip)%r = particles_mm(qm_atom_index(ip))%r
519      END DO
520
521      subsys_section => section_vals_get_subs_vals(qmmm_env%qs_env%input, "SUBSYS")
522
523      CALL get_qs_env(qs_env=qmmm_env%qs_env, qs_kind_set=qs_kind_set)
524      CALL write_qs_particle_coordinates(particles_qm, qs_kind_set, subsys_section, "QM/MM first QM, then MM (0 charges)")
525      ALLOCATE (charges(SIZE(particles_mm)))
526      charges = 0.0_dp
527      CALL write_fist_particle_coordinates(particles_mm, subsys_section, charges)
528      DEALLOCATE (charges)
529
530   END SUBROUTINE apply_qmmm_translate
531
532! **************************************************************************************************
533!> \brief pbc-aware mean QM atom position
534!> \param particles_mm ...
535!> \param mm_cell ...
536!> \param qm_atom_index ...
537!> \return ...
538! **************************************************************************************************
539   FUNCTION qmmm_pbc_aware_mean(particles_mm, mm_cell, qm_atom_index)
540      TYPE(particle_type), DIMENSION(:), POINTER         :: particles_mm
541      TYPE(cell_type), POINTER                           :: mm_cell
542      INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index
543      REAL(dp)                                           :: qmmm_pbc_aware_mean(3)
544
545      COMPLEX(dp)                                        :: mean_z(3)
546      INTEGER                                            :: ip
547
548      mean_z = 0.0_dp
549      DO ip = 1, SIZE(qm_atom_index)
550         mean_z = mean_z + EXP(CMPLX(0.0_dp, 1.0_dp, KIND=dp)*2.0*pi* &
551                               MATMUL(mm_cell%h_inv, particles_mm(qm_atom_index(ip))%r))
552      END DO
553      mean_z = mean_z/ABS(mean_z)
554      qmmm_pbc_aware_mean = MATMUL(mm_cell%hmat, &
555                                   REAL(LOG(mean_z)/(CMPLX(0.0_dp, 1.0_dp, KIND=dp)*2.0_dp*pi), dp))
556   END FUNCTION
557
558! **************************************************************************************************
559!> \brief minimum image lattice coordinates difference vector
560!> \param mm_cell ...
561!> \param p1 ...
562!> \param p2 ...
563!> \return ...
564! **************************************************************************************************
565   FUNCTION qmmm_lat_dv(mm_cell, p1, p2)
566      TYPE(cell_type), POINTER                           :: mm_cell
567      REAL(dp)                                           :: p1(3), p2(3), qmmm_lat_dv(3)
568
569      REAL(dp)                                           :: lat_v1(3), lat_v2(3)
570
571      lat_v1 = MATMUL(mm_cell%h_inv, p1)
572      lat_v2 = MATMUL(mm_cell%h_inv, p2)
573
574      qmmm_lat_dv = lat_v2 - lat_v1
575      qmmm_lat_dv = qmmm_lat_dv - FLOOR(qmmm_lat_dv)
576   END FUNCTION qmmm_lat_dv
577
578! **************************************************************************************************
579!> \brief find closest QM particle, in position/negative direction
580!>        if dir is 1 or -1, respectively
581!> \param particles_mm ...
582!> \param mm_cell ...
583!> \param qm_atom_index ...
584!> \param avoid ...
585!> \param p ...
586!> \param i_dim ...
587!> \param dir ...
588!> \param closest_dv ...
589!> \return ...
590! **************************************************************************************************
591   FUNCTION qmmm_find_closest(particles_mm, mm_cell, qm_atom_index, avoid, p, i_dim, dir, closest_dv) RESULT(closest_ip)
592      TYPE(particle_type), DIMENSION(:), POINTER         :: particles_mm
593      TYPE(cell_type), POINTER                           :: mm_cell
594      INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index
595      LOGICAL                                            :: avoid(:)
596      REAL(dp)                                           :: p(3)
597      INTEGER                                            :: i_dim, dir
598      REAL(dp), OPTIONAL                                 :: closest_dv
599      INTEGER                                            :: closest_ip
600
601      INTEGER                                            :: ip, shift
602      REAL(dp)                                           :: lat_dv3(3), lat_dv_shifted, my_closest_dv
603
604      closest_ip = -1
605      my_closest_dv = HUGE(0.0)
606      DO ip = 1, SIZE(qm_atom_index)
607         IF (avoid(ip)) CYCLE
608         lat_dv3 = qmmm_lat_dv(mm_cell, p, particles_mm(qm_atom_index(ip))%r)
609         DO shift = -1, 1
610            lat_dv_shifted = lat_dv3(i_dim) + shift*1.0_dp
611            IF (ABS(lat_dv_shifted) < ABS(my_closest_dv) .AND. (dir*lat_dv_shifted >= 0.0)) THEN
612               my_closest_dv = lat_dv_shifted
613               closest_ip = ip
614            ENDIF
615         END DO
616      END DO
617
618      IF (PRESENT(closest_dv)) THEN
619         closest_dv = my_closest_dv
620      ENDIF
621
622   END FUNCTION qmmm_find_closest
623
624! **************************************************************************************************
625!> \brief Computes a spherical cutoff factor for the QMMM interactions
626!> \param spherical_cutoff ...
627!> \param rij ...
628!> \param factor ...
629!> \par History
630!>      08.2008 created
631!> \author Teodoro Laino
632! **************************************************************************************************
633   SUBROUTINE spherical_cutoff_factor(spherical_cutoff, rij, factor)
634      REAL(KIND=dp), DIMENSION(2), INTENT(IN)            :: spherical_cutoff
635      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rij
636      REAL(KIND=dp), INTENT(OUT)                         :: factor
637
638      CHARACTER(len=*), PARAMETER :: routineN = 'spherical_cutoff_factor', &
639         routineP = moduleN//':'//routineN
640
641      REAL(KIND=dp)                                      :: r, r0
642
643      r = SQRT(DOT_PRODUCT(rij, rij))
644      r0 = spherical_cutoff(1) - 20.0_dp*spherical_cutoff(2)
645      factor = 0.5_dp*(1.0_dp - TANH((r - r0)/spherical_cutoff(2)))
646
647   END SUBROUTINE spherical_cutoff_factor
648
649END MODULE qmmm_util
650