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 2270 by gezelter, Tue Aug 9 22:33:37 2005 UTC vs.
Revision 2402 by chrisfen, Tue Nov 1 19:09:30 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.28 2005-08-09 22:33:37 gezelter Exp $, $Date: 2005-08-09 22:33:37 $, $Name: not supported by cvs2svn $, $Revision: 1.28 $
48 > !! @version $Id: doForces.F90,v 1.64 2005-11-01 19:09:23 chrisfen Exp $, $Date: 2005-11-01 19:09:23 $, $Name: not supported by cvs2svn $, $Revision: 1.64 $
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 73 | Line 72 | module doForces
72  
73   #define __FORTRAN90
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
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82  
# Line 83 | 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
91 –  logical, save :: FF_uses_RF
95  
96    logical, save :: SIM_uses_DirectionalAtoms
97    logical, save :: SIM_uses_EAM
95 –  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 :: electrostaticSummationMethod
103 +
104    public :: init_FF
105 +  public :: setDefaultCutoffs
106    public :: do_force_loop
107    public :: createInteractionHash
108    public :: createGtypeCutoffMap
109 +  public :: getStickyCut
110 +  public :: getStickyPowerCut
111 +  public :: getGayBerneCut
112 +  public :: getEAMCut
113 +  public :: getShapeCut
114  
115   #ifdef PROFILE
116    public :: getforcetime
# Line 113 | 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
137       real(kind=dp) :: rlistsq
138    end type gtypeCutoffs
139    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
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 146 | Line 165 | contains
165      logical :: j_is_GB
166      logical :: j_is_EAM
167      logical :: j_is_Shape
168 <    
168 >    real(kind=dp) :: myRcut
169 >
170      status = 0  
171  
172      if (.not. associated(atypes)) then
# Line 164 | 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 178 | Line 204 | contains
204         call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
205         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
206         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
181 –
182 –       if (i_is_LJ) then
183 –          thisCut = getDefaultLJCutoff(i)
184 –          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
185 –       endif
207  
187 –
188 –
208         do j = i, nAtypes
209  
210            iHash = 0
# Line 238 | Line 257 | contains
257      haveInteractionHash = .true.
258    end subroutine createInteractionHash
259  
260 <  subroutine createGtypeCutoffMap(defaultRcut, defaultSkinThickness, stat)
260 >  subroutine createGtypeCutoffMap(stat)
261  
262 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
263 <    integer, intent(out) :: stat
262 >    integer, intent(out), optional :: stat
263 >    logical :: i_is_LJ
264 >    logical :: i_is_Elect
265 >    logical :: i_is_Sticky
266 >    logical :: i_is_StickyP
267 >    logical :: i_is_GB
268 >    logical :: i_is_EAM
269 >    logical :: i_is_Shape
270 >    logical :: GtypeFound
271  
272 <    integer :: myStatus, nAtypes
272 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
273 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
274 >    integer :: nGroupsInRow
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
281      if (.not. haveInteractionHash) then
# Line 254 | Line 286 | contains
286            return
287         endif
288      endif
289 <
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 <      
298 <       atypeMaxCutoff(i) =
297 >       if (SimHasAtype(i)) then    
298 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
299 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
300 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
301 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
302 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
303 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
304 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
305 >          
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 +
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(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 <     haveGtypeCutoffMap = .true.
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 >    nGroupTypesRow = 0
419 >
420 >    do i = istart, iend      
421 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
422 >       groupMaxCutoffRow(i) = 0.0_dp
423 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
424 >          atom1 = groupListRow(ia)
425 > #ifdef IS_MPI
426 >          me_i = atid_row(atom1)
427 > #else
428 >          me_i = atid(atom1)
429 > #endif          
430 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
431 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
432 >          endif          
433 >       enddo
434 >
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, 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 >             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(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 >    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 = tradRcut
513 >          case(MIX_CUTOFF_POLICY)
514 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
515 >          case(MAX_CUTOFF_POLICY)
516 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
517 >          case default
518 >             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
519 >             return
520 >          end select
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
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
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()
274 –    SIM_uses_RF = SimUsesRF()
572      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574      SIM_uses_PBC = SimUsesPBC()
# Line 312 | Line 609 | contains
609         call setSimVariables()
610      endif
611  
612 <    if (.not. haveRlist) then
613 <       write(default_error, *) 'rList has not been set in doForces!'
614 <       error = -1
615 <       return
616 <    endif
612 >  !  if (.not. haveRlist) then
613 >  !     write(default_error, *) 'rList has not been set in doForces!'
614 >  !     error = -1
615 >  !     return
616 >  !  endif
617  
618      if (.not. haveNeighborList) then
619         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 341 | Line 638 | contains
638    end subroutine doReadyCheck
639  
640  
641 <  subroutine init_FF(use_RF_c, thisStat)
641 >  subroutine init_FF(thisESM, thisStat)
642  
643 <    logical, intent(in) :: use_RF_c
347 <
643 >    integer, intent(in) :: thisESM
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
646      integer, pointer :: MatchList(:) => null()
351 –    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:
357 <    FF_uses_RF = use_RF_c
651 >    electrostaticSummationMethod = thisESM
652  
653      !! init_FF is called *after* all of the atom types have been
654      !! defined in atype_module using the new_atype subroutine.
# Line 385 | Line 679 | contains
679  
680      haveSaneForceField = .true.
681  
388 –    !! check to make sure the FF_uses_RF setting makes sense
389 –
390 –    if (FF_uses_Dipoles) then
391 –       if (FF_uses_RF) then
392 –          dielect = getDielect()
393 –          call initialize_rf(dielect)
394 –       endif
395 –    else
396 –       if (FF_uses_RF) then          
397 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
398 –          thisStat = -1
399 –          haveSaneForceField = .false.
400 –          return
401 –       endif
402 –    endif
403 –
682      if (FF_uses_EAM) then
683         call init_EAM_FF(my_status)
684         if (my_status /= 0) then
# Line 411 | Line 689 | contains
689         end if
690      endif
691  
414 –    if (FF_uses_GayBerne) then
415 –       call check_gb_pair_FF(my_status)
416 –       if (my_status .ne. 0) then
417 –          thisStat = -1
418 –          haveSaneForceField = .false.
419 –          return
420 –       endif
421 –    endif
422 –
692      if (.not. haveNeighborList) then
693         !! Create neighbor lists
694         call expandNeighborList(nLocal, my_status)
# Line 453 | 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 474 | 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 484 | 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 575 | Line 846 | contains
846   #endif
847         outer: do i = istart, iend
848  
578 – #ifdef IS_MPI
579 –             me_i = atid_row(i)
580 – #else
581 –             me_i = atid(i)
582 – #endif
583 –
849            if (update_nlist) point(i) = nlist + 1
850  
851            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 615 | 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 < InteractionHash(me_i,me_j)%rListsq) then
885 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
886                  if (update_nlist) then
887                     nlist = nlist + 1
888  
# Line 655 | 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 684 | 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 738 | Line 1003 | contains
1003                  endif
1004               end if
1005            enddo
1006 +
1007         enddo outer
1008  
1009         if (update_nlist) then
# Line 797 | 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
821 <
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)
825 <          do i = 1,nlocal
826 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
827 <          end do
828 < #endif
829 <
830 <          do i = 1, nLocal
831 <
832 <             rfpot = 0.0_DP
833 < #ifdef IS_MPI
834 <             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 <
861 <
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
867 <       !! 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 > !!$  subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1167 > !!$       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, felec)
1168  
1169 <    real( kind = dp ) :: pot, vpair, sw
1169 >    real( kind = dp ) :: vpair, sw
1170 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1171      real( kind = dp ), dimension(3) :: fpair
1172 +    real( kind = dp ), dimension(3) :: felec
1173      real( kind = dp ), dimension(nLocal)   :: mfact
1174      real( kind = dp ), dimension(9,nLocal) :: eFrame
1175      real( kind = dp ), dimension(9,nLocal) :: A
# Line 899 | Line 1179 | contains
1179      logical, intent(inout) :: do_pot
1180      integer, intent(in) :: i, j
1181      real ( kind = dp ), intent(inout) :: rijsq
1182 <    real ( kind = dp )                :: r
1182 >    real ( kind = dp ), intent(inout) :: r_grp
1183      real ( kind = dp ), intent(inout) :: d(3)
1184 <    real ( kind = dp ) :: ebalance
1184 >    real ( kind = dp ), intent(inout) :: d_grp(3)
1185 >    real ( kind = dp ) :: r
1186      integer :: me_i, me_j
1187  
1188      integer :: iHash
# Line 919 | Line 1200 | contains
1200   #endif
1201  
1202      iHash = InteractionHash(me_i, me_j)
1203 <
1203 >    
1204      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1205 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1205 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1206 >            pot(VDW_POT), f, do_pot)
1207      endif
1208 <
1208 >    
1209      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1210         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1211 <            pot, eFrame, f, t, do_pot)
1212 <
1213 <       if (FF_uses_RF .and. SIM_uses_RF) then
932 <
933 <          ! CHECK ME (RF needs to know about all electrostatic types)
934 <          call accumulate_rf(i, j, r, eFrame, sw)
935 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
936 <       endif
937 <
1211 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1212 > !!$       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213 > !!$            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot, felec)
1214      endif
1215 <
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)
1218 >            pot(HB_POT), A, f, t, do_pot)
1219      endif
1220 <
1220 >    
1221      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1222         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1223 <            pot, A, f, t, do_pot)
1223 >            pot(HB_POT), A, f, t, do_pot)
1224      endif
1225 <
1225 >    
1226      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1227         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1228 <            pot, A, f, t, do_pot)
1228 >            pot(VDW_POT), A, f, t, do_pot)
1229      endif
1230      
1231      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1232 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1233 < !           pot, A, f, t, do_pot)
1232 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1233 >            pot(VDW_POT), A, f, t, do_pot)
1234      endif
1235 <
1235 >    
1236      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1237 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1238 <            do_pot)
1237 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1238 >            pot(METALLIC_POT), f, do_pot)
1239      endif
1240 <
1240 >    
1241      if ( iand(iHash, SHAPE_PAIR).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      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1247         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1248 <            pot, A, f, t, do_pot)
1248 >            pot(VDW_POT), A, f, t, do_pot)
1249      endif
1250 <    
1250 >    
1251    end subroutine do_pair
1252  
1253    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1254         do_pot, do_stress, eFrame, A, f, t, pot)
1255  
1256 <    real( kind = dp ) :: pot, sw
1256 >    real( kind = dp ) :: sw
1257 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1258      real( kind = dp ), dimension(9,nLocal) :: eFrame
1259      real (kind=dp), dimension(9,nLocal) :: A
1260      real (kind=dp), dimension(3,nLocal) :: f
# Line 991 | Line 1268 | contains
1268  
1269      integer :: me_i, me_j, iHash
1270  
1271 +    r = sqrt(rijsq)
1272 +
1273   #ifdef IS_MPI  
1274      me_i = atid_row(i)
1275      me_j = atid_col(j)  
# Line 1010 | Line 1289 | contains
1289  
1290    subroutine do_preforce(nlocal,pot)
1291      integer :: nlocal
1292 <    real( kind = dp ) :: pot
1292 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1293  
1294      if (FF_uses_EAM .and. SIM_uses_EAM) then
1295 <       call calc_EAM_preforce_Frho(nlocal,pot)
1295 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1296      endif
1297  
1298  
# Line 1099 | Line 1378 | contains
1378      pot_Col = 0.0_dp
1379      pot_Temp = 0.0_dp
1380  
1102 –    rf_Row = 0.0_dp
1103 –    rf_Col = 0.0_dp
1104 –    rf_Temp = 0.0_dp
1105 –
1381   #endif
1382  
1383      if (FF_uses_EAM .and. SIM_uses_EAM) then
# Line 1206 | Line 1481 | contains
1481      doesit = FF_uses_EAM
1482    end function FF_RequiresPrepairCalc
1483  
1209 –  function FF_RequiresPostpairCalc() result(doesit)
1210 –    logical :: doesit
1211 –    doesit = FF_uses_RF
1212 –  end function FF_RequiresPostpairCalc
1213 –
1484   #ifdef PROFILE
1485    function getforcetime() result(totalforcetime)
1486      real(kind=dp) :: totalforcetime

Diff Legend

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