mdtools/libmdCode/forceGlobals.F90

module forceGlobals
  
  use definitions
  use simulation
  use atype_typedefs
  use vector_class
#ifdef IS_MPI
  use mpiSimulation
#endif
    
  logical, save :: atype_vector_initialized = .false.
  logical, save :: ff_initialized = .false.
  
#ifdef IS_MPI
  real( kind = dp ), allocatable, dimension(:,:) :: q_Row
  real( kind = dp ), allocatable, dimension(:,:) :: q_Col
  real( kind = dp ), allocatable, dimension(:,:) :: u_l_Row
  real( kind = dp ), allocatable, dimension(:,:) :: u_l_Col
  real( kind = dp ), allocatable, dimension(:,:) :: A_Row
  real( kind = dp ), allocatable, dimension(:,:) :: A_Col
  
  real( kind = dp ), allocatable, dimension(:) :: pot_Row
  real( kind = dp ), allocatable, dimension(:) :: pot_Col
  real( kind = dp ), allocatable, dimension(:) :: pot_Temp
  real( kind = dp ), allocatable, dimension(:,:) :: f_Row
  real( kind = dp ), allocatable, dimension(:,:) :: f_Col
  real( kind = dp ), allocatable, dimension(:,:) :: f_Temp
  real( kind = dp ), allocatable, dimension(:,:) :: t_Row
  real( kind = dp ), allocatable, dimension(:,:) :: t_Col
  real( kind = dp ), allocatable, dimension(:,:) :: t_Temp
  integer, allocatable, dimension(:) :: atid_Row
  integer, allocatable, dimension(:) :: atid_Col
#else
  integer, allocatable, dimension(:) :: atid
#endif
  
  real(kind = dp) :: pot = 0.0_dp
  real(kind = dp), dimension(9) :: tau_Temp = 0.0_dp
  real(kind = dp) :: virial_Temp = 0.0_dp
  
  public :: new_atype 
  public :: init_FF

contains
  
  subroutine new_atype(c_ident, is_LJ, is_Sticky, is_DP, is_GB, &
       lj_epsilon, lj_sigma, & 
       dipole_moment, &
       sticky_w0, sticky_v0, &
       GB_sigma, GB_l2b_ratio, GB_eps, GB_eps_ratio, GB_mu, GB_nu, &
       status)
    
    real( kind = dp ), intent(in) :: lj_epsilon
    real( kind = dp ), intent(in) :: lj_sigma
    real( kind = dp ), intent(in) :: sticky_w0
    real( kind = dp ), intent(in) :: sticky_v0
    real( kind = dp ), intent(in) :: dipole_moment
    real( kind = dp ), intent(in) :: GB_sigma, GB_l2b_ratio, GB_eps
    real( kind = dp ), intent(in) :: GB_eps_ratio, GB_mu, GB_nu

    integer, intent(in)  :: c_ident
    integer, intent(out) :: status
    integer, intent(in)  :: is_Sticky
    integer, intent(in)  :: is_DP
    integer, intent(in)  :: is_GB
    integer, intent(in)  :: is_LJ
    integer :: me
    logical :: l_is_LJ, l_is_DP, l_is_Sticky, l_is_GB

    status = 0

    if (.not.atype_vector_initialized) then
       !! There are 16 properties to worry about for now.  
       !! Fix this if needed for more atomic properties
       atypes = initialize(16)
       if (.not.associated(atypes)) then
          status = -1
          return
       endif
    endif

    me = addElement(atypes)

    call setElementProperty(atypes, me, "c_ident", c_ident)

    l_is_LJ = (is_LJ .ne. 0)
    l_is_DP = (is_DP .ne. 0)
    l_is_Sticky = (is_Sticky .ne. 0)
    l_is_GB = (is_GB .ne. 0)

    call setElementProperty(atypes, me, "is_LJ", l_is_LJ)
    call setElementProperty(atypes, me, "is_DP", l_is_DP)
    call setElementProperty(atypes, me, "is_Sticky", l_is_Sticky)
    call setElementProperty(atypes, me, "is_GB", l_is_GB)

    if (l_is_LJ) then
       call setElementProperty(atypes, me, "lj_sigma", lj_sigma)
       call setElementProperty(atypes, me, "lj_epsilon", lj_epsilon)
    endif
    if (l_is_DP) then
       call setElementProperty(atypes, me, "dipole_moment", dipole_moment)
    endif
    if (l_is_Sticky) then
       call setElementProperty(atypes, me, "sticky_w0", sticky_w0)
       call setElementProperty(atypes, me, "sticky_v0", sticky_v0)
    endif
    if (l_is_GB) then
       call setElementProperty(atypes, me, "GB_sigma", GB_sigma)
       call setElementProperty(atypes, me, "GB_l2b_ratio", GB_l2b_ratio)
       call setElementProperty(atypes, me, "GB_eps", GB_eps)
       call setElementProperty(atypes, me, "GB_eps_ratio", GB_eps_ratio)
       call setElementProperty(atypes, me, "GB_mu", GB_mu)
       call setElementProperty(atypes, me, "GB_nu", GB_nu)
    endif

  end subroutine new_atype

  
  subroutine init_FF(nComponents, c_idents, nExcludes, excludesLocal, status)
    !! Number of components in c_idents array
    integer, intent(inout) :: nComponents
    !! Array of identities nComponents long corresponding to
    !! ljatype ident.
    integer, dimension(nComponents),intent(inout) :: c_idents
    integer :: nExcludes
    integer, dimension(nExcludes),intent(inout) :: excludesLocal
    !!  Result status, success = 0, error = -1
    integer, intent(out) :: status    
    integer :: i, me, thisStat
    integer :: myNode

#ifdef IS_MPI
    integer, allocatable, dimension(:) :: c_idents_Row
    integer, allocatable, dimension(:) :: c_idents_Col
    integer :: nrow
    integer :: ncol
#endif
 
    status = 0


#ifdef IS_MPI
    ! We can only set up forces if mpiSimulation has been setup.
    if (.not. isMPISimSet()) then
       write(default_error,*) "MPI is not set"
       status = -1
       return
    endif
    nrow = getNrow(plan_row)
    ncol = getNcol(plan_col)
    mynode = getMyNode()
    
    allocate(c_idents_Row(nrow),stat=alloc_stat)
    if (alloc_stat /= 0 ) then
       status = -1
       return
    endif

    allocate(c_idents_Col(ncol),stat=alloc_stat)
    if (alloc_stat /= 0 ) then
       status = -1
       return
    endif

    call gather(c_idents, c_idents_Row, plan_row)
    call gather(c_idents, c_idents_Col, plan_col)

    allocate(atid_Row(nrow),stat=alloc_stat)
    if (alloc_stat /= 0 ) then
       status = -1
       return
    endif

    allocate(atid_Col(ncol),stat=alloc_stat)
    if (alloc_stat /= 0 ) then
       status = -1
       return
    endif

    do i = 1, nrow
       me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
       atid_Row(i) = me
    enddo

    do i = 1, ncol
       me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
       atid_Col(i) = me
    enddo

    !! free temporary ident arrays
    if (allocated(identCol)) then
       deallocate(identCol)
    end if
    if (allocated(identCol)) then
       deallocate(identRow)
    endif
    
#else
    do i = 1, nComponents

       me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
       atid(i) = me

    enddo
#endif

    call initForce_Modules(thisStat)
    if (thisStat /= 0) then
       status = -1
       return
    endif
    
    !! Create neighbor lists
    call expandList(thisStat)
    if (thisStat /= 0) then
       status = -1
       return
    endif
    
    call setupGlobals(thisStat)
    if (thisStat /= 0) then
       status = -1
       return
    endif
    
    ff_initialized = .true.
    
    
  end subroutine init_FF
  
  subroutine setupGlobals(thisStat)
    integer, intent(out) :: thisStat
    
  end subroutine setupGlobals
  
  
  subroutine initForce_Modules(thisStat)
    integer, intent(out) :: thisStat
    integer :: my_status
    
    thisStat = 0
    call init_lj_FF(my_status)
    if (my_status /= 0) then
       thisStat = -1
       return
    end if
    
  end subroutine initForce_Modules
   
end module forceGlobals
Revision:	317
Committed:	Tue Mar 11 23:13:06 2003 UTC (21 years, 6 months ago) by gezelter
File size:	7071 byte(s)
Log Message:	Bug fixes
#	User	Rev	Content
1	chuckv	306	module forceGlobals
2	gezelter	317
3	gezelter	309	use definitions
4			use simulation
5			use atype_typedefs
6	gezelter	317	use vector_class
7	gezelter	309	#ifdef IS_MPI
8			use mpiSimulation
9			#endif
10	gezelter	317
11			logical, save :: atype_vector_initialized = .false.
12			logical, save :: ff_initialized = .false.
13
14	chuckv	306	#ifdef IS_MPI
15	gezelter	317	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	chuckv	306	#else
34	gezelter	317	integer, allocatable, dimension(:) :: atid
35	chuckv	306	#endif
36	gezelter	317
37	gezelter	309	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	gezelter	317
41	chuckv	306	public :: new_atype
42			public :: init_FF
43
44	gezelter	317	contains
45	chuckv	306
46	gezelter	317	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	chuckv	306
61	gezelter	317	integer, intent(in) :: c_ident
62	chuckv	306	integer, intent(out) :: status
63	gezelter	317	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	chuckv	306
70			status = 0
71
72	gezelter	317	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	chuckv	306
82	gezelter	317	me = addElement(atypes)
83	chuckv	306
84	gezelter	317	call setElementProperty(atypes, me, "c_ident", c_ident)
85	chuckv	306
86	gezelter	317	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	gezelter	309
91	gezelter	317	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	chuckv	306
96	gezelter	317	if (l_is_LJ) then
97			call setElementProperty(atypes, me, "lj_sigma", lj_sigma)
98			call setElementProperty(atypes, me, "lj_epsilon", lj_epsilon)
99	chuckv	306	endif
100	gezelter	317	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
116	chuckv	306	end subroutine new_atype
117
118	gezelter	317
119			subroutine init_FF(nComponents, c_idents, nExcludes, excludesLocal, status)
120			!! Number of components in c_idents array
121	chuckv	306	integer, intent(inout) :: nComponents
122	gezelter	317	!! Array of identities nComponents long corresponding to
123			!! ljatype ident.
124			integer, dimension(nComponents),intent(inout) :: c_idents
125	chuckv	306	integer :: nExcludes
126	gezelter	309	integer, dimension(nExcludes),intent(inout) :: excludesLocal
127	gezelter	317	!! Result status, success = 0, error = -1
128			integer, intent(out) :: status
129			integer :: i, me, thisStat
130			integer :: myNode
131	chuckv	306
132			#ifdef IS_MPI
133	gezelter	317	integer, allocatable, dimension(:) :: c_idents_Row
134			integer, allocatable, dimension(:) :: c_idents_Col
135	chuckv	306	integer :: nrow
136			integer :: ncol
137			#endif
138	gezelter	317
139	chuckv	306	status = 0
140
141
142	gezelter	317	#ifdef IS_MPI
143			! We can only set up forces if mpiSimulation has been setup.
144	chuckv	306	if (.not. isMPISimSet()) then
145			write(default_error,*) "MPI is not set"
146			status = -1
147			return
148			endif
149			nrow = getNrow(plan_row)
150			ncol = getNcol(plan_col)
151			mynode = getMyNode()
152	gezelter	317
153			allocate(c_idents_Row(nrow),stat=alloc_stat)
154	chuckv	306	if (alloc_stat /= 0 ) then
155			status = -1
156			return
157			endif
158
159	gezelter	317	allocate(c_idents_Col(ncol),stat=alloc_stat)
160	chuckv	306	if (alloc_stat /= 0 ) then
161			status = -1
162			return
163			endif
164
165	gezelter	317	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	chuckv	306	status = -1
171			return
172			endif
173	gezelter	317
174			allocate(atid_Col(ncol),stat=alloc_stat)
175			if (alloc_stat /= 0 ) then
176	chuckv	306	status = -1
177			return
178			endif
179
180	gezelter	317	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	chuckv	306	if (allocated(identCol)) then
192			deallocate(identCol)
193			end if
194			if (allocated(identCol)) then
195			deallocate(identRow)
196			endif
197	gezelter	317
198			#else
199			do i = 1, nComponents
200	chuckv	306
201	gezelter	317	me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
202			atid(i) = me
203
204			enddo
205	chuckv	306	#endif
206	gezelter	317
207	chuckv	306	call initForce_Modules(thisStat)
208			if (thisStat /= 0) then
209			status = -1
210			return
211			endif
212	gezelter	317
213			!! Create neighbor lists
214	chuckv	306	call expandList(thisStat)
215			if (thisStat /= 0) then
216			status = -1
217			return
218			endif
219	gezelter	317
220	chuckv	315	call setupGlobals(thisStat)
221			if (thisStat /= 0) then
222			status = -1
223			return
224			endif
225	gezelter	309
226	gezelter	317	ff_initialized = .true.
227	gezelter	309
228	gezelter	317
229	chuckv	306	end subroutine init_FF
230	gezelter	317
231	chuckv	315	subroutine setupGlobals(thisStat)
232			integer, intent(out) :: thisStat
233	gezelter	317
234	chuckv	315	end subroutine setupGlobals
235	gezelter	317
236
237	chuckv	306	subroutine initForce_Modules(thisStat)
238			integer, intent(out) :: thisStat
239			integer :: my_status
240
241			thisStat = 0
242	gezelter	317	call init_lj_FF(my_status)
243	chuckv	306	if (my_status /= 0) then
244			thisStat = -1
245			return
246			end if
247	gezelter	317
248	chuckv	306	end subroutine initForce_Modules
249	gezelter	317
250	chuckv	306	end module forceGlobals