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 2282 by chuckv, Tue Sep 6 17:32:42 2005 UTC vs.
Revision 2354 by chuckv, Tue Oct 11 22:00:30 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.36 2005-09-06 17:32:42 chuckv Exp $, $Date: 2005-09-06 17:32:42 $, $Name: not supported by cvs2svn $, $Revision: 1.36 $
48 > !! @version $Id: doForces.F90,v 1.53 2005-10-11 22:00:30 chuckv Exp $, $Date: 2005-10-11 22:00:30 $, $Name: not supported by cvs2svn $, $Revision: 1.53 $
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 75 | Line 75 | module doForces
75   #include "UseTheForce/fSwitchingFunction.h"
76   #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79  
80  
81    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
# Line 85 | Line 86 | module doForces
86    logical, save :: haveSaneForceField = .false.
87    logical, save :: haveInteractionHash = .false.
88    logical, save :: haveGtypeCutoffMap = .false.
89 +  logical, save :: haveDefaultCutoffs = .false.
90    logical, save :: haveRlist = .false.
91  
92    logical, save :: FF_uses_DirectionalAtoms
93    logical, save :: FF_uses_Dipoles
94    logical, save :: FF_uses_GayBerne
95    logical, save :: FF_uses_EAM
94 –  logical, save :: FF_uses_RF
96  
97    logical, save :: SIM_uses_DirectionalAtoms
98    logical, save :: SIM_uses_EAM
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
102  
103 <  integer, save :: corrMethod
103 >  integer, save :: electrostaticSummationMethod
104  
105    public :: init_FF
106    public :: setDefaultCutoffs
# Line 124 | Line 124 | module doForces
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
127 >  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
128 >  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
129 >
130 >  integer, dimension(:), allocatable, target :: groupToGtypeRow
131 >  integer, dimension(:), pointer :: groupToGtypeCol => null()
132 >
133 >  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
134 >  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
135    type ::gtypeCutoffs
136       real(kind=dp) :: rcut
137       real(kind=dp) :: rcutsq
# Line 136 | Line 141 | module doForces
141  
142    integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
143    real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
144 +  real(kind=dp),save :: listSkin
145    
146   contains
147  
# Line 179 | Line 185 | contains
185  
186      if (.not. allocated(InteractionHash)) then
187         allocate(InteractionHash(nAtypes,nAtypes))
188 +    else
189 +       deallocate(InteractionHash)
190 +       allocate(InteractionHash(nAtypes,nAtypes))
191      endif
192  
193      if (.not. allocated(atypeMaxCutoff)) then
194         allocate(atypeMaxCutoff(nAtypes))
195 +    else
196 +       deallocate(atypeMaxCutoff)
197 +       allocate(atypeMaxCutoff(nAtypes))
198      endif
199          
200      do i = 1, nAtypes
# Line 256 | Line 268 | contains
268      logical :: i_is_GB
269      logical :: i_is_EAM
270      logical :: i_is_Shape
271 +    logical :: GtypeFound
272  
273      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
274 <    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
275 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
274 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
275 >    integer :: nGroupsInRow
276 >    integer :: nGroupsInCol
277 >    integer :: nGroupTypesRow,nGroupTypesCol
278 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
279      real(kind=dp) :: biggestAtypeCutoff
280  
281      stat = 0
# Line 271 | Line 287 | contains
287            return
288         endif
289      endif
290 <
290 > #ifdef IS_MPI
291 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
292 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
293 > #endif
294      nAtypes = getSize(atypes)
295 <    
295 > ! Set all of the initial cutoffs to zero.
296 >    atypeMaxCutoff = 0.0_dp
297      do i = 1, nAtypes
298         if (SimHasAtype(i)) then    
299            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 284 | Line 304 | contains
304            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
305            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
306            
307 <          atypeMaxCutoff(i) = 0.0_dp
308 <          if (i_is_LJ) then
309 <             thisRcut = getSigma(i) * 2.5_dp
310 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
311 <          endif
312 <          if (i_is_Elect) then
313 <             thisRcut = defaultRcut
314 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
315 <          endif
316 <          if (i_is_Sticky) then
317 <             thisRcut = getStickyCut(i)
318 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
319 <          endif
320 <          if (i_is_StickyP) then
321 <             thisRcut = getStickyPowerCut(i)
322 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
323 <          endif
324 <          if (i_is_GB) then
325 <             thisRcut = getGayBerneCut(i)
326 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
327 <          endif
328 <          if (i_is_EAM) then
329 <             thisRcut = getEAMCut(i)
330 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
331 <          endif
332 <          if (i_is_Shape) then
333 <             thisRcut = getShapeCut(i)
334 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
307 >
308 >          if (haveDefaultCutoffs) then
309 >             atypeMaxCutoff(i) = defaultRcut
310 >          else
311 >             if (i_is_LJ) then          
312 >                thisRcut = getSigma(i) * 2.5_dp
313 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314 >             endif
315 >             if (i_is_Elect) then
316 >                thisRcut = defaultRcut
317 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
318 >             endif
319 >             if (i_is_Sticky) then
320 >                thisRcut = getStickyCut(i)
321 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322 >             endif
323 >             if (i_is_StickyP) then
324 >                thisRcut = getStickyPowerCut(i)
325 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326 >             endif
327 >             if (i_is_GB) then
328 >                thisRcut = getGayBerneCut(i)
329 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330 >             endif
331 >             if (i_is_EAM) then
332 >                thisRcut = getEAMCut(i)
333 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334 >             endif
335 >             if (i_is_Shape) then
336 >                thisRcut = getShapeCut(i)
337 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
338 >             endif
339            endif
340            
341 +          
342            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
343               biggestAtypeCutoff = atypeMaxCutoff(i)
344            endif
345 +
346         endif
347      enddo
348    
349 <    nGroupTypes = 0
349 >
350      
351      istart = 1
352 +    jstart = 1
353   #ifdef IS_MPI
354      iend = nGroupsInRow
355 +    jend = nGroupsInCol
356   #else
357      iend = nGroups
358 +    jend = nGroups
359   #endif
360      
361      !! allocate the groupToGtype and gtypeMaxCutoff here.
362 <    if(.not.allocated(groupToGtype)) then
363 <       allocate(groupToGtype(iend))
364 <       allocate(groupMaxCutoff(iend))
365 <       allocate(gtypeMaxCutoff(iend))
362 >    if(.not.allocated(groupToGtypeRow)) then
363 >     !  allocate(groupToGtype(iend))
364 >       allocate(groupToGtypeRow(iend))
365 >    else
366 >       deallocate(groupToGtypeRow)
367 >       allocate(groupToGtypeRow(iend))
368      endif
369 +    if(.not.allocated(groupMaxCutoffRow)) then
370 +       allocate(groupMaxCutoffRow(iend))
371 +    else
372 +       deallocate(groupMaxCutoffRow)
373 +       allocate(groupMaxCutoffRow(iend))
374 +    end if
375 +
376 +    if(.not.allocated(gtypeMaxCutoffRow)) then
377 +       allocate(gtypeMaxCutoffRow(iend))
378 +    else
379 +       deallocate(gtypeMaxCutoffRow)
380 +       allocate(gtypeMaxCutoffRow(iend))
381 +    endif
382 +
383 +
384 + #ifdef IS_MPI
385 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
386 +    if(.not.associated(groupToGtypeCol)) then
387 +       allocate(groupToGtypeCol(jend))
388 +    else
389 +       deallocate(groupToGtypeCol)
390 +       allocate(groupToGtypeCol(jend))
391 +    end if
392 +
393 +    if(.not.associated(groupToGtypeCol)) then
394 +       allocate(groupToGtypeCol(jend))
395 +    else
396 +       deallocate(groupToGtypeCol)
397 +       allocate(groupToGtypeCol(jend))
398 +    end if
399 +    if(.not.associated(gtypeMaxCutoffCol)) then
400 +       allocate(gtypeMaxCutoffCol(jend))
401 +    else
402 +       deallocate(gtypeMaxCutoffCol)      
403 +       allocate(gtypeMaxCutoffCol(jend))
404 +    end if
405 +
406 +       groupMaxCutoffCol = 0.0_dp
407 +       gtypeMaxCutoffCol = 0.0_dp
408 +
409 + #endif
410 +       groupMaxCutoffRow = 0.0_dp
411 +       gtypeMaxCutoffRow = 0.0_dp
412 +
413 +
414      !! first we do a single loop over the cutoff groups to find the
415      !! largest cutoff for any atypes present in this group.  We also
416      !! create gtypes at this point.
417      
418      tol = 1.0d-6
419 <    
419 >    nGroupTypesRow = 0
420 >
421      do i = istart, iend      
422         n_in_i = groupStartRow(i+1) - groupStartRow(i)
423 <       groupMaxCutoff(i) = 0.0_dp
423 >       groupMaxCutoffRow(i) = 0.0_dp
424         do ia = groupStartRow(i), groupStartRow(i+1)-1
425            atom1 = groupListRow(ia)
426   #ifdef IS_MPI
# Line 351 | Line 428 | contains
428   #else
429            me_i = atid(atom1)
430   #endif          
431 <          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
432 <             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
433 <          endif
431 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
432 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
433 >          endif          
434         enddo
435 <       if (nGroupTypes.eq.0) then
436 <          nGroupTypes = nGroupTypes + 1
437 <          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
438 <          groupToGtype(i) = nGroupTypes
435 >
436 >       if (nGroupTypesRow.eq.0) then
437 >          nGroupTypesRow = nGroupTypesRow + 1
438 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
439 >          groupToGtypeRow(i) = nGroupTypesRow
440 >       else
441 >          GtypeFound = .false.
442 >          do g = 1, nGroupTypesRow
443 >             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
444 >                groupToGtypeRow(i) = g
445 >                GtypeFound = .true.
446 >             endif
447 >          enddo
448 >          if (.not.GtypeFound) then            
449 >             nGroupTypesRow = nGroupTypesRow + 1
450 >             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
451 >             groupToGtypeRow(i) = nGroupTypesRow
452 >          endif
453 >       endif
454 >    enddo    
455 >
456 > #ifdef IS_MPI
457 >    do j = jstart, jend      
458 >       n_in_j = groupStartCol(j+1) - groupStartCol(j)
459 >       groupMaxCutoffCol(j) = 0.0_dp
460 >       do ja = groupStartCol(j), groupStartCol(j+1)-1
461 >          atom1 = groupListCol(ja)
462 >
463 >          me_j = atid_col(atom1)
464 >
465 >          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
466 >             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
467 >          endif          
468 >       enddo
469 >
470 >       if (nGroupTypesCol.eq.0) then
471 >          nGroupTypesCol = nGroupTypesCol + 1
472 >          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
473 >          groupToGtypeCol(j) = nGroupTypesCol
474         else
475 <          do g = 1, nGroupTypes
476 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
477 <                nGroupTypes = nGroupTypes + 1
478 <                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
479 <                groupToGtype(i) = nGroupTypes
368 <             else
369 <                groupToGtype(i) = g
475 >          GtypeFound = .false.
476 >          do g = 1, nGroupTypesCol
477 >             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
478 >                groupToGtypeCol(j) = g
479 >                GtypeFound = .true.
480               endif
481            enddo
482 +          if (.not.GtypeFound) then            
483 +             nGroupTypesCol = nGroupTypesCol + 1
484 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
485 +             groupToGtypeCol(j) = nGroupTypesCol
486 +          endif
487         endif
488 <    enddo
489 <    
488 >    enddo    
489 >
490 > #else
491 > ! Set pointers to information we just found
492 >    nGroupTypesCol = nGroupTypesRow
493 >    groupToGtypeCol => groupToGtypeRow
494 >    gtypeMaxCutoffCol => gtypeMaxCutoffRow
495 >    groupMaxCutoffCol => groupMaxCutoffRow
496 > #endif
497 >
498 >
499 >
500 >
501 >
502      !! allocate the gtypeCutoffMap here.
503 <    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
503 >    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
504      !! then we do a double loop over all the group TYPES to find the cutoff
505      !! map between groups of two types
506 <    
507 <    do i = 1, nGroupTypes
508 <       do j = 1, nGroupTypes
506 >    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507 >
508 >    do i = 1, nGroupTypesRow
509 >       do j = 1, nGroupTypesCol
510        
511            select case(cutoffPolicy)
512            case(TRADITIONAL_CUTOFF_POLICY)
513 <             thisRcut = maxval(gtypeMaxCutoff)
513 >             thisRcut = tradRcut
514            case(MIX_CUTOFF_POLICY)
515 <             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
515 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516            case(MAX_CUTOFF_POLICY)
517 <             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
517 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
518            case default
519               call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520               return
# Line 394 | Line 522 | contains
522            gtypeCutoffMap(i,j)%rcut = thisRcut
523            gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
524            skin = defaultRlist - defaultRcut
525 +          listSkin = skin ! set neighbor list skin thickness
526            gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
527 +
528 +          ! sanity check
529 +
530 +          if (haveDefaultCutoffs) then
531 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
532 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
533 +             endif
534 +          endif
535         enddo
536      enddo
537 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
538 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
539 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
540 + #ifdef IS_MPI
541 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
542 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
543 + #endif
544 +    groupMaxCutoffCol => null()
545 +    gtypeMaxCutoffCol => null()
546      
547      haveGtypeCutoffMap = .true.
548 <    
403 <  end subroutine createGtypeCutoffMap
404 <  
405 <  subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
406 <    real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
407 <    integer, intent(in) :: cutPolicy
408 <    
409 <    defaultRcut = defRcut
410 <    defaultRsw = defRsw
411 <    defaultRlist = defRlist
412 <    cutoffPolicy = cutPolicy
413 <  end subroutine setDefaultCutoffs
414 <  
415 <  subroutine setCutoffPolicy(cutPolicy)
548 >   end subroutine createGtypeCutoffMap
549  
550 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
551 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
552       integer, intent(in) :: cutPolicy
553 +
554 +     defaultRcut = defRcut
555 +     defaultRsw = defRsw
556 +     defaultRlist = defRlist
557       cutoffPolicy = cutPolicy
419 –     call createGtypeCutoffMap()
558  
559 +     haveDefaultCutoffs = .true.
560 +   end subroutine setDefaultCutoffs
561 +
562 +   subroutine setCutoffPolicy(cutPolicy)
563 +
564 +     integer, intent(in) :: cutPolicy
565 +     cutoffPolicy = cutPolicy
566 +     call createGtypeCutoffMap()
567     end subroutine setCutoffPolicy
568      
569      
570    subroutine setSimVariables()
571      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
572      SIM_uses_EAM = SimUsesEAM()
427 –    SIM_uses_RF = SimUsesRF()
573      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
574      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
575      SIM_uses_PBC = SimUsesPBC()
# Line 494 | Line 639 | contains
639    end subroutine doReadyCheck
640  
641  
642 <  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
642 >  subroutine init_FF(thisESM, thisStat)
643  
644 <    logical, intent(in) :: use_RF
500 <    logical, intent(in) :: use_UW
501 <    logical, intent(in) :: use_DW
644 >    integer, intent(in) :: thisESM
645      integer, intent(out) :: thisStat  
646      integer :: my_status, nMatches
504 –    integer :: corrMethod
647      integer, pointer :: MatchList(:) => null()
648      real(kind=dp) :: rcut, rrf, rt, dielect
649  
650      !! assume things are copacetic, unless they aren't
651      thisStat = 0
652  
653 <    !! Fortran's version of a cast:
512 <    FF_uses_RF = use_RF
653 >    electrostaticSummationMethod = thisESM
654  
514 –    !! set the electrostatic correction method
515 –    if (use_UW) then
516 –       corrMethod = 1
517 –    elseif (use_DW) then
518 –       corrMethod = 2
519 –    else
520 –       corrMethod = 0
521 –    endif
522 –    
655      !! init_FF is called *after* all of the atom types have been
656      !! defined in atype_module using the new_atype subroutine.
657      !!
# Line 549 | Line 681 | contains
681  
682      haveSaneForceField = .true.
683  
684 <    !! check to make sure the FF_uses_RF setting makes sense
684 >    !! check to make sure the reaction field setting makes sense
685  
686      if (FF_uses_Dipoles) then
687 <       if (FF_uses_RF) then
687 >       if (electrostaticSummationMethod == REACTION_FIELD) then
688            dielect = getDielect()
689            call initialize_rf(dielect)
690         endif
691      else
692 <       if (FF_uses_RF) then          
692 >       if (electrostaticSummationMethod == REACTION_FIELD) then
693            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
694            thisStat = -1
695            haveSaneForceField = .false.
# Line 648 | Line 780 | contains
780      integer :: propPack_i, propPack_j
781      integer :: loopStart, loopEnd, loop
782      integer :: iHash
783 <    real(kind=dp) :: listSkin = 1.0  
783 >  
784  
785      !! initialize local variables  
786  
# Line 773 | Line 905 | contains
905               me_j = atid(j)
906               call get_interatomic_vector(q_group(:,i), &
907                    q_group(:,j), d_grp, rgrpsq)
908 < #endif
908 > #endif      
909  
910 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
910 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
911                  if (update_nlist) then
912                     nlist = nlist + 1
913  
# Line 896 | Line 1028 | contains
1028                  endif
1029               end if
1030            enddo
1031 +
1032         enddo outer
1033  
1034         if (update_nlist) then
# Line 975 | Line 1108 | contains
1108  
1109      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1110  
1111 <       if (FF_uses_RF .and. SIM_uses_RF) then
1111 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1112  
1113   #ifdef IS_MPI
1114            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 1059 | Line 1192 | contains
1192      real ( kind = dp ), intent(inout) :: rijsq
1193      real ( kind = dp )                :: r
1194      real ( kind = dp ), intent(inout) :: d(3)
1062 –    real ( kind = dp ) :: ebalance
1195      integer :: me_i, me_j
1196  
1197      integer :: iHash
# Line 1084 | Line 1216 | contains
1216  
1217      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1218         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1219 <            pot, eFrame, f, t, do_pot, corrMethod)
1219 >            pot, eFrame, f, t, do_pot)
1220  
1221 <       if (FF_uses_RF .and. SIM_uses_RF) then
1221 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1222  
1223            ! CHECK ME (RF needs to know about all electrostatic types)
1224            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1149 | Line 1281 | contains
1281  
1282      integer :: me_i, me_j, iHash
1283  
1284 +    r = sqrt(rijsq)
1285 +
1286   #ifdef IS_MPI  
1287      me_i = atid_row(i)
1288      me_j = atid_col(j)  
# Line 1366 | Line 1500 | contains
1500  
1501    function FF_RequiresPostpairCalc() result(doesit)
1502      logical :: doesit
1503 <    doesit = FF_uses_RF
1503 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1504    end function FF_RequiresPostpairCalc
1505  
1506   #ifdef PROFILE

Diff Legend

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