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 306 by chuckv, Mon Mar 10 19:26:45 2003 UTC

# Line 1 | Line 1
1 + !! do_Forces.F90
2 + !! module do_Forces
3   !! Calculates Long Range forces.
4 +
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @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 $
7 > !! @version $Id: do_Forces.F90,v 1.8 2003-03-10 19:26:45 chuckv Exp $, $Date: 2003-03-10 19:26:45 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $
8  
9  
10  
11   module do_Forces
12    use simulation
13    use definitions
14 <  use generic_atypes
14 >  use forceGlobals
15 >  use atype_typedefs
16    use neighborLists
17 +
18    
19    use lj_FF
20    use sticky_FF
# Line 22 | Line 27 | module do_Forces
27    implicit none
28    PRIVATE
29  
30 < !! Number of lj_atypes in lj_atype_list
26 <  integer, save :: n_atypes = 0
30 > public :: do_force_loop
31  
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
32   contains
33  
34 < !! 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 <
304 < !! FORCE routine Calculates Lennard Jones forces.
34 > !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
35   !------------------------------------------------------------->
36    subroutine do_force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
37   !! Position array provided by C, dimensioned by getNlocal
# Line 356 | Line 86 | contains
86    real( kind = DP ) ::  rx_ij, ry_ij, rz_ij, rijsq
87    real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut
88  
89 <  real( kind = DP ) :: dielectric = 0.0_dp
89 >  
90  
91   ! a rig that need to be fixed.
92   #ifdef IS_MPI
# Line 735 | Line 465 | contains
465      real( kind = dp ) ::  drdz = 0.0_dp
466      
467  
468 + #ifdef IS_MPI
469 +
470      if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
471         call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
472      endif
473  
474      if (Atype_i%is_dp .and. Atype_j%is_dp) then
475  
744 #ifdef IS_MPI
476         call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
477              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
478  
479         if (do_reaction_field) then
753 #ifdef IS_MPI
480            call accumulate_rf(i, j, r_ij, rflRow(:,i), rflCol(:j), &
481                 ulRow(:i), ulCol(:,j), rt, rrf)
482 +       endif
483 +
484 +    endif
485 +
486 +    if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
487 +       call getstickyforce(r, pot, dudr, Atype_i, Atype_j)
488 +    endif
489 +
490   #else
491 +
492 +    if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
493 +       call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
494 +    endif
495 +
496 +    if (Atype_i%is_dp .and. Atype_j%is_dp) then
497 +       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
498 +            ul(:,i), ul(:,j), rt, rrf, pot)
499 +
500 +       if (do_reaction_field) then
501            call accumulate_rf(i, j, r_ij, rfl(:,i), rfl(:j), &
502                 ul(:,i), ul(:,j), rt, rrf)
759 #endif
503         endif
504  
762
505      endif
506  
507      if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
508 <       call getstickyforce(r,pot,dudr,ljAtype_i,ljAtype_j)
508 >       call getstickyforce(r,pot,dudr, Atype_i, Atype_j)
509      endif
510  
511 + #endif
512 +
513        
514   #ifdef IS_MPI
515                  eRow(i) = eRow(i) + pot*0.5
# Line 781 | Line 525 | contains
525                  fx = dudr * drdx
526                  fy = dudr * drdy
527                  fz = dudr * drdz
784
785
786
787
788
789
528                  
529   #ifdef IS_MPI
530                  fCol(1,j) = fCol(1,j) - fx
# Line 920 | Line 658 | contains
658  
659    end subroutine zero_module_variables
660  
661 < #ifdef IS_MPI
661 >
662   !! Function to properly build neighbor lists in MPI using newtons 3rd law.
663   !! We don't want 2 processors doing the same i j pair twice.
664   !! Also checks to see if i and j are the same particle.
665 <  function mpi_cycle_jLoop(i,j) result(do_cycle)
665 >  function checkExcludes(atom1,atom2) result(do_cycle)
666   !--------------- Arguments--------------------------
667   ! Index i
668 <    integer,intent(in) :: i
668 >    integer,intent(in) :: atom1
669   ! Index j
670 <    integer,intent(in) :: j
670 >    integer,intent(in), optional :: atom2
671   ! Result do_cycle
672      logical :: do_cycle
673   !--------------- Local variables--------------------
674      integer :: tag_i
675      integer :: tag_j
676 < !--------------- END DECLARATIONS------------------    
677 <    tag_i = tagRow(i)
676 >    integer :: i
677 > !--------------- END DECLARATIONS------------------  
678 >    do_cycle = .false.
679 >
680 > #ifdef IS_MPI
681 >    tag_i = tagRow(atom1)
682 > #else
683 >    tag_i = tag(atom1)
684 > #endif
685 >
686 > !! Check global excludes first
687 >    if (.not. present(atom2)) then
688 >       do i = 1,nGlobalExcludes
689 >          if (excludeGlobal(i) == tag_i) then
690 >             do_cycle = .true.
691 >             return
692 >          end if
693 >       end do
694 >       return !! return after checking globals
695 >    end if
696 >
697 > !! we return if j not present here.
698      tag_j = tagColumn(j)
699  
700 <    do_cycle = .false.
700 >
701  
702      if (tag_i == tag_j) then
703         do_cycle = .true.
# Line 952 | Line 710 | contains
710      else                
711         if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
712      endif
955  end function mpi_cycle_jLoop
956 #endif
713  
714 +
715 +
716 +    do i = 1, nLocalExcludes
717 +       if (tag_i = excludes(1,i) .and. excludes(2,i) < 0) then
718 +          do_cycle = .true.
719 +          return
720 +       end if
721 +    end do
722 +      
723 +
724 +  end function checkExcludes
725 +
726 +
727   end module do_Forces

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines