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 2282 by chuckv, Tue Sep 6 17:32:42 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.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.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    
293            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# 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 <          atypeMaxCutoff(i) = 0.0_dp
302 <          if (i_is_LJ) then
303 <             thisRcut = getSigma(i) * 2.5_dp
304 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
305 <          endif
306 <          if (i_is_Elect) then
307 <             thisRcut = defaultRcut
308 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
309 <          endif
310 <          if (i_is_Sticky) then
311 <             thisRcut = getStickyCut(i)
312 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313 <          endif
314 <          if (i_is_StickyP) then
315 <             thisRcut = getStickyPowerCut(i)
316 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317 <          endif
318 <          if (i_is_GB) then
319 <             thisRcut = getGayBerneCut(i)
320 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321 <          endif
322 <          if (i_is_EAM) then
323 <             thisRcut = getEAMCut(i)
324 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325 <          endif
326 <          if (i_is_Shape) then
327 <             thisRcut = getShapeCut(i)
328 <             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      
355      !! allocate the groupToGtype and gtypeMaxCutoff here.
356 <    if(.not.allocated(groupToGtype)) then
357 <       allocate(groupToGtype(iend))
358 <       allocate(groupMaxCutoff(iend))
359 <       allocate(gtypeMaxCutoff(iend))
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      
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 351 | 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 <          do g = 1, nGroupTypes
436 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
437 <                nGroupTypes = nGroupTypes + 1
438 <                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
439 <                groupToGtype(i) = nGroupTypes
368 <             else
369 <                groupToGtype(i) = g
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 <    
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 >          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    
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(nGroupTypes,nGroupTypes))
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)
507 >             thisRcut = tradRcut
508            case(MIX_CUTOFF_POLICY)
509 <             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
509 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
510            case(MAX_CUTOFF_POLICY)
511 <             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
511 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
512            case default
513               call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
514               return
# Line 394 | Line 516 | contains
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 <    
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)
542 >   end subroutine createGtypeCutoffMap
543  
544 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
545 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
546       integer, intent(in) :: cutPolicy
547 +
548 +     defaultRcut = defRcut
549 +     defaultRsw = defRsw
550 +     defaultRlist = defRlist
551       cutoffPolicy = cutPolicy
419 –     call createGtypeCutoffMap()
552  
553 +     haveDefaultCutoffs = .true.
554 +   end subroutine setDefaultCutoffs
555 +
556 +   subroutine setCutoffPolicy(cutPolicy)
557 +
558 +     integer, intent(in) :: cutPolicy
559 +     cutoffPolicy = cutPolicy
560 +     call createGtypeCutoffMap()
561     end subroutine setCutoffPolicy
562      
563      
564    subroutine setSimVariables()
565      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
566      SIM_uses_EAM = SimUsesEAM()
427 –    SIM_uses_RF = SimUsesRF()
567      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
568      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
569      SIM_uses_PBC = SimUsesPBC()
# Line 494 | 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
500 <    logical, intent(in) :: use_UW
501 <    logical, intent(in) :: use_DW
638 >    integer, intent(in) :: thisESM
639      integer, intent(out) :: thisStat  
640      integer :: my_status, nMatches
504 –    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:
512 <    FF_uses_RF = use_RF
647 >    electrostaticSummationMethod = thisESM
648  
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 –    
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 549 | 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 648 | 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 773 | 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 896 | Line 1022 | contains
1022                  endif
1023               end if
1024            enddo
1025 +
1026         enddo outer
1027  
1028         if (update_nlist) then
# Line 975 | 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 1059 | Line 1186 | contains
1186      real ( kind = dp ), intent(inout) :: rijsq
1187      real ( kind = dp )                :: r
1188      real ( kind = dp ), intent(inout) :: d(3)
1062 –    real ( kind = dp ) :: ebalance
1189      integer :: me_i, me_j
1190  
1191      integer :: iHash
# Line 1084 | 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 1148 | Line 1274 | contains
1274      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1275  
1276      integer :: me_i, me_j, iHash
1277 +
1278 +    r = sqrt(rijsq)
1279  
1280   #ifdef IS_MPI  
1281      me_i = atid_row(i)
# Line 1366 | 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