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 2309 by chrisfen, Sun Sep 18 20:45:38 2005 UTC vs.
Revision 2398 by chrisfen, Wed Oct 26 23:31:18 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.46 2005-09-18 20:45:38 chrisfen Exp $, $Date: 2005-09-18 20:45:38 $, $Name: not supported by cvs2svn $, $Revision: 1.46 $
48 > !! @version $Id: doForces.F90,v 1.63 2005-10-26 23:31:18 chrisfen Exp $, $Date: 2005-10-26 23:31:18 $, $Name: not supported by cvs2svn $, $Revision: 1.63 $
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_module
62 <  use gb_pair
61 >  use gayberne
62    use shapes
63    use vector_class
64    use eam
# Line 75 | Line 74 | module doForces
74   #include "UseTheForce/fSwitchingFunction.h"
75   #include "UseTheForce/fCutoffPolicy.h"
76   #include "UseTheForce/DarkSide/fInteractionMap.h"
77 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78  
79  
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
# Line 85 | Line 85 | module doForces
85    logical, save :: haveSaneForceField = .false.
86    logical, save :: haveInteractionHash = .false.
87    logical, save :: haveGtypeCutoffMap = .false.
88 +  logical, save :: haveDefaultCutoffs = .false.
89    logical, save :: haveRlist = .false.
90  
91    logical, save :: FF_uses_DirectionalAtoms
# Line 122 | Line 123 | module doForces
123    ! Bit hash to determine pair-pair interactions.
124    integer, dimension(:,:), allocatable :: InteractionHash
125    real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
126 <  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
127 <  integer, dimension(:), allocatable :: groupToGtype
128 <  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
126 >  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
127 >  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
128 >
129 >  integer, dimension(:), allocatable, target :: groupToGtypeRow
130 >  integer, dimension(:), pointer :: groupToGtypeCol => null()
131 >
132 >  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
133 >  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
134    type ::gtypeCutoffs
135       real(kind=dp) :: rcut
136       real(kind=dp) :: rcutsq
# Line 134 | Line 140 | module doForces
140  
141    integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
142    real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
143 +  real(kind=dp),save :: listSkin
144    
145   contains
146  
# Line 177 | Line 184 | contains
184  
185      if (.not. allocated(InteractionHash)) then
186         allocate(InteractionHash(nAtypes,nAtypes))
187 +    else
188 +       deallocate(InteractionHash)
189 +       allocate(InteractionHash(nAtypes,nAtypes))
190      endif
191  
192      if (.not. allocated(atypeMaxCutoff)) then
193         allocate(atypeMaxCutoff(nAtypes))
194 +    else
195 +       deallocate(atypeMaxCutoff)
196 +       allocate(atypeMaxCutoff(nAtypes))
197      endif
198          
199      do i = 1, nAtypes
# Line 257 | Line 270 | contains
270      logical :: GtypeFound
271  
272      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
273 <    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
273 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
274      integer :: nGroupsInRow
275 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
275 >    integer :: nGroupsInCol
276 >    integer :: nGroupTypesRow,nGroupTypesCol
277 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
278      real(kind=dp) :: biggestAtypeCutoff
279  
280      stat = 0
# Line 273 | Line 288 | contains
288      endif
289   #ifdef IS_MPI
290      nGroupsInRow = getNgroupsInRow(plan_group_row)
291 +    nGroupsInCol = getNgroupsInCol(plan_group_col)
292   #endif
293      nAtypes = getSize(atypes)
294   ! Set all of the initial cutoffs to zero.
# Line 288 | Line 304 | contains
304            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
305            
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 >          if (haveDefaultCutoffs) then
308 >             atypeMaxCutoff(i) = defaultRcut
309 >          else
310 >             if (i_is_LJ) then          
311 >                thisRcut = getSigma(i) * 2.5_dp
312 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313 >             endif
314 >             if (i_is_Elect) then
315 >                thisRcut = defaultRcut
316 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317 >             endif
318 >             if (i_is_Sticky) then
319 >                thisRcut = getStickyCut(i)
320 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321 >             endif
322 >             if (i_is_StickyP) then
323 >                thisRcut = getStickyPowerCut(i)
324 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325 >             endif
326 >             if (i_is_GB) then
327 >                thisRcut = getGayBerneCut(i)
328 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329 >             endif
330 >             if (i_is_EAM) then
331 >                thisRcut = getEAMCut(i)
332 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333 >             endif
334 >             if (i_is_Shape) then
335 >                thisRcut = getShapeCut(i)
336 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337 >             endif
338            endif
339            
340 +          
341            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
342               biggestAtypeCutoff = atypeMaxCutoff(i)
343            endif
344 +
345         endif
346      enddo
347    
348 <    nGroupTypes = 0
348 >
349      
350      istart = 1
351 +    jstart = 1
352   #ifdef IS_MPI
353      iend = nGroupsInRow
354 +    jend = nGroupsInCol
355   #else
356      iend = nGroups
357 +    jend = nGroups
358   #endif
359      
360      !! allocate the groupToGtype and gtypeMaxCutoff here.
361 <    if(.not.allocated(groupToGtype)) then
362 <       allocate(groupToGtype(iend))
363 <       allocate(groupMaxCutoff(iend))
364 <       allocate(gtypeMaxCutoff(iend))
365 <       groupMaxCutoff = 0.0_dp
366 <       gtypeMaxCutoff = 0.0_dp
361 >    if(.not.allocated(groupToGtypeRow)) then
362 >     !  allocate(groupToGtype(iend))
363 >       allocate(groupToGtypeRow(iend))
364 >    else
365 >       deallocate(groupToGtypeRow)
366 >       allocate(groupToGtypeRow(iend))
367      endif
368 +    if(.not.allocated(groupMaxCutoffRow)) then
369 +       allocate(groupMaxCutoffRow(iend))
370 +    else
371 +       deallocate(groupMaxCutoffRow)
372 +       allocate(groupMaxCutoffRow(iend))
373 +    end if
374 +
375 +    if(.not.allocated(gtypeMaxCutoffRow)) then
376 +       allocate(gtypeMaxCutoffRow(iend))
377 +    else
378 +       deallocate(gtypeMaxCutoffRow)
379 +       allocate(gtypeMaxCutoffRow(iend))
380 +    endif
381 +
382 +
383 + #ifdef IS_MPI
384 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
385 +    if(.not.associated(groupToGtypeCol)) then
386 +       allocate(groupToGtypeCol(jend))
387 +    else
388 +       deallocate(groupToGtypeCol)
389 +       allocate(groupToGtypeCol(jend))
390 +    end if
391 +
392 +    if(.not.associated(groupToGtypeCol)) then
393 +       allocate(groupToGtypeCol(jend))
394 +    else
395 +       deallocate(groupToGtypeCol)
396 +       allocate(groupToGtypeCol(jend))
397 +    end if
398 +    if(.not.associated(gtypeMaxCutoffCol)) then
399 +       allocate(gtypeMaxCutoffCol(jend))
400 +    else
401 +       deallocate(gtypeMaxCutoffCol)      
402 +       allocate(gtypeMaxCutoffCol(jend))
403 +    end if
404 +
405 +       groupMaxCutoffCol = 0.0_dp
406 +       gtypeMaxCutoffCol = 0.0_dp
407 +
408 + #endif
409 +       groupMaxCutoffRow = 0.0_dp
410 +       gtypeMaxCutoffRow = 0.0_dp
411 +
412 +
413      !! first we do a single loop over the cutoff groups to find the
414      !! largest cutoff for any atypes present in this group.  We also
415      !! create gtypes at this point.
416      
417      tol = 1.0d-6
418 <    
418 >    nGroupTypesRow = 0
419 >
420      do i = istart, iend      
421         n_in_i = groupStartRow(i+1) - groupStartRow(i)
422 <       groupMaxCutoff(i) = 0.0_dp
422 >       groupMaxCutoffRow(i) = 0.0_dp
423         do ia = groupStartRow(i), groupStartRow(i+1)-1
424            atom1 = groupListRow(ia)
425   #ifdef IS_MPI
# Line 356 | Line 427 | contains
427   #else
428            me_i = atid(atom1)
429   #endif          
430 <          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
431 <             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
430 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
431 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
432            endif          
433         enddo
434  
435 <       if (nGroupTypes.eq.0) then
436 <          nGroupTypes = nGroupTypes + 1
437 <          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
438 <          groupToGtype(i) = nGroupTypes
435 >       if (nGroupTypesRow.eq.0) then
436 >          nGroupTypesRow = nGroupTypesRow + 1
437 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
438 >          groupToGtypeRow(i) = nGroupTypesRow
439         else
440            GtypeFound = .false.
441 <          do g = 1, nGroupTypes
442 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
443 <                groupToGtype(i) = g
441 >          do g = 1, nGroupTypesRow
442 >             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
443 >                groupToGtypeRow(i) = g
444                  GtypeFound = .true.
445               endif
446            enddo
447            if (.not.GtypeFound) then            
448 <             nGroupTypes = nGroupTypes + 1
449 <             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
450 <             groupToGtype(i) = nGroupTypes
448 >             nGroupTypesRow = nGroupTypesRow + 1
449 >             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
450 >             groupToGtypeRow(i) = nGroupTypesRow
451 >          endif
452 >       endif
453 >    enddo    
454 >
455 > #ifdef IS_MPI
456 >    do j = jstart, jend      
457 >       n_in_j = groupStartCol(j+1) - groupStartCol(j)
458 >       groupMaxCutoffCol(j) = 0.0_dp
459 >       do ja = groupStartCol(j), groupStartCol(j+1)-1
460 >          atom1 = groupListCol(ja)
461 >
462 >          me_j = atid_col(atom1)
463 >
464 >          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
465 >             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
466 >          endif          
467 >       enddo
468 >
469 >       if (nGroupTypesCol.eq.0) then
470 >          nGroupTypesCol = nGroupTypesCol + 1
471 >          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
472 >          groupToGtypeCol(j) = nGroupTypesCol
473 >       else
474 >          GtypeFound = .false.
475 >          do g = 1, nGroupTypesCol
476 >             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
477 >                groupToGtypeCol(j) = g
478 >                GtypeFound = .true.
479 >             endif
480 >          enddo
481 >          if (.not.GtypeFound) then            
482 >             nGroupTypesCol = nGroupTypesCol + 1
483 >             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
484 >             groupToGtypeCol(j) = nGroupTypesCol
485            endif
486         endif
487      enddo    
488  
489 + #else
490 + ! Set pointers to information we just found
491 +    nGroupTypesCol = nGroupTypesRow
492 +    groupToGtypeCol => groupToGtypeRow
493 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
494 +    groupMaxCutoffCol => groupMaxCutoffRow
495 + #endif
496 +
497 +
498 +
499 +
500 +
501      !! allocate the gtypeCutoffMap here.
502 <    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
502 >    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
503      !! then we do a double loop over all the group TYPES to find the cutoff
504      !! map between groups of two types
505 <    
506 <    do i = 1, nGroupTypes
507 <       do j = 1, nGroupTypes
505 >    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
506 >
507 >    do i = 1, nGroupTypesRow
508 >       do j = 1, nGroupTypesCol
509        
510            select case(cutoffPolicy)
511            case(TRADITIONAL_CUTOFF_POLICY)
512 <             thisRcut = maxval(gtypeMaxCutoff)
512 >             thisRcut = tradRcut
513            case(MIX_CUTOFF_POLICY)
514 <             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
514 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
515            case(MAX_CUTOFF_POLICY)
516 <             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
516 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
517            case default
518               call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
519               return
# Line 403 | Line 521 | contains
521            gtypeCutoffMap(i,j)%rcut = thisRcut
522            gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
523            skin = defaultRlist - defaultRcut
524 +          listSkin = skin ! set neighbor list skin thickness
525            gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
526  
527 +          ! sanity check
528 +
529 +          if (haveDefaultCutoffs) then
530 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
531 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
532 +             endif
533 +          endif
534         enddo
535      enddo
536 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
537 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
538 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
539 + #ifdef IS_MPI
540 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
541 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
542 + #endif
543 +    groupMaxCutoffCol => null()
544 +    gtypeMaxCutoffCol => null()
545      
546      haveGtypeCutoffMap = .true.
547     end subroutine createGtypeCutoffMap
# Line 419 | Line 554 | contains
554       defaultRsw = defRsw
555       defaultRlist = defRlist
556       cutoffPolicy = cutPolicy
557 +
558 +     haveDefaultCutoffs = .true.
559     end subroutine setDefaultCutoffs
560  
561     subroutine setCutoffPolicy(cutPolicy)
# Line 504 | Line 641 | contains
641    subroutine init_FF(thisESM, thisStat)
642  
643      integer, intent(in) :: thisESM
507    real(kind=dp), intent(in) :: dampingAlpha
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
646      integer, pointer :: MatchList(:) => null()
# Line 544 | Line 680 | contains
680  
681      haveSaneForceField = .true.
682  
547    !! check to make sure the reaction field setting makes sense
548
549    if (FF_uses_Dipoles) then
550       if (electrostaticSummationMethod == 3) then
551          dielect = getDielect()
552          call initialize_rf(dielect)
553       endif
554    else
555       if (electrostaticSummationMethod == 3) then
556          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
557          thisStat = -1
558          haveSaneForceField = .false.
559          return
560       endif
561    endif
562
683      if (FF_uses_EAM) then
684         call init_EAM_FF(my_status)
685         if (my_status /= 0) then
# Line 568 | Line 688 | contains
688            haveSaneForceField = .false.
689            return
690         end if
571    endif
572
573    if (FF_uses_GayBerne) then
574       call check_gb_pair_FF(my_status)
575       if (my_status .ne. 0) then
576          thisStat = -1
577          haveSaneForceField = .false.
578          return
579       endif
691      endif
692  
693      if (.not. haveNeighborList) then
# Line 612 | Line 723 | contains
723  
724      !! Stress Tensor
725      real( kind = dp), dimension(9) :: tau  
726 <    real ( kind = dp ) :: pot
726 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
727      logical ( kind = 2) :: do_pot_c, do_stress_c
728      logical :: do_pot
729      logical :: do_stress
730      logical :: in_switching_region
731   #ifdef IS_MPI
732 <    real( kind = DP ) :: pot_local
732 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
733      integer :: nAtomsInRow
734      integer :: nAtomsInCol
735      integer :: nprocs
# Line 633 | Line 744 | contains
744      integer :: nlist
745      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
746      real( kind = DP ) :: sw, dswdr, swderiv, mf
747 +    real( kind = DP ) :: rVal
748      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
749      real(kind=dp) :: rfpot, mu_i, virial
750      integer :: me_i, me_j, n_in_i, n_in_j
# Line 643 | Line 755 | contains
755      integer :: propPack_i, propPack_j
756      integer :: loopStart, loopEnd, loop
757      integer :: iHash
758 <    real(kind=dp) :: listSkin = 1.0  
758 >    integer :: i1
759 >    logical :: indirect_only
760 >  
761  
762      !! initialize local variables  
763  
# Line 768 | Line 882 | contains
882               me_j = atid(j)
883               call get_interatomic_vector(q_group(:,i), &
884                    q_group(:,j), d_grp, rgrpsq)
885 < #endif
885 > #endif      
886  
887 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
887 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
888                  if (update_nlist) then
889                     nlist = nlist + 1
890  
# Line 809 | Line 923 | contains
923  
924                        atom2 = groupListCol(jb)
925  
926 <                      if (skipThisPair(atom1, atom2)) cycle inner
926 >                      indirect_only = .false.
927 >    
928 >                      if (skipThisPair(atom1, atom2)) then
929 >                         if (electrostaticSummationMethod.ne.REACTION_FIELD) then
930 >                            cycle inner
931 >                         else
932 >                            indirect_only = .true.
933 >                         endif
934 >                      endif
935 >    
936  
937                        if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
938                           d_atm(1:3) = d_grp(1:3)
# Line 837 | Line 960 | contains
960                        else
961   #ifdef IS_MPI                      
962                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
963 <                              do_pot, &
964 <                              eFrame, A, f, t, pot_local, vpair, fpair)
963 >                              do_pot, eFrame, A, f, t, pot_local, vpair, &
964 >                              fpair, d_grp, rgrp, indirect_only)
965   #else
966                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
967 <                              do_pot,  &
968 <                              eFrame, A, f, t, pot, vpair, fpair)
967 >                              do_pot, eFrame, A, f, t, pot, vpair, fpair, &
968 >                              d_grp, rgrp, indirect_only)
969   #endif
970  
971                           vij = vij + vpair
# Line 891 | Line 1014 | contains
1014                  endif
1015               end if
1016            enddo
1017 +
1018         enddo outer
1019  
1020         if (update_nlist) then
# Line 950 | Line 1074 | contains
1074  
1075      if (do_pot) then
1076         ! scatter/gather pot_row into the members of my column
1077 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1078 <
1077 >       do i = 1,LR_POT_TYPES
1078 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1079 >       end do
1080         ! scatter/gather pot_local into all other procs
1081         ! add resultant to get total pot
1082         do i = 1, nlocal
1083 <          pot_local = pot_local + pot_Temp(i)
1083 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1084 >               + pot_Temp(1:LR_POT_TYPES,i)
1085         enddo
1086  
1087         pot_Temp = 0.0_DP
1088 <
1089 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1088 >       do i = 1,LR_POT_TYPES
1089 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1090 >       end do
1091         do i = 1, nlocal
1092 <          pot_local = pot_local + pot_Temp(i)
1092 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1093 >               + pot_Temp(1:LR_POT_TYPES,i)
1094         enddo
1095  
1096      endif
1097   #endif
1098  
1099 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1099 >    if (SIM_requires_postpair_calc) then
1100 >       do i = 1, nlocal            
1101 >          
1102 >          ! we loop only over the local atoms, so we don't need row and column
1103 >          ! lookups for the types
1104 >          
1105 >          me_i = atid(i)
1106 >          
1107 >          ! is the atom electrostatic?  See if it would have an
1108 >          ! electrostatic interaction with itself
1109 >          iHash = InteractionHash(me_i,me_i)
1110  
1111 <       if (electrostaticSummationMethod == 3) then
974 <
1111 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1112   #ifdef IS_MPI
1113 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1114 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
978 <          do i = 1,nlocal
979 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
980 <          end do
981 < #endif
982 <
983 <          do i = 1, nLocal
984 <
985 <             rfpot = 0.0_DP
986 < #ifdef IS_MPI
987 <             me_i = atid_row(i)
1113 >             call rf_self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1114 >                  t, do_pot)
1115   #else
1116 <             me_i = atid(i)
1116 >             call rf_self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1117 >                  t, do_pot)
1118   #endif
1119 <             iHash = InteractionHash(me_i,me_j)
1120 <            
1121 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1119 >          endif
1120 >  
1121 >          ! loop over the excludes to accumulate any additional RF components
1122  
1123 <                mu_i = getDipoleMoment(me_i)
1124 <
1125 <                !! The reaction field needs to include a self contribution
1126 <                !! to the field:
1127 <                call accumulate_self_rf(i, mu_i, eFrame)
1128 <                !! Get the reaction field contribution to the
1129 <                !! potential and torques:
1130 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1123 >          do i1 = 1, nSkipsForAtom(i)
1124 >             j = skipsForAtom(i, i1)
1125 >            
1126 >             ! prevent overcounting of the skips
1127 >             if (i.lt.j) then
1128 >                call get_interatomic_vector(q(:,i), &
1129 >                     q(:,j), d_atm, ratmsq)
1130 >                rVal = dsqrt(ratmsq)
1131 >                call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1132 >                     in_switching_region)
1133   #ifdef IS_MPI
1134 <                pot_local = pot_local + rfpot
1134 >                call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, vpair, &
1135 >                     pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1136   #else
1137 <                pot = pot + rfpot
1138 <
1137 >                call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, vpair, &
1138 >                     pot(ELECTROSTATIC_POT), f, t, do_pot)
1139   #endif
1140               endif
1141            enddo
1142 <       endif
1142 >       enddo      
1143      endif
1144 <
1014 <
1144 >    
1145   #ifdef IS_MPI
1146 <
1146 >    
1147      if (do_pot) then
1148 <       pot = pot + pot_local
1149 <       !! we assume the c code will do the allreduce to get the total potential
1020 <       !! we could do it right here if we needed to...
1148 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1149 >            mpi_comm_world,mpi_err)            
1150      endif
1151 <
1151 >    
1152      if (do_stress) then
1153         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1154              mpi_comm_world,mpi_err)
1155         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1156              mpi_comm_world,mpi_err)
1157      endif
1158 <
1158 >    
1159   #else
1160 <
1160 >    
1161      if (do_stress) then
1162         tau = tau_Temp
1163         virial = virial_Temp
1164      endif
1165 <
1165 >    
1166   #endif
1167 <
1167 >    
1168    end subroutine do_force_loop
1169  
1170    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1171 <       eFrame, A, f, t, pot, vpair, fpair)
1171 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, indirect_only)
1172  
1173 <    real( kind = dp ) :: pot, vpair, sw
1173 >    real( kind = dp ) :: vpair, sw
1174 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1175      real( kind = dp ), dimension(3) :: fpair
1176      real( kind = dp ), dimension(nLocal)   :: mfact
1177      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1050 | Line 1180 | contains
1180      real( kind = dp ), dimension(3,nLocal) :: t
1181  
1182      logical, intent(inout) :: do_pot
1183 +    logical, intent(inout) :: indirect_only
1184      integer, intent(in) :: i, j
1185      real ( kind = dp ), intent(inout) :: rijsq
1186 <    real ( kind = dp )                :: r
1186 >    real ( kind = dp ), intent(inout) :: r_grp
1187      real ( kind = dp ), intent(inout) :: d(3)
1188 +    real ( kind = dp ), intent(inout) :: d_grp(3)
1189 +    real ( kind = dp ) :: r
1190      integer :: me_i, me_j
1191  
1192      integer :: iHash
# Line 1072 | Line 1205 | contains
1205  
1206      iHash = InteractionHash(me_i, me_j)
1207  
1208 <    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1209 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1210 <    endif
1211 <
1079 <    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1080 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1081 <            pot, eFrame, f, t, do_pot)
1082 <
1083 <       if (electrostaticSummationMethod == 3) then
1084 <
1085 <          ! CHECK ME (RF needs to know about all electrostatic types)
1086 <          call accumulate_rf(i, j, r, eFrame, sw)
1087 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1208 >    if (indirect_only) then
1209 >       if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1210 >          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1211 >               pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot, indirect_only)
1212         endif
1213 <
1214 <    endif
1215 <
1216 <    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1217 <       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1218 <            pot, A, f, t, do_pot)
1095 <    endif
1096 <
1097 <    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1098 <       call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1099 <            pot, A, f, t, do_pot)
1100 <    endif
1101 <
1102 <    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1103 <       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1104 <            pot, A, f, t, do_pot)
1105 <    endif
1106 <    
1107 <    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1108 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1109 < !           pot, A, f, t, do_pot)
1110 <    endif
1111 <
1112 <    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1113 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1114 <            do_pot)
1115 <    endif
1213 >    else
1214 >          
1215 >       if ( iand(iHash, LJ_PAIR).ne.0 ) then
1216 >          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1217 >               pot(VDW_POT), f, do_pot)
1218 >       endif
1219  
1220 <    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1221 <       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1222 <            pot, A, f, t, do_pot)
1223 <    endif
1224 <
1225 <    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1226 <       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1227 <            pot, A, f, t, do_pot)
1220 >       if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1221 >          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1222 >               pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot, indirect_only)
1223 >       endif
1224 >      
1225 >       if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1226 >          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1227 >               pot(HB_POT), A, f, t, do_pot)
1228 >       endif
1229 >      
1230 >       if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1231 >          call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1232 >               pot(HB_POT), A, f, t, do_pot)
1233 >       endif
1234 >      
1235 >       if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1236 >          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1237 >               pot(VDW_POT), A, f, t, do_pot)
1238 >       endif
1239 >      
1240 >       if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1241 >          call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1242 >               pot(VDW_POT), A, f, t, do_pot)
1243 >       endif
1244 >      
1245 >       if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1246 >          call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1247 >               pot(METALLIC_POT), f, do_pot)
1248 >       endif
1249 >      
1250 >       if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1251 >          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1252 >               pot(VDW_POT), A, f, t, do_pot)
1253 >       endif
1254 >      
1255 >       if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1256 >          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1257 >               pot(VDW_POT), A, f, t, do_pot)
1258 >       endif
1259      endif
1260      
1261    end subroutine do_pair
# Line 1129 | Line 1263 | contains
1263    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1264         do_pot, do_stress, eFrame, A, f, t, pot)
1265  
1266 <    real( kind = dp ) :: pot, sw
1266 >    real( kind = dp ) :: sw
1267 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1268      real( kind = dp ), dimension(9,nLocal) :: eFrame
1269      real (kind=dp), dimension(9,nLocal) :: A
1270      real (kind=dp), dimension(3,nLocal) :: f
# Line 1164 | Line 1299 | contains
1299  
1300    subroutine do_preforce(nlocal,pot)
1301      integer :: nlocal
1302 <    real( kind = dp ) :: pot
1302 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1303  
1304      if (FF_uses_EAM .and. SIM_uses_EAM) then
1305 <       call calc_EAM_preforce_Frho(nlocal,pot)
1305 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1306      endif
1307  
1308  
# Line 1360 | Line 1495 | contains
1495      doesit = FF_uses_EAM
1496    end function FF_RequiresPrepairCalc
1497  
1363  function FF_RequiresPostpairCalc() result(doesit)
1364    logical :: doesit
1365    if (electrostaticSummationMethod == 3) doesit = .true.
1366  end function FF_RequiresPostpairCalc
1367
1498   #ifdef PROFILE
1499    function getforcetime() result(totalforcetime)
1500      real(kind=dp) :: totalforcetime

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines