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

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2268 by gezelter, Fri Jul 29 19:38:27 2005 UTC vs.
Revision 2306 by chrisfen, Fri Sep 16 20:37:05 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.45 2005-09-16 20:36:55 chrisfen Exp $, $Date: 2005-09-16 20:36:55 $, $Name: not supported by cvs2svn $, $Revision: 1.45 $
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
99  logical, save :: FF_uses_RF
94  
95    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
96    logical, save :: SIM_uses_EAM
111  logical, save :: SIM_uses_Shapes
112  logical, save :: SIM_uses_FLARB
113  logical, save :: SIM_uses_RF
97    logical, save :: SIM_requires_postpair_calc
98    logical, save :: SIM_requires_prepair_calc
99    logical, save :: SIM_uses_PBC
117  logical, save :: SIM_uses_molecular_cutoffs
100  
101 <  !!!GO AWAY---------
120 <  !!!!!real(kind=dp), save :: rlist, rlistsq
101 >  integer, save :: electrostaticSummationMethod
102  
103    public :: init_FF
104 +  public :: setDefaultCutoffs
105    public :: do_force_loop
106 < !  public :: setRlistDF
107 <  !public :: addInteraction
108 <  !public :: setInteractionHash
109 <  !public :: getInteractionHash
110 <  public :: createInteractionMap
111 <  public :: createRcuts
106 >  public :: createInteractionHash
107 >  public :: createGtypeCutoffMap
108 >  public :: getStickyCut
109 >  public :: getStickyPowerCut
110 >  public :: getGayBerneCut
111 >  public :: getEAMCut
112 >  public :: getShapeCut
113  
114   #ifdef PROFILE
115    public :: getforcetime
# Line 134 | Line 117 | module doForces
117    real :: forceTimeInitial, forceTimeFinal
118    integer :: nLoops
119   #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
120    
121 <  type(Interaction), dimension(:,:),allocatable :: InteractionMap
122 <  
121 >  !! Variables for cutoff mapping and interaction mapping
122 >  ! Bit hash to determine pair-pair interactions.
123 >  integer, dimension(:,:), allocatable :: InteractionHash
124 >  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
125 >  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
126 >  integer, dimension(:), allocatable :: groupToGtype
127 >  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
128 >  type ::gtypeCutoffs
129 >     real(kind=dp) :: rcut
130 >     real(kind=dp) :: rcutsq
131 >     real(kind=dp) :: rlistsq
132 >  end type gtypeCutoffs
133 >  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
134  
135 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
136 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
137 +  real(kind=dp),save :: rcuti
138    
139   contains
140  
141 <
152 <  subroutine createInteractionMap(status)
141 >  subroutine createInteractionHash(status)
142      integer :: nAtypes
143      integer, intent(out) :: status
144      integer :: i
145      integer :: j
146 <    integer :: ihash
158 <    real(kind=dp) :: myRcut
146 >    integer :: iHash
147      !! Test Types
148      logical :: i_is_LJ
149      logical :: i_is_Elect
# Line 171 | 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
167 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
167 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
168         status = -1
169         return
170      endif
# Line 187 | Line 176 | contains
176         return
177      end if
178  
179 <    if (.not. allocated(InteractionMap)) then
180 <       allocate(InteractionMap(nAtypes,nAtypes))
179 >    if (.not. allocated(InteractionHash)) then
180 >       allocate(InteractionHash(nAtypes,nAtypes))
181      endif
182 +
183 +    if (.not. allocated(atypeMaxCutoff)) then
184 +       allocate(atypeMaxCutoff(nAtypes))
185 +    endif
186          
187      do i = 1, nAtypes
188         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 242 | Line 235 | contains
235            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
236  
237  
238 <          InteractionMap(i,j)%InteractionHash = iHash
239 <          InteractionMap(j,i)%InteractionHash = iHash
238 >          InteractionHash(i,j) = iHash
239 >          InteractionHash(j,i) = iHash
240  
241         end do
242  
243      end do
244  
245 <    haveInteractionMap = .true.
246 <  end subroutine createInteractionMap
245 >    haveInteractionHash = .true.
246 >  end subroutine createInteractionHash
247  
248 <  ! 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)
248 >  subroutine createGtypeCutoffMap(stat)
249  
250 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
251 <    integer :: iMap
252 <    integer :: map_i,map_j
253 <    real(kind=dp) :: thisRCut = 0.0_dp
254 <    real(kind=dp) :: actualCutoff = 0.0_dp
255 <    integer, intent(out) :: stat
256 <    integer :: nAtypes
257 <    integer :: myStatus
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 >    logical :: GtypeFound
259  
260 <    stat = 0
261 <    if (.not. haveInteractionMap) then
260 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
261 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
262 >    integer :: nGroupsInRow
263 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
264 >    real(kind=dp) :: biggestAtypeCutoff
265  
266 <       call createInteractionMap(myStatus)
267 <
266 >    stat = 0
267 >    if (.not. haveInteractionHash) then
268 >       call createInteractionHash(myStatus)      
269         if (myStatus .ne. 0) then
270 <          write(default_error, *) 'createInteractionMap failed in doForces!'
270 >          write(default_error, *) 'createInteractionHash failed in doForces!'
271            stat = -1
272            return
273         endif
274      endif
275 <
276 <    nAtypes = getSize(atypes)
277 <    !! If we pass a default rcut, set all atypes to that cutoff distance
278 <    if(present(defaultRList)) then
279 <       InteractionMap(:,:)%rCut = defaultRCut
280 <       InteractionMap(:,:)%rCutSq = defaultRCut*defaultRCut
281 <       InteractionMap(:,:)%rListSq = (defaultRCut+defaultSkinThickness)**2
282 <       haveRlist = .true.
283 <       return
284 <    end if
285 <
286 <    do map_i = 1,nAtypes
287 <       do map_j = map_i,nAtypes
288 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
275 > #ifdef IS_MPI
276 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
277 > #endif
278 >    nAtypes = getSize(atypes)
279 > ! Set all of the initial cutoffs to zero.
280 >    atypeMaxCutoff = 0.0_dp
281 >    do i = 1, nAtypes
282 >       if (SimHasAtype(i)) then    
283 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
284 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
285 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
286 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
287 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
288 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
289 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
290            
291 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
292 <             ! thisRCut = getLJCutOff(map_i,map_j)
293 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
291 >
292 >          if (i_is_LJ) then
293 >             thisRcut = getSigma(i) * 2.5_dp
294 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
295            endif
296 <          
297 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
298 <             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
304 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
296 >          if (i_is_Elect) then
297 >             thisRcut = defaultRcut
298 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
299            endif
300 +          if (i_is_Sticky) then
301 +             thisRcut = getStickyCut(i)
302 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
303 +          endif
304 +          if (i_is_StickyP) then
305 +             thisRcut = getStickyPowerCut(i)
306 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
307 +          endif
308 +          if (i_is_GB) then
309 +             thisRcut = getGayBerneCut(i)
310 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
311 +          endif
312 +          if (i_is_EAM) then
313 +             thisRcut = getEAMCut(i)
314 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
315 +          endif
316 +          if (i_is_Shape) then
317 +             thisRcut = getShapeCut(i)
318 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
319 +          endif
320            
321 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
322 <             ! thisRCut = getStickyCutOff(map_i,map_j)
323 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
324 <           endif
325 <          
326 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
327 <              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
328 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
329 <           endif
330 <          
331 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
332 <              ! thisRCut = getGayberneCutOff(map_i,map_j)
333 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
334 <           endif
335 <          
336 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
337 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
338 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
339 <           endif
340 <          
341 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
342 < !              thisRCut = getEAMCutOff(map_i,map_j)
343 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
344 <           endif
345 <          
346 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
347 < !              thisRCut = getShapeCutOff(map_i,map_j)
348 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
349 <           endif
350 <          
351 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
352 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
353 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
354 <           endif
355 <           InteractionMap(map_i, map_j)%rCut = actualCutoff
356 <           InteractionMap(map_i, map_j)%rCutSq = actualCutoff * actualCutoff
357 <           InteractionMap(map_i, map_j)%rListSq = (actualCutoff + skinThickness)**2
321 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
322 >             biggestAtypeCutoff = atypeMaxCutoff(i)
323 >          endif
324 >       endif
325 >    enddo
326 >  
327 >    nGroupTypes = 0
328 >    
329 >    istart = 1
330 > #ifdef IS_MPI
331 >    iend = nGroupsInRow
332 > #else
333 >    iend = nGroups
334 > #endif
335 >    
336 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
337 >    if(.not.allocated(groupToGtype)) then
338 >       allocate(groupToGtype(iend))
339 >       allocate(groupMaxCutoff(iend))
340 >       allocate(gtypeMaxCutoff(iend))
341 >       groupMaxCutoff = 0.0_dp
342 >       gtypeMaxCutoff = 0.0_dp
343 >    endif
344 >    !! first we do a single loop over the cutoff groups to find the
345 >    !! largest cutoff for any atypes present in this group.  We also
346 >    !! create gtypes at this point.
347 >    
348 >    tol = 1.0d-6
349 >    
350 >    do i = istart, iend      
351 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
352 >       groupMaxCutoff(i) = 0.0_dp
353 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
354 >          atom1 = groupListRow(ia)
355 > #ifdef IS_MPI
356 >          me_i = atid_row(atom1)
357 > #else
358 >          me_i = atid(atom1)
359 > #endif          
360 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
361 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
362 >          endif          
363 >       enddo
364  
365 <           InteractionMap(map_j, map_i)%rCut = InteractionMap(map_i, map_j)%rCut
366 <           InteractionMap(map_j, map_i)%rCutSq = InteractionMap(map_i, map_j)%rCutSq
367 <           InteractionMap(map_j, map_i)%rListSq = InteractionMap(map_i, map_j)%rListSq
368 <        end do
369 <     end do
370 <     ! now the groups
365 >       if (nGroupTypes.eq.0) then
366 >          nGroupTypes = nGroupTypes + 1
367 >          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
368 >          groupToGtype(i) = nGroupTypes
369 >       else
370 >          GtypeFound = .false.
371 >          do g = 1, nGroupTypes
372 >             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
373 >                groupToGtype(i) = g
374 >                GtypeFound = .true.
375 >             endif
376 >          enddo
377 >          if (.not.GtypeFound) then            
378 >             nGroupTypes = nGroupTypes + 1
379 >             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
380 >             groupToGtype(i) = nGroupTypes
381 >          endif
382 >       endif
383 >    enddo    
384  
385 +    !! allocate the gtypeCutoffMap here.
386 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
387 +    !! then we do a double loop over all the group TYPES to find the cutoff
388 +    !! map between groups of two types
389 +    
390 +    do i = 1, nGroupTypes
391 +       do j = 1, nGroupTypes
392 +      
393 +          select case(cutoffPolicy)
394 +          case(TRADITIONAL_CUTOFF_POLICY)
395 +             thisRcut = maxval(gtypeMaxCutoff)
396 +          case(MIX_CUTOFF_POLICY)
397 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
398 +          case(MAX_CUTOFF_POLICY)
399 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
400 +          case default
401 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
402 +             return
403 +          end select
404 +          gtypeCutoffMap(i,j)%rcut = thisRcut
405 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
406 +          skin = defaultRlist - defaultRcut
407 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
408  
409 +       enddo
410 +    enddo
411 +    
412 +    haveGtypeCutoffMap = .true.
413 +   end subroutine createGtypeCutoffMap
414  
415 <     haveRlist = .true.
416 <  end subroutine createRcuts
415 >   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
416 >     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
417 >     integer, intent(in) :: cutPolicy
418  
419 +     defaultRcut = defRcut
420 +     defaultRsw = defRsw
421 +     defaultRlist = defRlist
422 +     cutoffPolicy = cutPolicy
423 +     rcuti = 1.0_dp / defaultRcut
424 +   end subroutine setDefaultCutoffs
425  
426 < !!! 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
426 >   subroutine setCutoffPolicy(cutPolicy)
427  
428 <
428 >     integer, intent(in) :: cutPolicy
429 >     cutoffPolicy = cutPolicy
430 >     call createGtypeCutoffMap()
431 >   end subroutine setCutoffPolicy
432 >    
433 >    
434    subroutine setSimVariables()
435      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()
436      SIM_uses_EAM = SimUsesEAM()
381    SIM_uses_Shapes = SimUsesShapes()
382    SIM_uses_FLARB = SimUsesFLARB()
383    SIM_uses_RF = SimUsesRF()
437      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
438      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
439      SIM_uses_PBC = SimUsesPBC()
# Line 397 | Line 450 | contains
450  
451      error = 0
452  
453 <    if (.not. haveInteractionMap) then
401 <      
453 >    if (.not. haveInteractionHash) then      
454         myStatus = 0      
455 <       call createInteractionMap(myStatus)
404 <      
455 >       call createInteractionHash(myStatus)      
456         if (myStatus .ne. 0) then
457 <          write(default_error, *) 'createInteractionMap failed in doForces!'
457 >          write(default_error, *) 'createInteractionHash failed in doForces!'
458            error = -1
459            return
460         endif
461      endif
462  
463 +    if (.not. haveGtypeCutoffMap) then        
464 +       myStatus = 0      
465 +       call createGtypeCutoffMap(myStatus)      
466 +       if (myStatus .ne. 0) then
467 +          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
468 +          error = -1
469 +          return
470 +       endif
471 +    endif
472 +
473      if (.not. haveSIMvariables) then
474         call setSimVariables()
475      endif
476  
477 <    if (.not. haveRlist) then
478 <       write(default_error, *) 'rList has not been set in doForces!'
479 <       error = -1
480 <       return
481 <    endif
477 >  !  if (.not. haveRlist) then
478 >  !     write(default_error, *) 'rList has not been set in doForces!'
479 >  !     error = -1
480 >  !     return
481 >  !  endif
482  
483      if (.not. haveNeighborList) then
484         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 442 | Line 503 | contains
503    end subroutine doReadyCheck
504  
505  
506 <  subroutine init_FF(use_RF_c, thisStat)
506 >  subroutine init_FF(thisESM, dampingAlpha, thisStat)
507  
508 <    logical, intent(in) :: use_RF_c
509 <
508 >    integer, intent(in) :: thisESM
509 >    real(kind=dp), intent(in) :: dampingAlpha
510      integer, intent(out) :: thisStat  
511      integer :: my_status, nMatches
512      integer, pointer :: MatchList(:) => null()
# Line 454 | Line 515 | contains
515      !! assume things are copacetic, unless they aren't
516      thisStat = 0
517  
518 <    !! Fortran's version of a cast:
458 <    FF_uses_RF = use_RF_c
518 >    electrostaticSummationMethod = thisESM
519  
520      !! init_FF is called *after* all of the atom types have been
521      !! defined in atype_module using the new_atype subroutine.
# Line 464 | Line 524 | contains
524      !! interactions are used by the force field.    
525  
526      FF_uses_DirectionalAtoms = .false.
467    FF_uses_LennardJones = .false.
468    FF_uses_Electrostatics = .false.
469    FF_uses_Charges = .false.    
527      FF_uses_Dipoles = .false.
471    FF_uses_Sticky = .false.
472    FF_uses_StickyPower = .false.
528      FF_uses_GayBerne = .false.
529      FF_uses_EAM = .false.
475    FF_uses_Shapes = .false.
476    FF_uses_FLARB = .false.
530  
531      call getMatchingElementList(atypes, "is_Directional", .true., &
532           nMatches, MatchList)
533      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
481
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
534  
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
535      call getMatchingElementList(atypes, "is_Dipole", .true., &
536           nMatches, MatchList)
537 <    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
537 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
538      
539      call getMatchingElementList(atypes, "is_GayBerne", .true., &
540           nMatches, MatchList)
541 <    if (nMatches .gt. 0) then
532 <       FF_uses_GayBerne = .true.
533 <       FF_uses_DirectionalAtoms = .true.
534 <    endif
541 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
542  
543      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
544      if (nMatches .gt. 0) FF_uses_EAM = .true.
545  
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
546  
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)
547      haveSaneForceField = .true.
548  
549 <    !! check to make sure the FF_uses_RF setting makes sense
549 >    !! check to make sure the reaction field setting makes sense
550  
551 <    if (FF_uses_dipoles) then
552 <       if (FF_uses_RF) then
551 >    if (FF_uses_Dipoles) then
552 >       if (electrostaticSummationMethod == 3) then
553            dielect = getDielect()
554            call initialize_rf(dielect)
555         endif
556      else
557 <       if (FF_uses_RF) then          
557 >       if (electrostaticSummationMethod == 3) then
558            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
559            thisStat = -1
560            haveSaneForceField = .false.
561            return
562         endif
563      endif
568
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
564  
565      if (FF_uses_EAM) then
566         call init_EAM_FF(my_status)
# Line 595 | Line 581 | contains
581         endif
582      endif
583  
598    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
599    endif
600
584      if (.not. haveNeighborList) then
585         !! Create neighbor lists
586         call expandNeighborList(nLocal, my_status)
# Line 661 | Line 644 | contains
644      integer :: localError
645      integer :: propPack_i, propPack_j
646      integer :: loopStart, loopEnd, loop
647 <    integer :: iMap
647 >    integer :: iHash
648      real(kind=dp) :: listSkin = 1.0  
649  
650      !! initialize local variables  
# Line 753 | Line 736 | contains
736   #endif
737         outer: do i = istart, iend
738  
756 #ifdef IS_MPI
757             me_i = atid_row(i)
758 #else
759             me_i = atid(i)
760 #endif
761
739            if (update_nlist) point(i) = nlist + 1
740  
741            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 795 | Line 772 | contains
772                    q_group(:,j), d_grp, rgrpsq)
773   #endif
774  
775 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
775 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
776                  if (update_nlist) then
777                     nlist = nlist + 1
778  
# Line 995 | Line 972 | contains
972  
973      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
974  
975 <       if (FF_uses_RF .and. SIM_uses_RF) then
975 >       if (electrostaticSummationMethod == 3) then
976  
977   #ifdef IS_MPI
978            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 1013 | Line 990 | contains
990   #else
991               me_i = atid(i)
992   #endif
993 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
993 >             iHash = InteractionHash(me_i,me_j)
994              
995 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
995 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
996  
997                  mu_i = getDipoleMoment(me_i)
998  
# Line 1079 | Line 1056 | contains
1056      real ( kind = dp ), intent(inout) :: rijsq
1057      real ( kind = dp )                :: r
1058      real ( kind = dp ), intent(inout) :: d(3)
1082    real ( kind = dp ) :: ebalance
1059      integer :: me_i, me_j
1060  
1061 <    integer :: iMap
1061 >    integer :: iHash
1062  
1063      r = sqrt(rijsq)
1064      vpair = 0.0d0
# Line 1096 | Line 1072 | contains
1072      me_j = atid(j)
1073   #endif
1074  
1075 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1075 >    iHash = InteractionHash(me_i, me_j)
1076  
1077 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1077 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1078         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1079      endif
1080  
1081 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1081 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1082         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1083              pot, eFrame, f, t, do_pot)
1084  
1085 <       if (FF_uses_RF .and. SIM_uses_RF) then
1085 >       if (electrostaticSummationMethod == 3) then
1086  
1087            ! CHECK ME (RF needs to know about all electrostatic types)
1088            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1115 | Line 1091 | contains
1091  
1092      endif
1093  
1094 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1094 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1095         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1096              pot, A, f, t, do_pot)
1097      endif
1098  
1099 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1099 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1100         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1101              pot, A, f, t, do_pot)
1102      endif
1103  
1104 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1104 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1105         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1106              pot, A, f, t, do_pot)
1107      endif
1108      
1109 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1109 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1110   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1111   !           pot, A, f, t, do_pot)
1112      endif
1113  
1114 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1114 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1115         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1116              do_pot)
1117      endif
1118  
1119 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1119 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1120         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1121              pot, A, f, t, do_pot)
1122      endif
1123  
1124 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1124 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1125         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1126              pot, A, f, t, do_pot)
1127      endif
# Line 1167 | Line 1143 | contains
1143      real ( kind = dp )                :: r, rc
1144      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1145  
1146 <    integer :: me_i, me_j, iMap
1146 >    integer :: me_i, me_j, iHash
1147  
1148 +    r = sqrt(rijsq)
1149 +
1150   #ifdef IS_MPI  
1151      me_i = atid_row(i)
1152      me_j = atid_col(j)  
# Line 1177 | Line 1155 | contains
1155      me_j = atid(j)  
1156   #endif
1157  
1158 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1158 >    iHash = InteractionHash(me_i, me_j)
1159  
1160 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1160 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1161              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1162      endif
1163      
# Line 1376 | Line 1354 | contains
1354  
1355    function FF_UsesDirectionalAtoms() result(doesit)
1356      logical :: doesit
1357 <    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
1357 >    doesit = FF_uses_DirectionalAtoms
1358    end function FF_UsesDirectionalAtoms
1359  
1360    function FF_RequiresPrepairCalc() result(doesit)
# Line 1388 | Line 1364 | contains
1364  
1365    function FF_RequiresPostpairCalc() result(doesit)
1366      logical :: doesit
1367 <    doesit = FF_uses_RF
1367 >    if (electrostaticSummationMethod == 3) doesit = .true.
1368    end function FF_RequiresPostpairCalc
1369  
1370   #ifdef PROFILE

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines