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 2288 by chuckv, Wed Sep 7 22:44:48 2005 UTC vs.
Revision 2405 by chrisfen, Tue Nov 1 19:24:57 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.40 2005-09-07 22:44:48 chuckv Exp $, $Date: 2005-09-07 22:44:48 $, $Name: not supported by cvs2svn $, $Revision: 1.40 $
48 > !! @version $Id: doForces.F90,v 1.65 2005-11-01 19:24:54 chrisfen Exp $, $Date: 2005-11-01 19:24:54 $, $Name: not supported by cvs2svn $, $Revision: 1.65 $
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
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
92    logical, save :: FF_uses_Dipoles
93    logical, save :: FF_uses_GayBerne
94    logical, save :: FF_uses_EAM
94 –  logical, save :: FF_uses_RF
95  
96    logical, save :: SIM_uses_DirectionalAtoms
97    logical, save :: SIM_uses_EAM
98 –  logical, save :: SIM_uses_RF
98    logical, save :: SIM_requires_postpair_calc
99    logical, save :: SIM_requires_prepair_calc
100    logical, save :: SIM_uses_PBC
101  
102 <  integer, save :: corrMethod
102 >  integer, save :: electrostaticSummationMethod
103  
104    public :: init_FF
105    public :: setDefaultCutoffs
# Line 124 | 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 136 | 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 179 | 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 259 | 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 275 | 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 <    
294 > ! Set all of the initial cutoffs to zero.
295 >    atypeMaxCutoff = 0.0_dp
296      do i = 1, nAtypes
297         if (SimHasAtype(i)) then    
298            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 288 | Line 303 | contains
303            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
304            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
305            
306 <          atypeMaxCutoff(i) = 0.0_dp
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
306 >
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))
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 355 | 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 402 | 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 <    
412 <  end subroutine createGtypeCutoffMap
413 <  
414 <  subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
415 <    real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
416 <    integer, intent(in) :: cutPolicy
417 <    
418 <    defaultRcut = defRcut
419 <    defaultRsw = defRsw
420 <    defaultRlist = defRlist
421 <    cutoffPolicy = cutPolicy
422 <  end subroutine setDefaultCutoffs
423 <  
424 <  subroutine setCutoffPolicy(cutPolicy)
547 >   end subroutine createGtypeCutoffMap
548  
549 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
550 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
551       integer, intent(in) :: cutPolicy
552 +
553 +     defaultRcut = defRcut
554 +     defaultRsw = defRsw
555 +     defaultRlist = defRlist
556       cutoffPolicy = cutPolicy
428 –     call createGtypeCutoffMap()
557  
558 +     haveDefaultCutoffs = .true.
559 +   end subroutine setDefaultCutoffs
560 +
561 +   subroutine setCutoffPolicy(cutPolicy)
562 +
563 +     integer, intent(in) :: cutPolicy
564 +     cutoffPolicy = cutPolicy
565 +     call createGtypeCutoffMap()
566     end subroutine setCutoffPolicy
567      
568      
569    subroutine setSimVariables()
570      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
571      SIM_uses_EAM = SimUsesEAM()
436 –    SIM_uses_RF = SimUsesRF()
572      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574      SIM_uses_PBC = SimUsesPBC()
# Line 503 | Line 638 | contains
638    end subroutine doReadyCheck
639  
640  
641 <  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
641 >  subroutine init_FF(thisESM, thisStat)
642  
643 <    logical, intent(in) :: use_RF
509 <    logical, intent(in) :: use_UW
510 <    logical, intent(in) :: use_DW
643 >    integer, intent(in) :: thisESM
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
513 –    integer :: corrMethod
646      integer, pointer :: MatchList(:) => null()
515 –    real(kind=dp) :: rcut, rrf, rt, dielect
647  
648      !! assume things are copacetic, unless they aren't
649      thisStat = 0
650  
651 <    !! Fortran's version of a cast:
521 <    FF_uses_RF = use_RF
651 >    electrostaticSummationMethod = thisESM
652  
523 –    !! set the electrostatic correction method
524 –    if (use_UW) then
525 –       corrMethod = 1
526 –    elseif (use_DW) then
527 –       corrMethod = 2
528 –    else
529 –       corrMethod = 0
530 –    endif
531 –    
653      !! init_FF is called *after* all of the atom types have been
654      !! defined in atype_module using the new_atype subroutine.
655      !!
# Line 557 | Line 678 | contains
678  
679  
680      haveSaneForceField = .true.
560 –
561 –    !! check to make sure the FF_uses_RF setting makes sense
681  
563 –    if (FF_uses_Dipoles) then
564 –       if (FF_uses_RF) then
565 –          dielect = getDielect()
566 –          call initialize_rf(dielect)
567 –       endif
568 –    else
569 –       if (FF_uses_RF) then          
570 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
571 –          thisStat = -1
572 –          haveSaneForceField = .false.
573 –          return
574 –       endif
575 –    endif
576 –
682      if (FF_uses_EAM) then
683         call init_EAM_FF(my_status)
684         if (my_status /= 0) then
# Line 582 | Line 687 | contains
687            haveSaneForceField = .false.
688            return
689         end if
585 –    endif
586 –
587 –    if (FF_uses_GayBerne) then
588 –       call check_gb_pair_FF(my_status)
589 –       if (my_status .ne. 0) then
590 –          thisStat = -1
591 –          haveSaneForceField = .false.
592 –          return
593 –       endif
690      endif
691  
692      if (.not. haveNeighborList) then
# Line 626 | Line 722 | contains
722  
723      !! Stress Tensor
724      real( kind = dp), dimension(9) :: tau  
725 <    real ( kind = dp ) :: pot
725 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
726      logical ( kind = 2) :: do_pot_c, do_stress_c
727      logical :: do_pot
728      logical :: do_stress
729      logical :: in_switching_region
730   #ifdef IS_MPI
731 <    real( kind = DP ) :: pot_local
731 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
732      integer :: nAtomsInRow
733      integer :: nAtomsInCol
734      integer :: nprocs
# Line 647 | Line 743 | contains
743      integer :: nlist
744      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
745      real( kind = DP ) :: sw, dswdr, swderiv, mf
746 +    real( kind = DP ) :: rVal
747      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
748      real(kind=dp) :: rfpot, mu_i, virial
749      integer :: me_i, me_j, n_in_i, n_in_j
# Line 657 | Line 754 | contains
754      integer :: propPack_i, propPack_j
755      integer :: loopStart, loopEnd, loop
756      integer :: iHash
757 <    real(kind=dp) :: listSkin = 1.0  
757 >    integer :: i1
758 >  
759  
760      !! initialize local variables  
761  
# Line 782 | Line 880 | contains
880               me_j = atid(j)
881               call get_interatomic_vector(q_group(:,i), &
882                    q_group(:,j), d_grp, rgrpsq)
883 < #endif
883 > #endif      
884  
885 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
885 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
886                  if (update_nlist) then
887                     nlist = nlist + 1
888  
# Line 822 | Line 920 | contains
920                     inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
921  
922                        atom2 = groupListCol(jb)
923 <
924 <                      if (skipThisPair(atom1, atom2)) cycle inner
925 <
923 >    
924 >                      if (skipThisPair(atom1, atom2))  cycle inner
925 >                      
926                        if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
927                           d_atm(1:3) = d_grp(1:3)
928                           ratmsq = rgrpsq
# Line 851 | Line 949 | contains
949                        else
950   #ifdef IS_MPI                      
951                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
952 <                              do_pot, &
953 <                              eFrame, A, f, t, pot_local, vpair, fpair)
952 >                              do_pot, eFrame, A, f, t, pot_local, vpair, &
953 >                              fpair, d_grp, rgrp)
954   #else
955                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
956 <                              do_pot,  &
957 <                              eFrame, A, f, t, pot, vpair, fpair)
956 >                              do_pot, eFrame, A, f, t, pot, vpair, fpair, &
957 >                              d_grp, rgrp)
958   #endif
959  
960                           vij = vij + vpair
# Line 905 | Line 1003 | contains
1003                  endif
1004               end if
1005            enddo
1006 +
1007         enddo outer
1008  
1009         if (update_nlist) then
# Line 964 | Line 1063 | contains
1063  
1064      if (do_pot) then
1065         ! scatter/gather pot_row into the members of my column
1066 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1067 <
1066 >       do i = 1,LR_POT_TYPES
1067 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1068 >       end do
1069         ! scatter/gather pot_local into all other procs
1070         ! add resultant to get total pot
1071         do i = 1, nlocal
1072 <          pot_local = pot_local + pot_Temp(i)
1072 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1073 >               + pot_Temp(1:LR_POT_TYPES,i)
1074         enddo
1075  
1076         pot_Temp = 0.0_DP
1077 <
1078 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1077 >       do i = 1,LR_POT_TYPES
1078 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1079 >       end do
1080         do i = 1, nlocal
1081 <          pot_local = pot_local + pot_Temp(i)
1081 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1082 >               + pot_Temp(1:LR_POT_TYPES,i)
1083         enddo
1084  
1085      endif
1086   #endif
1087  
1088 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1088 >    if (SIM_requires_postpair_calc) then
1089 >       do i = 1, nlocal            
1090 >          
1091 >          ! we loop only over the local atoms, so we don't need row and column
1092 >          ! lookups for the types
1093 >          
1094 >          me_i = atid(i)
1095 >          
1096 >          ! is the atom electrostatic?  See if it would have an
1097 >          ! electrostatic interaction with itself
1098 >          iHash = InteractionHash(me_i,me_i)
1099  
1100 <       if (FF_uses_RF .and. SIM_uses_RF) then
988 <
1100 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1101   #ifdef IS_MPI
1102 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1103 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
992 <          do i = 1,nlocal
993 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
994 <          end do
995 < #endif
996 <
997 <          do i = 1, nLocal
998 <
999 <             rfpot = 0.0_DP
1000 < #ifdef IS_MPI
1001 <             me_i = atid_row(i)
1102 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1103 >                  t, do_pot)
1104   #else
1105 <             me_i = atid(i)
1105 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1106 >                  t, do_pot)
1107   #endif
1108 <             iHash = InteractionHash(me_i,me_j)
1108 >          endif
1109 >  
1110 >          
1111 > !!$          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1112              
1113 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1114 <
1115 <                mu_i = getDipoleMoment(me_i)
1116 <
1117 <                !! The reaction field needs to include a self contribution
1118 <                !! to the field:
1119 <                call accumulate_self_rf(i, mu_i, eFrame)
1120 <                !! Get the reaction field contribution to the
1121 <                !! potential and torques:
1122 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1113 >             ! loop over the excludes to accumulate RF stuff we've
1114 >             ! left out of the normal pair loop
1115 >            
1116 >             do i1 = 1, nSkipsForAtom(i)
1117 >                j = skipsForAtom(i, i1)
1118 >                
1119 >                ! prevent overcounting of the skips
1120 >                if (i.lt.j) then
1121 >                   call get_interatomic_vector(q(:,i), &
1122 >                        q(:,j), d_atm, ratmsq)
1123 >                   rVal = dsqrt(ratmsq)
1124 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1125 >                        in_switching_region)
1126   #ifdef IS_MPI
1127 <                pot_local = pot_local + rfpot
1127 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1128 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1129   #else
1130 <                pot = pot + rfpot
1131 <
1130 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1131 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1132   #endif
1133 <             endif
1134 <          enddo
1135 <       endif
1133 >                endif
1134 >             enddo
1135 > !!$          endif
1136 >       enddo
1137      endif
1138 <
1028 <
1138 >    
1139   #ifdef IS_MPI
1140 <
1140 >    
1141      if (do_pot) then
1142 <       pot = pot + pot_local
1143 <       !! we assume the c code will do the allreduce to get the total potential
1034 <       !! we could do it right here if we needed to...
1142 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1143 >            mpi_comm_world,mpi_err)            
1144      endif
1145 <
1145 >    
1146      if (do_stress) then
1147         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1148              mpi_comm_world,mpi_err)
1149         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1150              mpi_comm_world,mpi_err)
1151      endif
1152 <
1152 >    
1153   #else
1154 <
1154 >    
1155      if (do_stress) then
1156         tau = tau_Temp
1157         virial = virial_Temp
1158      endif
1159 <
1159 >    
1160   #endif
1161 <
1161 >    
1162    end subroutine do_force_loop
1163  
1164    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1165 <       eFrame, A, f, t, pot, vpair, fpair)
1165 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1166  
1167 <    real( kind = dp ) :: pot, vpair, sw
1167 >    real( kind = dp ) :: vpair, sw
1168 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1169      real( kind = dp ), dimension(3) :: fpair
1170      real( kind = dp ), dimension(nLocal)   :: mfact
1171      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1066 | Line 1176 | contains
1176      logical, intent(inout) :: do_pot
1177      integer, intent(in) :: i, j
1178      real ( kind = dp ), intent(inout) :: rijsq
1179 <    real ( kind = dp )                :: r
1179 >    real ( kind = dp ), intent(inout) :: r_grp
1180      real ( kind = dp ), intent(inout) :: d(3)
1181 <    real ( kind = dp ) :: ebalance
1181 >    real ( kind = dp ), intent(inout) :: d_grp(3)
1182 >    real ( kind = dp ) :: r
1183      integer :: me_i, me_j
1184  
1185      integer :: iHash
# Line 1086 | Line 1197 | contains
1197   #endif
1198  
1199      iHash = InteractionHash(me_i, me_j)
1200 <
1200 >    
1201      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1202 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1202 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1203 >            pot(VDW_POT), f, do_pot)
1204      endif
1205 <
1205 >    
1206      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1207         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1208 <            pot, eFrame, f, t, do_pot, corrMethod)
1097 <
1098 <       if (FF_uses_RF .and. SIM_uses_RF) then
1099 <
1100 <          ! CHECK ME (RF needs to know about all electrostatic types)
1101 <          call accumulate_rf(i, j, r, eFrame, sw)
1102 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1103 <       endif
1104 <
1208 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1209      endif
1210 <
1210 >    
1211      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1212         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213 <            pot, A, f, t, do_pot)
1213 >            pot(HB_POT), A, f, t, do_pot)
1214      endif
1215 <
1215 >    
1216      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1217         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1218 <            pot, A, f, t, do_pot)
1218 >            pot(HB_POT), A, f, t, do_pot)
1219      endif
1220 <
1220 >    
1221      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1222         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1223 <            pot, A, f, t, do_pot)
1223 >            pot(VDW_POT), A, f, t, do_pot)
1224      endif
1225      
1226      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1227 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1228 < !           pot, A, f, t, do_pot)
1227 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1228 >            pot(VDW_POT), A, f, t, do_pot)
1229      endif
1230 <
1230 >    
1231      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1232 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1233 <            do_pot)
1232 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1233 >            pot(METALLIC_POT), f, do_pot)
1234      endif
1235 <
1235 >    
1236      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1237         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1238 <            pot, A, f, t, do_pot)
1238 >            pot(VDW_POT), A, f, t, do_pot)
1239      endif
1240 <
1240 >    
1241      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1242         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1243 <            pot, A, f, t, do_pot)
1243 >            pot(VDW_POT), A, f, t, do_pot)
1244      endif
1245 <    
1245 >    
1246    end subroutine do_pair
1247  
1248    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1249         do_pot, do_stress, eFrame, A, f, t, pot)
1250  
1251 <    real( kind = dp ) :: pot, sw
1251 >    real( kind = dp ) :: sw
1252 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1253      real( kind = dp ), dimension(9,nLocal) :: eFrame
1254      real (kind=dp), dimension(9,nLocal) :: A
1255      real (kind=dp), dimension(3,nLocal) :: f
# Line 1158 | Line 1263 | contains
1263  
1264      integer :: me_i, me_j, iHash
1265  
1266 +    r = sqrt(rijsq)
1267 +
1268   #ifdef IS_MPI  
1269      me_i = atid_row(i)
1270      me_j = atid_col(j)  
# Line 1177 | Line 1284 | contains
1284  
1285    subroutine do_preforce(nlocal,pot)
1286      integer :: nlocal
1287 <    real( kind = dp ) :: pot
1287 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1288  
1289      if (FF_uses_EAM .and. SIM_uses_EAM) then
1290 <       call calc_EAM_preforce_Frho(nlocal,pot)
1290 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1291      endif
1292  
1293  
# Line 1266 | Line 1373 | contains
1373      pot_Col = 0.0_dp
1374      pot_Temp = 0.0_dp
1375  
1269 –    rf_Row = 0.0_dp
1270 –    rf_Col = 0.0_dp
1271 –    rf_Temp = 0.0_dp
1272 –
1376   #endif
1377  
1378      if (FF_uses_EAM .and. SIM_uses_EAM) then
# Line 1373 | Line 1476 | contains
1476      doesit = FF_uses_EAM
1477    end function FF_RequiresPrepairCalc
1478  
1376 –  function FF_RequiresPostpairCalc() result(doesit)
1377 –    logical :: doesit
1378 –    doesit = FF_uses_RF
1379 –  end function FF_RequiresPostpairCalc
1380 –
1479   #ifdef PROFILE
1480    function getforcetime() result(totalforcetime)
1481      real(kind=dp) :: totalforcetime

Diff Legend

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