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 2282 by chuckv, Tue Sep 6 17:32:42 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.36 2005-09-06 17:32:42 chuckv Exp $, $Date: 2005-09-06 17:32:42 $, $Name: not supported by cvs2svn $, $Revision: 1.36 $
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 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 256 | Line 267 | contains
267      logical :: i_is_GB
268      logical :: i_is_EAM
269      logical :: i_is_Shape
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
274 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
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
# Line 271 | 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         if (SimHasAtype(i)) then    
298            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 284 | 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 351 | 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 (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 <       if (nGroupTypes.eq.0) then
469 <          nGroupTypes = nGroupTypes + 1
470 <          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
471 <          groupToGtype(i) = nGroupTypes
468 >
469 >       if (nGroupTypesCol.eq.0) then
470 >          nGroupTypesCol = nGroupTypesCol + 1
471 >          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
472 >          groupToGtypeCol(j) = nGroupTypesCol
473         else
474 <          do g = 1, nGroupTypes
475 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
476 <                nGroupTypes = nGroupTypes + 1
477 <                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
478 <                groupToGtype(i) = nGroupTypes
368 <             else
369 <                groupToGtype(i) = g
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 <    
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 394 | 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 <    
403 <  end subroutine createGtypeCutoffMap
404 <  
405 <  subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
406 <    real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
407 <    integer, intent(in) :: cutPolicy
408 <    
409 <    defaultRcut = defRcut
410 <    defaultRsw = defRsw
411 <    defaultRlist = defRlist
412 <    cutoffPolicy = cutPolicy
413 <  end subroutine setDefaultCutoffs
414 <  
415 <  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
419 –     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()
427 –    SIM_uses_RF = SimUsesRF()
572      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574      SIM_uses_PBC = SimUsesPBC()
# Line 494 | 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
500 <    logical, intent(in) :: use_UW
501 <    logical, intent(in) :: use_DW
643 >    integer, intent(in) :: thisESM
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
504 –    integer :: corrMethod
646      integer, pointer :: MatchList(:) => null()
506 –    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:
512 <    FF_uses_RF = use_RF
651 >    electrostaticSummationMethod = thisESM
652  
514 –    !! set the electrostatic correction method
515 –    if (use_UW) then
516 –       corrMethod = 1
517 –    elseif (use_DW) then
518 –       corrMethod = 2
519 –    else
520 –       corrMethod = 0
521 –    endif
522 –    
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 549 | Line 679 | contains
679  
680      haveSaneForceField = .true.
681  
552 –    !! check to make sure the FF_uses_RF setting makes sense
553 –
554 –    if (FF_uses_Dipoles) then
555 –       if (FF_uses_RF) then
556 –          dielect = getDielect()
557 –          call initialize_rf(dielect)
558 –       endif
559 –    else
560 –       if (FF_uses_RF) then          
561 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
562 –          thisStat = -1
563 –          haveSaneForceField = .false.
564 –          return
565 –       endif
566 –    endif
567 –
682      if (FF_uses_EAM) then
683         call init_EAM_FF(my_status)
684         if (my_status /= 0) then
# Line 575 | Line 689 | contains
689         end if
690      endif
691  
578 –    if (FF_uses_GayBerne) then
579 –       call check_gb_pair_FF(my_status)
580 –       if (my_status .ne. 0) then
581 –          thisStat = -1
582 –          haveSaneForceField = .false.
583 –          return
584 –       endif
585 –    endif
586 –
692      if (.not. haveNeighborList) then
693         !! Create neighbor lists
694         call expandNeighborList(nLocal, my_status)
# Line 617 | 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 638 | 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 648 | 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 773 | 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 813 | 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 842 | 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 896 | Line 1003 | contains
1003                  endif
1004               end if
1005            enddo
1006 +
1007         enddo outer
1008  
1009         if (update_nlist) then
# Line 955 | 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
979 <
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)
983 <          do i = 1,nlocal
984 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
985 <          end do
986 < #endif
987 <
988 <          do i = 1, nLocal
989 <
990 <             rfpot = 0.0_DP
991 < #ifdef IS_MPI
992 <             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 <
1019 <
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
1025 <       !! 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 1057 | 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 1077 | 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, corrMethod)
1212 <
1213 <       if (FF_uses_RF .and. SIM_uses_RF) then
1090 <
1091 <          ! CHECK ME (RF needs to know about all electrostatic types)
1092 <          call accumulate_rf(i, j, r, eFrame, sw)
1093 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1094 <       endif
1095 <
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 1149 | 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 1168 | 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 1257 | Line 1378 | contains
1378      pot_Col = 0.0_dp
1379      pot_Temp = 0.0_dp
1380  
1260 –    rf_Row = 0.0_dp
1261 –    rf_Col = 0.0_dp
1262 –    rf_Temp = 0.0_dp
1263 –
1381   #endif
1382  
1383      if (FF_uses_EAM .and. SIM_uses_EAM) then
# Line 1364 | Line 1481 | contains
1481      doesit = FF_uses_EAM
1482    end function FF_RequiresPrepairCalc
1483  
1367 –  function FF_RequiresPostpairCalc() result(doesit)
1368 –    logical :: doesit
1369 –    doesit = FF_uses_RF
1370 –  end function FF_RequiresPostpairCalc
1371 –
1484   #ifdef PROFILE
1485    function getforcetime() result(totalforcetime)
1486      real(kind=dp) :: totalforcetime

Diff Legend

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