1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Multipole structure: for multipole (fixed and induced) in FF based MD
8!> \author Teodoro Laino [tlaino] - University of Zurich - 12.2007
9! **************************************************************************************************
10MODULE multipole_types
11   USE atomic_kind_types,               ONLY: get_atomic_kind
12   USE external_potential_types,        ONLY: fist_potential_type,&
13                                              get_potential
14   USE input_section_types,             ONLY: section_vals_get,&
15                                              section_vals_get_subs_vals,&
16                                              section_vals_type,&
17                                              section_vals_val_get
18   USE kinds,                           ONLY: dp
19   USE particle_types,                  ONLY: particle_type
20#include "../base/base_uses.f90"
21
22   IMPLICIT NONE
23
24   PRIVATE
25   PUBLIC :: multipole_type, &
26             create_multipole_type, &
27             release_multipole_type, &
28             retain_multipole_type
29
30   INTEGER, PARAMETER, PUBLIC               :: do_multipole_none = -1, &
31                                               do_multipole_charge = 0, &
32                                               do_multipole_dipole = 1, &
33                                               do_multipole_quadrupole = 2
34
35! **************************************************************************************************
36!> \brief Define multipole type
37!> \param error variable to control error logging, stopping,...
38!>        see module cp_error_handling
39!> \par History
40!>      12.2007 created [tlaino] - Teodoro Laino - University of Zurich
41!> \author Teodoro Laino
42! **************************************************************************************************
43   TYPE multipole_type
44      INTEGER :: id_nr, ref_count
45      LOGICAL, DIMENSION(3)                    :: task
46      REAL(KIND=dp), DIMENSION(:), POINTER     :: charges
47      REAL(KIND=dp), DIMENSION(:), POINTER     :: radii
48      REAL(KIND=dp), DIMENSION(:, :), POINTER   :: dipoles
49      REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: quadrupoles
50   END TYPE multipole_type
51
52   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'multipole_types'
53   INTEGER, PRIVATE, SAVE :: last_multipole_id_nr = 0
54
55CONTAINS
56
57! **************************************************************************************************
58!> \brief Create a multipole type
59!> \param multipoles ...
60!> \param particle_set ...
61!> \param subsys_section ...
62!> \param max_multipole ...
63!> \par History
64!>      12.2007 created [tlaino] - Teodoro Laino - University of Zurich
65!> \author Teodoro Laino
66! **************************************************************************************************
67   SUBROUTINE create_multipole_type(multipoles, particle_set, subsys_section, max_multipole)
68      TYPE(multipole_type), POINTER                      :: multipoles
69      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
70      TYPE(section_vals_type), POINTER                   :: subsys_section
71      INTEGER, INTENT(IN)                                :: max_multipole
72
73      CHARACTER(len=*), PARAMETER :: routineN = 'create_multipole_type', &
74         routineP = moduleN//':'//routineN
75
76      INTEGER                                            :: i, ind2, iparticle, j, n_rep, nparticles
77      LOGICAL                                            :: explicit
78      REAL(KIND=dp), DIMENSION(:), POINTER               :: work
79      TYPE(fist_potential_type), POINTER                 :: fist_potential
80      TYPE(section_vals_type), POINTER                   :: work_section
81
82      ALLOCATE (multipoles)
83
84      last_multipole_id_nr = last_multipole_id_nr + 1
85      multipoles%id_nr = last_multipole_id_nr
86      multipoles%ref_count = 1
87      multipoles%task = .FALSE.
88      NULLIFY (multipoles%charges)
89      NULLIFY (multipoles%radii)
90      NULLIFY (multipoles%dipoles)
91      NULLIFY (multipoles%quadrupoles)
92      SELECT CASE (max_multipole)
93      CASE (do_multipole_none)
94         ! Do nothing..
95      CASE (do_multipole_charge)
96         multipoles%task(1:1) = .TRUE.
97      CASE (do_multipole_dipole)
98         multipoles%task(1:2) = .TRUE.
99      CASE (do_multipole_quadrupole)
100         multipoles%task(1:3) = .TRUE.
101      CASE DEFAULT
102         CPABORT("")
103      END SELECT
104      nparticles = SIZE(particle_set)
105      IF (multipoles%task(1)) THEN
106         ALLOCATE (multipoles%charges(nparticles))
107         ALLOCATE (multipoles%radii(nparticles))
108         ! Fill in charge array
109         DO iparticle = 1, nparticles
110            !atomic_kind =>
111            CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, &
112                                 fist_potential=fist_potential)
113            CALL get_potential(fist_potential, qeff=multipoles%charges(iparticle), &
114                               mm_radius=multipoles%radii(iparticle))
115         END DO
116      END IF
117      IF (multipoles%task(2)) THEN
118         ALLOCATE (multipoles%dipoles(3, nparticles))
119         ! Fill in dipole array (if specified)
120         work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES%DIPOLES")
121         CALL section_vals_get(work_section, explicit=explicit)
122         IF (explicit) THEN
123            CALL section_vals_val_get(work_section, "_DEFAULT_KEYWORD_", n_rep_val=n_rep)
124            CPASSERT(n_rep == nparticles)
125            DO iparticle = 1, n_rep
126               CALL section_vals_val_get(work_section, "_DEFAULT_KEYWORD_", i_rep_val=iparticle, r_vals=work)
127               multipoles%dipoles(1:3, iparticle) = work
128            END DO
129         ELSE
130            multipoles%dipoles = 0.0_dp
131         END IF
132      END IF
133      IF (multipoles%task(3)) THEN
134         ALLOCATE (multipoles%quadrupoles(3, 3, nparticles))
135         ! Fill in quadrupole array (if specified)
136         work_section => section_vals_get_subs_vals(subsys_section, "MULTIPOLES%QUADRUPOLES")
137         CALL section_vals_get(work_section, explicit=explicit)
138         IF (explicit) THEN
139            CALL section_vals_val_get(work_section, "_DEFAULT_KEYWORD_", n_rep_val=n_rep)
140            CPASSERT(n_rep == nparticles)
141            DO iparticle = 1, n_rep
142               CALL section_vals_val_get(work_section, "_DEFAULT_KEYWORD_", i_rep_val=iparticle, r_vals=work)
143               DO i = 1, 3
144                  DO j = 1, 3
145                     ind2 = 3*(MIN(i, j) - 1) - (MIN(i, j)*(MIN(i, j) - 1))/2 + MAX(i, j)
146                     multipoles%quadrupoles(i, j, iparticle) = work(ind2)
147                  END DO
148               END DO
149            END DO
150         ELSE
151            multipoles%quadrupoles = 0.0_dp
152         END IF
153      END IF
154   END SUBROUTINE create_multipole_type
155
156! **************************************************************************************************
157!> \brief ...
158!> \param multipoles ...
159!> \par History
160!>      12.2007 created [tlaino] - Teodoro Laino - University of Zurich
161!> \author Teodoro Laino
162! **************************************************************************************************
163   SUBROUTINE release_multipole_type(multipoles)
164      TYPE(multipole_type), POINTER                      :: multipoles
165
166      CHARACTER(len=*), PARAMETER :: routineN = 'release_multipole_type', &
167         routineP = moduleN//':'//routineN
168
169      IF (ASSOCIATED(multipoles)) THEN
170         CPASSERT(multipoles%ref_count > 0)
171         multipoles%ref_count = multipoles%ref_count - 1
172         IF (multipoles%ref_count == 0) THEN
173            IF (ASSOCIATED(multipoles%charges)) THEN
174               DEALLOCATE (multipoles%charges)
175            END IF
176            IF (ASSOCIATED(multipoles%radii)) THEN
177               DEALLOCATE (multipoles%radii)
178            END IF
179            IF (ASSOCIATED(multipoles%dipoles)) THEN
180               DEALLOCATE (multipoles%dipoles)
181            END IF
182            IF (ASSOCIATED(multipoles%quadrupoles)) THEN
183               DEALLOCATE (multipoles%quadrupoles)
184            END IF
185            DEALLOCATE (multipoles)
186         END IF
187      END IF
188   END SUBROUTINE release_multipole_type
189
190! **************************************************************************************************
191!> \brief ...
192!> \param multipoles ...
193!> \par History
194!>      12.2007 created [tlaino] - Teodoro Laino - University of Zurich
195!> \author Teodoro Laino
196! **************************************************************************************************
197   SUBROUTINE retain_multipole_type(multipoles)
198      TYPE(multipole_type), POINTER                      :: multipoles
199
200      CHARACTER(len=*), PARAMETER :: routineN = 'retain_multipole_type', &
201         routineP = moduleN//':'//routineN
202
203      IF (ASSOCIATED(multipoles)) THEN
204         CPASSERT(multipoles%ref_count > 0)
205         multipoles%ref_count = multipoles%ref_count + 1
206      END IF
207   END SUBROUTINE retain_multipole_type
208
209END MODULE multipole_types
210