ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/UseTheForce/DarkSide/simulation.F90
Revision: 1528
Committed: Fri Dec 17 20:11:05 2010 UTC (14 years, 7 months ago) by gezelter
File size: 23654 byte(s)
Log Message:
Doesn't build yet, but getting much closer...


File Contents

# User Rev Content
1 gezelter 246 !!
2     !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     !!
4     !! The University of Notre Dame grants you ("Licensee") a
5     !! non-exclusive, royalty free, license to use, modify and
6     !! redistribute this software in source and binary code form, provided
7     !! that the following conditions are met:
8     !!
9 gezelter 1390 !! 1. Redistributions of source code must retain the above copyright
10 gezelter 246 !! notice, this list of conditions and the following disclaimer.
11     !!
12 gezelter 1390 !! 2. Redistributions in binary form must reproduce the above copyright
13 gezelter 246 !! notice, this list of conditions and the following disclaimer in the
14     !! documentation and/or other materials provided with the
15     !! distribution.
16     !!
17     !! This software is provided "AS IS," without a warranty of any
18     !! kind. All express or implied conditions, representations and
19     !! warranties, including any implied warranty of merchantability,
20     !! fitness for a particular purpose or non-infringement, are hereby
21     !! excluded. The University of Notre Dame and its licensors shall not
22     !! be liable for any damages suffered by licensee as a result of
23     !! using, modifying or distributing the software or its
24     !! derivatives. In no event will the University of Notre Dame or its
25     !! licensors be liable for any lost revenue, profit or data, or for
26     !! direct, indirect, special, consequential, incidental or punitive
27     !! damages, however caused and regardless of the theory of liability,
28     !! arising out of the use of or inability to use software, even if the
29     !! University of Notre Dame has been advised of the possibility of
30     !! such damages.
31     !!
32 gezelter 1390 !! SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33     !! research, please cite the appropriate papers when you publish your
34     !! work. Good starting points are:
35     !!
36     !! [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37     !! [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38     !! [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39     !! [4] Vardeman & Gezelter, in progress (2009).
40     !!
41 gezelter 246
42 gezelter 115 !! Fortran interface to C entry plug.
43    
44     module simulation
45     use definitions
46 gezelter 889 use status
47     use linearAlgebra
48 gezelter 115 use neighborLists
49     use force_globals
50     use vector_class
51     use atype_module
52     use switcheroo
53     #ifdef IS_MPI
54     use mpiSimulation
55     #endif
56    
57     implicit none
58     PRIVATE
59    
60     #define __FORTRAN90
61     #include "brains/fSimulation.h"
62     #include "UseTheForce/fSwitchingFunction.h"
63    
64     type (simtype), public, save :: thisSim
65    
66     logical, save :: simulation_setup_complete = .false.
67    
68     integer, public, save :: nLocal, nGlobal
69     integer, public, save :: nGroups, nGroupGlobal
70 gezelter 1286 integer, public, save :: nExcludes = 0
71     integer, public, save :: nOneTwo = 0
72     integer, public, save :: nOneThree = 0
73     integer, public, save :: nOneFour = 0
74    
75     integer, allocatable, dimension(:,:), public :: excludes
76 gezelter 115 integer, allocatable, dimension(:), public :: molMembershipList
77     integer, allocatable, dimension(:), public :: groupListRow
78     integer, allocatable, dimension(:), public :: groupStartRow
79     integer, allocatable, dimension(:), public :: groupListCol
80     integer, allocatable, dimension(:), public :: groupStartCol
81     integer, allocatable, dimension(:), public :: groupListLocal
82     integer, allocatable, dimension(:), public :: groupStartLocal
83 gezelter 1346 integer, allocatable, dimension(:), public :: nSkipsForLocalAtom
84     integer, allocatable, dimension(:,:), public :: skipsForLocalAtom
85     integer, allocatable, dimension(:), public :: nSkipsForRowAtom
86     integer, allocatable, dimension(:,:), public :: skipsForRowAtom
87 gezelter 1286 integer, allocatable, dimension(:), public :: nTopoPairsForAtom
88     integer, allocatable, dimension(:,:), public :: toposForAtom
89     integer, allocatable, dimension(:,:), public :: topoDistance
90 gezelter 115 real(kind=dp), allocatable, dimension(:), public :: mfactRow
91     real(kind=dp), allocatable, dimension(:), public :: mfactCol
92     real(kind=dp), allocatable, dimension(:), public :: mfactLocal
93    
94     real(kind=dp), public, dimension(3,3), save :: Hmat, HmatInv
95 gezelter 889 real(kind=dp), save :: DangerRcut
96 gezelter 115 logical, public, save :: boxIsOrthorhombic
97 gezelter 507
98 gezelter 115 public :: SimulationSetup
99     public :: getNlocal
100     public :: setBox
101 gezelter 889 public :: checkBox
102 gezelter 115 public :: SimUsesPBC
103 gezelter 1126 public :: SimUsesAtomicVirial
104 gezelter 140 public :: SimUsesDirectionalAtoms
105 gezelter 1528 public :: SimUsesMetallicAtoms
106     public :: SimRequiresSkipCorrection
107     public :: SimRequiresSelfCorrection
108 gezelter 889 public :: setHmatDangerousRcutValue
109 gezelter 140
110 gezelter 115 contains
111 gezelter 507
112 gezelter 115 subroutine SimulationSetup(setThisSim, CnGlobal, CnLocal, c_idents, &
113 gezelter 1286 CnExcludes, Cexcludes, CnOneTwo, ConeTwo, CnOneThree, ConeThree, &
114     CnOneFour, ConeFour, CmolMembership, Cmfact, CnGroups, &
115     CglobalGroupMembership, status)
116 gezelter 115
117     type (simtype) :: setThisSim
118     integer, intent(inout) :: CnGlobal, CnLocal
119 gezelter 1473 integer, dimension(CnLocal), intent(inout) :: c_idents
120 gezelter 115
121 gezelter 1286 integer :: CnExcludes
122     integer, dimension(2,CnExcludes), intent(in) :: Cexcludes
123     integer :: CnOneTwo
124     integer, dimension(2,CnOneTwo), intent(in) :: ConeTwo
125     integer :: CnOneThree
126     integer, dimension(2,CnOneThree), intent(in) :: ConeThree
127     integer :: CnOneFour
128     integer, dimension(2,CnOneFour), intent(in) :: ConeFour
129    
130 gezelter 115 integer, dimension(CnGlobal),intent(in) :: CmolMembership
131     !! Result status, success = 0, status = -1
132     integer, intent(out) :: status
133     integer :: i, j, me, thisStat, alloc_stat, myNode, id1, id2
134 gezelter 1286 integer :: ia, jend
135 gezelter 115
136     !! mass factors used for molecular cutoffs
137     real ( kind = dp ), dimension(CnLocal) :: Cmfact
138     integer, intent(in):: CnGroups
139     integer, dimension(CnGlobal), intent(in):: CglobalGroupMembership
140 gezelter 1346 integer :: maxSkipsForLocalAtom, maxToposForAtom, glPointer
141     integer :: maxSkipsForRowAtom
142 gezelter 115
143     #ifdef IS_MPI
144     integer, allocatable, dimension(:) :: c_idents_Row
145     integer, allocatable, dimension(:) :: c_idents_Col
146     integer :: nAtomsInRow, nGroupsInRow, aid
147     integer :: nAtomsInCol, nGroupsInCol, gid
148     #endif
149    
150     simulation_setup_complete = .false.
151     status = 0
152    
153     ! copy C struct into fortran type
154    
155     nLocal = CnLocal
156     nGlobal = CnGlobal
157     nGroups = CnGroups
158    
159     thisSim = setThisSim
160    
161 gezelter 1286 nExcludes = CnExcludes
162 gezelter 115
163     call InitializeForceGlobals(nLocal, thisStat)
164     if (thisStat /= 0) then
165     write(default_error,*) "SimSetup: InitializeForceGlobals error"
166     status = -1
167     return
168     endif
169    
170     call InitializeSimGlobals(thisStat)
171     if (thisStat /= 0) then
172     write(default_error,*) "SimSetup: InitializeSimGlobals error"
173     status = -1
174     return
175     endif
176    
177     #ifdef IS_MPI
178     ! We can only set up forces if mpiSimulation has been setup.
179     if (.not. isMPISimSet()) then
180     write(default_error,*) "MPI is not set"
181     status = -1
182     return
183     endif
184     nAtomsInRow = getNatomsInRow(plan_atom_row)
185     nAtomsInCol = getNatomsInCol(plan_atom_col)
186     nGroupsInRow = getNgroupsInRow(plan_group_row)
187     nGroupsInCol = getNgroupsInCol(plan_group_col)
188     mynode = getMyNode()
189 gezelter 507
190 gezelter 115 call gather(c_idents, c_idents_Row, plan_atom_row)
191     call gather(c_idents, c_idents_Col, plan_atom_col)
192    
193     do i = 1, nAtomsInRow
194     me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
195     atid_Row(i) = me
196     enddo
197    
198     do i = 1, nAtomsInCol
199     me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
200     atid_Col(i) = me
201     enddo
202    
203     #endif
204    
205     #ifdef IS_MPI
206     allocate(groupStartRow(nGroupsInRow+1),stat=alloc_stat)
207     if (alloc_stat /= 0 ) then
208     status = -1
209     return
210     endif
211     allocate(groupStartCol(nGroupsInCol+1),stat=alloc_stat)
212     if (alloc_stat /= 0 ) then
213     status = -1
214     return
215     endif
216     allocate(groupListRow(nAtomsInRow),stat=alloc_stat)
217     if (alloc_stat /= 0 ) then
218     status = -1
219     return
220     endif
221     allocate(groupListCol(nAtomsInCol),stat=alloc_stat)
222     if (alloc_stat /= 0 ) then
223     status = -1
224     return
225     endif
226     allocate(mfactRow(nAtomsInRow),stat=alloc_stat)
227     if (alloc_stat /= 0 ) then
228     status = -1
229     return
230     endif
231     allocate(mfactCol(nAtomsInCol),stat=alloc_stat)
232     if (alloc_stat /= 0 ) then
233     status = -1
234     return
235     endif
236     allocate(mfactLocal(nLocal),stat=alloc_stat)
237     if (alloc_stat /= 0 ) then
238     status = -1
239     return
240     endif
241 gezelter 1286
242 gezelter 115 glPointer = 1
243 gezelter 1286
244 gezelter 115 do i = 1, nGroupsInRow
245 gezelter 1286
246 gezelter 115 gid = GroupRowToGlobal(i)
247     groupStartRow(i) = glPointer
248 gezelter 1286
249 gezelter 115 do j = 1, nAtomsInRow
250     aid = AtomRowToGlobal(j)
251     if (CglobalGroupMembership(aid) .eq. gid) then
252     groupListRow(glPointer) = j
253     glPointer = glPointer + 1
254     endif
255     enddo
256     enddo
257     groupStartRow(nGroupsInRow+1) = nAtomsInRow + 1
258 gezelter 1286
259 gezelter 115 glPointer = 1
260 gezelter 1286
261 gezelter 115 do i = 1, nGroupsInCol
262 gezelter 1286
263 gezelter 115 gid = GroupColToGlobal(i)
264     groupStartCol(i) = glPointer
265 gezelter 1286
266 gezelter 115 do j = 1, nAtomsInCol
267     aid = AtomColToGlobal(j)
268     if (CglobalGroupMembership(aid) .eq. gid) then
269     groupListCol(glPointer) = j
270     glPointer = glPointer + 1
271     endif
272     enddo
273     enddo
274     groupStartCol(nGroupsInCol+1) = nAtomsInCol + 1
275    
276     mfactLocal = Cmfact
277    
278     call gather(mfactLocal, mfactRow, plan_atom_row)
279     call gather(mfactLocal, mfactCol, plan_atom_col)
280 gezelter 507
281 gezelter 115 if (allocated(mfactLocal)) then
282     deallocate(mfactLocal)
283     end if
284     #else
285     allocate(groupStartRow(nGroups+1),stat=alloc_stat)
286     if (alloc_stat /= 0 ) then
287     status = -1
288     return
289     endif
290     allocate(groupStartCol(nGroups+1),stat=alloc_stat)
291     if (alloc_stat /= 0 ) then
292     status = -1
293     return
294     endif
295     allocate(groupListRow(nLocal),stat=alloc_stat)
296     if (alloc_stat /= 0 ) then
297     status = -1
298     return
299     endif
300     allocate(groupListCol(nLocal),stat=alloc_stat)
301     if (alloc_stat /= 0 ) then
302     status = -1
303     return
304     endif
305     allocate(mfactRow(nLocal),stat=alloc_stat)
306     if (alloc_stat /= 0 ) then
307     status = -1
308     return
309     endif
310     allocate(mfactCol(nLocal),stat=alloc_stat)
311     if (alloc_stat /= 0 ) then
312     status = -1
313     return
314     endif
315     allocate(mfactLocal(nLocal),stat=alloc_stat)
316     if (alloc_stat /= 0 ) then
317     status = -1
318     return
319     endif
320    
321     glPointer = 1
322     do i = 1, nGroups
323     groupStartRow(i) = glPointer
324     groupStartCol(i) = glPointer
325     do j = 1, nLocal
326     if (CglobalGroupMembership(j) .eq. i) then
327     groupListRow(glPointer) = j
328     groupListCol(glPointer) = j
329     glPointer = glPointer + 1
330     endif
331     enddo
332     enddo
333     groupStartRow(nGroups+1) = nLocal + 1
334     groupStartCol(nGroups+1) = nLocal + 1
335    
336     do i = 1, nLocal
337     mfactRow(i) = Cmfact(i)
338     mfactCol(i) = Cmfact(i)
339     end do
340 gezelter 507
341 gezelter 115 #endif
342    
343 gezelter 507 ! We build the local atid's for both mpi and nonmpi
344 gezelter 115 do i = 1, nLocal
345     me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
346     atid(i) = me
347 gezelter 1473 c_idents_local(i) = c_idents(i)
348 gezelter 115 enddo
349    
350 gezelter 1286 do i = 1, nExcludes
351     excludes(1,i) = Cexcludes(1,i)
352     excludes(2,i) = Cexcludes(2,i)
353 gezelter 115 enddo
354    
355     #ifdef IS_MPI
356 gezelter 1346 allocate(nSkipsForRowAtom(nAtomsInRow), stat=alloc_stat)
357 gezelter 115 #endif
358 gezelter 1346
359     allocate(nSkipsForLocalAtom(nLocal), stat=alloc_stat)
360    
361 gezelter 115 if (alloc_stat /= 0 ) then
362     thisStat = -1
363     write(*,*) 'Could not allocate nSkipsForAtom array'
364     return
365     endif
366    
367     #ifdef IS_MPI
368 gezelter 1346 maxSkipsForRowAtom = 0
369     do j = 1, nAtomsInRow
370     nSkipsForRowAtom(j) = 0
371     id1 = AtomRowToGlobal(j)
372     do i = 1, nExcludes
373     if (excludes(1,i) .eq. id1 ) then
374     nSkipsForRowAtom(j) = nSkipsForRowAtom(j) + 1
375     if (nSkipsForRowAtom(j) .gt. maxSkipsForRowAtom) then
376     maxSkipsForRowAtom = nSkipsForRowAtom(j)
377     endif
378     endif
379     if (excludes(2,i) .eq. id1 ) then
380     nSkipsForRowAtom(j) = nSkipsForRowAtom(j) + 1
381     if (nSkipsForRowAtom(j) .gt. maxSkipsForRowAtom) then
382     maxSkipsForRowAtom = nSkipsForRowAtom(j)
383     endif
384     endif
385     end do
386     enddo
387 gezelter 115 #endif
388 gezelter 1346 maxSkipsForLocalAtom = 0
389     do j = 1, nLocal
390     nSkipsForLocalAtom(j) = 0
391 gezelter 115 #ifdef IS_MPI
392 gezelter 1346 id1 = AtomLocalToGlobal(j)
393     #else
394 gezelter 1286 id1 = j
395 gezelter 115 #endif
396 gezelter 1286 do i = 1, nExcludes
397     if (excludes(1,i) .eq. id1 ) then
398 gezelter 1346 nSkipsForLocalAtom(j) = nSkipsForLocalAtom(j) + 1
399     if (nSkipsForLocalAtom(j) .gt. maxSkipsForLocalAtom) then
400     maxSkipsForLocalAtom = nSkipsForLocalAtom(j)
401 gezelter 115 endif
402 gezelter 1286 endif
403     if (excludes(2,i) .eq. id1 ) then
404 gezelter 1346 nSkipsForLocalAtom(j) = nSkipsForLocalAtom(j) + 1
405     if (nSkipsForLocalAtom(j) .gt. maxSkipsForLocalAtom) then
406     maxSkipsForLocalAtom = nSkipsForLocalAtom(j)
407 gezelter 115 endif
408 gezelter 1286 endif
409     end do
410     enddo
411    
412     #ifdef IS_MPI
413 gezelter 1346 allocate(skipsForRowAtom(nAtomsInRow, maxSkipsForRowAtom), stat=alloc_stat)
414 gezelter 1286 #endif
415 gezelter 1346 allocate(skipsForLocalAtom(nLocal, maxSkipsForLocalAtom), stat=alloc_stat)
416    
417 gezelter 1286 if (alloc_stat /= 0 ) then
418 gezelter 1346 write(*,*) 'Could not allocate skipsForAtom arrays'
419 gezelter 1286 return
420     endif
421    
422     #ifdef IS_MPI
423 gezelter 1346 do j = 1, nAtomsInRow
424     nSkipsForRowAtom(j) = 0
425     id1 = AtomRowToGlobal(j)
426     do i = 1, nExcludes
427     if (excludes(1,i) .eq. id1 ) then
428     nSkipsForRowAtom(j) = nSkipsForRowAtom(j) + 1
429     ! exclude lists have global ID's
430     id2 = excludes(2,i)
431     skipsForRowAtom(j, nSkipsForRowAtom(j)) = id2
432     endif
433     if (excludes(2, i) .eq. id1 ) then
434     nSkipsForRowAtom(j) = nSkipsForRowAtom(j) + 1
435     ! exclude lists have global ID's
436     id2 = excludes(1,i)
437     skipsForRowAtom(j, nSkipsForRowAtom(j)) = id2
438     endif
439     end do
440     enddo
441 gezelter 1286 #endif
442 gezelter 1346 do j = 1, nLocal
443     nSkipsForLocalAtom(j) = 0
444 gezelter 1286 #ifdef IS_MPI
445 gezelter 1346 id1 = AtomLocalToGlobal(j)
446     #else
447 gezelter 1286 id1 = j
448     #endif
449     do i = 1, nExcludes
450     if (excludes(1,i) .eq. id1 ) then
451 gezelter 1346 nSkipsForLocalAtom(j) = nSkipsForLocalAtom(j) + 1
452     ! exclude lists have global ID's
453     #ifdef IS_MPI
454     id2 = AtomGlobalToLocal(excludes(2,i))
455     #else
456 gezelter 1286 id2 = excludes(2,i)
457 gezelter 1346 #endif
458     skipsForLocalAtom(j, nSkipsForLocalAtom(j)) = id2
459 gezelter 1286 endif
460     if (excludes(2, i) .eq. id1 ) then
461 gezelter 1346 nSkipsForLocalAtom(j) = nSkipsForLocalAtom(j) + 1
462     ! exclude lists have global ID's
463     #ifdef IS_MPI
464     id2 = AtomGlobalToLocal(excludes(1,i))
465     #else
466 gezelter 1286 id2 = excludes(1,i)
467 gezelter 1346 #endif
468     skipsForLocalAtom(j, nSkipsForLocalAtom(j)) = id2
469 gezelter 1286 endif
470     end do
471     enddo
472    
473     do i = 1, nGlobal
474     molMemberShipList(i) = CmolMembership(i)
475     enddo
476 gezelter 115
477     #ifdef IS_MPI
478 gezelter 1286 allocate(nTopoPairsForAtom(nAtomsInRow), stat=alloc_stat)
479 gezelter 115 #else
480 gezelter 1286 allocate(nTopoPairsForAtom(nLocal), stat=alloc_stat)
481 gezelter 115 #endif
482 gezelter 1286 if (alloc_stat /= 0 ) then
483     thisStat = -1
484     write(*,*) 'Could not allocate nTopoPairsForAtom array'
485     return
486     endif
487 gezelter 115
488     #ifdef IS_MPI
489 gezelter 1286 jend = nAtomsInRow
490 gezelter 115 #else
491 gezelter 1286 jend = nLocal
492 gezelter 115 #endif
493 gezelter 1286
494     do j = 1, jend
495     nTopoPairsForAtom(j) = 0
496 gezelter 115 #ifdef IS_MPI
497 gezelter 1286 id1 = AtomRowToGlobal(j)
498 gezelter 115 #else
499 gezelter 1286 id1 = j
500 gezelter 115 #endif
501 gezelter 1286 do i = 1, CnOneTwo
502     if (ConeTwo(1,i) .eq. id1 ) then
503     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
504     endif
505     if (ConeTwo(2,i) .eq. id1 ) then
506     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
507     endif
508     end do
509 gezelter 507
510 gezelter 1286 do i = 1, CnOneThree
511     if (ConeThree(1,i) .eq. id1 ) then
512     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
513     endif
514     if (ConeThree(2,i) .eq. id1 ) then
515     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
516     endif
517     end do
518 gezelter 507
519 gezelter 1286 do i = 1, CnOneFour
520     if (ConeFour(1,i) .eq. id1 ) then
521     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
522     endif
523     if (ConeFour(2,i) .eq. id1 ) then
524     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
525     endif
526     end do
527     enddo
528    
529     maxToposForAtom = maxval(nTopoPairsForAtom)
530     #ifdef IS_MPI
531     allocate(toposForAtom(nAtomsInRow, maxToposForAtom), stat=alloc_stat)
532     allocate(topoDistance(nAtomsInRow, maxToposForAtom), stat=alloc_stat)
533     #else
534     allocate(toposForAtom(nLocal, maxToposForAtom), stat=alloc_stat)
535     allocate(topoDistance(nLocal, maxToposForAtom), stat=alloc_stat)
536     #endif
537     if (alloc_stat /= 0 ) then
538     write(*,*) 'Could not allocate topoDistance array'
539     return
540     endif
541    
542     #ifdef IS_MPI
543     jend = nAtomsInRow
544     #else
545     jend = nLocal
546     #endif
547     do j = 1, jend
548     nTopoPairsForAtom(j) = 0
549     #ifdef IS_MPI
550     id1 = AtomRowToGlobal(j)
551     #else
552     id1 = j
553     #endif
554     do i = 1, CnOneTwo
555     if (ConeTwo(1,i) .eq. id1 ) then
556     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
557     id2 = ConeTwo(2,i)
558     toposForAtom(j, nTopoPairsForAtom(j)) = id2
559     topoDistance(j, nTopoPairsForAtom(j)) = 1
560     endif
561     if (ConeTwo(2, i) .eq. id1 ) then
562     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
563     id2 = ConeTwo(1,i)
564     toposForAtom(j, nTopoPairsForAtom(j)) = id2
565     topoDistance(j, nTopoPairsForAtom(j)) = 1
566     endif
567     end do
568 gezelter 507
569 gezelter 1286 do i = 1, CnOneThree
570     if (ConeThree(1,i) .eq. id1 ) then
571     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
572     id2 = ConeThree(2,i)
573     toposForAtom(j, nTopoPairsForAtom(j)) = id2
574     topoDistance(j, nTopoPairsForAtom(j)) = 2
575     endif
576     if (ConeThree(2, i) .eq. id1 ) then
577     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
578     id2 = ConeThree(1,i)
579     toposForAtom(j, nTopoPairsForAtom(j)) = id2
580     topoDistance(j, nTopoPairsForAtom(j)) = 2
581     endif
582     end do
583 gezelter 507
584 gezelter 1286 do i = 1, CnOneFour
585     if (ConeFour(1,i) .eq. id1 ) then
586     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
587     id2 = ConeFour(2,i)
588     toposForAtom(j, nTopoPairsForAtom(j)) = id2
589     topoDistance(j, nTopoPairsForAtom(j)) = 3
590     endif
591     if (ConeFour(2, i) .eq. id1 ) then
592     nTopoPairsForAtom(j) = nTopoPairsForAtom(j) + 1
593     id2 = ConeFour(1,i)
594     toposForAtom(j, nTopoPairsForAtom(j)) = id2
595     topoDistance(j, nTopoPairsForAtom(j)) = 3
596     endif
597     end do
598     enddo
599    
600     if (status == 0) simulation_setup_complete = .true.
601    
602     end subroutine SimulationSetup
603 gezelter 115
604 gezelter 1286 subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
605     real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
606     integer :: cBoxIsOrthorhombic
607     integer :: smallest, status
608    
609     Hmat = cHmat
610     HmatInv = cHmatInv
611     if (cBoxIsOrthorhombic .eq. 0 ) then
612     boxIsOrthorhombic = .false.
613     else
614     boxIsOrthorhombic = .true.
615     endif
616    
617     call checkBox()
618     return
619     end subroutine setBox
620    
621 gezelter 889 subroutine checkBox()
622    
623     integer :: i
624     real(kind=dp), dimension(3) :: hx, hy, hz, ax, ay, az, piped
625     character(len = statusMsgSize) :: errMsg
626    
627     hx = Hmat(1,:)
628     hy = Hmat(2,:)
629     hz = Hmat(3,:)
630    
631     ax = cross_product(hy, hz)
632     ay = cross_product(hx, hz)
633     az = cross_product(hx, hy)
634    
635     ax = ax / length(ax)
636     ay = ay / length(ay)
637     az = az / length(az)
638    
639     piped(1) = abs(dot_product(ax, hx))
640     piped(2) = abs(dot_product(ay, hy))
641     piped(3) = abs(dot_product(az, hz))
642    
643     do i = 1, 3
644     if ((0.5_dp * piped(i)).lt.DangerRcut) then
645     write(errMsg, '(a94,f9.4,a1)') 'One of the dimensions of the Periodic ' &
646     // 'Box is smaller than ' // newline // tab // &
647     'the largest cutoff radius' // &
648     ' (rCut = ', DangerRcut, ')'
649     call handleError("checkBox", errMsg)
650 chuckv 1135
651 gezelter 889 end if
652     enddo
653     return
654     end subroutine checkBox
655 gezelter 1528
656 gezelter 507 function SimUsesPBC() result(doesit)
657     logical :: doesit
658     doesit = thisSim%SIM_uses_PBC
659     end function SimUsesPBC
660 gezelter 1528
661 gezelter 1126 function SimUsesAtomicVirial() result(doesit)
662     logical :: doesit
663     doesit = thisSim%SIM_uses_AtomicVirial
664     end function SimUsesAtomicVirial
665 gezelter 1528
666 gezelter 507 function SimUsesDirectionalAtoms() result(doesit)
667     logical :: doesit
668 gezelter 1528 doesit = thisSim%SIM_uses_DirectionalAtoms
669 gezelter 507 end function SimUsesDirectionalAtoms
670 gezelter 1528
671     function SimUsesMetallicAtoms() result(doesit)
672 gezelter 507 logical :: doesit
673 gezelter 1528 doesit = thisSim%SIM_uses_MetallicAtoms
674     end function SimUsesMetallicAtoms
675    
676     function SimRequiresSkipCorrection() result(doesit)
677 gezelter 507 logical :: doesit
678 gezelter 1528 doesit = thisSim%SIM_requires_SkipCorrection
679     end function SimRequiresSkipCorrection
680 chrisfen 691
681 gezelter 1528 function SimRequiresSelfCorrection() result(doesit)
682 gezelter 507 logical :: doesit
683 gezelter 1528 doesit = thisSim%SIM_requires_SelfCorrection
684     end function SimRequiresSelfCorrection
685 gezelter 575
686 gezelter 1528 subroutine InitializeSimGlobals(thisStat)
687 gezelter 507 integer, intent(out) :: thisStat
688     integer :: alloc_stat
689 gezelter 1528
690 gezelter 507 thisStat = 0
691 gezelter 1528
692 gezelter 507 call FreeSimGlobals()
693 gezelter 1528
694 gezelter 1286 allocate(excludes(2,nExcludes), stat=alloc_stat)
695 gezelter 507 if (alloc_stat /= 0 ) then
696     thisStat = -1
697     return
698     endif
699    
700     allocate(molMembershipList(nGlobal), stat=alloc_stat)
701     if (alloc_stat /= 0 ) then
702     thisStat = -1
703     return
704     endif
705    
706     end subroutine InitializeSimGlobals
707    
708     subroutine FreeSimGlobals()
709    
710     !We free in the opposite order in which we allocate in.
711 gezelter 1286 if (allocated(topoDistance)) deallocate(topoDistance)
712     if (allocated(toposForAtom)) deallocate(toposForAtom)
713     if (allocated(nTopoPairsForAtom)) deallocate(nTopoPairsForAtom)
714 gezelter 1346 if (allocated(skipsForLocalAtom)) deallocate(skipsForLocalAtom)
715     if (allocated(nSkipsForLocalAtom)) deallocate(nSkipsForLocalAtom)
716     if (allocated(skipsForRowAtom)) deallocate(skipsForRowAtom)
717     if (allocated(nSkipsForRowAtom)) deallocate(nSkipsForRowAtom)
718 gezelter 507 if (allocated(mfactLocal)) deallocate(mfactLocal)
719     if (allocated(mfactCol)) deallocate(mfactCol)
720     if (allocated(mfactRow)) deallocate(mfactRow)
721     if (allocated(groupListCol)) deallocate(groupListCol)
722     if (allocated(groupListRow)) deallocate(groupListRow)
723     if (allocated(groupStartCol)) deallocate(groupStartCol)
724     if (allocated(groupStartRow)) deallocate(groupStartRow)
725 gezelter 1286 if (allocated(molMembershipList)) deallocate(molMembershipList)
726     if (allocated(excludes)) deallocate(excludes)
727 gezelter 507
728     end subroutine FreeSimGlobals
729    
730     pure function getNlocal() result(n)
731     integer :: n
732     n = nLocal
733     end function getNlocal
734    
735 gezelter 889 subroutine setHmatDangerousRcutValue(dangerWillRobinson)
736     real(kind=dp), intent(in) :: dangerWillRobinson
737     DangerRcut = dangerWillRobinson
738 chuckv 563
739 gezelter 889 call checkBox()
740 chuckv 563
741 gezelter 889 return
742     end subroutine setHmatDangerousRcutValue
743 gezelter 1528
744 gezelter 507 end module simulation
745 gezelter 1528

Properties

Name Value
svn:keywords Author Id Revision Date