1 |
module force_utilities |
2 |
use definitions, ONLY : DP |
3 |
use parameter |
4 |
use simulation, ONLY : q,q0,natoms,use_pbc,box,nlocal |
5 |
use status, ONLY: error,info |
6 |
#ifdef MPI |
7 |
use mpi_module |
8 |
#endif |
9 |
implicit none |
10 |
PRIVATE |
11 |
|
12 |
public :: check |
13 |
public :: save_nlist |
14 |
public :: wrap |
15 |
|
16 |
contains |
17 |
|
18 |
subroutine check(update_nlist) |
19 |
|
20 |
integer :: i |
21 |
real( kind = DP ) :: dispmx |
22 |
logical, intent(out) :: update_nlist |
23 |
real( kind = DP ) :: dispmx_tmp |
24 |
|
25 |
dispmx = 0.0E0_DP |
26 |
|
27 |
!! calculate the largest displacement of any atom in any direction |
28 |
|
29 |
#ifdef MPI |
30 |
dispmx_tmp = 0.0E0_DP |
31 |
do i = 1, nlocal |
32 |
dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx ) |
33 |
dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx ) |
34 |
dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx ) |
35 |
end do |
36 |
call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, & |
37 |
mpi_max,mpi_comm_world,mpi_err) |
38 |
#else |
39 |
do i = 1, natoms |
40 |
dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx ) |
41 |
dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx ) |
42 |
dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx ) |
43 |
end do |
44 |
#endif |
45 |
|
46 |
!! a conservative test of list skin crossings |
47 |
dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx) |
48 |
|
49 |
update_nlist = (dispmx.gt.(skin_thickness)) |
50 |
|
51 |
end subroutine check |
52 |
|
53 |
subroutine save_nlist() |
54 |
integer :: i |
55 |
#ifdef MPI |
56 |
do i = 1, nlocal |
57 |
#else |
58 |
do i = 1, natoms |
59 |
#endif |
60 |
q0(1,i) = q(1,i) |
61 |
q0(2,i) = q(2,i) |
62 |
q0(3,i) = q(3,i) |
63 |
end do |
64 |
|
65 |
end subroutine save_nlist |
66 |
|
67 |
function wrap(r,dim) result(this_wrap) |
68 |
|
69 |
|
70 |
real( kind = DP ) :: r |
71 |
real( kind = DP ) :: this_wrap |
72 |
integer :: dim |
73 |
|
74 |
if (use_pbc) then |
75 |
! this_wrap = r - box(dim)*dsign(1.0E0_DP,r)*int(abs(r/box(dim)) + 0.5E0_DP) |
76 |
this_wrap = r - box(dim)*nint(r/box(dim)) |
77 |
else |
78 |
this_wrap = r |
79 |
endif |
80 |
|
81 |
return |
82 |
end function wrap |
83 |
|
84 |
|
85 |
end module force_utilities |