Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of dispersion using pair potentials
10 : !> \author Johann Pototschnig
11 : ! **************************************************************************************************
12 : MODULE qs_dispersion_d4
13 :
14 : USE atomic_kind_types, ONLY: get_atomic_kind
15 : USE cell_types, ONLY: cell_type, &
16 : get_cell
17 : #if defined(__DFTD4)
18 : USE dftd4, ONLY: d4_model, &
19 : damping_param, &
20 : get_dispersion, &
21 : get_rational_damping, &
22 : new, &
23 : new_d4_model, &
24 : realspace_cutoff, &
25 : structure_type
26 : #endif
27 : USE kinds, ONLY: dp
28 : USE particle_types, ONLY: particle_type
29 : USE qs_dispersion_types, ONLY: qs_dispersion_type
30 : USE qs_force_types, ONLY: qs_force_type
31 : USE message_passing, ONLY: mp_para_env_type
32 : #include "./base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 :
36 : PRIVATE
37 :
38 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d4'
39 :
40 : PUBLIC :: calculate_dispersion_d4_pairpot
41 :
42 : ! **************************************************************************************************
43 :
44 : CONTAINS
45 :
46 : #if defined(__DFTD4)
47 : ! **************************************************************************************************
48 : !> \brief ...
49 : !> \param dispersion_env ...
50 : !> \param particle_set ...
51 : !> \param cell ...
52 : !> \param para_env ...
53 : !> \param energy ...
54 : !> \param force ...
55 : !> \param atom_of_kind ...
56 : !> \param virial ...
57 : ! **************************************************************************************************
58 6 : SUBROUTINE calculate_dispersion_d4_pairpot(dispersion_env, particle_set, cell, para_env, energy, force, atom_of_kind, virial)
59 :
60 : TYPE(qs_dispersion_type), INTENT(IN), POINTER :: dispersion_env
61 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
62 : POINTER :: particle_set
63 : TYPE(cell_type), INTENT(IN), POINTER :: cell
64 : TYPE(mp_para_env_type), POINTER :: para_env
65 : REAL(KIND=dp), INTENT(OUT) :: energy
66 : TYPE(qs_force_type), DIMENSION(:), INTENT(INOUT), &
67 : OPTIONAL, POINTER :: force
68 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: atom_of_kind
69 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
70 : OPTIONAL :: virial
71 :
72 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d4_pairpot'
73 :
74 12 : class(damping_param), ALLOCATABLE :: param
75 6 : TYPE(d4_model) :: disp
76 6 : TYPE(structure_type) :: mol
77 : TYPE(realspace_cutoff) :: cutoff
78 :
79 : INTEGER :: iatom, natom, ind_atom
80 : INTEGER, ALLOCATABLE, DIMENSION(:) :: el_num
81 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gradient, xyz
82 : REAL(KIND=dp), DIMENSION(3, 3) :: stress
83 : INTEGER :: handle, ikind
84 : !REAL(KIND=dp), DIMENSION(3, 3) :: h
85 : INTEGER, DIMENSION(3) :: periodic
86 : LOGICAL, DIMENSION(3) :: lperiod
87 :
88 6 : CALL timeset(routineN, handle)
89 :
90 : !get information about particles
91 6 : natom = SIZE(particle_set)
92 18 : ALLOCATE (xyz(3, natom))
93 18 : ALLOCATE (el_num(natom))
94 54 : DO iatom = 1, natom
95 192 : xyz(:, iatom) = particle_set(iatom)%r(:)
96 48 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
97 54 : el_num(iatom) = ikind
98 : END DO
99 :
100 : !get information about cell / lattice
101 6 : CALL get_cell(cell=cell, periodic=periodic)
102 6 : lperiod(1) = periodic(1) == 1
103 6 : lperiod(2) = periodic(2) == 1
104 6 : lperiod(3) = periodic(3) == 1
105 :
106 : !prepare for the call to the dispersion function
107 6 : CALL new(mol, el_num, xyz, lattice=cell%hmat, periodic=lperiod)
108 6 : CALL new_d4_model(disp, mol)
109 6 : CALL get_rational_damping(dispersion_env%ref_functional, param, s9=dispersion_env%s9)
110 :
111 : ! Coordination number cutoff
112 6 : cutoff%cn = dispersion_env%rc_cn
113 : ! Two-body interaction cutoff
114 6 : cutoff%disp2 = dispersion_env%rc_d4*2._dp
115 : ! Three-body interaction cutoff
116 6 : cutoff%disp3 = dispersion_env%rc_disp*2._dp
117 :
118 : !> Wrapper to handle the evaluation of dispersion energy and derivatives
119 6 : IF (PRESENT(force)) THEN
120 4 : IF (PRESENT(atom_of_kind)) THEN
121 12 : ALLOCATE (gradient(3, mol%nat))
122 4 : CALL get_dispersion(mol, disp, param, cutoff, energy, gradient, stress)
123 4 : IF (PRESENT(virial)) THEN
124 26 : virial = -1.00_dp*stress
125 2 : IF (para_env%num_pe > 1 .AND. para_env%mepos > 0) virial = 0.00_dp
126 : END IF
127 32 : DO iatom = 1, natom
128 28 : ikind = el_num(iatom)
129 28 : ind_atom = atom_of_kind(iatom)
130 116 : force(ikind)%dispersion(:, ind_atom) = force(ikind)%dispersion(:, ind_atom) + gradient(:, iatom)
131 : END DO
132 4 : DEALLOCATE (gradient)
133 : ELSE
134 0 : CPABORT("missing atom_of_kind")
135 : END IF
136 : ELSE
137 2 : IF (PRESENT(virial)) THEN
138 0 : CPABORT("missing force for virial term")
139 : END IF
140 2 : CALL get_dispersion(mol, disp, param, cutoff, energy)
141 : END IF
142 :
143 : !dispersion energy is computed by every MPI process
144 6 : IF (para_env%num_pe > 1 .AND. para_env%mepos > 0) energy = 0.00_dp
145 :
146 6 : DEALLOCATE (el_num, xyz)
147 :
148 6 : CALL timestop(handle)
149 :
150 18 : END SUBROUTINE calculate_dispersion_d4_pairpot
151 :
152 : #else
153 :
154 : ! **************************************************************************************************
155 : !> \brief ...
156 : !> \param dispersion_env ...
157 : !> \param particle_set ...
158 : !> \param cell ...
159 : !> \param para_env ...
160 : !> \param energy ...
161 : !> \param force ...
162 : !> \param atom_of_kind ...
163 : !> \param virial ...
164 : ! **************************************************************************************************
165 : SUBROUTINE calculate_dispersion_d4_pairpot(dispersion_env, particle_set, cell, para_env, energy, force, atom_of_kind, virial)
166 :
167 : TYPE(qs_dispersion_type), INTENT(IN), POINTER :: dispersion_env
168 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
169 : POINTER :: particle_set
170 : TYPE(cell_type), INTENT(IN), POINTER :: cell
171 : TYPE(mp_para_env_type), POINTER :: para_env
172 : REAL(KIND=dp), INTENT(OUT) :: energy
173 : TYPE(qs_force_type), DIMENSION(:), INTENT(INOUT), &
174 : OPTIONAL, POINTER :: force
175 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: atom_of_kind
176 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
177 : OPTIONAL :: virial
178 :
179 : MARK_USED(dispersion_env)
180 : MARK_USED(particle_set)
181 : MARK_USED(cell)
182 : MARK_USED(para_env)
183 : MARK_USED(energy)
184 : MARK_USED(force)
185 : MARK_USED(atom_of_kind)
186 : MARK_USED(virial)
187 :
188 : CPABORT("Build without DFTD4")
189 :
190 : END SUBROUTINE calculate_dispersion_d4_pairpot
191 :
192 : #endif
193 :
194 : END MODULE qs_dispersion_d4
|