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 : !> \par History
10 : !> EAM potential
11 : !> \author CJM, I-Feng W. Kuo, Teodoro Laino
12 : ! **************************************************************************************************
13 : MODULE manybody_eam
14 :
15 : USE bibliography, ONLY: Foiles1986,&
16 : cite_reference
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: eam_type,&
21 : fist_nonbond_env_get,&
22 : fist_nonbond_env_set,&
23 : fist_nonbond_env_type,&
24 : pos_type
25 : USE kinds, ONLY: dp
26 : USE message_passing, ONLY: mp_para_env_type
27 : USE pair_potential_types, ONLY: ea_type,&
28 : eam_pot_type,&
29 : pair_potential_pp_type
30 : USE particle_types, ONLY: particle_type
31 : #include "./base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 : PUBLIC :: get_force_eam, density_nonbond
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_eam'
38 :
39 : CONTAINS
40 :
41 : ! **************************************************************************************************
42 : !> \brief ...
43 : !> \param fist_nonbond_env ...
44 : !> \param particle_set ...
45 : !> \param cell ...
46 : !> \param para_env ...
47 : !> \author CJM
48 : ! **************************************************************************************************
49 83143 : SUBROUTINE density_nonbond(fist_nonbond_env, particle_set, cell, para_env)
50 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
51 : TYPE(particle_type), DIMENSION(:), INTENT(INOUT) :: particle_set
52 : TYPE(cell_type), POINTER :: cell
53 : TYPE(mp_para_env_type), POINTER :: para_env
54 :
55 : CHARACTER(LEN=*), PARAMETER :: routineN = 'density_nonbond'
56 :
57 : INTEGER :: atom_a, atom_b, handle, i, i_a, i_b, &
58 : iend, igrp, ikind, ilist, ipair, &
59 : iparticle, istart, jkind, kind_a, &
60 : kind_b, nkinds, nparticle
61 83143 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eam_kinds_index
62 : LOGICAL :: do_eam
63 : REAL(KIND=dp) :: fac, rab2, rab2_max
64 : REAL(KIND=dp), DIMENSION(3) :: cell_v, rab
65 83143 : REAL(KIND=dp), DIMENSION(:), POINTER :: rho
66 : TYPE(eam_pot_type), POINTER :: eam_a, eam_b
67 83143 : TYPE(eam_type), DIMENSION(:), POINTER :: eam_data
68 : TYPE(fist_neighbor_type), POINTER :: nonbonded
69 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
70 : TYPE(pair_potential_pp_type), POINTER :: potparm
71 83143 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update, r_last_update_pbc
72 :
73 83143 : CALL timeset(routineN, handle)
74 83143 : do_eam = .FALSE.
75 : CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, &
76 : potparm=potparm, r_last_update=r_last_update, &
77 83143 : r_last_update_pbc=r_last_update_pbc, eam_data=eam_data)
78 83143 : nkinds = SIZE(potparm%pot, 1)
79 332572 : ALLOCATE (eam_kinds_index(nkinds, nkinds))
80 3033069 : eam_kinds_index = -1
81 316642 : DO ikind = 1, nkinds
82 1791605 : DO jkind = ikind, nkinds
83 3183473 : DO i = 1, SIZE(potparm%pot(ikind, jkind)%pot%type)
84 2949974 : IF (potparm%pot(ikind, jkind)%pot%type(i) == ea_type) THEN
85 : ! At the moment we allow only 1 EAM per each kinds pair..
86 692 : CPASSERT(eam_kinds_index(ikind, jkind) == -1)
87 692 : CPASSERT(eam_kinds_index(jkind, ikind) == -1)
88 692 : eam_kinds_index(ikind, jkind) = i
89 692 : eam_kinds_index(jkind, ikind) = i
90 692 : do_eam = .TRUE.
91 : END IF
92 : END DO
93 : END DO
94 : END DO
95 :
96 83143 : nparticle = SIZE(particle_set)
97 :
98 83143 : IF (do_eam) THEN
99 284 : IF (.NOT. ASSOCIATED(eam_data)) THEN
100 246 : ALLOCATE (eam_data(nparticle))
101 12 : CALL fist_nonbond_env_set(fist_nonbond_env, eam_data=eam_data)
102 : END IF
103 7586 : DO i = 1, nparticle
104 7302 : eam_data(i)%rho = 0.0_dp
105 7586 : eam_data(i)%f_embed = 0.0_dp
106 : END DO
107 : END IF
108 :
109 : ! Only if EAM potential are present
110 83143 : IF (do_eam) THEN
111 : ! Add citation
112 284 : CALL cite_reference(Foiles1986)
113 284 : NULLIFY (eam_a, eam_b)
114 852 : ALLOCATE (rho(nparticle))
115 7586 : rho = 0._dp
116 : ! Starting the force loop
117 46056 : DO ilist = 1, nonbonded%nlists
118 45772 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
119 45772 : IF (neighbor_kind_pair%npairs == 0) CYCLE
120 21866 : Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
121 13581 : istart = neighbor_kind_pair%grp_kind_start(igrp)
122 13581 : iend = neighbor_kind_pair%grp_kind_end(igrp)
123 13581 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
124 13581 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
125 :
126 13581 : i = eam_kinds_index(ikind, jkind)
127 13581 : IF (i == -1) CYCLE
128 13535 : rab2_max = potparm%pot(ikind, jkind)%pot%rcutsq
129 216560 : cell_v = MATMUL(cell%hmat, REAL(neighbor_kind_pair%cell_vector, KIND=dp))
130 217208 : DO ipair = istart, iend
131 157901 : atom_a = neighbor_kind_pair%list(1, ipair)
132 157901 : atom_b = neighbor_kind_pair%list(2, ipair)
133 157901 : fac = 1.0_dp
134 157901 : IF (atom_a == atom_b) fac = 0.5_dp
135 631604 : rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
136 631604 : rab = rab + cell_v
137 157901 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
138 171436 : IF (rab2 <= rab2_max) THEN
139 97493 : kind_a = particle_set(atom_a)%atomic_kind%kind_number
140 97493 : kind_b = particle_set(atom_b)%atomic_kind%kind_number
141 97493 : i_a = eam_kinds_index(kind_a, kind_a)
142 97493 : i_b = eam_kinds_index(kind_b, kind_b)
143 97493 : eam_a => potparm%pot(kind_a, kind_a)%pot%set(i_a)%eam
144 97493 : eam_b => potparm%pot(kind_b, kind_b)%pot%set(i_b)%eam
145 97493 : CALL get_rho_eam(eam_a, eam_b, rab2, atom_a, atom_b, rho, fac)
146 : END IF
147 : END DO
148 : END DO Kind_Group_Loop
149 : END DO
150 14888 : CALL para_env%sum(rho)
151 7586 : DO iparticle = 1, nparticle
152 7586 : eam_data(iparticle)%rho = rho(iparticle)
153 : END DO
154 :
155 284 : DEALLOCATE (rho)
156 : END IF
157 83143 : DEALLOCATE (eam_kinds_index)
158 83143 : CALL timestop(handle)
159 :
160 83143 : END SUBROUTINE density_nonbond
161 :
162 : ! **************************************************************************************************
163 : !> \brief ...
164 : !> \param eam_a ...
165 : !> \param eam_b ...
166 : !> \param rab2 ...
167 : !> \param atom_a ...
168 : !> \param atom_b ...
169 : !> \param rho ...
170 : !> \param fac ...
171 : !> \author CJM
172 : ! **************************************************************************************************
173 97493 : SUBROUTINE get_rho_eam(eam_a, eam_b, rab2, atom_a, atom_b, rho, fac)
174 : TYPE(eam_pot_type), POINTER :: eam_a, eam_b
175 : REAL(dp), INTENT(IN) :: rab2
176 : INTEGER, INTENT(IN) :: atom_a, atom_b
177 : REAL(dp), INTENT(INOUT) :: rho(:)
178 : REAL(dp), INTENT(IN) :: fac
179 :
180 : INTEGER :: index
181 : REAL(dp) :: qq, rab, rhoi, rhoj
182 :
183 97493 : rab = SQRT(rab2)
184 :
185 : ! Particle A
186 97493 : index = INT(rab/eam_b%drar) + 1
187 97493 : IF (index > eam_b%npoints) THEN
188 : index = eam_b%npoints
189 : ELSEIF (index < 1) THEN
190 : index = 1
191 : END IF
192 97493 : qq = rab - eam_b%rval(index)
193 97493 : rhoi = eam_b%rho(index) + qq*eam_b%rhop(index)
194 :
195 : ! Particle B
196 97493 : index = INT(rab/eam_a%drar) + 1
197 97493 : IF (index > eam_a%npoints) THEN
198 : index = eam_a%npoints
199 : ELSEIF (index < 1) THEN
200 : index = 1
201 : END IF
202 97493 : qq = rab - eam_a%rval(index)
203 97493 : rhoj = eam_a%rho(index) + qq*eam_a%rhop(index)
204 :
205 97493 : rho(atom_a) = rho(atom_a) + rhoi*fac
206 97493 : rho(atom_b) = rho(atom_b) + rhoj*fac
207 97493 : END SUBROUTINE get_rho_eam
208 :
209 : ! **************************************************************************************************
210 : !> \brief ...
211 : !> \param rab2 ...
212 : !> \param eam_a ...
213 : !> \param eam_b ...
214 : !> \param eam_data ...
215 : !> \param atom_a ...
216 : !> \param atom_b ...
217 : !> \param f_eam ...
218 : !> \author CJM
219 : ! **************************************************************************************************
220 97493 : SUBROUTINE get_force_eam(rab2, eam_a, eam_b, eam_data, atom_a, atom_b, f_eam)
221 : REAL(dp), INTENT(IN) :: rab2
222 : TYPE(eam_pot_type), POINTER :: eam_a, eam_b
223 : TYPE(eam_type), INTENT(IN) :: eam_data(:)
224 : INTEGER, INTENT(IN) :: atom_a, atom_b
225 : REAL(dp), INTENT(OUT) :: f_eam
226 :
227 : INTEGER :: index
228 : REAL(KIND=dp) :: denspi, denspj, fcp, qq, rab
229 :
230 97493 : rab = SQRT(rab2)
231 :
232 : ! Particle A
233 97493 : index = INT(rab/eam_a%drar) + 1
234 97493 : IF (index > eam_a%npoints) THEN
235 : index = eam_a%npoints
236 : ELSEIF (index < 1) THEN
237 : index = 1
238 : END IF
239 97493 : qq = rab - eam_a%rval(index)
240 97493 : IF (index == eam_a%npoints) THEN
241 1 : denspi = eam_a%rhop(index) + qq*(eam_a%rhop(index) - eam_a%rhop(index - 1))/eam_a%drar
242 : ELSE
243 97492 : denspi = eam_a%rhop(index) + qq*(eam_a%rhop(index + 1) - eam_a%rhop(index))/eam_a%drar
244 : END IF
245 :
246 : ! Particle B
247 97493 : index = INT(rab/eam_b%drar) + 1
248 97493 : IF (index > eam_b%npoints) THEN
249 : index = eam_b%npoints
250 : ELSEIF (index < 1) THEN
251 : index = 1
252 : END IF
253 97493 : qq = rab - eam_b%rval(index)
254 97493 : IF (index == eam_b%npoints) THEN
255 1 : denspj = eam_b%rhop(index) + qq*(eam_b%rhop(index) - eam_b%rhop(index - 1))/eam_b%drar
256 : ELSE
257 97492 : denspj = eam_b%rhop(index) + qq*(eam_b%rhop(index + 1) - eam_b%rhop(index))/eam_b%drar
258 : END IF
259 :
260 97493 : fcp = denspj*eam_data(atom_a)%f_embed + denspi*eam_data(atom_b)%f_embed
261 97493 : f_eam = fcp/rab
262 97493 : END SUBROUTINE get_force_eam
263 :
264 : END MODULE manybody_eam
265 :
|
