ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/forceGlobals.F90
(Generate patch)

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/forceGlobals.F90 (file contents):
Revision 306 by chuckv, Mon Mar 10 19:26:45 2003 UTC vs.
Revision 317 by gezelter, Tue Mar 11 23:13:06 2003 UTC

# Line 1 | Line 1
1   module forceGlobals
2 <
3 <
4 <
5 < !! Number of lj_atypes in lj_atype_list
6 <  integer, save :: n_atypes = 0
7 <
8 < !! Global list of lj atypes in simulation
9 <  type (atype), pointer :: ListHead => null()
10 <  type (atype), pointer :: ListTail => null()
11 <
12 <
13 <
14 <
15 <  logical, save :: firstTime = .True.
16 <
17 < !! Atype identity pointer lists
2 >  
3 >  use definitions
4 >  use simulation
5 >  use atype_typedefs
6 >  use vector_class
7   #ifdef IS_MPI
8 < !! Row lj_atype pointer list
20 <  type (identPtrList), dimension(:), pointer :: identPtrListRow => null()
21 < !! Column lj_atype pointer list
22 <  type (identPtrList), dimension(:), pointer :: identPtrListColumn => null()
23 < #else
24 <  type(identPtrList ), dimension(:), pointer :: identPtrList => null()
8 >  use mpiSimulation
9   #endif
10 <
11 <
12 < !! Logical has lj force field module been initialized?
13 <  logical, save :: isFFinit = .false.
30 <
31 < !! Use periodic boundry conditions
32 <  logical :: wrap = .false.
33 <
34 < !! Potential energy global module variables
10 >    
11 >  logical, save :: atype_vector_initialized = .false.
12 >  logical, save :: ff_initialized = .false.
13 >  
14   #ifdef IS_MPI
15 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
16 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
17 <
18 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: muRow = 0.0_dp
19 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: muCol = 0.0_dp
20 <
21 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
22 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp
23 <
24 <  real(kind = dp), dimension(9,getNrow(plan_row)) :: ARow = 0.0_dp
25 <  real(kind = dp), dimension(9,getNcol(plan_col)) :: ACol = 0.0_dp  
26 <
27 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
28 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
29 <  real(kind = dp), dimension(3,getNlocal()) :: fTemp1 = 0.0_dp
30 <  real(kind = dp), dimension(3,getNlocal()) :: tTemp1 = 0.0_dp
31 <  real(kind = dp), dimension(3,getNlocal()) :: fTemp2 = 0.0_dp
32 <  real(kind = dp), dimension(3,getNlocal()) :: tTemp2 = 0.0_dp
33 <  real(kind = dp), dimension(3,getNlocal()) :: fTemp = 0.0_dp
34 <  real(kind = dp), dimension(3,getNlocal()) :: tTemp = 0.0_dp
56 <
57 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
58 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp
59 <
60 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: rflRow = 0.0_dp
61 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: rflCol = 0.0_dp
62 <  real(kind = dp), dimension(3,getNlocal()) :: rflTemp = 0.0_dp
63 <
64 <  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
65 <  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
66 <
67 <  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
15 >  real( kind = dp ), allocatable, dimension(:,:) :: q_Row
16 >  real( kind = dp ), allocatable, dimension(:,:) :: q_Col
17 >  real( kind = dp ), allocatable, dimension(:,:) :: u_l_Row
18 >  real( kind = dp ), allocatable, dimension(:,:) :: u_l_Col
19 >  real( kind = dp ), allocatable, dimension(:,:) :: A_Row
20 >  real( kind = dp ), allocatable, dimension(:,:) :: A_Col
21 >  
22 >  real( kind = dp ), allocatable, dimension(:) :: pot_Row
23 >  real( kind = dp ), allocatable, dimension(:) :: pot_Col
24 >  real( kind = dp ), allocatable, dimension(:) :: pot_Temp
25 >  real( kind = dp ), allocatable, dimension(:,:) :: f_Row
26 >  real( kind = dp ), allocatable, dimension(:,:) :: f_Col
27 >  real( kind = dp ), allocatable, dimension(:,:) :: f_Temp
28 >  real( kind = dp ), allocatable, dimension(:,:) :: t_Row
29 >  real( kind = dp ), allocatable, dimension(:,:) :: t_Col
30 >  real( kind = dp ), allocatable, dimension(:,:) :: t_Temp
31 >  integer, allocatable, dimension(:) :: atid_Row
32 >  integer, allocatable, dimension(:) :: atid_Col
33 > #else
34 >  integer, allocatable, dimension(:) :: atid
35   #endif
36 <  real(kind = dp) :: pe = 0.0_dp
37 <  real(kind = dp), dimension(3,natoms) :: fTemp = 0.0_dp
38 <  real(kind = dp), dimension(3,natoms) :: tTemp = 0.0_dp
39 <  real(kind = dp), dimension(3,natoms) :: rflTemp = 0.0_dp
40 <  real(kind = dp), dimension(9) :: tauTemp = 0.0_dp
74 <
75 <  logical :: do_preForce  = .false.
76 <  logical :: do_postForce = .false.
77 <
78 <
79 <
80 < !! Public methods and data
36 >  
37 >  real(kind = dp) :: pot = 0.0_dp
38 >  real(kind = dp), dimension(9) :: tau_Temp = 0.0_dp
39 >  real(kind = dp) :: virial_Temp = 0.0_dp
40 >  
41    public :: new_atype
82  public :: do_forceLoop
42    public :: init_FF
43  
44 + contains
45    
46 +  subroutine new_atype(c_ident, is_LJ, is_Sticky, is_DP, is_GB, &
47 +       lj_epsilon, lj_sigma, &
48 +       dipole_moment, &
49 +       sticky_w0, sticky_v0, &
50 +       GB_sigma, GB_l2b_ratio, GB_eps, GB_eps_ratio, GB_mu, GB_nu, &
51 +       status)
52 +    
53 +    real( kind = dp ), intent(in) :: lj_epsilon
54 +    real( kind = dp ), intent(in) :: lj_sigma
55 +    real( kind = dp ), intent(in) :: sticky_w0
56 +    real( kind = dp ), intent(in) :: sticky_v0
57 +    real( kind = dp ), intent(in) :: dipole_moment
58 +    real( kind = dp ), intent(in) :: GB_sigma, GB_l2b_ratio, GB_eps
59 +    real( kind = dp ), intent(in) :: GB_eps_ratio, GB_mu, GB_nu
60  
61 <
88 < contains
89 <
90 < !! Adds a new lj_atype to the list.
91 <  subroutine new_atype(ident,mass,epsilon,sigma, &
92 <       is_LJ,is_Sticky,is_DP,is_GB,w0,v0,dipoleMoment,status)
93 <    real( kind = dp ), intent(in) :: mass
94 <    real( kind = dp ), intent(in) :: epsilon
95 <    real( kind = dp ), intent(in) :: sigma
96 <    real( kind = dp ), intent(in) :: w0
97 <    real( kind = dp ), intent(in) :: v0
98 <    real( kind = dp ), intent(in) :: dipoleMoment
99 <
100 <    integer, intent(in) :: ident
61 >    integer, intent(in)  :: c_ident
62      integer, intent(out) :: status
63 <    integer, intent(in) :: is_Sticky
64 <    integer, intent(in) :: is_DP
65 <    integer, intent(in) :: is_GB
66 <    integer, intent(in) :: is_LJ
63 >    integer, intent(in)  :: is_Sticky
64 >    integer, intent(in)  :: is_DP
65 >    integer, intent(in)  :: is_GB
66 >    integer, intent(in)  :: is_LJ
67 >    integer :: me
68 >    logical :: l_is_LJ, l_is_DP, l_is_Sticky, l_is_GB
69  
107
108    type (atype), pointer :: the_new_atype
109    integer :: alloc_error
110    integer :: atype_counter = 0
111    integer :: alloc_size
112    integer :: err_stat
70      status = 0
71  
72 +    if (.not.atype_vector_initialized) then
73 +       !! There are 16 properties to worry about for now.  
74 +       !! Fix this if needed for more atomic properties
75 +       atypes = initialize(16)
76 +       if (.not.associated(atypes)) then
77 +          status = -1
78 +          return
79 +       endif
80 +    endif
81  
82 +    me = addElement(atypes)
83  
84 < ! allocate a new atype    
118 <    allocate(the_new_atype,stat=alloc_error)
119 <    if (alloc_error /= 0 ) then
120 <       status = -1
121 <       return
122 <    end if
84 >    call setElementProperty(atypes, me, "c_ident", c_ident)
85  
86 < ! assign our new atype information
87 <    the_new_atype%mass        = mass
88 <    the_new_atype%epsilon     = epsilon
89 <    the_new_atype%sigma       = sigma
128 <    the_new_atype%sigma2      = sigma * sigma
129 <    the_new_atype%sigma6      = the_new_atype%sigma2 * the_new_atype%sigma2 &
130 <         * the_new_atype%sigma2
131 <    the_new_atype%w0       = w0
132 <    the_new_atype%v0       = v0
133 <    the_new_atype%dipoleMoment       = dipoleMoment
134 <
135 <    
136 < ! assume that this atype will be successfully added
137 <    the_new_atype%atype_ident = ident
138 <    the_new_atype%atype_number = n_lj_atypes + 1
139 <    
140 <    if ( is_Sticky /= 0 )    the_new_atype%is_Sticky   = .true.
141 <    if ( is_GB /= 0 )        the_new_atype%is_GB       = .true.
142 <    if ( is_LJ /= 0 )        the_new_atype%is_LJ       = .true.
143 <    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.
86 >    l_is_LJ = (is_LJ .ne. 0)
87 >    l_is_DP = (is_DP .ne. 0)
88 >    l_is_Sticky = (is_Sticky .ne. 0)
89 >    l_is_GB = (is_GB .ne. 0)
90  
91 <    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
92 <    if (err_stat /= 0 ) then
93 <       status = -1
94 <       return
91 >    call setElementProperty(atypes, me, "is_LJ", l_is_LJ)
92 >    call setElementProperty(atypes, me, "is_DP", l_is_DP)
93 >    call setElementProperty(atypes, me, "is_Sticky", l_is_Sticky)
94 >    call setElementProperty(atypes, me, "is_GB", l_is_GB)
95 >
96 >    if (l_is_LJ) then
97 >       call setElementProperty(atypes, me, "lj_sigma", lj_sigma)
98 >       call setElementProperty(atypes, me, "lj_epsilon", lj_epsilon)
99      endif
100 +    if (l_is_DP) then
101 +       call setElementProperty(atypes, me, "dipole_moment", dipole_moment)
102 +    endif
103 +    if (l_is_Sticky) then
104 +       call setElementProperty(atypes, me, "sticky_w0", sticky_w0)
105 +       call setElementProperty(atypes, me, "sticky_v0", sticky_v0)
106 +    endif
107 +    if (l_is_GB) then
108 +       call setElementProperty(atypes, me, "GB_sigma", GB_sigma)
109 +       call setElementProperty(atypes, me, "GB_l2b_ratio", GB_l2b_ratio)
110 +       call setElementProperty(atypes, me, "GB_eps", GB_eps)
111 +       call setElementProperty(atypes, me, "GB_eps_ratio", GB_eps_ratio)
112 +       call setElementProperty(atypes, me, "GB_mu", GB_mu)
113 +       call setElementProperty(atypes, me, "GB_nu", GB_nu)
114 +    endif
115  
151    n_atypes = n_atypes + 1
152
153
116    end subroutine new_atype
117  
118 <
119 <  subroutine init_FF(nComponents,ident, nExcludes,excludesLocal,status)
120 < !! Number of components in ident array
118 >  
119 >  subroutine init_FF(nComponents, c_idents, nExcludes, excludesLocal, status)
120 >    !! Number of components in c_idents array
121      integer, intent(inout) :: nComponents
122 < !! Array of identities nComponents long corresponding to
123 < !! ljatype ident.
124 <    integer, dimension(nComponents),intent(inout) :: ident
122 >    !! Array of identities nComponents long corresponding to
123 >    !! ljatype ident.
124 >    integer, dimension(nComponents),intent(inout) :: c_idents
125      integer :: nExcludes
126 <    integer, dimension(n),intent(inout) :: excludesLocal
127 < !!  Result status, success = 0, error = -1
128 <    integer, intent(out) :: Status
129 <
168 <    integer :: alloc_stat
169 <
170 <    integer :: thisStat
171 <    integer :: i
172 <
126 >    integer, dimension(nExcludes),intent(inout) :: excludesLocal
127 >    !!  Result status, success = 0, error = -1
128 >    integer, intent(out) :: status    
129 >    integer :: i, me, thisStat
130      integer :: myNode
131 +
132   #ifdef IS_MPI
133 <    integer, allocatable, dimension(:) :: identRow
134 <    integer, allocatable, dimension(:) :: identCol
133 >    integer, allocatable, dimension(:) :: c_idents_Row
134 >    integer, allocatable, dimension(:) :: c_idents_Col
135      integer :: nrow
136      integer :: ncol
137   #endif
138 +
139      status = 0
181  
140  
183    
141  
142 < !! if were're not in MPI, we just update ljatypePtrList
143 < #ifndef IS_MPI
187 <    call create_IdentPtrlst(ident,ListHead,identPtrList,thisStat)
188 <    if ( thisStat /= 0 ) then
189 <       status = -1
190 <       return
191 <    endif
192 <
193 <    
194 < ! if were're in MPI, we also have to worry about row and col lists    
195 < #else
196 <  
197 < ! We can only set up forces if mpiSimulation has been setup.
142 > #ifdef IS_MPI
143 >    ! We can only set up forces if mpiSimulation has been setup.
144      if (.not. isMPISimSet()) then
145         write(default_error,*) "MPI is not set"
146         status = -1
# Line 203 | Line 149 | contains
149      nrow = getNrow(plan_row)
150      ncol = getNcol(plan_col)
151      mynode = getMyNode()
152 < !! Allocate temperary arrays to hold gather information
153 <    allocate(identRow(nrow),stat=alloc_stat)
152 >    
153 >    allocate(c_idents_Row(nrow),stat=alloc_stat)
154      if (alloc_stat /= 0 ) then
155         status = -1
156         return
157      endif
158  
159 <    allocate(identCol(ncol),stat=alloc_stat)
159 >    allocate(c_idents_Col(ncol),stat=alloc_stat)
160      if (alloc_stat /= 0 ) then
161         status = -1
162         return
163      endif
164  
165 < !! Gather idents into row and column idents
166 <
167 <    call gather(ident,identRow,plan_row)
168 <    call gather(ident,identCol,plan_col)
169 <    
224 <  
225 < !! Create row and col pointer lists
226 <  
227 <    call create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
228 <    if (thisStat /= 0 ) then
165 >    call gather(c_idents, c_idents_Row, plan_row)
166 >    call gather(c_idents, c_idents_Col, plan_col)
167 >
168 >    allocate(atid_Row(nrow),stat=alloc_stat)
169 >    if (alloc_stat /= 0 ) then
170         status = -1
171         return
172      endif
173 <  
174 <    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
175 <    if (thisStat /= 0 ) then
173 >
174 >    allocate(atid_Col(ncol),stat=alloc_stat)
175 >    if (alloc_stat /= 0 ) then
176         status = -1
177         return
178      endif
179  
180 < !! free temporary ident arrays
180 >    do i = 1, nrow
181 >       me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
182 >       atid_Row(i) = me
183 >    enddo
184 >
185 >    do i = 1, ncol
186 >       me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
187 >       atid_Col(i) = me
188 >    enddo
189 >
190 >    !! free temporary ident arrays
191      if (allocated(identCol)) then
192         deallocate(identCol)
193      end if
194      if (allocated(identCol)) then
195         deallocate(identRow)
196      endif
197 +    
198 + #else
199 +    do i = 1, nComponents
200  
201 +       me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
202 +       atid(i) = me
203 +
204 +    enddo
205   #endif
206 <    
206 >
207      call initForce_Modules(thisStat)
208      if (thisStat /= 0) then
209         status = -1
210         return
211      endif
212 <
213 < !! Create neighbor lists
212 >    
213 >    !! Create neighbor lists
214      call expandList(thisStat)
215      if (thisStat /= 0) then
216         status = -1
217         return
218      endif
219 <
220 <    isFFinit = .true.
221 <
222 <
219 >    
220 >    call setupGlobals(thisStat)
221 >    if (thisStat /= 0) then
222 >       status = -1
223 >       return
224 >    endif
225 >    
226 >    ff_initialized = .true.
227 >    
228 >    
229    end subroutine init_FF
230 <
231 <
232 <
233 <
230 >  
231 >  subroutine setupGlobals(thisStat)
232 >    integer, intent(out) :: thisStat
233 >    
234 >  end subroutine setupGlobals
235 >  
236 >  
237    subroutine initForce_Modules(thisStat)
238      integer, intent(out) :: thisStat
239      integer :: my_status
240      
241      thisStat = 0
242 <    call init_lj_FF(ListHead,my_status)
242 >    call init_lj_FF(my_status)
243      if (my_status /= 0) then
244         thisStat = -1
245         return
246      end if
247 <
247 >    
248    end subroutine initForce_Modules
249 <
283 <
284 <
285 <
286 <
287 <
288 <
249 >  
250   end module forceGlobals

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines