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

Comparing trunk/OOPSE-3.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2280 by gezelter, Thu Sep 1 20:17:55 2005 UTC vs.
Revision 2350 by chuckv, Mon Oct 10 21:20:46 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.34 2005-09-01 20:17:55 gezelter Exp $, $Date: 2005-09-01 20:17:55 $, $Name: not supported by cvs2svn $, $Revision: 1.34 $
48 > !! @version $Id: doForces.F90,v 1.51 2005-10-10 21:20:46 chuckv Exp $, $Date: 2005-10-10 21:20:46 $, $Name: not supported by cvs2svn $, $Revision: 1.51 $
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 256 | Line 262 | contains
262      logical :: i_is_GB
263      logical :: i_is_EAM
264      logical :: i_is_Shape
265 +    logical :: GtypeFound
266  
267      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
268 <    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
269 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
268 >    integer :: n_in_i, me_i, ia, g, atom1
269 >    integer :: nGroupsInRow
270 >    integer :: nGroupsInCol
271 >    integer :: nGroupTypesRow,nGroupTypesCol
272 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
273      real(kind=dp) :: biggestAtypeCutoff
274  
275      stat = 0
# Line 271 | Line 281 | contains
281            return
282         endif
283      endif
284 <
284 > #ifdef IS_MPI
285 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
286 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
287 > #endif
288      nAtypes = getSize(atypes)
289 <    
289 > ! Set all of the initial cutoffs to zero.
290 >    atypeMaxCutoff = 0.0_dp
291      do i = 1, nAtypes
292 <       if (SimHasAtype(i)) then          
292 >       if (SimHasAtype(i)) then    
293            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
294            call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
295            call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
# Line 284 | Line 298 | contains
298            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
299            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
300            
301 <          if (i_is_LJ) then
302 <             thisRcut = getSigma(i) * 2.5_dp
303 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
304 <          endif
305 <          if (i_is_Elect) then
306 <             thisRcut = defaultRcut
307 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 <          endif
309 <          if (i_is_Sticky) then
310 <             thisRcut = getStickyCut(i)
311 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312 <          endif
313 <          if (i_is_StickyP) then
314 <             thisRcut = getStickyPowerCut(i)
315 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 <          endif
317 <          if (i_is_GB) then
318 <             thisRcut = getGayBerneCut(i)
319 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
320 <          endif
321 <          if (i_is_EAM) then
322 <             thisRcut = getEAMCut(i)
323 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
324 <          endif
325 <          if (i_is_Shape) then
326 <             thisRcut = getShapeCut(i)
327 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
301 >
302 >          if (haveDefaultCutoffs) then
303 >             atypeMaxCutoff(i) = defaultRcut
304 >          else
305 >             if (i_is_LJ) then          
306 >                thisRcut = getSigma(i) * 2.5_dp
307 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 >             endif
309 >             if (i_is_Elect) then
310 >                thisRcut = defaultRcut
311 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312 >             endif
313 >             if (i_is_Sticky) then
314 >                thisRcut = getStickyCut(i)
315 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 >             endif
317 >             if (i_is_StickyP) then
318 >                thisRcut = getStickyPowerCut(i)
319 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
320 >             endif
321 >             if (i_is_GB) then
322 >                thisRcut = getGayBerneCut(i)
323 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
324 >             endif
325 >             if (i_is_EAM) then
326 >                thisRcut = getEAMCut(i)
327 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
328 >             endif
329 >             if (i_is_Shape) then
330 >                thisRcut = getShapeCut(i)
331 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
332 >             endif
333            endif
334            
335 +          
336            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
337               biggestAtypeCutoff = atypeMaxCutoff(i)
338            endif
339 +
340         endif
341      enddo
342    
343 <    nGroupTypes = 0
343 >
344      
345      istart = 1
346 +    jstart = 1
347   #ifdef IS_MPI
348      iend = nGroupsInRow
349 +    jend = nGroupsInCol
350   #else
351      iend = nGroups
352 +    jend = nGroups
353   #endif
354 <
354 >    
355      !! allocate the groupToGtype and gtypeMaxCutoff here.
356 +    if(.not.allocated(groupToGtypeRow)) then
357 +     !  allocate(groupToGtype(iend))
358 +       allocate(groupToGtypeRow(iend))
359 +    else
360 +       deallocate(groupToGtypeRow)
361 +       allocate(groupToGtypeRow(iend))
362 +    endif
363 +    if(.not.allocated(groupMaxCutoffRow)) then
364 +       allocate(groupMaxCutoffRow(iend))
365 +    else
366 +       deallocate(groupMaxCutoffRow)
367 +       allocate(groupMaxCutoffRow(iend))
368 +    end if
369 +
370 +    if(.not.allocated(gtypeMaxCutoffRow)) then
371 +       allocate(gtypeMaxCutoffRow(iend))
372 +    else
373 +       deallocate(gtypeMaxCutoffRow)
374 +       allocate(gtypeMaxCutoffRow(iend))
375 +    endif
376 +
377 +
378 + #ifdef IS_MPI
379 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
380 +    if(.not.allocated(groupToGtypeCol)) then
381 +       allocate(groupToGtypeCol(jend))
382 +    else
383 +       deallocate(groupToGtypeCol)
384 +       allocate(groupToGtypeCol(jend))
385 +    end if
386 +
387 +    if(.not.allocated(groupToGtypeCol)) then
388 +       allocate(groupToGtypeCol(jend))
389 +    else
390 +       deallocate(groupToGtypeCol)
391 +       allocate(groupToGtypeCol(jend))
392 +    end if
393 +    if(.not.allocated(gtypeMaxCutoffCol)) then
394 +       allocate(gtypeMaxCutoffCol(jend))
395 +    else
396 +       deallocate(gtypeMaxCutoffCol)      
397 +       allocate(gtypeMaxCutoffCol(jend))
398 +    end if
399 +
400 +       groupMaxCutoffCol = 0.0_dp
401 +       gtypeMaxCutoffCol = 0.0_dp
402 +
403 + #endif
404 +       groupMaxCutoffRow = 0.0_dp
405 +       gtypeMaxCutoffRow = 0.0_dp
406 +
407 +
408 +    !! first we do a single loop over the cutoff groups to find the
409 +    !! largest cutoff for any atypes present in this group.  We also
410 +    !! create gtypes at this point.
411      
333 –    !! first we do a single loop over the cutoff groups to find the largest cutoff for any atypes
334 –    !! present in this group.   We also create gtypes at this point.
412      tol = 1.0d-6
413 <    
413 >    nGroupTypesRow = 0
414 >
415      do i = istart, iend      
416         n_in_i = groupStartRow(i+1) - groupStartRow(i)
417 <       groupMaxCutoff(i) = 0.0_dp
417 >       groupMaxCutoffRow(i) = 0.0_dp
418         do ia = groupStartRow(i), groupStartRow(i+1)-1
419            atom1 = groupListRow(ia)
420   #ifdef IS_MPI
# Line 344 | Line 422 | contains
422   #else
423            me_i = atid(atom1)
424   #endif          
425 <          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
426 <             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
427 <          endif
425 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
426 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
427 >          endif          
428         enddo
429 <       if (nGroupTypes.eq.0) then
430 <          nGroupTypes = nGroupTypes + 1
431 <          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
432 <          groupToGtype(i) = nGroupTypes
429 >
430 >       if (nGroupTypesRow.eq.0) then
431 >          nGroupTypesRow = nGroupTypesRow + 1
432 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
433 >          groupToGtypeRow(i) = nGroupTypesRow
434 >       else
435 >          GtypeFound = .false.
436 >          do g = 1, nGroupTypesRow
437 >             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
438 >                groupToGtypeRow(i) = g
439 >                GtypeFound = .true.
440 >             endif
441 >          enddo
442 >          if (.not.GtypeFound) then            
443 >             nGroupTypesRow = nGroupTypesRow + 1
444 >             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
445 >             groupToGtypeRow(i) = nGroupTypesRow
446 >          endif
447 >       endif
448 >    enddo    
449 >
450 > #ifdef IS_MPI
451 >    do j = jstart, jend      
452 >       n_in_j = groupStartCol(j+1) - groupStartCol(j)
453 >       groupMaxCutoffCol(j) = 0.0_dp
454 >       do ja = groupStartCol(j), groupStartCol(j+1)-1
455 >          atom1 = groupListCol(ja)
456 >
457 >          me_j = atid_col(atom1)
458 >
459 >          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
460 >             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
461 >          endif          
462 >       enddo
463 >
464 >       if (nGroupTypesCol.eq.0) then
465 >          nGroupTypesCol = nGroupTypesCol + 1
466 >          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
467 >          groupToGtypeCol(j) = nGroupTypesCol
468         else
469 <          do g = 1, nGroupTypes
470 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
471 <                nGroupTypes = nGroupTypes + 1
472 <                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
473 <                groupToGtype(i) = nGroupTypes
361 <             else
362 <                groupToGtype(i) = g
469 >          GtypeFound = .false.
470 >          do g = 1, nGroupTypesCol
471 >             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
472 >                groupToGtypeCol(j) = g
473 >                GtypeFound = .true.
474               endif
475            enddo
476 +          if (.not.GtypeFound) then            
477 +             nGroupTypesCol = nGroupTypesCol + 1
478 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
479 +             groupToGtypeCol(j) = nGroupTypesCol
480 +          endif
481         endif
482 <    enddo
367 <    
368 <    !! allocate the gtypeCutoffMap here.
482 >    enddo    
483  
484 + #else
485 + ! Set pointers to information we just found
486 +    nGroupTypesCol = nGroupTypesRow
487 +    groupToGtypeCol => groupToGtypeRow
488 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
489 +    groupMaxCutoffCol => groupMaxCutoffRow
490 + #endif
491 +
492 +
493 +
494 +
495 +
496 +    !! allocate the gtypeCutoffMap here.
497 +    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
498      !! then we do a double loop over all the group TYPES to find the cutoff
499      !! map between groups of two types
500 <    
501 <    do i = 1, nGroupTypes
502 <       do j = 1, nGroupTypes
500 >    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
501 >
502 >    do i = 1, nGroupTypesRow
503 >       do j = 1, nGroupTypesCol
504        
505            select case(cutoffPolicy)
506 <             case(TRADITIONAL_CUTOFF_POLICY)
507 <                thisRcut = maxval(gtypeMaxCutoff)
508 <             case(MIX_CUTOFF_POLICY)
509 <                thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
510 <             case(MAX_CUTOFF_POLICY)
511 <                thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
512 <             case default
513 <                call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
514 <                return
515 <          end select      
516 <         gtypeCutoffMap(i,j)%rcut = thisRcut
517 <         gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
518 <         skin = defaultRlist - defaultRcut
519 <         gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
506 >          case(TRADITIONAL_CUTOFF_POLICY)
507 >             thisRcut = tradRcut
508 >          case(MIX_CUTOFF_POLICY)
509 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
510 >          case(MAX_CUTOFF_POLICY)
511 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
512 >          case default
513 >             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
514 >             return
515 >          end select
516 >          gtypeCutoffMap(i,j)%rcut = thisRcut
517 >          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
518 >          skin = defaultRlist - defaultRcut
519 >          listSkin = skin ! set neighbor list skin thickness
520 >          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
521 >
522 >          ! sanity check
523 >
524 >          if (haveDefaultCutoffs) then
525 >             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
526 >                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
527 >             endif
528 >          endif
529         enddo
530      enddo
531 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
532 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
533 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
534 + #ifdef IS_MPI
535 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
536 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
537 + #endif
538 +    groupMaxCutoffCol => null()
539 +    gtypeMaxCutoffCol => null()
540      
541      haveGtypeCutoffMap = .true.
542     end subroutine createGtypeCutoffMap
# Line 402 | Line 549 | contains
549       defaultRsw = defRsw
550       defaultRlist = defRlist
551       cutoffPolicy = cutPolicy
552 +
553 +     haveDefaultCutoffs = .true.
554     end subroutine setDefaultCutoffs
555  
556     subroutine setCutoffPolicy(cutPolicy)
# Line 409 | Line 558 | contains
558       integer, intent(in) :: cutPolicy
559       cutoffPolicy = cutPolicy
560       call createGtypeCutoffMap()
412 –
561     end subroutine setCutoffPolicy
562      
563      
564    subroutine setSimVariables()
565      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
566      SIM_uses_EAM = SimUsesEAM()
419 –    SIM_uses_RF = SimUsesRF()
567      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
568      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
569      SIM_uses_PBC = SimUsesPBC()
# Line 457 | Line 604 | contains
604         call setSimVariables()
605      endif
606  
607 <    if (.not. haveRlist) then
608 <       write(default_error, *) 'rList has not been set in doForces!'
609 <       error = -1
610 <       return
611 <    endif
607 >  !  if (.not. haveRlist) then
608 >  !     write(default_error, *) 'rList has not been set in doForces!'
609 >  !     error = -1
610 >  !     return
611 >  !  endif
612  
613      if (.not. haveNeighborList) then
614         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 486 | Line 633 | contains
633    end subroutine doReadyCheck
634  
635  
636 <  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
636 >  subroutine init_FF(thisESM, thisStat)
637  
638 <    logical, intent(in) :: use_RF
492 <    logical, intent(in) :: use_UW
493 <    logical, intent(in) :: use_DW
638 >    integer, intent(in) :: thisESM
639      integer, intent(out) :: thisStat  
640      integer :: my_status, nMatches
496 –    integer :: corrMethod
641      integer, pointer :: MatchList(:) => null()
642      real(kind=dp) :: rcut, rrf, rt, dielect
643  
644      !! assume things are copacetic, unless they aren't
645      thisStat = 0
646  
647 <    !! Fortran's version of a cast:
504 <    FF_uses_RF = use_RF
647 >    electrostaticSummationMethod = thisESM
648  
506 –    !! set the electrostatic correction method
507 –    if (use_UW) then
508 –       corrMethod = 1
509 –    elseif (use_DW) then
510 –       corrMethod = 2
511 –    else
512 –       corrMethod = 0
513 –    endif
514 –    
649      !! init_FF is called *after* all of the atom types have been
650      !! defined in atype_module using the new_atype subroutine.
651      !!
# Line 541 | Line 675 | contains
675  
676      haveSaneForceField = .true.
677  
678 <    !! check to make sure the FF_uses_RF setting makes sense
678 >    !! check to make sure the reaction field setting makes sense
679  
680      if (FF_uses_Dipoles) then
681 <       if (FF_uses_RF) then
681 >       if (electrostaticSummationMethod == REACTION_FIELD) then
682            dielect = getDielect()
683            call initialize_rf(dielect)
684         endif
685      else
686 <       if (FF_uses_RF) then          
686 >       if (electrostaticSummationMethod == REACTION_FIELD) then
687            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
688            thisStat = -1
689            haveSaneForceField = .false.
# Line 640 | Line 774 | contains
774      integer :: propPack_i, propPack_j
775      integer :: loopStart, loopEnd, loop
776      integer :: iHash
777 <    real(kind=dp) :: listSkin = 1.0  
777 >  
778  
779      !! initialize local variables  
780  
# Line 765 | Line 899 | contains
899               me_j = atid(j)
900               call get_interatomic_vector(q_group(:,i), &
901                    q_group(:,j), d_grp, rgrpsq)
902 < #endif
902 > #endif      
903  
904 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
904 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
905                  if (update_nlist) then
906                     nlist = nlist + 1
907  
# Line 888 | Line 1022 | contains
1022                  endif
1023               end if
1024            enddo
1025 +
1026         enddo outer
1027  
1028         if (update_nlist) then
# Line 967 | Line 1102 | contains
1102  
1103      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1104  
1105 <       if (FF_uses_RF .and. SIM_uses_RF) then
1105 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1106  
1107   #ifdef IS_MPI
1108            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 1051 | Line 1186 | contains
1186      real ( kind = dp ), intent(inout) :: rijsq
1187      real ( kind = dp )                :: r
1188      real ( kind = dp ), intent(inout) :: d(3)
1054 –    real ( kind = dp ) :: ebalance
1189      integer :: me_i, me_j
1190  
1191      integer :: iHash
# Line 1076 | Line 1210 | contains
1210  
1211      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1212         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213 <            pot, eFrame, f, t, do_pot, corrMethod)
1213 >            pot, eFrame, f, t, do_pot)
1214  
1215 <       if (FF_uses_RF .and. SIM_uses_RF) then
1215 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1216  
1217            ! CHECK ME (RF needs to know about all electrostatic types)
1218            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1141 | Line 1275 | contains
1275  
1276      integer :: me_i, me_j, iHash
1277  
1278 +    r = sqrt(rijsq)
1279 +
1280   #ifdef IS_MPI  
1281      me_i = atid_row(i)
1282      me_j = atid_col(j)  
# Line 1358 | Line 1494 | contains
1494  
1495    function FF_RequiresPostpairCalc() result(doesit)
1496      logical :: doesit
1497 <    doesit = FF_uses_RF
1497 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1498    end function FF_RequiresPostpairCalc
1499  
1500   #ifdef PROFILE

Diff Legend

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