ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2270 by gezelter, Tue Aug 9 22:33:37 2005 UTC vs.
Revision 2279 by chrisfen, Tue Aug 30 18:23:50 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.28 2005-08-09 22:33:37 gezelter Exp $, $Date: 2005-08-09 22:33:37 $, $Name: not supported by cvs2svn $, $Revision: 1.28 $
48 > !! @version $Id: doForces.F90,v 1.33 2005-08-30 18:23:29 chrisfen Exp $, $Date: 2005-08-30 18:23:29 $, $Name: not supported by cvs2svn $, $Revision: 1.33 $
49  
50  
51   module doForces
# Line 73 | Line 73 | module doForces
73  
74   #define __FORTRAN90
75   #include "UseTheForce/fSwitchingFunction.h"
76 + #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78  
79 +
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82  
# Line 97 | Line 99 | module doForces
99    logical, save :: SIM_requires_prepair_calc
100    logical, save :: SIM_uses_PBC
101  
102 +  integer, save :: corrMethod
103 +
104    public :: init_FF
105 +  public :: setDefaultCutoffs
106    public :: do_force_loop
107    public :: createInteractionHash
108    public :: createGtypeCutoffMap
109 +  public :: getStickyCut
110 +  public :: getStickyPowerCut
111 +  public :: getGayBerneCut
112 +  public :: getEAMCut
113 +  public :: getShapeCut
114  
115   #ifdef PROFILE
116    public :: getforcetime
# Line 122 | Line 132 | module doForces
132       real(kind=dp) :: rlistsq
133    end type gtypeCutoffs
134    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
135 +
136 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
137 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
138    
139   contains
140  
# Line 146 | Line 159 | contains
159      logical :: j_is_GB
160      logical :: j_is_EAM
161      logical :: j_is_Shape
162 <    
162 >    real(kind=dp) :: myRcut
163 >
164      status = 0  
165  
166      if (.not. associated(atypes)) then
# Line 179 | Line 193 | contains
193         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
194         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
195  
182 –       if (i_is_LJ) then
183 –          thisCut = getDefaultLJCutoff(i)
184 –          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
185 –       endif
186 –
187 –
188 –
196         do j = i, nAtypes
197  
198            iHash = 0
# Line 238 | Line 245 | contains
245      haveInteractionHash = .true.
246    end subroutine createInteractionHash
247  
248 <  subroutine createGtypeCutoffMap(defaultRcut, defaultSkinThickness, stat)
248 >  subroutine createGtypeCutoffMap(stat)
249  
250 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
251 <    integer, intent(out) :: stat
250 >    integer, intent(out), optional :: stat
251 >    logical :: i_is_LJ
252 >    logical :: i_is_Elect
253 >    logical :: i_is_Sticky
254 >    logical :: i_is_StickyP
255 >    logical :: i_is_GB
256 >    logical :: i_is_EAM
257 >    logical :: i_is_Shape
258  
259 <    integer :: myStatus, nAtypes
259 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
260 >    integer :: n_in_i
261 >    real(kind=dp):: thisSigma, bigSigma, thisRcut
262 >    real(kind=dp) :: biggestAtypeCutoff
263  
264      stat = 0
265      if (.not. haveInteractionHash) then
# Line 256 | Line 272 | contains
272      endif
273  
274      nAtypes = getSize(atypes)
275 <
275 >    
276      do i = 1, nAtypes
277 <      
278 <       atypeMaxCutoff(i) =
277 >       if (SimHasAtype(i)) then          
278 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
279 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
280 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
281 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
282 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
283 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
284 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
285 >          
286 >          if (i_is_LJ) then
287 >             thisRcut = getSigma(i) * 2.5_dp
288 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
289 >          endif
290 >          if (i_is_Elect) then
291 >             thisRcut = defaultRcut
292 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
293 >          endif
294 >          if (i_is_Sticky) then
295 >             thisRcut = getStickyCut(i)
296 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
297 >          endif
298 >          if (i_is_StickyP) then
299 >             thisRcut = getStickyPowerCut(i)
300 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
301 >          endif
302 >          if (i_is_GB) then
303 >             thisRcut = getGayBerneCut(i)
304 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
305 >          endif
306 >          if (i_is_EAM) then
307 >             thisRcut = getEAMCut(i)
308 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
309 >          endif
310 >          if (i_is_Shape) then
311 >             thisRcut = getShapeCut(i)
312 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313 >          endif
314 >          
315 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
316 >             biggestAtypeCutoff = atypeMaxCutoff(i)
317 >          endif
318 >       endif
319 >    enddo
320  
321 +    istart = 1
322 + #ifdef IS_MPI
323 +    iend = nGroupsInRow
324 + #else
325 +    iend = nGroups
326 + #endif
327 +    outer: do i = istart, iend
328 +      
329 +       n_in_i = groupStartRow(i+1) - groupStartRow(i)
330 +      
331 + #ifdef IS_MPI
332 +       jstart = 1
333 +       jend = nGroupsInCol
334 + #else
335 +       jstart = i+1
336 +       jend = nGroups
337 + #endif
338 +      
339 +      
340 +      
341 +      
342 +      
343 +      
344 +    enddo outer        
345      
265 –
266 –
267 –
346       haveGtypeCutoffMap = .true.
347     end subroutine createGtypeCutoffMap
348 +
349 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
350 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
351 +     integer, intent(in) :: cutPolicy
352  
353 +     defaultRcut = defRcut
354 +     defaultRsw = defRsw
355 +     defaultRlist = defRlist
356 +     cutoffPolicy = cutPolicy
357 +   end subroutine setDefaultCutoffs
358 +
359 +   subroutine setCutoffPolicy(cutPolicy)
360 +
361 +     integer, intent(in) :: cutPolicy
362 +     cutoffPolicy = cutPolicy
363 +     call createGtypeCutoffMap()
364 +
365 +   end subroutine setCutoffPolicy
366 +    
367 +    
368    subroutine setSimVariables()
369      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
370      SIM_uses_EAM = SimUsesEAM()
# Line 341 | Line 438 | contains
438    end subroutine doReadyCheck
439  
440  
441 <  subroutine init_FF(use_RF_c, thisStat)
441 >  subroutine init_FF(use_RF_c, use_UW_c, use_DW_c, thisStat)
442  
443      logical, intent(in) :: use_RF_c
444 <
444 >    logical, intent(in) :: use_UW_c
445 >    logical, intent(in) :: use_DW_c
446      integer, intent(out) :: thisStat  
447      integer :: my_status, nMatches
448 +    integer :: corrMethod
449      integer, pointer :: MatchList(:) => null()
450      real(kind=dp) :: rcut, rrf, rt, dielect
451  
# Line 356 | Line 455 | contains
455      !! Fortran's version of a cast:
456      FF_uses_RF = use_RF_c
457  
458 +    !! set the electrostatic correction method
459 +    if (use_UW_c .eq. .true.) then
460 +       corrMethod = 1
461 +    elseif (use_DW_c .eq. .true.) then
462 +       corrMethod = 2
463 +    else
464 +       corrMethod = 0
465 +    endif
466 +    
467      !! init_FF is called *after* all of the atom types have been
468      !! defined in atype_module using the new_atype subroutine.
469      !!
# Line 574 | Line 682 | contains
682         iend = nGroups - 1
683   #endif
684         outer: do i = istart, iend
577 –
578 – #ifdef IS_MPI
579 –             me_i = atid_row(i)
580 – #else
581 –             me_i = atid(i)
582 – #endif
685  
686            if (update_nlist) point(i) = nlist + 1
687  
# Line 926 | Line 1028 | contains
1028  
1029      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1030         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1031 <            pot, eFrame, f, t, do_pot)
1031 >            pot, eFrame, f, t, do_pot, corrMethod)
1032  
1033         if (FF_uses_RF .and. SIM_uses_RF) then
1034  

Diff Legend

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