ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/md_code/simulation_module.F90
(Generate patch)

Comparing trunk/mdtools/md_code/simulation_module.F90 (file contents):
Revision 215 by chuckv, Thu Dec 19 21:59:51 2002 UTC vs.
Revision 241 by mmeineke, Thu Jan 23 16:19:59 2003 UTC

# Line 1 | Line 1
1   module simulation
2    use definitions, ONLY :dp
3 <  use force_wrappers, ONLY : alloc_force_wrappers
3 > #ifdef IS_MPI
4 >  use mpiSimulation
5 > #endif
6  
7    implicit none
8    PRIVATE
9  
10  
9 –  
10 –  real ( kind = dp ), public, dimension(3) :: box
11  
12 +  type, public :: simtype
13 +     PRIVATE
14 +     SEQUENCE
15 + !! Number of particles on this processor
16 +     integer :: nLRparticles
17 + !! Periodic Box    
18 +     real ( kind = dp ), dimension(3) :: box
19 + !! List Cutoff    
20 +     real ( kind = dp ) :: rlist = 0.0_dp
21 + !! Radial cutoff
22 +     real ( kind = dp ) :: rcut  = 0.0_dp
23 + !! List cutoff squared
24 +     real ( kind = dp ) :: rlistsq = 0.0_dp
25 + !! Radial Cutoff squared
26 +     real ( kind = dp ) :: rcutsq  = 0.0_dp
27 + !! Radial Cutoff^6
28 +     real ( kind = dp ) :: rcut6  = 0.0_dp
29  
30 <  real ( kind = dp ), public :: rlist
14 <  real ( kind = dp ), public :: rcut
15 <  real ( kind = dp ), public :: rlistsq
16 <  real ( kind = dp ), public :: rcutsq
30 >  end type simtype
31  
32 <  integer,public, allocatable, dimension(:)   :: point
33 <  integer,public, allocatable, dimension(:)   :: list
32 >  type (simtype), public :: thisSim
33 > !! Tag for MPI calculations  
34 >  integer, allocatable, dimension(:) :: tag
35  
36 + #ifdef IS_MPI
37 +  integer, allocatable, dimension(:) :: tag_row
38 +  integer, allocatable, dimension(:) :: tag_column
39 + #endif
40  
22 – !MPI dependent routines
41  
42 < #ifdef IS_MPI
43 < ! Universal routines: All types of force calculations will need these arrays
44 < ! Arrays specific to a type of force calculation should be declared in that module.
27 <  real( kind = dp ), allocatable, dimension(:,:) :: qRow
28 <  real( kind = dp ), allocatable, dimension(:,:) :: qColumn
42 >  logical :: setSim = .false.
43 >  
44 >  real( kind = dp ), allocatable(:,:),save :: q0
45  
46 <  real( kind = dp ), allocatable, dimension(:,:) :: fRow
47 <  real( kind = dp ), allocatable, dimension(:,:) :: fColumn
46 >  public :: check
47 >  public :: save_nlist
48 >  public :: wrap
49 >  public :: getBox
50 >  public :: getRcut
51 >  public :: setRcut
52  
53 <  real( kind = dp ), allocatable, dimension(:,:) :: tRow
54 <  real( kind = dp ), allocatable, dimension(:,:) :: tColumn
53 >  interface wrap
54 >     module procedure wrap_1d
55 >     module procedure wrap_3d
56 >  end interface
57  
58 < #endif
58 >  interface getBox
59 >     module procedure getBox_3d
60 >     module procedure getBox_dim
61 >  end interface
62  
63  
64  
65  
66  
67  
68 +
69 +
70 +
71 +
72   contains
73  
74 < #ifdef MPI
75 < ! Allocated work arrays for MPI
76 <  subroutine allocate_mpi_arrays(nDimensions,numComponents)
77 <    integer, intent(in) :: nDimensions
78 <    integer, intent(in) :: numComponents
74 >  subroutine setSimulation(nLRParticles,box,rlist,rcut)
75 >    integer, intent(in) :: nLRParticles
76 >    real(kind = dp ), intent(in), dimension(3) :: box
77 >    real(kind = dp ), intent(in) :: rlist
78 >    real(kind = dp ), intent(in) :: rcut
79  
80 +    if( setsim ) return  ! simulation is already initialized
81 +    setSim = .true.
82  
83 +    thisSim%nLRParticles = nLRParticles
84 +    thisSim%box          = box
85 +    thisSim%rlist        = rlist
86 +    thisSim%rcut         = rcut
87 +    thisSim%rcutsq       = rcut * rcut
88 +    thisSim%rcut6        = rcutsq * rcutsq * rcutsq
89  
90 +  end subroutine setSimulation
91  
92 +  function getNparticles() result(nparticles)
93 +    integer :: nparticles
94 +    nparticles = thisSim%nLRparticles
95 +  end function getNparticles
96  
97 <  end subroutine allocate_mpi_arrays
97 >
98 >  subroutine change_box_size(new_box_size)
99 >    real(kind=dp), dimension(3) :: new_box_size
100 >
101 >    thisSim%box = new_box_size
102 >
103 >  end subroutine change_box_size
104 >
105 >
106 >  elemental function getBox_3d() result(thisBox)
107 >    real( kind = dp ), dimension(3) :: thisBox
108 >    thisBox = thisSim%box
109 >  end function getBox_3d
110 >
111 >  function getBox_dim(dim) result(thisBox)
112 >    integer, intent(in) :: dim
113 >    real( kind = dp ) :: thisBox
114 >
115 >    thisBox = thisSim%box(dim)
116 >  end function getBox_dim
117 >    
118 >  subroutine check(update_nlist)
119 >  
120 >    integer :: i
121 >    real( kind = DP ) :: dispmx
122 >    logical, intent(out) :: update_nlist
123 >    real( kind = DP ) :: dispmx_tmp
124 >    
125 >    dispmx = 0.0E0_DP
126 >    
127 >    !! calculate the largest displacement of any atom in any direction
128 >    
129 > #ifdef MPI
130 >    dispmx_tmp = 0.0E0_DP
131 >    do i = 1, nlocal
132 >       dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx )
133 >       dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx )
134 >       dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx )
135 >    end do
136 >    call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
137 >       mpi_max,mpi_comm_world,mpi_err)
138 > #else
139 >    do i = 1, natoms
140 >       dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
141 >       dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
142 >       dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
143 >    end do
144   #endif
145 +  
146 +    !! a conservative test of list skin crossings
147 +    dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
148 +    
149 +    update_nlist = (dispmx.gt.(skin_thickness))
150 +    
151 +  end subroutine check
152 +  
153 +  subroutine save_nlist(q)
154 +    real(kind = dp ), dimension(:,:), intent(in)  :: q
155 +    integer :: list_size
156 +    
157  
158 <  subroutine set_simulation(box,rlist,rcut)
158 >    list_size = size(q)
159 >
160 >    if (.not. allocated(q0)) then
161 >       allocate(q0(3,list_size))
162 >    else if( list_size > size(q0))
163 >       deallocate(q0)
164 >       allocate(q0(3,list_size))
165 >    endif
166 >
167 >    q0 = q
168 >
169 >  end subroutine save_nlist
170 >  
171 >
172 >  function wrap_1d(r,dim) result(this_wrap)
173      
174 +    
175 +    real( kind = DP ) :: r
176 +    real( kind = DP ) :: this_wrap
177 +    integer           :: dim
178 +    
179 +    if (use_pbc) then
180 +       !     this_wrap = r - box(dim)*dsign(1.0E0_DP,r)*int(abs(r/box(dim)) + 0.5E0_DP)
181 +       this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim))
182 +    else
183 +       this_wrap = r
184 +    endif
185 +    
186 +    return
187 +  end function wrap_1d
188  
189 +  elemental function wrap_3d(r) result(this_wrap)
190 +    real( kind = dp ), dimension(3), intent(in) :: r
191 +    real( kind = dp ), dimension(3) :: this_wrap
192  
193 <  end subroutine set_simulation
193 >    
194 >    if (use_pbc) then
195 >       !     this_wrap = r - box(dim)*dsign(1.0E0_DP,r)*int(abs(r/box(dim)) + 0.5E0_DP)
196 >       this_wrap(1:3) = r(1:3) - box(1:3)*nint(r(1:3)/box(1:3))
197 >    else
198 >       this_wrap = r
199 >    endif
200 >  end function wrap_3d
201  
202  
203  
204 <  subroutine change_box_size(new_box_size)
205 <    real(kind=dp), dimension(3) :: new_box_size
204 >  subroutine getRcut(thisrcut,rcut2,rcut6,status)
205 >    real( kind = dp ), intent(out) :: thisrcut
206 >    real( kind = dp ), intent(out), optional :: rcut2
207 >    real( kind = dp ), intent(out), optional :: thisrcut6
208 >    integer, optional :: status
209  
210 <    box = new_box_size
210 >    if (present(status)) status = 0
211  
212 <  end subroutine change_box_size
212 >    if (.not.setSim ) then
213 >       if (present(status)) status = -1
214 >       return
215 >    end if
216 >    
217 >    thisrcut = rcut
218 >    if(present(rcut2)) rcut2 = rcutsq
219 >    if(present(rcut2)) rcut6 = rcut_6
220  
221 +  end subroutine getRcut
222    
223 +
224 +
225 +
226   end module simulation

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines