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

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/do_Forces.F90 (file contents):
Revision 297 by gezelter, Thu Mar 6 22:08:29 2003 UTC vs.
Revision 302 by gezelter, Mon Mar 10 14:53:36 2003 UTC

# Line 1 | Line 1
1   !! Calculates Long Range forces.
2   !! @author Charles F. Vardeman II
3   !! @author Matthew Meineke
4 < !! @version $Id: do_Forces.F90,v 1.4 2003-03-06 22:08:29 gezelter Exp $, $Date: 2003-03-06 22:08:29 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $
4 > !! @version $Id: do_Forces.F90,v 1.6 2003-03-10 14:53:36 gezelter Exp $, $Date: 2003-03-10 14:53:36 $, $Name: not supported by cvs2svn $, $Revision: 1.6 $
5  
6  
7  
8   module do_Forces
9    use simulation
10    use definitions
11 <  use generic_atypes
11 >  use forceGlobals
12 >  use atype_typedefs
13    use neighborLists
14 +
15    
16    use lj_FF
17    use sticky_FF
# Line 22 | Line 24 | module do_Forces
24    implicit none
25    PRIVATE
26  
27 < !! Number of lj_atypes in lj_atype_list
26 <  integer, save :: n_atypes = 0
27 > public :: do_force_loop
28  
28 – !! Global list of lj atypes in simulation
29 –  type (atype), pointer :: ListHead => null()
30 –  type (atype), pointer :: ListTail => null()
31 –
32 –
33 –
34 –
35 –  logical, save :: firstTime = .True.
36 –
37 – !! Atype identity pointer lists
38 – #ifdef IS_MPI
39 – !! Row lj_atype pointer list
40 –  type (identPtrList), dimension(:), pointer :: identPtrListRow => null()
41 – !! Column lj_atype pointer list
42 –  type (identPtrList), dimension(:), pointer :: identPtrListColumn => null()
43 – #else
44 –  type(identPtrList ), dimension(:), pointer :: identPtrList => null()
45 – #endif
46 –
47 –
48 – !! Logical has lj force field module been initialized?
49 –  logical, save :: isFFinit = .false.
50 –
51 – !! Use periodic boundry conditions
52 –  logical :: wrap = .false.
53 –
54 – !! Potential energy global module variables
55 – #ifdef IS_MPI
56 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
57 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
58 –
59 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: muRow = 0.0_dp
60 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: muCol = 0.0_dp
61 –
62 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
63 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp
64 –
65 –  real(kind = dp), dimension(9,getNrow(plan_row)) :: ARow = 0.0_dp
66 –  real(kind = dp), dimension(9,getNcol(plan_col)) :: ACol = 0.0_dp  
67 –
68 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
69 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
70 –  real(kind = dp), dimension(3,getNlocal()) :: fTemp1 = 0.0_dp
71 –  real(kind = dp), dimension(3,getNlocal()) :: tTemp1 = 0.0_dp
72 –  real(kind = dp), dimension(3,getNlocal()) :: fTemp2 = 0.0_dp
73 –  real(kind = dp), dimension(3,getNlocal()) :: tTemp2 = 0.0_dp
74 –  real(kind = dp), dimension(3,getNlocal()) :: fTemp = 0.0_dp
75 –  real(kind = dp), dimension(3,getNlocal()) :: tTemp = 0.0_dp
76 –
77 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
78 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp
79 –
80 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: rflRow = 0.0_dp
81 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: rflCol = 0.0_dp
82 –  real(kind = dp), dimension(3,getNlocal()) :: rflTemp = 0.0_dp
83 –
84 –  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
85 –  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
86 –
87 –  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
88 – #endif
89 –  real(kind = dp) :: pe = 0.0_dp
90 –  real(kind = dp), dimension(3,natoms) :: fTemp = 0.0_dp
91 –  real(kind = dp), dimension(3,natoms) :: tTemp = 0.0_dp
92 –  real(kind = dp), dimension(3,natoms) :: rflTemp = 0.0_dp
93 –  real(kind = dp), dimension(9) :: tauTemp = 0.0_dp
94 –
95 –  logical :: do_preForce  = .false.
96 –  logical :: do_postForce = .false.
97 –
98 –
99 –
100 – !! Public methods and data
101 –  public :: new_atype
102 –  public :: do_forceLoop
103 –  public :: init_FF
104 –
105 –  
106 –
107 –
29   contains
30  
110 – !! Adds a new lj_atype to the list.
111 –  subroutine new_atype(ident,mass,epsilon,sigma, &
112 –       is_LJ,is_Sticky,is_DP,is_GB,w0,v0,dipoleMoment,status)
113 –    real( kind = dp ), intent(in) :: mass
114 –    real( kind = dp ), intent(in) :: epsilon
115 –    real( kind = dp ), intent(in) :: sigma
116 –    real( kind = dp ), intent(in) :: w0
117 –    real( kind = dp ), intent(in) :: v0
118 –    real( kind = dp ), intent(in) :: dipoleMoment
119 –
120 –    integer, intent(in) :: ident
121 –    integer, intent(out) :: status
122 –    integer, intent(in) :: is_Sticky
123 –    integer, intent(in) :: is_DP
124 –    integer, intent(in) :: is_GB
125 –    integer, intent(in) :: is_LJ
126 –
127 –
128 –    type (atype), pointer :: the_new_atype
129 –    integer :: alloc_error
130 –    integer :: atype_counter = 0
131 –    integer :: alloc_size
132 –    integer :: err_stat
133 –    status = 0
134 –
135 –
136 –
137 – ! allocate a new atype    
138 –    allocate(the_new_atype,stat=alloc_error)
139 –    if (alloc_error /= 0 ) then
140 –       status = -1
141 –       return
142 –    end if
143 –
144 – ! assign our new atype information
145 –    the_new_atype%mass        = mass
146 –    the_new_atype%epsilon     = epsilon
147 –    the_new_atype%sigma       = sigma
148 –    the_new_atype%sigma2      = sigma * sigma
149 –    the_new_atype%sigma6      = the_new_atype%sigma2 * the_new_atype%sigma2 &
150 –         * the_new_atype%sigma2
151 –    the_new_atype%w0       = w0
152 –    the_new_atype%v0       = v0
153 –    the_new_atype%dipoleMoment       = dipoleMoment
154 –
155 –    
156 – ! assume that this atype will be successfully added
157 –    the_new_atype%atype_ident = ident
158 –    the_new_atype%atype_number = n_lj_atypes + 1
159 –    
160 –    if ( is_Sticky /= 0 )    the_new_atype%is_Sticky   = .true.
161 –    if ( is_GB /= 0 )        the_new_atype%is_GB       = .true.
162 –    if ( is_LJ /= 0 )        the_new_atype%is_LJ       = .true.
163 –    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.
164 –
165 –    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
166 –    if (err_stat /= 0 ) then
167 –       status = -1
168 –       return
169 –    endif
170 –
171 –    n_atypes = n_atypes + 1
172 –
173 –
174 –  end subroutine new_atype
175 –
176 –
177 –  subroutine init_FF(nComponents,ident, status)
178 – !! Number of components in ident array
179 –    integer, intent(inout) :: nComponents
180 – !! Array of identities nComponents long corresponding to
181 – !! ljatype ident.
182 –    integer, dimension(nComponents),intent(inout) :: ident
183 – !!  Result status, success = 0, error = -1
184 –    integer, intent(out) :: Status
185 –
186 –    integer :: alloc_stat
187 –
188 –    integer :: thisStat
189 –    integer :: i
190 –
191 –    integer :: myNode
192 – #ifdef IS_MPI
193 –    integer, allocatable, dimension(:) :: identRow
194 –    integer, allocatable, dimension(:) :: identCol
195 –    integer :: nrow
196 –    integer :: ncol
197 – #endif
198 –    status = 0
199 –  
200 –
201 –    
202 –
203 – !! if were're not in MPI, we just update ljatypePtrList
204 – #ifndef IS_MPI
205 –    call create_IdentPtrlst(ident,ListHead,identPtrList,thisStat)
206 –    if ( thisStat /= 0 ) then
207 –       status = -1
208 –       return
209 –    endif
210 –
211 –    
212 – ! if were're in MPI, we also have to worry about row and col lists    
213 – #else
214 –  
215 – ! We can only set up forces if mpiSimulation has been setup.
216 –    if (.not. isMPISimSet()) then
217 –       write(default_error,*) "MPI is not set"
218 –       status = -1
219 –       return
220 –    endif
221 –    nrow = getNrow(plan_row)
222 –    ncol = getNcol(plan_col)
223 –    mynode = getMyNode()
224 – !! Allocate temperary arrays to hold gather information
225 –    allocate(identRow(nrow),stat=alloc_stat)
226 –    if (alloc_stat /= 0 ) then
227 –       status = -1
228 –       return
229 –    endif
230 –
231 –    allocate(identCol(ncol),stat=alloc_stat)
232 –    if (alloc_stat /= 0 ) then
233 –       status = -1
234 –       return
235 –    endif
236 –
237 – !! Gather idents into row and column idents
238 –
239 –    call gather(ident,identRow,plan_row)
240 –    call gather(ident,identCol,plan_col)
241 –    
242 –  
243 – !! Create row and col pointer lists
244 –  
245 –    call create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
246 –    if (thisStat /= 0 ) then
247 –       status = -1
248 –       return
249 –    endif
250 –  
251 –    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
252 –    if (thisStat /= 0 ) then
253 –       status = -1
254 –       return
255 –    endif
256 –
257 – !! free temporary ident arrays
258 –    if (allocated(identCol)) then
259 –       deallocate(identCol)
260 –    end if
261 –    if (allocated(identCol)) then
262 –       deallocate(identRow)
263 –    endif
264 –
265 – #endif
266 –    
267 –    call initForce_Modules(thisStat)
268 –    if (thisStat /= 0) then
269 –       status = -1
270 –       return
271 –    endif
272 –
273 – !! Create neighbor lists
274 –    call expandList(thisStat)
275 –    if (thisStat /= 0) then
276 –       status = -1
277 –       return
278 –    endif
279 –
280 –    isFFinit = .true.
281 –
282 –
283 –  end subroutine init_FF
284 –
285 –
286 –
287 –
288 –  subroutine initForce_Modules(thisStat)
289 –    integer, intent(out) :: thisStat
290 –    integer :: my_status
291 –    
292 –    thisStat = 0
293 –    call init_lj_FF(ListHead,my_status)
294 –    if (my_status /= 0) then
295 –       thisStat = -1
296 –       return
297 –    end if
298 –
299 –  end subroutine initForce_Modules
300 –
301 –
302 –
303 –
31   !! FORCE routine Calculates Lennard Jones forces.
32   !------------------------------------------------------------->
33    subroutine do_force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
# Line 734 | Line 461 | contains
461      real( kind = dp ) ::  drdy = 0.0_dp
462      real( kind = dp ) ::  drdz = 0.0_dp
463      
464 +
465 + #ifdef IS_MPI
466  
467      if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
468         call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
# Line 741 | Line 470 | contains
470  
471      if (Atype_i%is_dp .and. Atype_j%is_dp) then
472  
744 – #ifdef IS_MPI
473         call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
474              ulRow(:,i), ulCol(:,j), rt, rrf, pot)
747 – #else
748 –       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
749 –            ul(:,i), ul(:,j), rt, rrf, pot)
750 – #endif
475  
476         if (do_reaction_field) then
753 – #ifdef IS_MPI
477            call accumulate_rf(i, j, r_ij, rflRow(:,i), rflCol(:j), &
478                 ulRow(:i), ulCol(:,j), rt, rrf)
479 +       endif
480 +
481 +    endif
482 +
483 +    if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
484 +       call getstickyforce(r,pot,dudr,ljAtype_i,ljAtype_j)
485 +    endif
486 +
487   #else
488 +
489 +    if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
490 +       call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
491 +    endif
492 +
493 +    if (Atype_i%is_dp .and. Atype_j%is_dp) then
494 +       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
495 +            ul(:,i), ul(:,j), rt, rrf, pot)
496 +
497 +       if (do_reaction_field) then
498            call accumulate_rf(i, j, r_ij, rfl(:,i), rfl(:j), &
499                 ul(:,i), ul(:,j), rt, rrf)
759 – #endif
500         endif
501  
762 –
502      endif
503  
504      if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
505         call getstickyforce(r,pot,dudr,ljAtype_i,ljAtype_j)
506      endif
507  
508 + #endif
509 +
510        
511   #ifdef IS_MPI
512                  eRow(i) = eRow(i) + pot*0.5
# Line 781 | Line 522 | contains
522                  fx = dudr * drdx
523                  fy = dudr * drdy
524                  fz = dudr * drdz
784 –
785 –
786 –
787 –
788 –
789 –
525                  
526   #ifdef IS_MPI
527                  fCol(1,j) = fCol(1,j) - fx

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines