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