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 2277 by chrisfen, Fri Aug 26 21:30:41 2005 UTC vs.
Revision 2355 by chuckv, Wed Oct 12 18:59:16 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.32 2005-08-26 21:30:30 chrisfen Exp $, $Date: 2005-08-26 21:30:30 $, $Name: not supported by cvs2svn $, $Revision: 1.32 $
48 > !! @version $Id: doForces.F90,v 1.54 2005-10-12 18:59:16 chuckv Exp $, $Date: 2005-10-12 18:59:16 $, $Name: not supported by cvs2svn $, $Revision: 1.54 $
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
93 –  logical, save :: FF_uses_RF
96  
97    logical, save :: SIM_uses_DirectionalAtoms
98    logical, save :: SIM_uses_EAM
97 –  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 :: electrostaticSummationMethod
104 +
105    public :: init_FF
106    public :: setDefaultCutoffs
107    public :: do_force_loop
# Line 121 | 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 133 | 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 176 | 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 253 | 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
275 <    real(kind=dp):: thisSigma, bigSigma, thisRcut
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 268 | 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 281 | 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  
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 <    outer: do i = istart, iend
361 <      
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 >    
418 >    tol = 1.0d-6
419 >    nGroupTypesRow = 0
420 >
421 >    do i = istart, iend      
422         n_in_i = groupStartRow(i+1) - groupStartRow(i)
423 <      
423 >       groupMaxCutoffRow(i) = 0.0_dp
424 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
425 >          atom1 = groupListRow(ia)
426   #ifdef IS_MPI
427 <       jstart = 1
331 <       jend = nGroupsInCol
427 >          me_i = atid_row(atom1)
428   #else
429 <       jstart = i+1
430 <       jend = nGroups
429 >          me_i = atid(atom1)
430 > #endif          
431 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
432 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
433 >          endif          
434 >       enddo
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 >          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 +    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507 +
508 +    do i = 1, nGroupTypesRow
509 +       do j = 1, nGroupTypesCol
510        
511 <      
512 <      
513 <      
514 <      
515 <      
516 <    enddo outer        
511 >          select case(cutoffPolicy)
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.
547 >    haveGtypeCutoffMap = .true.
548     end subroutine createGtypeCutoffMap
549  
550     subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
# Line 352 | 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 359 | Line 564 | contains
564       integer, intent(in) :: cutPolicy
565       cutoffPolicy = cutPolicy
566       call createGtypeCutoffMap()
362 –
567     end subroutine setCutoffPolicy
568      
569      
570    subroutine setSimVariables()
571      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
572      SIM_uses_EAM = SimUsesEAM()
369 –    SIM_uses_RF = SimUsesRF()
573      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
574      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
575      SIM_uses_PBC = SimUsesPBC()
# Line 407 | 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 436 | Line 639 | contains
639    end subroutine doReadyCheck
640  
641  
642 <  subroutine init_FF(use_RF_c, thisStat)
642 >  subroutine init_FF(thisESM, thisStat)
643  
644 <    logical, intent(in) :: use_RF_c
442 <
644 >    integer, intent(in) :: thisESM
645      integer, intent(out) :: thisStat  
646      integer :: my_status, nMatches
647      integer, pointer :: MatchList(:) => null()
# Line 448 | Line 650 | contains
650      !! assume things are copacetic, unless they aren't
651      thisStat = 0
652  
653 <    !! Fortran's version of a cast:
452 <    FF_uses_RF = use_RF_c
653 >    electrostaticSummationMethod = thisESM
654  
655      !! init_FF is called *after* all of the atom types have been
656      !! defined in atype_module using the new_atype subroutine.
# Line 480 | 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 548 | Line 749 | contains
749  
750      !! Stress Tensor
751      real( kind = dp), dimension(9) :: tau  
752 <    real ( kind = dp ) :: pot
752 >    real ( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
753      logical ( kind = 2) :: do_pot_c, do_stress_c
754      logical :: do_pot
755      logical :: do_stress
756      logical :: in_switching_region
757   #ifdef IS_MPI
758 <    real( kind = DP ) :: pot_local
758 >    real( kind = DP ), dimension(POT_ARRAY_SIZE) :: pot_local
759      integer :: nAtomsInRow
760      integer :: nAtomsInCol
761      integer :: nprocs
# Line 579 | 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 704 | 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 < InteractionHash(me_i,me_j)%rListsq) then
910 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
911                  if (update_nlist) then
912                     nlist = nlist + 1
913  
# Line 827 | Line 1028 | contains
1028                  endif
1029               end if
1030            enddo
1031 +
1032         enddo outer
1033  
1034         if (update_nlist) then
# Line 906 | 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 937 | Line 1139 | contains
1139                  !! potential and torques:
1140                  call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1141   #ifdef IS_MPI
1142 <                pot_local = pot_local + rfpot
1142 >                pot_local(RF_POT) = pot_local(RF_POT) + rfpot
1143   #else
1144 <                pot = pot + rfpot
1144 >                pot(RF_POT) = pot(RF_POT) + rfpot
1145  
1146   #endif
1147               endif
# Line 951 | Line 1153 | contains
1153   #ifdef IS_MPI
1154  
1155      if (do_pot) then
1156 <       pot = pot + pot_local
1156 >       pot(1:SIZE_POT_ARRAY) = pot(1:SIZE_POT_ARRAY) &
1157 >            + pot_local(1:SIZE_POT_ARRAY)
1158         !! we assume the c code will do the allreduce to get the total potential
1159         !! we could do it right here if we needed to...
1160      endif
# Line 977 | Line 1180 | contains
1180    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1181         eFrame, A, f, t, pot, vpair, fpair)
1182  
1183 <    real( kind = dp ) :: pot, vpair, sw
1183 >    real( kind = dp ) :: vpair, sw
1184 >    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1185      real( kind = dp ), dimension(3) :: fpair
1186      real( kind = dp ), dimension(nLocal)   :: mfact
1187      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 990 | Line 1194 | contains
1194      real ( kind = dp ), intent(inout) :: rijsq
1195      real ( kind = dp )                :: r
1196      real ( kind = dp ), intent(inout) :: d(3)
993 –    real ( kind = dp ) :: ebalance
1197      integer :: me_i, me_j
1198  
1199      integer :: iHash
# Line 1010 | Line 1213 | contains
1213      iHash = InteractionHash(me_i, me_j)
1214  
1215      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1216 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1216 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(LJ_POT), f, do_pot)
1217      endif
1218  
1219      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1220         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1221 <            pot, eFrame, f, t, do_pot)
1221 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1222  
1223 <       if (FF_uses_RF .and. SIM_uses_RF) then
1223 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1224  
1225            ! CHECK ME (RF needs to know about all electrostatic types)
1226            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1028 | Line 1231 | contains
1231  
1232      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1233         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1234 <            pot, A, f, t, do_pot)
1234 >            pot(STICKY_POT), A, f, t, do_pot)
1235      endif
1236  
1237      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1238         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1239 <            pot, A, f, t, do_pot)
1239 >            pot(STICKYPOWER_POT), A, f, t, do_pot)
1240      endif
1241  
1242      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1243         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1244 <            pot, A, f, t, do_pot)
1244 >            pot(GAYBERNE_POT), A, f, t, do_pot)
1245      endif
1246      
1247      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1248   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1249 < !           pot, A, f, t, do_pot)
1249 > !           pot(GAYBERNE_LJ_POT), A, f, t, do_pot)
1250      endif
1251  
1252      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1253 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1253 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(EAM_POT), f, &
1254              do_pot)
1255      endif
1256  
1257      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1258         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1259 <            pot, A, f, t, do_pot)
1259 >            pot(SHAPE_POT), A, f, t, do_pot)
1260      endif
1261  
1262      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1263         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1264 <            pot, A, f, t, do_pot)
1264 >            pot(SHAPE_LJ_POT), A, f, t, do_pot)
1265      endif
1266      
1267    end subroutine do_pair
# Line 1066 | Line 1269 | contains
1269    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1270         do_pot, do_stress, eFrame, A, f, t, pot)
1271  
1272 <    real( kind = dp ) :: pot, sw
1272 >    real( kind = dp ) :: sw
1273 >    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1274      real( kind = dp ), dimension(9,nLocal) :: eFrame
1275      real (kind=dp), dimension(9,nLocal) :: A
1276      real (kind=dp), dimension(3,nLocal) :: f
# Line 1080 | Line 1284 | contains
1284  
1285      integer :: me_i, me_j, iHash
1286  
1287 +    r = sqrt(rijsq)
1288 +
1289   #ifdef IS_MPI  
1290      me_i = atid_row(i)
1291      me_j = atid_col(j)  
# Line 1099 | Line 1305 | contains
1305  
1306    subroutine do_preforce(nlocal,pot)
1307      integer :: nlocal
1308 <    real( kind = dp ) :: pot
1308 >    real( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
1309  
1310      if (FF_uses_EAM .and. SIM_uses_EAM) then
1311 <       call calc_EAM_preforce_Frho(nlocal,pot)
1311 >       call calc_EAM_preforce_Frho(nlocal,pot(EAM_POT))
1312      endif
1313  
1314  
# Line 1297 | Line 1503 | contains
1503  
1504    function FF_RequiresPostpairCalc() result(doesit)
1505      logical :: doesit
1506 <    doesit = FF_uses_RF
1506 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1507    end function FF_RequiresPostpairCalc
1508  
1509   #ifdef PROFILE

Diff Legend

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