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 2268 by gezelter, Fri Jul 29 19:38:27 2005 UTC vs.
Revision 2301 by gezelter, Thu Sep 15 22:05:21 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.26 2005-07-29 19:38:27 gezelter Exp $, $Date: 2005-07-29 19:38:27 $, $Name: not supported by cvs2svn $, $Revision: 1.26 $
48 > !! @version $Id: doForces.F90,v 1.44 2005-09-15 22:05:17 gezelter Exp $, $Date: 2005-09-15 22:05:17 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $
49  
50  
51   module doForces
# Line 58 | Line 58 | module doForces
58    use lj
59    use sticky
60    use electrostatic_module
61 <  use reaction_field
61 >  use reaction_field_module
62    use gb_pair
63    use shapes
64    use vector_class
# 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  
81  logical, save :: haveRlist = .false.
83    logical, save :: haveNeighborList = .false.
84    logical, save :: haveSIMvariables = .false.
85    logical, save :: haveSaneForceField = .false.
86 <  logical, save :: haveInteractionMap = .false.
86 >  logical, save :: haveInteractionHash = .false.
87 >  logical, save :: haveGtypeCutoffMap = .false.
88 >  logical, save :: haveRlist = .false.
89  
90    logical, save :: FF_uses_DirectionalAtoms
88  logical, save :: FF_uses_LennardJones
89  logical, save :: FF_uses_Electrostatics
90  logical, save :: FF_uses_Charges
91    logical, save :: FF_uses_Dipoles
92  logical, save :: FF_uses_Quadrupoles
93  logical, save :: FF_uses_Sticky
94  logical, save :: FF_uses_StickyPower
92    logical, save :: FF_uses_GayBerne
93    logical, save :: FF_uses_EAM
97  logical, save :: FF_uses_Shapes
98  logical, save :: FF_uses_FLARB
94    logical, save :: FF_uses_RF
95  
96    logical, save :: SIM_uses_DirectionalAtoms
102  logical, save :: SIM_uses_LennardJones
103  logical, save :: SIM_uses_Electrostatics
104  logical, save :: SIM_uses_Charges
105  logical, save :: SIM_uses_Dipoles
106  logical, save :: SIM_uses_Quadrupoles
107  logical, save :: SIM_uses_Sticky
108  logical, save :: SIM_uses_StickyPower
109  logical, save :: SIM_uses_GayBerne
97    logical, save :: SIM_uses_EAM
111  logical, save :: SIM_uses_Shapes
112  logical, save :: SIM_uses_FLARB
98    logical, save :: SIM_uses_RF
99    logical, save :: SIM_requires_postpair_calc
100    logical, save :: SIM_requires_prepair_calc
101    logical, save :: SIM_uses_PBC
117  logical, save :: SIM_uses_molecular_cutoffs
102  
103 <  !!!GO AWAY---------
120 <  !!!!!real(kind=dp), save :: rlist, rlistsq
103 >  integer, save :: corrMethod
104  
105    public :: init_FF
106 +  public :: setDefaultCutoffs
107    public :: do_force_loop
108 < !  public :: setRlistDF
109 <  !public :: addInteraction
110 <  !public :: setInteractionHash
111 <  !public :: getInteractionHash
112 <  public :: createInteractionMap
113 <  public :: createRcuts
108 >  public :: createInteractionHash
109 >  public :: createGtypeCutoffMap
110 >  public :: getStickyCut
111 >  public :: getStickyPowerCut
112 >  public :: getGayBerneCut
113 >  public :: getEAMCut
114 >  public :: getShapeCut
115  
116   #ifdef PROFILE
117    public :: getforcetime
# Line 134 | Line 119 | module doForces
119    real :: forceTimeInitial, forceTimeFinal
120    integer :: nLoops
121   #endif
137
138  type, public :: Interaction
139     integer :: InteractionHash
140     real(kind=dp) :: rCut = 0.0_dp
141     real(kind=dp) :: rCutSq = 0.0_dp    
142     real(kind=dp) :: rListSq = 0.0_dp
143  end type Interaction
122    
123 <  type(Interaction), dimension(:,:),allocatable :: InteractionMap
124 <  
123 >  !! Variables for cutoff mapping and interaction mapping
124 >  ! Bit hash to determine pair-pair interactions.
125 >  integer, dimension(:,:), allocatable :: InteractionHash
126 >  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
127 >  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
128 >  integer, dimension(:), allocatable :: groupToGtype
129 >  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
130 >  type ::gtypeCutoffs
131 >     real(kind=dp) :: rcut
132 >     real(kind=dp) :: rcutsq
133 >     real(kind=dp) :: rlistsq
134 >  end type gtypeCutoffs
135 >  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
136  
137 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
138 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
139 +  real(kind=dp),save :: rcuti
140    
141   contains
142  
143 <
152 <  subroutine createInteractionMap(status)
143 >  subroutine createInteractionHash(status)
144      integer :: nAtypes
145      integer, intent(out) :: status
146      integer :: i
147      integer :: j
148 <    integer :: ihash
158 <    real(kind=dp) :: myRcut
148 >    integer :: iHash
149      !! Test Types
150      logical :: i_is_LJ
151      logical :: i_is_Elect
# Line 171 | Line 161 | contains
161      logical :: j_is_GB
162      logical :: j_is_EAM
163      logical :: j_is_Shape
164 <    
164 >    real(kind=dp) :: myRcut
165 >
166      status = 0  
167  
168      if (.not. associated(atypes)) then
169 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
169 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
170         status = -1
171         return
172      endif
# Line 187 | Line 178 | contains
178         return
179      end if
180  
181 <    if (.not. allocated(InteractionMap)) then
182 <       allocate(InteractionMap(nAtypes,nAtypes))
181 >    if (.not. allocated(InteractionHash)) then
182 >       allocate(InteractionHash(nAtypes,nAtypes))
183      endif
184 +
185 +    if (.not. allocated(atypeMaxCutoff)) then
186 +       allocate(atypeMaxCutoff(nAtypes))
187 +    endif
188          
189      do i = 1, nAtypes
190         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 242 | Line 237 | contains
237            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
238  
239  
240 <          InteractionMap(i,j)%InteractionHash = iHash
241 <          InteractionMap(j,i)%InteractionHash = iHash
240 >          InteractionHash(i,j) = iHash
241 >          InteractionHash(j,i) = iHash
242  
243         end do
244  
245      end do
246  
247 <    haveInteractionMap = .true.
248 <  end subroutine createInteractionMap
247 >    haveInteractionHash = .true.
248 >  end subroutine createInteractionHash
249  
250 <  ! Query each potential and return the cutoff for that potential. We
256 <  ! build the neighbor list based on the largest cutoff value for that
257 <  ! atype. Each potential can decide whether to calculate the force for
258 <  ! that atype based upon it's own cutoff.
259 <  
260 <  subroutine createRcuts(defaultRcut, defaultSkinThickness, stat)
250 >  subroutine createGtypeCutoffMap(stat)
251  
252 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
253 <    integer :: iMap
254 <    integer :: map_i,map_j
255 <    real(kind=dp) :: thisRCut = 0.0_dp
256 <    real(kind=dp) :: actualCutoff = 0.0_dp
257 <    integer, intent(out) :: stat
258 <    integer :: nAtypes
259 <    integer :: myStatus
252 >    integer, intent(out), optional :: stat
253 >    logical :: i_is_LJ
254 >    logical :: i_is_Elect
255 >    logical :: i_is_Sticky
256 >    logical :: i_is_StickyP
257 >    logical :: i_is_GB
258 >    logical :: i_is_EAM
259 >    logical :: i_is_Shape
260 >    logical :: GtypeFound
261  
262 <    stat = 0
263 <    if (.not. haveInteractionMap) then
262 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
263 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
264 >    integer :: nGroupsInRow
265 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
266 >    real(kind=dp) :: biggestAtypeCutoff
267  
268 <       call createInteractionMap(myStatus)
269 <
268 >    stat = 0
269 >    if (.not. haveInteractionHash) then
270 >       call createInteractionHash(myStatus)      
271         if (myStatus .ne. 0) then
272 <          write(default_error, *) 'createInteractionMap failed in doForces!'
272 >          write(default_error, *) 'createInteractionHash failed in doForces!'
273            stat = -1
274            return
275         endif
276      endif
277 <
278 <    nAtypes = getSize(atypes)
279 <    !! If we pass a default rcut, set all atypes to that cutoff distance
280 <    if(present(defaultRList)) then
281 <       InteractionMap(:,:)%rCut = defaultRCut
282 <       InteractionMap(:,:)%rCutSq = defaultRCut*defaultRCut
283 <       InteractionMap(:,:)%rListSq = (defaultRCut+defaultSkinThickness)**2
284 <       haveRlist = .true.
285 <       return
286 <    end if
287 <
288 <    do map_i = 1,nAtypes
289 <       do map_j = map_i,nAtypes
290 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
277 > #ifdef IS_MPI
278 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
279 > #endif
280 >    nAtypes = getSize(atypes)
281 > ! Set all of the initial cutoffs to zero.
282 >    atypeMaxCutoff = 0.0_dp
283 >    do i = 1, nAtypes
284 >       if (SimHasAtype(i)) then    
285 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
286 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
287 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
288 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
289 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
290 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
291 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
292            
293 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
294 <             ! thisRCut = getLJCutOff(map_i,map_j)
295 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
293 >
294 >          if (i_is_LJ) then
295 >             thisRcut = getSigma(i) * 2.5_dp
296 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
297            endif
298 <          
299 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
300 <             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
304 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
298 >          if (i_is_Elect) then
299 >             thisRcut = defaultRcut
300 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
301            endif
302 +          if (i_is_Sticky) then
303 +             thisRcut = getStickyCut(i)
304 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
305 +          endif
306 +          if (i_is_StickyP) then
307 +             thisRcut = getStickyPowerCut(i)
308 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
309 +          endif
310 +          if (i_is_GB) then
311 +             thisRcut = getGayBerneCut(i)
312 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313 +          endif
314 +          if (i_is_EAM) then
315 +             thisRcut = getEAMCut(i)
316 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317 +          endif
318 +          if (i_is_Shape) then
319 +             thisRcut = getShapeCut(i)
320 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321 +          endif
322            
323 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
324 <             ! thisRCut = getStickyCutOff(map_i,map_j)
325 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
326 <           endif
327 <          
328 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
329 <              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
330 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
331 <           endif
332 <          
333 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
334 <              ! thisRCut = getGayberneCutOff(map_i,map_j)
335 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
336 <           endif
337 <          
338 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
339 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
340 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
341 <           endif
342 <          
343 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
344 < !              thisRCut = getEAMCutOff(map_i,map_j)
345 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
346 <           endif
347 <          
348 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
349 < !              thisRCut = getShapeCutOff(map_i,map_j)
350 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
351 <           endif
352 <          
353 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
354 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
355 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
356 <           endif
357 <           InteractionMap(map_i, map_j)%rCut = actualCutoff
358 <           InteractionMap(map_i, map_j)%rCutSq = actualCutoff * actualCutoff
359 <           InteractionMap(map_i, map_j)%rListSq = (actualCutoff + skinThickness)**2
323 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
324 >             biggestAtypeCutoff = atypeMaxCutoff(i)
325 >          endif
326 >       endif
327 >    enddo
328 >  
329 >    nGroupTypes = 0
330 >    
331 >    istart = 1
332 > #ifdef IS_MPI
333 >    iend = nGroupsInRow
334 > #else
335 >    iend = nGroups
336 > #endif
337 >    
338 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
339 >    if(.not.allocated(groupToGtype)) then
340 >       allocate(groupToGtype(iend))
341 >       allocate(groupMaxCutoff(iend))
342 >       allocate(gtypeMaxCutoff(iend))
343 >       groupMaxCutoff = 0.0_dp
344 >       gtypeMaxCutoff = 0.0_dp
345 >    endif
346 >    !! first we do a single loop over the cutoff groups to find the
347 >    !! largest cutoff for any atypes present in this group.  We also
348 >    !! create gtypes at this point.
349 >    
350 >    tol = 1.0d-6
351 >    
352 >    do i = istart, iend      
353 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
354 >       groupMaxCutoff(i) = 0.0_dp
355 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
356 >          atom1 = groupListRow(ia)
357 > #ifdef IS_MPI
358 >          me_i = atid_row(atom1)
359 > #else
360 >          me_i = atid(atom1)
361 > #endif          
362 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
363 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
364 >          endif          
365 >       enddo
366  
367 <           InteractionMap(map_j, map_i)%rCut = InteractionMap(map_i, map_j)%rCut
368 <           InteractionMap(map_j, map_i)%rCutSq = InteractionMap(map_i, map_j)%rCutSq
369 <           InteractionMap(map_j, map_i)%rListSq = InteractionMap(map_i, map_j)%rListSq
370 <        end do
371 <     end do
372 <     ! now the groups
367 >       if (nGroupTypes.eq.0) then
368 >          nGroupTypes = nGroupTypes + 1
369 >          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
370 >          groupToGtype(i) = nGroupTypes
371 >       else
372 >          GtypeFound = .false.
373 >          do g = 1, nGroupTypes
374 >             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
375 >                groupToGtype(i) = g
376 >                GtypeFound = .true.
377 >             endif
378 >          enddo
379 >          if (.not.GtypeFound) then            
380 >             nGroupTypes = nGroupTypes + 1
381 >             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
382 >             groupToGtype(i) = nGroupTypes
383 >          endif
384 >       endif
385 >    enddo    
386  
387 +    !! allocate the gtypeCutoffMap here.
388 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
389 +    !! then we do a double loop over all the group TYPES to find the cutoff
390 +    !! map between groups of two types
391 +    
392 +    do i = 1, nGroupTypes
393 +       do j = 1, nGroupTypes
394 +      
395 +          select case(cutoffPolicy)
396 +          case(TRADITIONAL_CUTOFF_POLICY)
397 +             thisRcut = maxval(gtypeMaxCutoff)
398 +          case(MIX_CUTOFF_POLICY)
399 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
400 +          case(MAX_CUTOFF_POLICY)
401 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
402 +          case default
403 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
404 +             return
405 +          end select
406 +          gtypeCutoffMap(i,j)%rcut = thisRcut
407 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
408 +          skin = defaultRlist - defaultRcut
409 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
410  
411 +       enddo
412 +    enddo
413 +    
414 +    haveGtypeCutoffMap = .true.
415 +   end subroutine createGtypeCutoffMap
416  
417 <     haveRlist = .true.
418 <  end subroutine createRcuts
417 >   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
418 >     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
419 >     integer, intent(in) :: cutPolicy
420  
421 +     defaultRcut = defRcut
422 +     defaultRsw = defRsw
423 +     defaultRlist = defRlist
424 +     cutoffPolicy = cutPolicy
425 +     rcuti = 1.0_dp / defaultRcut
426 +   end subroutine setDefaultCutoffs
427  
428 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
359 < !!$  subroutine setRlistDF( this_rlist )
360 < !!$
361 < !!$   real(kind=dp) :: this_rlist
362 < !!$
363 < !!$    rlist = this_rlist
364 < !!$    rlistsq = rlist * rlist
365 < !!$
366 < !!$    haveRlist = .true.
367 < !!$
368 < !!$  end subroutine setRlistDF
428 >   subroutine setCutoffPolicy(cutPolicy)
429  
430 <
430 >     integer, intent(in) :: cutPolicy
431 >     cutoffPolicy = cutPolicy
432 >     call createGtypeCutoffMap()
433 >   end subroutine setCutoffPolicy
434 >    
435 >    
436    subroutine setSimVariables()
437      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
373    SIM_uses_LennardJones = SimUsesLennardJones()
374    SIM_uses_Electrostatics = SimUsesElectrostatics()
375    SIM_uses_Charges = SimUsesCharges()
376    SIM_uses_Dipoles = SimUsesDipoles()
377    SIM_uses_Sticky = SimUsesSticky()
378    SIM_uses_StickyPower = SimUsesStickyPower()
379    SIM_uses_GayBerne = SimUsesGayBerne()
438      SIM_uses_EAM = SimUsesEAM()
381    SIM_uses_Shapes = SimUsesShapes()
382    SIM_uses_FLARB = SimUsesFLARB()
383    SIM_uses_RF = SimUsesRF()
439      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
440      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
441      SIM_uses_PBC = SimUsesPBC()
442 +    SIM_uses_RF = SimUsesRF()
443  
444      haveSIMvariables = .true.
445  
# Line 397 | Line 453 | contains
453  
454      error = 0
455  
456 <    if (.not. haveInteractionMap) then
401 <      
456 >    if (.not. haveInteractionHash) then      
457         myStatus = 0      
458 <       call createInteractionMap(myStatus)
404 <      
458 >       call createInteractionHash(myStatus)      
459         if (myStatus .ne. 0) then
460 <          write(default_error, *) 'createInteractionMap failed in doForces!'
460 >          write(default_error, *) 'createInteractionHash failed in doForces!'
461            error = -1
462            return
463         endif
464      endif
465  
466 +    if (.not. haveGtypeCutoffMap) then        
467 +       myStatus = 0      
468 +       call createGtypeCutoffMap(myStatus)      
469 +       if (myStatus .ne. 0) then
470 +          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
471 +          error = -1
472 +          return
473 +       endif
474 +    endif
475 +
476      if (.not. haveSIMvariables) then
477         call setSimVariables()
478      endif
479  
480 <    if (.not. haveRlist) then
481 <       write(default_error, *) 'rList has not been set in doForces!'
482 <       error = -1
483 <       return
484 <    endif
480 >  !  if (.not. haveRlist) then
481 >  !     write(default_error, *) 'rList has not been set in doForces!'
482 >  !     error = -1
483 >  !     return
484 >  !  endif
485  
486      if (.not. haveNeighborList) then
487         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 442 | Line 506 | contains
506    end subroutine doReadyCheck
507  
508  
509 <  subroutine init_FF(use_RF_c, thisStat)
509 >  subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat)
510  
511 <    logical, intent(in) :: use_RF_c
512 <
511 >    logical, intent(in) :: use_RF
512 >    integer, intent(in) :: correctionMethod
513 >    real(kind=dp), intent(in) :: dampingAlpha
514      integer, intent(out) :: thisStat  
515      integer :: my_status, nMatches
516      integer, pointer :: MatchList(:) => null()
# Line 455 | Line 520 | contains
520      thisStat = 0
521  
522      !! Fortran's version of a cast:
523 <    FF_uses_RF = use_RF_c
523 >    FF_uses_RF = use_RF
524  
525 +        
526      !! init_FF is called *after* all of the atom types have been
527      !! defined in atype_module using the new_atype subroutine.
528      !!
# Line 464 | Line 530 | contains
530      !! interactions are used by the force field.    
531  
532      FF_uses_DirectionalAtoms = .false.
467    FF_uses_LennardJones = .false.
468    FF_uses_Electrostatics = .false.
469    FF_uses_Charges = .false.    
533      FF_uses_Dipoles = .false.
471    FF_uses_Sticky = .false.
472    FF_uses_StickyPower = .false.
534      FF_uses_GayBerne = .false.
535      FF_uses_EAM = .false.
475    FF_uses_Shapes = .false.
476    FF_uses_FLARB = .false.
536  
537      call getMatchingElementList(atypes, "is_Directional", .true., &
538           nMatches, MatchList)
539      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
540  
482    call getMatchingElementList(atypes, "is_LennardJones", .true., &
483         nMatches, MatchList)
484    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
485
486    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
487         nMatches, MatchList)
488    if (nMatches .gt. 0) then
489       FF_uses_Electrostatics = .true.
490    endif
491
492    call getMatchingElementList(atypes, "is_Charge", .true., &
493         nMatches, MatchList)
494    if (nMatches .gt. 0) then
495       FF_uses_Charges = .true.  
496       FF_uses_Electrostatics = .true.
497    endif
498
541      call getMatchingElementList(atypes, "is_Dipole", .true., &
542           nMatches, MatchList)
543 <    if (nMatches .gt. 0) then
502 <       FF_uses_Dipoles = .true.
503 <       FF_uses_Electrostatics = .true.
504 <       FF_uses_DirectionalAtoms = .true.
505 <    endif
506 <
507 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
508 <         nMatches, MatchList)
509 <    if (nMatches .gt. 0) then
510 <       FF_uses_Quadrupoles = .true.
511 <       FF_uses_Electrostatics = .true.
512 <       FF_uses_DirectionalAtoms = .true.
513 <    endif
514 <
515 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
516 <         MatchList)
517 <    if (nMatches .gt. 0) then
518 <       FF_uses_Sticky = .true.
519 <       FF_uses_DirectionalAtoms = .true.
520 <    endif
521 <
522 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
523 <         MatchList)
524 <    if (nMatches .gt. 0) then
525 <       FF_uses_StickyPower = .true.
526 <       FF_uses_DirectionalAtoms = .true.
527 <    endif
543 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
544      
545      call getMatchingElementList(atypes, "is_GayBerne", .true., &
546           nMatches, MatchList)
547 <    if (nMatches .gt. 0) then
532 <       FF_uses_GayBerne = .true.
533 <       FF_uses_DirectionalAtoms = .true.
534 <    endif
547 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
548  
549      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
550      if (nMatches .gt. 0) FF_uses_EAM = .true.
551  
539    call getMatchingElementList(atypes, "is_Shape", .true., &
540         nMatches, MatchList)
541    if (nMatches .gt. 0) then
542       FF_uses_Shapes = .true.
543       FF_uses_DirectionalAtoms = .true.
544    endif
552  
546    call getMatchingElementList(atypes, "is_FLARB", .true., &
547         nMatches, MatchList)
548    if (nMatches .gt. 0) FF_uses_FLARB = .true.
549
550    !! Assume sanity (for the sake of argument)
553      haveSaneForceField = .true.
554  
555      !! check to make sure the FF_uses_RF setting makes sense
556  
557 <    if (FF_uses_dipoles) then
557 >    if (FF_uses_Dipoles) then
558         if (FF_uses_RF) then
559            dielect = getDielect()
560            call initialize_rf(dielect)
561         endif
562      else
563 <       if (FF_uses_RF) then          
563 >       if ((corrMethod == 3) .or. FF_uses_RF) then
564            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
565            thisStat = -1
566            haveSaneForceField = .false.
# Line 566 | Line 568 | contains
568         endif
569      endif
570  
569    !sticky module does not contain check_sticky_FF anymore
570    !if (FF_uses_sticky) then
571    !   call check_sticky_FF(my_status)
572    !   if (my_status /= 0) then
573    !      thisStat = -1
574    !      haveSaneForceField = .false.
575    !      return
576    !   end if
577    !endif
578
571      if (FF_uses_EAM) then
572         call init_EAM_FF(my_status)
573         if (my_status /= 0) then
# Line 593 | Line 585 | contains
585            haveSaneForceField = .false.
586            return
587         endif
596    endif
597
598    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
588      endif
589  
590      if (.not. haveNeighborList) then
# Line 661 | Line 650 | contains
650      integer :: localError
651      integer :: propPack_i, propPack_j
652      integer :: loopStart, loopEnd, loop
653 <    integer :: iMap
653 >    integer :: iHash
654      real(kind=dp) :: listSkin = 1.0  
655  
656      !! initialize local variables  
# Line 753 | Line 742 | contains
742   #endif
743         outer: do i = istart, iend
744  
756 #ifdef IS_MPI
757             me_i = atid_row(i)
758 #else
759             me_i = atid(i)
760 #endif
761
745            if (update_nlist) point(i) = nlist + 1
746  
747            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 795 | Line 778 | contains
778                    q_group(:,j), d_grp, rgrpsq)
779   #endif
780  
781 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
781 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
782                  if (update_nlist) then
783                     nlist = nlist + 1
784  
# Line 995 | Line 978 | contains
978  
979      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
980  
981 <       if (FF_uses_RF .and. SIM_uses_RF) then
981 >       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
982  
983   #ifdef IS_MPI
984            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 1013 | Line 996 | contains
996   #else
997               me_i = atid(i)
998   #endif
999 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
999 >             iHash = InteractionHash(me_i,me_j)
1000              
1001 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1001 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1002  
1003                  mu_i = getDipoleMoment(me_i)
1004  
# Line 1079 | Line 1062 | contains
1062      real ( kind = dp ), intent(inout) :: rijsq
1063      real ( kind = dp )                :: r
1064      real ( kind = dp ), intent(inout) :: d(3)
1082    real ( kind = dp ) :: ebalance
1065      integer :: me_i, me_j
1066  
1067 <    integer :: iMap
1067 >    integer :: iHash
1068  
1069      r = sqrt(rijsq)
1070      vpair = 0.0d0
# Line 1096 | Line 1078 | contains
1078      me_j = atid(j)
1079   #endif
1080  
1081 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1081 >    iHash = InteractionHash(me_i, me_j)
1082  
1083 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1083 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1084         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1085      endif
1086  
1087 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1087 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1088         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1089 <            pot, eFrame, f, t, do_pot)
1089 >            pot, eFrame, f, t, do_pot, corrMethod, rcuti)
1090  
1091 <       if (FF_uses_RF .and. SIM_uses_RF) then
1091 >       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
1092  
1093            ! CHECK ME (RF needs to know about all electrostatic types)
1094            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1115 | Line 1097 | contains
1097  
1098      endif
1099  
1100 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1100 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1101         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1102              pot, A, f, t, do_pot)
1103      endif
1104  
1105 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1105 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1106         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1107              pot, A, f, t, do_pot)
1108      endif
1109  
1110 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1110 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1111         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1112              pot, A, f, t, do_pot)
1113      endif
1114      
1115 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1115 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1116   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1117   !           pot, A, f, t, do_pot)
1118      endif
1119  
1120 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1120 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1121         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1122              do_pot)
1123      endif
1124  
1125 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1125 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1126         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1127              pot, A, f, t, do_pot)
1128      endif
1129  
1130 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1130 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1131         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1132              pot, A, f, t, do_pot)
1133      endif
# Line 1167 | Line 1149 | contains
1149      real ( kind = dp )                :: r, rc
1150      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1151  
1152 <    integer :: me_i, me_j, iMap
1152 >    integer :: me_i, me_j, iHash
1153  
1154 +    r = sqrt(rijsq)
1155 +
1156   #ifdef IS_MPI  
1157      me_i = atid_row(i)
1158      me_j = atid_col(j)  
# Line 1177 | Line 1161 | contains
1161      me_j = atid(j)  
1162   #endif
1163  
1164 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1164 >    iHash = InteractionHash(me_i, me_j)
1165  
1166 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1166 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1167              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1168      endif
1169      
# Line 1376 | Line 1360 | contains
1360  
1361    function FF_UsesDirectionalAtoms() result(doesit)
1362      logical :: doesit
1363 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1380 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1381 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1363 >    doesit = FF_uses_DirectionalAtoms
1364    end function FF_UsesDirectionalAtoms
1365  
1366    function FF_RequiresPrepairCalc() result(doesit)
# Line 1389 | Line 1371 | contains
1371    function FF_RequiresPostpairCalc() result(doesit)
1372      logical :: doesit
1373      doesit = FF_uses_RF
1374 +    if (corrMethod == 3) doesit = .true.
1375    end function FF_RequiresPostpairCalc
1376  
1377   #ifdef PROFILE

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines