ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/UseTheForce/DarkSide/simulation.F90
Revision: 1532
Committed: Wed Dec 29 19:59:21 2010 UTC (14 years, 6 months ago) by gezelter
File size: 23201 byte(s)
Log Message:
Added MAW to the C++ side, removed a whole bunch more fortran. 

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

Properties

Name Value
svn:keywords Author Id Revision Date