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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2280 by gezelter, Thu Sep 1 20:17:55 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.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.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          
# 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          
298 >       if (SimHasAtype(i)) then    
299            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
300            call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
301            call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
# 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 <          if (i_is_LJ) then
308 <             thisRcut = getSigma(i) * 2.5_dp
309 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
310 <          endif
311 <          if (i_is_Elect) then
312 <             thisRcut = defaultRcut
313 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314 <          endif
315 <          if (i_is_Sticky) then
316 <             thisRcut = getStickyCut(i)
317 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
318 <          endif
319 <          if (i_is_StickyP) then
320 <             thisRcut = getStickyPowerCut(i)
321 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322 <          endif
323 <          if (i_is_GB) then
324 <             thisRcut = getGayBerneCut(i)
325 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326 <          endif
327 <          if (i_is_EAM) then
328 <             thisRcut = getEAMCut(i)
329 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330 <          endif
331 <          if (i_is_Shape) then
332 <             thisRcut = getShapeCut(i)
333 <             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 <
360 >    
361      !! allocate the groupToGtype and gtypeMaxCutoff here.
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      
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.
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 344 | 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 <          do g = 1, nGroupTypes
442 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
443 <                nGroupTypes = nGroupTypes + 1
444 <                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
445 <                groupToGtype(i) = nGroupTypes
361 <             else
362 <                groupToGtype(i) = g
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
367 <    
368 <    !! allocate the gtypeCutoffMap here.
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 +          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 +
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(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)
514 <             case(MIX_CUTOFF_POLICY)
515 <                thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
516 <             case(MAX_CUTOFF_POLICY)
517 <                thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
518 <             case default
519 <                call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520 <                return
521 <          end select      
522 <         gtypeCutoffMap(i,j)%rcut = thisRcut
523 <         gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
524 <         skin = defaultRlist - defaultRcut
525 <         gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
512 >          case(TRADITIONAL_CUTOFF_POLICY)
513 >             thisRcut = tradRcut
514 >          case(MIX_CUTOFF_POLICY)
515 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516 >          case(MAX_CUTOFF_POLICY)
517 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
518 >          case default
519 >             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520 >             return
521 >          end select
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     end subroutine createGtypeCutoffMap
# Line 402 | Line 555 | contains
555       defaultRsw = defRsw
556       defaultRlist = defRlist
557       cutoffPolicy = cutPolicy
558 +
559 +     haveDefaultCutoffs = .true.
560     end subroutine setDefaultCutoffs
561  
562     subroutine setCutoffPolicy(cutPolicy)
# Line 409 | Line 564 | contains
564       integer, intent(in) :: cutPolicy
565       cutoffPolicy = cutPolicy
566       call createGtypeCutoffMap()
412 –
567     end subroutine setCutoffPolicy
568      
569      
570    subroutine setSimVariables()
571      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
572      SIM_uses_EAM = SimUsesEAM()
419 –    SIM_uses_RF = SimUsesRF()
573      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
574      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
575      SIM_uses_PBC = SimUsesPBC()
# Line 457 | Line 610 | contains
610         call setSimVariables()
611      endif
612  
613 <    if (.not. haveRlist) then
614 <       write(default_error, *) 'rList has not been set in doForces!'
615 <       error = -1
616 <       return
617 <    endif
613 >  !  if (.not. haveRlist) then
614 >  !     write(default_error, *) 'rList has not been set in doForces!'
615 >  !     error = -1
616 >  !     return
617 >  !  endif
618  
619      if (.not. haveNeighborList) then
620         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 486 | 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
492 <    logical, intent(in) :: use_UW
493 <    logical, intent(in) :: use_DW
644 >    integer, intent(in) :: thisESM
645      integer, intent(out) :: thisStat  
646      integer :: my_status, nMatches
496 –    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:
504 <    FF_uses_RF = use_RF
653 >    electrostaticSummationMethod = thisESM
654  
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 –    
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 541 | 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 640 | 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 765 | 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 888 | Line 1028 | contains
1028                  endif
1029               end if
1030            enddo
1031 +
1032         enddo outer
1033  
1034         if (update_nlist) then
# Line 967 | 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 1051 | Line 1192 | contains
1192      real ( kind = dp ), intent(inout) :: rijsq
1193      real ( kind = dp )                :: r
1194      real ( kind = dp ), intent(inout) :: d(3)
1054 –    real ( kind = dp ) :: ebalance
1195      integer :: me_i, me_j
1196  
1197      integer :: iHash
# Line 1076 | 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 1141 | 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 1358 | 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