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

# User Rev Content
1 mmeineke 377 !! 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 mmeineke 491 integer, public, save :: nLocal, nGlobal
24 mmeineke 377 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 chuckv 482 integer, allocatable, dimension(:), public :: molMembershipList
29 mmeineke 377
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 mmeineke 569 real(kind=dp), public, dimension(3,3), save :: Hmat, HmatInv
34 gezelter 570 logical, public, save :: boxIsOrthorhombic
35 mmeineke 377
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 mmeineke 491 subroutine SimulationSetup(setThisSim, CnGlobal, CnLocal, c_idents, &
59 gezelter 483 CnLocalExcludes, CexcludesLocal, CnGlobalExcludes, CexcludesGlobal, &
60     CmolMembership, &
61 mmeineke 377 status)
62    
63     type (simtype) :: setThisSim
64 mmeineke 491 integer, intent(inout) :: CnGlobal, CnLocal
65     integer, dimension(CnLocal),intent(inout) :: c_idents
66 mmeineke 377
67     integer :: CnLocalExcludes
68     integer, dimension(2,CnLocalExcludes), intent(in) :: CexcludesLocal
69     integer :: CnGlobalExcludes
70     integer, dimension(CnGlobalExcludes), intent(in) :: CexcludesGlobal
71 mmeineke 491 integer, dimension(CnGlobal),intent(in) :: CmolMembership
72 mmeineke 377 !! 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 mmeineke 491
87     nLocal = CnLocal
88     nGlobal = CnGlobal
89    
90 mmeineke 377 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 mmeineke 491 call InitializeForceGlobals(nLocal, thisStat)
99 mmeineke 377 if (thisStat /= 0) then
100 chuckv 480 write(default_error,*) "SimSetup: InitializeForceGlobals error"
101 mmeineke 377 status = -1
102     return
103     endif
104    
105     call InitializeSimGlobals(thisStat)
106     if (thisStat /= 0) then
107 chuckv 480 write(default_error,*) "SimSetup: InitializeSimGlobals error"
108 mmeineke 377 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 gezelter 490 do i = 1, nLocal
158 mmeineke 377
159     me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
160     atid(i) = me
161    
162     enddo
163     #endif
164    
165    
166 chuckv 388
167 mmeineke 377 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 mmeineke 435
176 gezelter 490 do i = 1, nGlobal
177 gezelter 483 molMemberShipList(i) = CmolMembership(i)
178 mmeineke 491 enddo
179 chuckv 482
180 mmeineke 377 if (status == 0) simulation_setup_complete = .true.
181    
182     end subroutine SimulationSetup
183    
184 mmeineke 569 subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
185 gezelter 570 real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
186 mmeineke 569 integer :: cBoxIsOrthorhombic
187 mmeineke 377 integer :: smallest, status, i
188 gezelter 570
189 mmeineke 569 Hmat = cHmat
190     HmatInv = cHmatInv
191 gezelter 570 if (cBoxIsOrthorhombic .eq. 0 ) then
192 mmeineke 569 boxIsOrthorhombic = .false.
193 gezelter 570 else
194     boxIsOrthorhombic = .true.
195 mmeineke 569 endif
196    
197 mmeineke 377 return
198 mmeineke 569 end subroutine setBox
199 mmeineke 377
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 mmeineke 619 end function getRrf
249 mmeineke 377
250     function getRt() result(rt)
251     real( kind = dp ) :: rt
252     rt = thisSim%rt
253 mmeineke 619 end function getRt
254 mmeineke 377
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 chuckv 482
331 mmeineke 491 allocate(molMembershipList(nGlobal), stat=alloc_stat)
332 chuckv 482 if (alloc_stat /= 0 ) then
333     thisStat = -1
334     return
335     endif
336 mmeineke 377
337     end subroutine InitializeSimGlobals
338    
339     subroutine FreeSimGlobals()
340    
341     !We free in the opposite order in which we allocate in.
342 gezelter 483
343     if (allocated(molMembershipList)) deallocate(molMembershipList)
344 mmeineke 377 if (allocated(excludesGlobal)) deallocate(excludesGlobal)
345     if (allocated(excludesLocal)) deallocate(excludesLocal)
346 gezelter 483
347 mmeineke 377 end subroutine FreeSimGlobals
348    
349 mmeineke 491 pure function getNlocal() result(n)
350     integer :: n
351     n = nLocal
352 mmeineke 377 end function getNlocal
353    
354    
355     end module simulation