ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/simulation_module.F90
Revision: 619
Committed: Tue Jul 15 22:22:41 2003 UTC (20 years, 11 months ago) by mmeineke
File size: 9122 byte(s)
Log Message:
fixed a long lived bug in do_forces. Rrf was not being used in the neighborlist correctly. rcut was conssistently being set lowere than Rrf causing the dipole cutoff region to be to small. Also led to the removal of the taper region to buffer the dipole cutoff.

File Contents

# Content
1 !! Fortran interface to C entry plug.
2
3 module simulation
4 use definitions
5 use neighborLists
6 use force_globals
7 use vector_class
8 use atype_module
9 #ifdef IS_MPI
10 use mpiSimulation
11 #endif
12
13 implicit none
14 PRIVATE
15
16 #define __FORTRAN90
17 #include "fSimulation.h"
18
19 type (simtype), public :: thisSim
20
21 logical, save :: simulation_setup_complete = .false.
22
23 integer, public, save :: nLocal, nGlobal
24 integer, public, save :: nExcludes_Global = 0
25 integer, public, save :: nExcludes_Local = 0
26 integer, allocatable, dimension(:,:), public :: excludesLocal
27 integer, allocatable, dimension(:), public :: excludesGlobal
28 integer, allocatable, dimension(:), public :: molMembershipList
29
30 real(kind=dp), save :: rcut2 = 0.0_DP
31 real(kind=dp), save :: rcut6 = 0.0_DP
32 real(kind=dp), save :: rlist2 = 0.0_DP
33 real(kind=dp), public, dimension(3,3), save :: Hmat, HmatInv
34 logical, public, save :: boxIsOrthorhombic
35
36 public :: SimulationSetup
37 public :: getNlocal
38 public :: setBox
39 public :: setRcut
40 public :: getRcut
41 public :: getRlist
42 public :: getRrf
43 public :: getRt
44 public :: getDielect
45 public :: SimUsesPBC
46 public :: SimUsesLJ
47 public :: SimUsesDipoles
48 public :: SimUsesSticky
49 public :: SimUsesRF
50 public :: SimUsesGB
51 public :: SimUsesEAM
52 public :: SimRequiresPrepairCalc
53 public :: SimRequiresPostpairCalc
54 public :: SimUsesDirectionalAtoms
55
56 contains
57
58 subroutine SimulationSetup(setThisSim, CnGlobal, CnLocal, c_idents, &
59 CnLocalExcludes, CexcludesLocal, CnGlobalExcludes, CexcludesGlobal, &
60 CmolMembership, &
61 status)
62
63 type (simtype) :: setThisSim
64 integer, intent(inout) :: CnGlobal, CnLocal
65 integer, dimension(CnLocal),intent(inout) :: c_idents
66
67 integer :: CnLocalExcludes
68 integer, dimension(2,CnLocalExcludes), intent(in) :: CexcludesLocal
69 integer :: CnGlobalExcludes
70 integer, dimension(CnGlobalExcludes), intent(in) :: CexcludesGlobal
71 integer, dimension(CnGlobal),intent(in) :: CmolMembership
72 !! Result status, success = 0, status = -1
73 integer, intent(out) :: status
74 integer :: i, me, thisStat, alloc_stat, myNode
75 #ifdef IS_MPI
76 integer, allocatable, dimension(:) :: c_idents_Row
77 integer, allocatable, dimension(:) :: c_idents_Col
78 integer :: nrow
79 integer :: ncol
80 #endif
81
82 simulation_setup_complete = .false.
83 status = 0
84
85 ! copy C struct into fortran type
86
87 nLocal = CnLocal
88 nGlobal = CnGlobal
89
90 thisSim = setThisSim
91 rcut2 = thisSim%rcut * thisSim%rcut
92 rcut6 = rcut2 * rcut2 * rcut2
93 rlist2 = thisSim%rlist * thisSim%rlist
94
95 nExcludes_Global = CnGlobalExcludes
96 nExcludes_Local = CnLocalExcludes
97
98 call InitializeForceGlobals(nLocal, thisStat)
99 if (thisStat /= 0) then
100 write(default_error,*) "SimSetup: InitializeForceGlobals error"
101 status = -1
102 return
103 endif
104
105 call InitializeSimGlobals(thisStat)
106 if (thisStat /= 0) then
107 write(default_error,*) "SimSetup: InitializeSimGlobals error"
108 status = -1
109 return
110 endif
111
112 #ifdef IS_MPI
113 ! We can only set up forces if mpiSimulation has been setup.
114 if (.not. isMPISimSet()) then
115 write(default_error,*) "MPI is not set"
116 status = -1
117 return
118 endif
119 nrow = getNrow(plan_row)
120 ncol = getNcol(plan_col)
121 mynode = getMyNode()
122
123 allocate(c_idents_Row(nrow),stat=alloc_stat)
124 if (alloc_stat /= 0 ) then
125 status = -1
126 return
127 endif
128
129 allocate(c_idents_Col(ncol),stat=alloc_stat)
130 if (alloc_stat /= 0 ) then
131 status = -1
132 return
133 endif
134
135 call gather(c_idents, c_idents_Row, plan_row)
136 call gather(c_idents, c_idents_Col, plan_col)
137
138 do i = 1, nrow
139 me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
140 atid_Row(i) = me
141 enddo
142
143 do i = 1, ncol
144 me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
145 atid_Col(i) = me
146 enddo
147
148 !! free temporary ident arrays
149 if (allocated(c_idents_Col)) then
150 deallocate(c_idents_Col)
151 end if
152 if (allocated(c_idents_Row)) then
153 deallocate(c_idents_Row)
154 endif
155
156 #else
157 do i = 1, nLocal
158
159 me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
160 atid(i) = me
161
162 enddo
163 #endif
164
165
166
167 do i = 1, nExcludes_Local
168 excludesLocal(1,i) = CexcludesLocal(1,i)
169 excludesLocal(2,i) = CexcludesLocal(2,i)
170 enddo
171
172 do i = 1, nExcludes_Global
173 excludesGlobal(i) = CexcludesGlobal(i)
174 enddo
175
176 do i = 1, nGlobal
177 molMemberShipList(i) = CmolMembership(i)
178 enddo
179
180 if (status == 0) simulation_setup_complete = .true.
181
182 end subroutine SimulationSetup
183
184 subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
185 real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
186 integer :: cBoxIsOrthorhombic
187 integer :: smallest, status, i
188
189 Hmat = cHmat
190 HmatInv = cHmatInv
191 if (cBoxIsOrthorhombic .eq. 0 ) then
192 boxIsOrthorhombic = .false.
193 else
194 boxIsOrthorhombic = .true.
195 endif
196
197 return
198 end subroutine setBox
199
200 subroutine setRcut(new_rcut, status)
201 real(kind = dp) :: new_rcut
202 integer :: myStatus, status
203 thisSim%rcut = new_rcut
204 rcut2 = thisSim%rcut * thisSim%rcut
205 rcut6 = rcut2 * rcut2 * rcut2
206 status = 0
207 return
208 end subroutine setRcut
209
210 subroutine getRcut(thisrcut,rc2,rc6,status)
211 real( kind = dp ), intent(out) :: thisrcut
212 real( kind = dp ), intent(out), optional :: rc2
213 real( kind = dp ), intent(out), optional :: rc6
214 integer, optional :: status
215
216 if (present(status)) status = 0
217
218 if (.not.simulation_setup_complete ) then
219 if (present(status)) status = -1
220 return
221 end if
222
223 thisrcut = thisSim%rcut
224 if(present(rc2)) rc2 = rcut2
225 if(present(rc6)) rc6 = rcut6
226 end subroutine getRcut
227
228 subroutine getRlist(thisrlist,rl2,status)
229 real( kind = dp ), intent(out) :: thisrlist
230 real( kind = dp ), intent(out), optional :: rl2
231
232 integer, optional :: status
233
234 if (present(status)) status = 0
235
236 if (.not.simulation_setup_complete ) then
237 if (present(status)) status = -1
238 return
239 end if
240
241 thisrlist = thisSim%rlist
242 if(present(rl2)) rl2 = rlist2
243 end subroutine getRlist
244
245 function getRrf() result(rrf)
246 real( kind = dp ) :: rrf
247 rrf = thisSim%rrf
248 end function getRrf
249
250 function getRt() result(rt)
251 real( kind = dp ) :: rt
252 rt = thisSim%rt
253 end function getRt
254
255 function getDielect() result(dielect)
256 real( kind = dp ) :: dielect
257 dielect = thisSim%dielect
258 end function getDielect
259
260 function SimUsesPBC() result(doesit)
261 logical :: doesit
262 doesit = thisSim%SIM_uses_PBC
263 end function SimUsesPBC
264
265 function SimUsesLJ() result(doesit)
266 logical :: doesit
267 doesit = thisSim%SIM_uses_LJ
268 end function SimUsesLJ
269
270 function SimUsesSticky() result(doesit)
271 logical :: doesit
272 doesit = thisSim%SIM_uses_sticky
273 end function SimUsesSticky
274
275 function SimUsesDipoles() result(doesit)
276 logical :: doesit
277 doesit = thisSim%SIM_uses_dipoles
278 end function SimUsesDipoles
279
280 function SimUsesRF() result(doesit)
281 logical :: doesit
282 doesit = thisSim%SIM_uses_RF
283 end function SimUsesRF
284
285 function SimUsesGB() result(doesit)
286 logical :: doesit
287 doesit = thisSim%SIM_uses_GB
288 end function SimUsesGB
289
290 function SimUsesEAM() result(doesit)
291 logical :: doesit
292 doesit = thisSim%SIM_uses_EAM
293 end function SimUsesEAM
294
295 function SimUsesDirectionalAtoms() result(doesit)
296 logical :: doesit
297 doesit = thisSim%SIM_uses_dipoles .or. thisSim%SIM_uses_sticky .or. &
298 thisSim%SIM_uses_GB .or. thisSim%SIM_uses_RF
299 end function SimUsesDirectionalAtoms
300
301 function SimRequiresPrepairCalc() result(doesit)
302 logical :: doesit
303 doesit = thisSim%SIM_uses_EAM
304 end function SimRequiresPrepairCalc
305
306 function SimRequiresPostpairCalc() result(doesit)
307 logical :: doesit
308 doesit = thisSim%SIM_uses_RF
309 end function SimRequiresPostpairCalc
310
311 subroutine InitializeSimGlobals(thisStat)
312 integer, intent(out) :: thisStat
313 integer :: alloc_stat
314
315 thisStat = 0
316
317 call FreeSimGlobals()
318
319 allocate(excludesLocal(2,nExcludes_Local), stat=alloc_stat)
320 if (alloc_stat /= 0 ) then
321 thisStat = -1
322 return
323 endif
324
325 allocate(excludesGlobal(nExcludes_Global), stat=alloc_stat)
326 if (alloc_stat /= 0 ) then
327 thisStat = -1
328 return
329 endif
330
331 allocate(molMembershipList(nGlobal), stat=alloc_stat)
332 if (alloc_stat /= 0 ) then
333 thisStat = -1
334 return
335 endif
336
337 end subroutine InitializeSimGlobals
338
339 subroutine FreeSimGlobals()
340
341 !We free in the opposite order in which we allocate in.
342
343 if (allocated(molMembershipList)) deallocate(molMembershipList)
344 if (allocated(excludesGlobal)) deallocate(excludesGlobal)
345 if (allocated(excludesLocal)) deallocate(excludesLocal)
346
347 end subroutine FreeSimGlobals
348
349 pure function getNlocal() result(n)
350 integer :: n
351 n = nLocal
352 end function getNlocal
353
354
355 end module simulation