1 |
module force_module |
2 |
use definitions, ONLY : DP,ndim |
3 |
use parameter |
4 |
use simulation |
5 |
use second_deriv |
6 |
use eam_module, ONLY: eam_check_type,eam_rcut,eam_atype_map |
7 |
use glue_module, ONLY: get_dpars |
8 |
use goddard_module, ONLY: rcutg |
9 |
use lj_module, ONLY: lj_sigma |
10 |
use status, ONLY: error,info,warning |
11 |
PRIVATE |
12 |
|
13 |
public :: initialize_forcefield,calc_forces |
14 |
public :: init_forcefield_parameters |
15 |
|
16 |
|
17 |
contains |
18 |
|
19 |
|
20 |
|
21 |
subroutine init_forcefield_parameters() |
22 |
|
23 |
use eam_module, ONLY: initialize_eam |
24 |
use glue_module, ONLY: initialize_glue |
25 |
use goddard_module, ONLY: initialize_goddard |
26 |
use lj_module, ONLY: initialize_lj |
27 |
|
28 |
! initialize sutton-chen potential parameters |
29 |
if (force_field.eq.'goddard') then |
30 |
call initialize_goddard() |
31 |
elseif(force_field.eq.'eam') then |
32 |
call initialize_eam() |
33 |
elseif(force_field.eq.'glue') then |
34 |
call initialize_glue() |
35 |
elseif(force_field.eq.'lj') then |
36 |
call initialize_lj() |
37 |
else |
38 |
call error('init_forcefield', 'UNKNOWN FORCE FIELD!') |
39 |
endif |
40 |
|
41 |
return |
42 |
end subroutine init_forcefield_parameters |
43 |
|
44 |
subroutine initialize_forcefield() |
45 |
use model_module |
46 |
#ifdef MPI |
47 |
use mpi_module |
48 |
#endif |
49 |
integer :: i, j |
50 |
|
51 |
integer :: atype1, atype2 |
52 |
real( kind = DP ) :: cutoffmax, maxcut, rc |
53 |
real( kind = DP ), dimension(13) :: dptmp |
54 |
logical :: gotit |
55 |
character(len=80) msg |
56 |
#ifdef MPI |
57 |
integer, allocatable, dimension(:) :: present_atypes_local |
58 |
integer, allocatable, dimension(:) :: global_present_atypes |
59 |
integer, allocatable, dimension(:) :: displs, counts |
60 |
integer :: max_present_local |
61 |
integer :: n_atypes_known |
62 |
integer :: n_present_local |
63 |
#endif |
64 |
n_atypes_present = 0 |
65 |
|
66 |
if (.not. allocated(present_atypes)) allocate(present_atypes(get_max_atype())) |
67 |
|
68 |
! This first section is to determine the number of atom types present |
69 |
! in the simulation. |
70 |
! |
71 |
!! we need to calculate the number of local atypes and then sync with |
72 |
!! the rest of the cluster, so we need to do this loop twice, one local |
73 |
!! and one global after a mpi_allreducev |
74 |
#ifdef MPI |
75 |
!! allocate temp arrays for allreducev |
76 |
allocate(present_atypes_local(get_max_atype())) |
77 |
allocate(displs(nprocs)) |
78 |
allocate(counts(nprocs)) |
79 |
n_present_local = 0 |
80 |
|
81 |
do i = 1, nlocal |
82 |
atype1 = ident(i) |
83 |
if (n_present_local.eq.0) then |
84 |
present_atypes_local(1) = atype1 |
85 |
n_present_local = 1 |
86 |
else |
87 |
gotit = .false. |
88 |
do j = 1, n_present_local |
89 |
if (present_atypes_local(j).eq.atype1) then |
90 |
gotit = .true. |
91 |
endif |
92 |
enddo |
93 |
|
94 |
if(.not.gotit) then |
95 |
present_atypes_local(n_present_local+1) = atype1 |
96 |
n_present_local = n_present_local + 1 |
97 |
endif |
98 |
endif |
99 |
enddo |
100 |
|
101 |
!! now gather the present_atypes to all nodes |
102 |
|
103 |
call mpi_allgather(n_present_local,1,mpi_integer,counts,1, & |
104 |
mpi_integer, mpi_comm_world,mpi_err) |
105 |
|
106 |
! call mpi_allreduce(n_present_local,max_present_local,1,mpi_integer, & |
107 |
! mpi_max,mpi_comm_world,mpi_err) |
108 |
!! figure out how large of an array to allocate for global_present_atypes |
109 |
max_present_local = 0 |
110 |
do i = 1, nprocs |
111 |
max_present_local = max_present_local + counts(i) |
112 |
end do |
113 |
|
114 |
allocate(global_present_atypes(max_present_local)) |
115 |
displs(1) = 0 |
116 |
do i = 2, nprocs |
117 |
displs(i) = displs(i-1) + counts(i-1) |
118 |
end do |
119 |
|
120 |
|
121 |
call mpi_allgatherv(present_atypes_local,n_present_local, & |
122 |
mpi_integer,global_present_atypes,counts,displs, & |
123 |
mpi_integer,mpi_comm_world,mpi_err) |
124 |
deallocate(present_atypes_local) |
125 |
deallocate(displs) |
126 |
deallocate(counts) |
127 |
|
128 |
|
129 |
do i = 1, max_present_local |
130 |
atype1 = global_present_atypes(i) |
131 |
if (n_atypes_present.eq.0) then |
132 |
present_atypes(1) = atype1 |
133 |
n_atypes_present = 1 |
134 |
else |
135 |
gotit = .false. |
136 |
do j = 1, n_atypes_present |
137 |
if (present_atypes(j).eq.atype1) then |
138 |
gotit = .true. |
139 |
endif |
140 |
enddo |
141 |
|
142 |
if(.not.gotit) then |
143 |
present_atypes(n_atypes_present+1) = atype1 |
144 |
n_atypes_present = n_atypes_present + 1 |
145 |
endif |
146 |
endif |
147 |
enddo |
148 |
deallocate(global_present_atypes) |
149 |
|
150 |
#else |
151 |
do i = 1, natoms |
152 |
atype1 = ident(i) |
153 |
if (n_atypes_present.eq.0) then |
154 |
present_atypes(1) = atype1 |
155 |
n_atypes_present = 1 |
156 |
else |
157 |
gotit = .false. |
158 |
do j = 1, n_atypes_present |
159 |
if (present_atypes(j).eq.atype1) then |
160 |
gotit = .true. |
161 |
endif |
162 |
enddo |
163 |
|
164 |
if(.not.gotit) then |
165 |
present_atypes(n_atypes_present+1) = atype1 |
166 |
n_atypes_present = n_atypes_present + 1 |
167 |
endif |
168 |
endif |
169 |
enddo |
170 |
#endif |
171 |
! loop over all possible pairs of atom types to determine the max |
172 |
! cutoff radius |
173 |
|
174 |
cutoffmax = 0.0E0_DP |
175 |
|
176 |
if (force_field.eq.'goddard') then |
177 |
|
178 |
do i = 1, n_atypes_present |
179 |
atype1 = present_atypes(i) |
180 |
do j = 1, n_atypes_present |
181 |
atype2 = present_atypes(j) |
182 |
rc = rcutg(atype1,atype2) |
183 |
|
184 |
if (rc.gt.cutoffmax) cutoffmax = rc |
185 |
enddo |
186 |
enddo |
187 |
|
188 |
elseif (force_field.eq.'glue') then |
189 |
|
190 |
do i = 1, n_atypes_present |
191 |
atype1 = present_atypes(i) |
192 |
call get_dpars(atype1,dptmp) |
193 |
|
194 |
rc = dptmp(3) |
195 |
|
196 |
if (rc.gt.cutoffmax) cutoffmax = rc |
197 |
enddo |
198 |
elseif (force_field.eq.'eam') then |
199 |
do i = 1, n_atypes_present |
200 |
atype1 = present_atypes(i) |
201 |
if(eam_check_type(atype1)) then |
202 |
rc = eam_rcut(eam_atype_map(atype1)) |
203 |
else |
204 |
write(msg,*) atype1, " Unknown eam atype" |
205 |
call error('INITIALIZE_FORCEFIELD',msg) |
206 |
end if |
207 |
if (rc.gt.cutoffmax) cutoffmax = rc |
208 |
enddo |
209 |
|
210 |
elseif (force_field.eq.'lj') then |
211 |
|
212 |
do i = 1, n_atypes_present |
213 |
|
214 |
atype1 = present_atypes(i) |
215 |
|
216 |
rc = 3.0E0_DP*lj_sigma(atype1) |
217 |
|
218 |
if (rc.gt.cutoffmax) cutoffmax = rc |
219 |
enddo |
220 |
|
221 |
endif |
222 |
|
223 |
|
224 |
if (cutoffmax.lt.rcut) then |
225 |
call info('INITIALIZE_FORCEFIELD',& |
226 |
'The force field defines a cutoff radius that is smaller') |
227 |
write(msg,'(a28,es12.3)') & |
228 |
'than rcut, so rcut is now: ', cutoffmax |
229 |
call info('INITIALIZE_FORCEFIELD', trim(msg)) |
230 |
rcut = cutoffmax |
231 |
endif |
232 |
|
233 |
if (use_pbc) then |
234 |
maxcut = rcut |
235 |
do i = 1, 3 |
236 |
if (maxcut.gt.box(i)/2.0E0_DP) then |
237 |
maxcut = box(i) / 2.0E0_DP |
238 |
endif |
239 |
enddo |
240 |
|
241 |
if (maxcut.ne.rcut) then |
242 |
call warning('INITIALIZE_FORCEFIELD',& |
243 |
'The size of the simulation box determined a maximum') |
244 |
write(msg,'(a18,es12.3)') & |
245 |
'cutoff radius of: ', maxcut |
246 |
call warning('INITIALIZE_FORCEFIELD',trim(msg)) |
247 |
rcut = maxcut |
248 |
endif |
249 |
endif |
250 |
|
251 |
rcutsq = rcut**2 |
252 |
|
253 |
rlstsq = (rcut+skin_thickness)**2 |
254 |
|
255 |
write(msg,'(a8,es12.3)') & |
256 |
'rlist = ', rcut + skin_thickness |
257 |
call info('INITIALIZE_FORCEFIELD',trim(msg)) |
258 |
write(msg,'(a8,es12.3)') & |
259 |
'rcut = ', rcut |
260 |
call info('INITIALIZE_FORCEFIELD',trim(msg)) |
261 |
return |
262 |
end subroutine initialize_forcefield |
263 |
|
264 |
|
265 |
subroutine calc_forces(update_nlist,nmflag,pe) |
266 |
use eam_module, ONLY: calc_eam_dens,calc_eam_forces |
267 |
use glue_module, ONLY: calc_glue_dens,calc_glue_forces |
268 |
use goddard_module, ONLY: calc_goddard_dens,calc_goddard_forces |
269 |
use lj_module, ONLY: calc_lj_forces |
270 |
|
271 |
logical, intent(inout) :: nmflag |
272 |
logical, intent(inout) :: update_nlist |
273 |
real( kind = DP ), intent(out), optional :: pe |
274 |
real( kind = DP ) :: dummy_pot |
275 |
|
276 |
call clean_forces(nmflag) |
277 |
if (present(pe)) then |
278 |
pe = 0.0 |
279 |
dummy_pot = 0.0E0_DP |
280 |
|
281 |
if (force_field.eq.'goddard') then |
282 |
call calc_goddard_dens(update_nlist) |
283 |
call calc_goddard_forces(nmflag,pot=dummy_pot) |
284 |
elseif(force_field.eq.'glue') then |
285 |
call calc_glue_dens(update_nlist) |
286 |
call calc_glue_forces(nmflag,pe) |
287 |
elseif(force_field.eq.'lj') then |
288 |
call calc_lj_forces(update_nlist, nmflag, pe=dummy_pot) |
289 |
elseif(force_field.eq.'eam') then |
290 |
call calc_eam_dens(update_nlist) |
291 |
call calc_eam_forces(nmflag,pot=dummy_pot) |
292 |
else |
293 |
call error('CALC_FORCES', 'UNKNOWN FORCE FIELD!') |
294 |
endif |
295 |
pe = dummy_pot |
296 |
|
297 |
else |
298 |
|
299 |
if (force_field.eq.'goddard') then |
300 |
call calc_goddard_dens(update_nlist) |
301 |
call calc_goddard_forces(nmflag) |
302 |
elseif(force_field.eq.'glue') then |
303 |
call calc_glue_dens(update_nlist) |
304 |
call calc_glue_forces(nmflag) |
305 |
elseif(force_field.eq.'lj') then |
306 |
call calc_lj_forces(update_nlist, nmflag) |
307 |
elseif(force_field.eq.'eam') then |
308 |
call calc_eam_dens(update_nlist) |
309 |
|
310 |
call calc_eam_forces(nmflag) |
311 |
else |
312 |
call error('CALC_FORCES', 'UNKNOWN FORCE FIELD!') |
313 |
endif |
314 |
endif |
315 |
end subroutine calc_forces |
316 |
|
317 |
subroutine clean_forces(nmflag) |
318 |
logical, intent(in) :: nmflag |
319 |
|
320 |
! clean out the regions of the force vector and Hessian that |
321 |
! we are going to use. You were using temporary force |
322 |
! and hessian arrays anyway, right? |
323 |
! |
324 |
|
325 |
! |
326 |
! 'if' is expensive. Put it outside the loop. |
327 |
|
328 |
if (nmflag) then |
329 |
f = 0.0E0_DP |
330 |
d = 0.0E0_DP |
331 |
else |
332 |
f = 0.0E0_DP |
333 |
end if |
334 |
|
335 |
end subroutine clean_forces |
336 |
end module force_module |