1 |
chuckv |
4 |
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 |