Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par History
10 : !> HAF (16-Apr-2025) : Import into CP2K
11 : !> \author HAF and yury-lysogorskiy and ralf-drautz
12 : ! **************************************************************************************************
13 :
14 : MODULE manybody_ace
15 :
16 : USE ISO_C_BINDING, ONLY: C_ASSOCIATED
17 : USE ace_nlist, ONLY: ace_interface
18 : USE atomic_kind_types, ONLY: atomic_kind_type,&
19 : get_atomic_kind
20 : USE bibliography, ONLY: Bochkarev2024,&
21 : Drautz2019,&
22 : Lysogorskiy2021,&
23 : cite_reference
24 : USE cell_types, ONLY: cell_type
25 : USE fist_nonbond_env_types, ONLY: ace_data_type,&
26 : fist_nonbond_env_get,&
27 : fist_nonbond_env_set,&
28 : fist_nonbond_env_type
29 : USE kinds, ONLY: dp
30 : USE pair_potential_types, ONLY: ace_type,&
31 : pair_potential_pp_type,&
32 : pair_potential_single_type
33 : USE particle_types, ONLY: particle_type
34 : USE physcon, ONLY: angstrom,&
35 : evolt
36 : #include "./base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 :
40 : PRIVATE
41 : PUBLIC ace_energy_store_force_virial, ace_add_force_virial
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_ace'
44 :
45 : CONTAINS
46 :
47 : ! **************************************************************************************************
48 : !> \brief ...
49 : !> \param particle_set ...
50 : !> \param atomic_kind_set ...
51 : !> \param potparm ...
52 : !> \param ace_data ...
53 : ! **************************************************************************************************
54 6 : SUBROUTINE init_ace_data(particle_set, atomic_kind_set, potparm, &
55 : ace_data)
56 :
57 : TYPE(particle_type), POINTER :: particle_set(:)
58 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
59 : TYPE(pair_potential_pp_type), POINTER :: potparm
60 : TYPE(ace_data_type), POINTER :: ace_data
61 :
62 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_ace_data'
63 :
64 : CHARACTER(2) :: element_symbol
65 : INTEGER :: ace_natom, handle, i, iat, iat_use, &
66 : ikind, jkind, lkind, n_atoms
67 6 : INTEGER, ALLOCATABLE :: use_atom_type(:)
68 6 : INTEGER, DIMENSION(:), POINTER :: ak_alist
69 6 : LOGICAL, ALLOCATABLE :: use_atom(:)
70 : TYPE(pair_potential_single_type), POINTER :: pot
71 :
72 6 : CALL timeset(routineN, handle)
73 :
74 : ! init ace_data
75 6 : IF (.NOT. ASSOCIATED(ace_data)) THEN
76 0 : ALLOCATE (ace_data)
77 : END IF
78 :
79 6 : n_atoms = SIZE(particle_set)
80 18 : ALLOCATE (use_atom(n_atoms))
81 12 : ALLOCATE (use_atom_type(n_atoms))
82 996 : use_atom = .FALSE.
83 996 : use_atom_type = 0
84 :
85 18 : DO ikind = 1, SIZE(atomic_kind_set)
86 12 : pot => potparm%pot(ikind, ikind)%pot
87 30 : DO i = 1, SIZE(pot%type)
88 12 : IF (pot%type(i) /= ace_type) CYCLE
89 : CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), &
90 : element_symbol=element_symbol, &
91 12 : natom=lkind, atom_list=ak_alist)
92 12 : IF (lkind < 1) CYCLE
93 12 : ace_data%model = pot%set(i)%ace%model
94 12 : jkind = 0
95 18 : DO iat = 1, SIZE(ace_data%model%symbolc)
96 18 : IF (element_symbol == ace_data%model%symbolc(iat)) THEN
97 : jkind = iat
98 : EXIT
99 : END IF
100 : END DO
101 12 : CPASSERT(jkind > 0)
102 1026 : DO iat = 1, lkind
103 990 : use_atom_type(ak_alist(iat)) = jkind
104 1002 : use_atom(ak_alist(iat)) = .TRUE.
105 : END DO
106 : END DO ! i
107 : END DO ! ikind
108 6 : CPASSERT(C_ASSOCIATED(ace_data%model%c_ptr))
109 :
110 996 : ace_natom = COUNT(use_atom)
111 :
112 6 : IF (.NOT. ALLOCATED(ace_data%uctype)) THEN
113 18 : ALLOCATE (ace_data%uctype(1:ace_natom))
114 : END IF
115 :
116 : iat_use = 0
117 996 : DO iat = 1, n_atoms
118 990 : IF (.NOT. use_atom(iat)) CYCLE
119 990 : iat_use = iat_use + 1
120 996 : ace_data%uctype(iat_use) = use_atom_type(iat)
121 : END DO
122 :
123 6 : IF (iat_use > 0) THEN
124 6 : CALL cite_reference(Drautz2019)
125 6 : CALL cite_reference(Lysogorskiy2021)
126 6 : CALL cite_reference(Bochkarev2024)
127 : END IF
128 :
129 6 : IF (.NOT. ALLOCATED(ace_data%force)) THEN
130 18 : ALLOCATE (ace_data%force(3, ace_natom))
131 18 : ALLOCATE (ace_data%use_indices(ace_natom))
132 12 : ALLOCATE (ace_data%inverse_index_map(n_atoms))
133 : END IF
134 6 : CPASSERT(SIZE(ace_data%force, 2) == ace_natom)
135 :
136 6 : iat_use = 0
137 996 : ace_data%inverse_index_map(:) = 0
138 996 : DO iat = 1, n_atoms
139 996 : IF (use_atom(iat)) THEN
140 990 : iat_use = iat_use + 1
141 990 : ace_data%use_indices(iat_use) = iat
142 990 : ace_data%inverse_index_map(iat) = iat_use
143 : END IF
144 : END DO
145 6 : ace_data%natom = ace_natom
146 6 : DEALLOCATE (use_atom, use_atom_type)
147 :
148 6 : CALL timestop(handle)
149 :
150 6 : END SUBROUTINE init_ace_data
151 :
152 : ! **************************************************************************************************
153 : !> \brief ...
154 : !> > \param particle_set ...
155 : !> \param particle_set ...
156 : !> \param cell ...
157 : !> \param atomic_kind_set ...
158 : !> \param potparm ...
159 : !> \param fist_nonbond_env ...
160 : !> \param pot_ace ...
161 : ! **************************************************************************************************
162 412 : SUBROUTINE ace_energy_store_force_virial(particle_set, cell, atomic_kind_set, potparm, &
163 : fist_nonbond_env, pot_ace)
164 :
165 : TYPE(particle_type), POINTER :: particle_set(:)
166 : TYPE(cell_type), POINTER :: cell
167 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
168 : TYPE(pair_potential_pp_type), POINTER :: potparm
169 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
170 : REAL(kind=dp), INTENT(OUT) :: pot_ace
171 :
172 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ace_energy_store_force_virial'
173 :
174 : INTEGER :: handle
175 : REAL(kind=dp) :: ace_virial(1:6)
176 : TYPE(ace_data_type), POINTER :: ace_data
177 :
178 206 : CALL timeset(routineN, handle)
179 :
180 : ! get ace_data to save force, virial info
181 206 : CALL fist_nonbond_env_get(fist_nonbond_env, ace_data=ace_data)
182 206 : IF (.NOT. ASSOCIATED(ace_data)) THEN
183 84 : ALLOCATE (ace_data)
184 : !initialize ace_data:
185 6 : CALL init_ace_data(particle_set, atomic_kind_set, potparm, ace_data)
186 6 : CALL fist_nonbond_env_set(fist_nonbond_env, ace_data=ace_data)
187 : END IF
188 :
189 : CALL ace_interface(ace_data%natom, ace_data%uctype, &
190 : pot_ace, ace_data%force, ace_virial, &
191 206 : fist_nonbond_env, cell, ace_data)
192 :
193 : ! convert units
194 206 : pot_ace = pot_ace/evolt
195 157766 : ace_data%force = ace_data%force/(evolt/angstrom)
196 1442 : ace_virial = ace_virial/evolt
197 :
198 : ! minus sign due to CP2K conventions
199 206 : ace_data%virial(1, 1) = -ace_virial(1)
200 206 : ace_data%virial(2, 2) = -ace_virial(2)
201 206 : ace_data%virial(3, 3) = -ace_virial(3)
202 206 : ace_data%virial(1, 2) = -ace_virial(4)
203 206 : ace_data%virial(2, 1) = -ace_virial(4)
204 206 : ace_data%virial(1, 3) = -ace_virial(5)
205 206 : ace_data%virial(3, 1) = -ace_virial(5)
206 206 : ace_data%virial(2, 3) = -ace_virial(6)
207 206 : ace_data%virial(3, 2) = -ace_virial(6)
208 :
209 206 : CALL timestop(handle)
210 206 : END SUBROUTINE ace_energy_store_force_virial
211 :
212 : ! **************************************************************************************************
213 : !> \brief ...
214 : !> \param fist_nonbond_env ...
215 : !> \param force ...
216 : !> \param pv_nonbond ...
217 : !> \param use_virial ...
218 : ! **************************************************************************************************
219 206 : SUBROUTINE ace_add_force_virial(fist_nonbond_env, force, pv_nonbond, use_virial)
220 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
221 : REAL(KIND=dp), INTENT(INOUT) :: force(:, :), pv_nonbond(3, 3)
222 : LOGICAL, OPTIONAL :: use_virial
223 :
224 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ace_add_force_virial'
225 :
226 : INTEGER :: handle, iat, iat_use
227 : TYPE(ace_data_type), POINTER :: ace_data
228 :
229 206 : CALL timeset(routineN, handle)
230 :
231 206 : CALL fist_nonbond_env_get(fist_nonbond_env, ace_data=ace_data)
232 :
233 206 : IF (.NOT. ASSOCIATED(ace_data)) RETURN
234 :
235 39596 : DO iat_use = 1, SIZE(ace_data%use_indices)
236 39390 : iat = ace_data%use_indices(iat_use)
237 39390 : CPASSERT(iat >= 1 .AND. iat <= SIZE(force, 2))
238 157766 : force(1:3, iat) = force(1:3, iat) + ace_data%force(1:3, iat_use)
239 : END DO
240 :
241 206 : IF (use_virial) THEN
242 26 : pv_nonbond = pv_nonbond + ace_data%virial
243 : END IF
244 :
245 206 : CALL timestop(handle)
246 : END SUBROUTINE ace_add_force_virial
247 :
248 : END MODULE manybody_ace
|
