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 ace_nlist
15 :
16 : USE ace_wrapper, ONLY: ace_model_compute
17 : USE cell_types, ONLY: cell_type
18 : USE fist_neighbor_list_types, ONLY: fist_neighbor_type,&
19 : neighbor_kind_pairs_type
20 : USE fist_nonbond_env_types, ONLY: ace_data_type,&
21 : fist_nonbond_env_get,&
22 : fist_nonbond_env_type,&
23 : pos_type
24 : USE kinds, ONLY: dp
25 : USE physcon, ONLY: angstrom
26 : #include "./base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 :
30 : PRIVATE
31 : PUBLIC ace_interface
32 :
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ace_nlist'
34 :
35 : CONTAINS
36 :
37 : !
38 : !-------------------------------------------------------------------------------------
39 :
40 : ! **************************************************************************************************
41 : !> \brief ...
42 : !> \param ace_natom ...
43 : !> \param ace_atype ...
44 : !> \param pot_ace ...
45 : !> \param ace_force ...
46 : !> \param ace_virial ...
47 : !> \param fist_nonbond_env ...
48 : !> \param cell ...
49 : !> \param ace_data ...
50 : ! **************************************************************************************************
51 206 : SUBROUTINE ace_interface(ace_natom, ace_atype, pot_ace, ace_force, ace_virial, &
52 : fist_nonbond_env, cell, ace_data)
53 :
54 : INTEGER, INTENT(IN) :: ace_natom, ace_atype(1:ace_natom)
55 : REAL(kind=dp), INTENT(OUT) :: pot_ace, ace_force(1:3, 1:ace_natom), &
56 : ace_virial(1:6)
57 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
58 : TYPE(cell_type), POINTER :: cell
59 : TYPE(ace_data_type), POINTER :: ace_data
60 :
61 : #if defined(__ACE)
62 : INTEGER :: atom_a, atom_b, counter, ilist, k, m, n, &
63 : natom, nghost, num_update
64 206 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ghostidx, listidx
65 206 : REAL(KIND=8), ALLOCATABLE :: forceunroll(:)
66 412 : REAL(kind=dp) :: cell_v(3), dv(1:3), energy(1:ace_natom)
67 : TYPE(fist_neighbor_type), POINTER :: nonbonded
68 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
69 206 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
70 :
71 206 : natom = ace_natom
72 :
73 : CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, &
74 : r_last_update_pbc=r_last_update_pbc, &
75 206 : num_update=num_update, counter=counter)
76 :
77 206 : IF ((counter == 1) .OR. (ace_data%refupdate /= num_update)) THEN
78 : ! fist neighborlist are new:
79 14 : ace_data%refupdate = num_update
80 :
81 14 : IF (.NOT. ALLOCATED(ace_data%neiat)) THEN
82 18 : ALLOCATE (ace_data%neiat(0:natom))
83 : ELSE
84 8 : CPASSERT(SIZE(ace_data%neiat) > natom)
85 : END IF
86 :
87 : !if more than one mpi task, some neighorlists might be empty
88 : !do we need to check for lone atoms?
89 56 : ALLOCATE (ghostidx(natom), listidx(natom))
90 14 : nghost = 0
91 2554 : ace_data%neiat(:) = 0
92 14 : ace_data%nei = 0
93 1796 : DO n = 1, SIZE(nonbonded%neighbor_kind_pairs)
94 1782 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(n)
95 1796 : IF (neighbor_kind_pair%npairs > 0) THEN
96 : IF ((neighbor_kind_pair%cell_vector(1) == 0) .AND. &
97 439 : (neighbor_kind_pair%cell_vector(2) == 0) .AND. &
98 : (neighbor_kind_pair%cell_vector(3) == 0)) THEN
99 2540 : DO ilist = 1, natom
100 2540 : ghostidx(ilist) = ilist
101 : END DO
102 : ELSE
103 66284 : ghostidx = 0
104 : END IF
105 282151 : DO ilist = 1, neighbor_kind_pair%npairs
106 281712 : atom_a = ace_data%inverse_index_map(neighbor_kind_pair%list(1, ilist))
107 281712 : atom_b = ace_data%inverse_index_map(neighbor_kind_pair%list(2, ilist))
108 281712 : IF ((atom_a == 0) .OR. (atom_b == 0)) CYCLE
109 281712 : ace_data%neiat(atom_a) = ace_data%neiat(atom_a) + 1
110 282151 : IF (ghostidx(atom_b) == 0) THEN
111 : ! new ghost
112 15930 : nghost = nghost + 1
113 15930 : ghostidx(atom_b) = nghost + natom
114 : END IF
115 : END DO
116 : END IF
117 : END DO
118 : ! build running sum
119 2540 : DO n = 1, natom
120 2540 : ace_data%neiat(n) = ace_data%neiat(n) + ace_data%neiat(n - 1)
121 : END DO
122 14 : ace_data%nei = ace_data%neiat(natom)
123 14 : ace_data%nghost = nghost
124 :
125 14 : IF (ALLOCATED(ace_data%nlist)) THEN
126 8 : IF (SIZE(ace_data%nlist) < ace_data%nei) THEN
127 2 : DEALLOCATE (ace_data%nlist)
128 6 : ALLOCATE (ace_data%nlist(1:ace_data%nei))
129 : END IF
130 : ELSE
131 18 : ALLOCATE (ace_data%nlist(1:ace_data%nei))
132 : END IF
133 :
134 14 : IF (ALLOCATED(ace_data%attype)) THEN
135 8 : IF (SIZE(ace_data%attype) < natom + nghost) THEN
136 3 : DEALLOCATE (ace_data%atpos)
137 3 : DEALLOCATE (ace_data%attype)
138 3 : DEALLOCATE (ace_data%origin)
139 3 : DEALLOCATE (ace_data%shift)
140 9 : ALLOCATE (ace_data%atpos(1:3, 1:natom + nghost))
141 9 : ALLOCATE (ace_data%attype(1:natom + nghost))
142 6 : ALLOCATE (ace_data%origin(1:natom + nghost))
143 9 : ALLOCATE (ace_data%shift(1:3, 1:natom + nghost))
144 : END IF
145 : ELSE
146 18 : ALLOCATE (ace_data%atpos(1:3, 1:natom + nghost))
147 18 : ALLOCATE (ace_data%attype(1:natom + nghost))
148 12 : ALLOCATE (ace_data%origin(1:natom + nghost))
149 18 : ALLOCATE (ace_data%shift(1:3, 1:natom + nghost))
150 : END IF
151 2540 : ace_data%attype(1:natom) = ace_atype(:)
152 :
153 2540 : DO n = 1, natom
154 20208 : ace_data%atpos(:, n) = r_last_update_pbc(ace_data%use_indices(n))%r*angstrom
155 2540 : ace_data%origin(n) = n
156 : END DO
157 74222 : ace_data%shift(:, :) = 0
158 :
159 14 : k = natom
160 2540 : listidx(1:natom) = ace_data%neiat(0:natom - 1)
161 1796 : DO n = 1, SIZE(nonbonded%neighbor_kind_pairs)
162 1782 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(n)
163 1796 : IF (neighbor_kind_pair%npairs > 0) THEN
164 : IF ((neighbor_kind_pair%cell_vector(1) == 0) .AND. &
165 439 : (neighbor_kind_pair%cell_vector(2) == 0) .AND. &
166 : (neighbor_kind_pair%cell_vector(3) == 0)) THEN
167 2540 : DO m = 1, natom
168 2540 : ghostidx(m) = m
169 : END DO
170 : ELSE
171 66284 : ghostidx = 0
172 : END IF
173 1756 : dv = REAL(neighbor_kind_pair%cell_vector, KIND=dp)
174 : ! dimensions it's not periodic with should be zero in cell_vector
175 : ! so should always be valid:
176 7463 : cell_v = MATMUL(cell%hmat, dv)*angstrom
177 282151 : DO ilist = 1, neighbor_kind_pair%npairs
178 281712 : atom_a = ace_data%inverse_index_map(neighbor_kind_pair%list(1, ilist))
179 281712 : atom_b = ace_data%inverse_index_map(neighbor_kind_pair%list(2, ilist))
180 281712 : IF ((atom_a == 0) .OR. (atom_b == 0)) CYCLE
181 281712 : IF (ghostidx(atom_b) == 0) THEN
182 15930 : k = k + 1
183 63720 : ace_data%atpos(:, k) = ace_data%atpos(:, atom_b) + cell_v
184 127440 : ace_data%shift(:, k) = neighbor_kind_pair%cell_vector
185 15930 : ace_data%origin(k) = atom_b
186 15930 : ace_data%attype(k) = ace_atype(atom_b)
187 15930 : ghostidx(atom_b) = k
188 : END IF
189 281712 : listidx(atom_a) = listidx(atom_a) + 1
190 282151 : ace_data%nlist(listidx(atom_a)) = ghostidx(atom_b)
191 : END DO
192 : END IF
193 : END DO
194 :
195 14 : DEALLOCATE (ghostidx)
196 14 : DEALLOCATE (listidx)
197 :
198 : ! transforming to c call
199 : ! -> atom index will change from 1...n to 0...n-1 -> subtract 1 from neighborlist
200 281726 : ace_data%nlist(1:ace_data%nei) = ace_data%nlist(1:ace_data%nei) - 1
201 18470 : ace_data%origin(1:natom + nghost) = ace_data%origin(1:natom + nghost) - 1
202 : !-----------------------------------------------
203 :
204 : ELSE ! no changes in neighbor list just update positions:
205 192 : nghost = ace_data%nghost
206 37056 : DO n = 1, natom
207 295104 : ace_data%atpos(:, n) = r_last_update_pbc(ace_data%use_indices(n))%r*angstrom
208 : END DO
209 236023 : DO n = natom + 1, natom + nghost
210 943324 : dv = REAL(ace_data%shift(:, n), KIND=dp)*angstrom
211 : !origin+1 since we already changed to C notation for origin:
212 3773488 : ace_data%atpos(:, n) = ace_data%atpos(:, ace_data%origin(n) + 1) + MATMUL(cell%hmat, dv)
213 : END DO
214 : END IF
215 :
216 : ! -> force array
217 618 : ALLOCATE (forceunroll(1:3*natom))
218 118376 : forceunroll = 0.0
219 206 : pot_ace = 0.0
220 206 : ace_virial = 0.0
221 :
222 : CALL ace_model_compute( &
223 : natomc=natom, &
224 : nghostc=nghost, &
225 : neic=ace_data%nei, &
226 : neiatc=ace_data%neiat, &
227 : originc=ace_data%origin, &
228 : nlistc=ace_data%nlist, &
229 : attypec=ace_data%attype, &
230 : atposc=RESHAPE(ace_data%atpos, (/3*(natom + nghost)/)), &
231 : forcec=forceunroll, &
232 : virialc=ace_virial, &
233 : energyc=energy, &
234 412 : model=ace_data%model)
235 :
236 618 : ace_force = RESHAPE(forceunroll, (/3, natom/))
237 39596 : pot_ace = SUM(energy(1:natom))
238 :
239 206 : DEALLOCATE (forceunroll)
240 :
241 : #else
242 : MARK_USED(ace_natom)
243 : MARK_USED(ace_atype)
244 : MARK_USED(pot_ace)
245 : MARK_USED(ace_force)
246 : MARK_USED(ace_virial)
247 : MARK_USED(fist_nonbond_env)
248 : MARK_USED(cell)
249 : MARK_USED(ace_data)
250 : CPABORT("CP2K was compiled without ACE library.")
251 : #endif
252 :
253 206 : END SUBROUTINE ace_interface
254 :
255 : !----------------------------------------------------------------------------------
256 :
257 : END MODULE ace_nlist
|
