ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2350 by chuckv, Mon Oct 10 21:20:46 2005 UTC vs.
Revision 2530 by chuckv, Fri Dec 30 00:18:28 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.51 2005-10-10 21:20:46 chuckv Exp $, $Date: 2005-10-10 21:20:46 $, $Name: not supported by cvs2svn $, $Revision: 1.51 $
48 > !! @version $Id: doForces.F90,v 1.72 2005-12-30 00:18:28 chuckv Exp $, $Date: 2005-12-30 00:18:28 $, $Name: not supported by cvs2svn $, $Revision: 1.72 $
49  
50  
51   module doForces
# Line 58 | Line 58 | module doForces
58    use lj
59    use sticky
60    use electrostatic_module
61 <  use reaction_field_module
62 <  use gb_pair
61 >  use gayberne
62    use shapes
63    use vector_class
64    use eam
65 +  use suttonchen
66    use status
67   #ifdef IS_MPI
68    use mpiSimulation
# Line 87 | Line 87 | module doForces
87    logical, save :: haveInteractionHash = .false.
88    logical, save :: haveGtypeCutoffMap = .false.
89    logical, save :: haveDefaultCutoffs = .false.
90 <  logical, save :: haveRlist = .false.
90 >  logical, save :: haveSkinThickness = .false.
91 >  logical, save :: haveElectrostaticSummationMethod = .false.
92 >  logical, save :: haveCutoffPolicy = .false.
93 >  logical, save :: VisitCutoffsAfterComputing = .false.
94  
95    logical, save :: FF_uses_DirectionalAtoms
96    logical, save :: FF_uses_Dipoles
97    logical, save :: FF_uses_GayBerne
98    logical, save :: FF_uses_EAM
99 +  logical, save :: FF_uses_SC
100 +  logical, save :: FF_uses_MEAM
101 +
102  
103    logical, save :: SIM_uses_DirectionalAtoms
104    logical, save :: SIM_uses_EAM
105 +  logical, save :: SIM_uses_SC
106 +  logical, save :: SIM_uses_MEAM
107    logical, save :: SIM_requires_postpair_calc
108    logical, save :: SIM_requires_prepair_calc
109    logical, save :: SIM_uses_PBC
110  
111    integer, save :: electrostaticSummationMethod
112 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
113  
114 +  real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
115 +  real(kind=dp), save :: skinThickness
116 +  logical, save :: defaultDoShift
117 +
118    public :: init_FF
119 <  public :: setDefaultCutoffs
119 >  public :: setCutoffs
120 >  public :: cWasLame
121 >  public :: setElectrostaticMethod
122 >  public :: setCutoffPolicy
123 >  public :: setSkinThickness
124    public :: do_force_loop
108 –  public :: createInteractionHash
109 –  public :: createGtypeCutoffMap
110 –  public :: getStickyCut
111 –  public :: getStickyPowerCut
112 –  public :: getGayBerneCut
113 –  public :: getEAMCut
114 –  public :: getShapeCut
125  
126   #ifdef PROFILE
127    public :: getforcetime
# Line 139 | Line 149 | module doForces
149    end type gtypeCutoffs
150    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
151  
142 –  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
143 –  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
144 –  real(kind=dp),save :: listSkin
145 –  
152   contains
153  
154 <  subroutine createInteractionHash(status)
154 >  subroutine createInteractionHash()
155      integer :: nAtypes
150 –    integer, intent(out) :: status
156      integer :: i
157      integer :: j
158      integer :: iHash
# Line 159 | Line 164 | contains
164      logical :: i_is_GB
165      logical :: i_is_EAM
166      logical :: i_is_Shape
167 +    logical :: i_is_SC
168 +    logical :: i_is_MEAM
169      logical :: j_is_LJ
170      logical :: j_is_Elect
171      logical :: j_is_Sticky
# Line 166 | Line 173 | contains
173      logical :: j_is_GB
174      logical :: j_is_EAM
175      logical :: j_is_Shape
176 +    logical :: j_is_SC
177 +    logical :: j_is_MEAM
178      real(kind=dp) :: myRcut
179  
171 –    status = 0  
172 –
180      if (.not. associated(atypes)) then
181 <       call handleError("atype", "atypes was not present before call of createInteractionHash!")
175 <       status = -1
181 >       call handleError("doForces", "atypes was not present before call of createInteractionHash!")
182         return
183      endif
184      
185      nAtypes = getSize(atypes)
186      
187      if (nAtypes == 0) then
188 <       status = -1
188 >       call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
189         return
190      end if
191  
192      if (.not. allocated(InteractionHash)) then
193         allocate(InteractionHash(nAtypes,nAtypes))
194 +    else
195 +       deallocate(InteractionHash)
196 +       allocate(InteractionHash(nAtypes,nAtypes))
197      endif
198  
199      if (.not. allocated(atypeMaxCutoff)) then
200         allocate(atypeMaxCutoff(nAtypes))
201 +    else
202 +       deallocate(atypeMaxCutoff)
203 +       allocate(atypeMaxCutoff(nAtypes))
204      endif
205          
206      do i = 1, nAtypes
# Line 199 | Line 211 | contains
211         call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
212         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
213         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
214 +       call getElementProperty(atypes, i, "is_SC", i_is_SC)
215 +       call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
216  
217         do j = i, nAtypes
218  
# Line 212 | Line 226 | contains
226            call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
227            call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
228            call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
229 +          call getElementProperty(atypes, j, "is_SC", j_is_SC)
230 +          call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
231  
232            if (i_is_LJ .and. j_is_LJ) then
233               iHash = ior(iHash, LJ_PAIR)            
# Line 233 | Line 249 | contains
249               iHash = ior(iHash, EAM_PAIR)
250            endif
251  
252 +          if (i_is_SC .and. j_is_SC) then
253 +             iHash = ior(iHash, SC_PAIR)
254 +          endif
255 +
256            if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
257            if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
258            if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
# Line 252 | Line 272 | contains
272      haveInteractionHash = .true.
273    end subroutine createInteractionHash
274  
275 <  subroutine createGtypeCutoffMap(stat)
275 >  subroutine createGtypeCutoffMap()
276  
257 –    integer, intent(out), optional :: stat
277      logical :: i_is_LJ
278      logical :: i_is_Elect
279      logical :: i_is_Sticky
# Line 262 | Line 281 | contains
281      logical :: i_is_GB
282      logical :: i_is_EAM
283      logical :: i_is_Shape
284 +    logical :: i_is_SC
285      logical :: GtypeFound
286  
287      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
288 <    integer :: n_in_i, me_i, ia, g, atom1
288 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
289      integer :: nGroupsInRow
290      integer :: nGroupsInCol
291      integer :: nGroupTypesRow,nGroupTypesCol
292 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
292 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
293      real(kind=dp) :: biggestAtypeCutoff
294  
275 –    stat = 0
295      if (.not. haveInteractionHash) then
296 <       call createInteractionHash(myStatus)      
278 <       if (myStatus .ne. 0) then
279 <          write(default_error, *) 'createInteractionHash failed in doForces!'
280 <          stat = -1
281 <          return
282 <       endif
296 >       call createInteractionHash()      
297      endif
298   #ifdef IS_MPI
299      nGroupsInRow = getNgroupsInRow(plan_group_row)
# Line 297 | Line 311 | contains
311            call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
312            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
313            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
314 <          
314 >          call getElementProperty(atypes, i, "is_SC", i_is_SC)
315  
316            if (haveDefaultCutoffs) then
317               atypeMaxCutoff(i) = defaultRcut
# Line 330 | Line 344 | contains
344                  thisRcut = getShapeCut(i)
345                  if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
346               endif
347 +             if (i_is_SC) then
348 +                thisRcut = getSCCut(i)
349 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
350 +             endif
351            endif
352 <          
335 <          
352 >                    
353            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
354               biggestAtypeCutoff = atypeMaxCutoff(i)
355            endif
356  
357         endif
358      enddo
342 –  
343 –
359      
360      istart = 1
361      jstart = 1
# Line 377 | Line 392 | contains
392  
393   #ifdef IS_MPI
394         ! We only allocate new storage if we are in MPI because Ncol /= Nrow
395 <    if(.not.allocated(groupToGtypeCol)) then
395 >    if(.not.associated(groupToGtypeCol)) then
396         allocate(groupToGtypeCol(jend))
397      else
398         deallocate(groupToGtypeCol)
399         allocate(groupToGtypeCol(jend))
400      end if
401  
402 <    if(.not.allocated(groupToGtypeCol)) then
402 >    if(.not.associated(groupToGtypeCol)) then
403         allocate(groupToGtypeCol(jend))
404      else
405         deallocate(groupToGtypeCol)
406         allocate(groupToGtypeCol(jend))
407      end if
408 <    if(.not.allocated(gtypeMaxCutoffCol)) then
408 >    if(.not.associated(gtypeMaxCutoffCol)) then
409         allocate(gtypeMaxCutoffCol(jend))
410      else
411         deallocate(gtypeMaxCutoffCol)      
# Line 426 | Line 441 | contains
441               groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
442            endif          
443         enddo
429 –
444         if (nGroupTypesRow.eq.0) then
445            nGroupTypesRow = nGroupTypesRow + 1
446            gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
# Line 489 | Line 503 | contains
503      groupMaxCutoffCol => groupMaxCutoffRow
504   #endif
505  
492 –
493 –
494 –
495 –
506      !! allocate the gtypeCutoffMap here.
507      allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
508      !! then we do a double loop over all the group TYPES to find the cutoff
509      !! map between groups of two types
510      tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
511  
512 <    do i = 1, nGroupTypesRow
512 >    do i = 1, nGroupTypesRow      
513         do j = 1, nGroupTypesCol
514        
515            select case(cutoffPolicy)
# Line 514 | Line 524 | contains
524               return
525            end select
526            gtypeCutoffMap(i,j)%rcut = thisRcut
527 +          
528 +          if (thisRcut.gt.largestRcut) largestRcut = thisRcut
529 +
530            gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
518 –          skin = defaultRlist - defaultRcut
519 –          listSkin = skin ! set neighbor list skin thickness
520 –          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
531  
532 +          if (.not.haveSkinThickness) then
533 +             skinThickness = 1.0_dp
534 +          endif
535 +
536 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
537 +
538            ! sanity check
539  
540            if (haveDefaultCutoffs) then
# Line 528 | Line 544 | contains
544            endif
545         enddo
546      enddo
547 +
548      if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
549      if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
550      if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
# Line 541 | Line 558 | contains
558      haveGtypeCutoffMap = .true.
559     end subroutine createGtypeCutoffMap
560  
561 <   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
545 <     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
546 <     integer, intent(in) :: cutPolicy
561 >   subroutine setCutoffs(defRcut, defRsw)
562  
563 +     real(kind=dp),intent(in) :: defRcut, defRsw
564 +     character(len = statusMsgSize) :: errMsg
565 +     integer :: localError
566 +
567       defaultRcut = defRcut
568       defaultRsw = defRsw
569 <     defaultRlist = defRlist
570 <     cutoffPolicy = cutPolicy
569 >    
570 >     defaultDoShift = .false.
571 >     if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
572 >        
573 >        write(errMsg, *) &
574 >             'cutoffRadius and switchingRadius are set to the same', newline &
575 >             // tab, 'value.  OOPSE will use shifted ', newline &
576 >             // tab, 'potentials instead of switching functions.'
577 >        
578 >        call handleInfo("setCutoffs", errMsg)
579 >        
580 >        defaultDoShift = .true.
581 >        
582 >     endif
583  
584 +     localError = 0
585 +     call setLJDefaultCutoff( defaultRcut, defaultDoShift )
586 +     call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
587 +     call setCutoffEAM( defaultRcut, localError)
588 +     if (localError /= 0) then
589 +       write(errMsg, *) 'An error has occured in setting the EAM cutoff'
590 +       call handleError("setCutoffs", errMsg)
591 +     end if
592 +     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
593 +
594       haveDefaultCutoffs = .true.
595 <   end subroutine setDefaultCutoffs
595 >     haveGtypeCutoffMap = .false.
596 >   end subroutine setCutoffs
597  
598 +   subroutine cWasLame()
599 +    
600 +     VisitCutoffsAfterComputing = .true.
601 +     return
602 +    
603 +   end subroutine cWasLame
604 +  
605     subroutine setCutoffPolicy(cutPolicy)
606 <
606 >    
607       integer, intent(in) :: cutPolicy
608 +    
609       cutoffPolicy = cutPolicy
610 <     call createGtypeCutoffMap()
611 <   end subroutine setCutoffPolicy
562 <    
610 >     haveCutoffPolicy = .true.
611 >     haveGtypeCutoffMap = .false.
612      
613 <  subroutine setSimVariables()
614 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
615 <    SIM_uses_EAM = SimUsesEAM()
567 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
568 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
569 <    SIM_uses_PBC = SimUsesPBC()
613 >   end subroutine setCutoffPolicy
614 >  
615 >   subroutine setElectrostaticMethod( thisESM )
616  
617 <    haveSIMvariables = .true.
617 >     integer, intent(in) :: thisESM
618  
619 <    return
620 <  end subroutine setSimVariables
619 >     electrostaticSummationMethod = thisESM
620 >     haveElectrostaticSummationMethod = .true.
621 >    
622 >   end subroutine setElectrostaticMethod
623 >
624 >   subroutine setSkinThickness( thisSkin )
625 >    
626 >     real(kind=dp), intent(in) :: thisSkin
627 >    
628 >     skinThickness = thisSkin
629 >     haveSkinThickness = .true.    
630 >     haveGtypeCutoffMap = .false.
631 >    
632 >   end subroutine setSkinThickness
633 >      
634 >   subroutine setSimVariables()
635 >     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
636 >     SIM_uses_EAM = SimUsesEAM()
637 >     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
638 >     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
639 >     SIM_uses_PBC = SimUsesPBC()
640 >    
641 >     haveSIMvariables = .true.
642 >    
643 >     return
644 >   end subroutine setSimVariables
645  
646    subroutine doReadyCheck(error)
647      integer, intent(out) :: error
# Line 581 | Line 651 | contains
651      error = 0
652  
653      if (.not. haveInteractionHash) then      
654 <       myStatus = 0      
585 <       call createInteractionHash(myStatus)      
586 <       if (myStatus .ne. 0) then
587 <          write(default_error, *) 'createInteractionHash failed in doForces!'
588 <          error = -1
589 <          return
590 <       endif
654 >       call createInteractionHash()      
655      endif
656  
657      if (.not. haveGtypeCutoffMap) then        
658 <       myStatus = 0      
595 <       call createGtypeCutoffMap(myStatus)      
596 <       if (myStatus .ne. 0) then
597 <          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
598 <          error = -1
599 <          return
600 <       endif
658 >       call createGtypeCutoffMap()      
659      endif
660  
661 +
662 +    if (VisitCutoffsAfterComputing) then
663 +       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
664 +    endif
665 +
666 +
667      if (.not. haveSIMvariables) then
668         call setSimVariables()
669      endif
# Line 633 | Line 697 | contains
697    end subroutine doReadyCheck
698  
699  
700 <  subroutine init_FF(thisESM, thisStat)
700 >  subroutine init_FF(thisStat)
701  
638 –    integer, intent(in) :: thisESM
702      integer, intent(out) :: thisStat  
703      integer :: my_status, nMatches
704      integer, pointer :: MatchList(:) => null()
642 –    real(kind=dp) :: rcut, rrf, rt, dielect
705  
706      !! assume things are copacetic, unless they aren't
707      thisStat = 0
708  
647 –    electrostaticSummationMethod = thisESM
648 –
709      !! init_FF is called *after* all of the atom types have been
710      !! defined in atype_module using the new_atype subroutine.
711      !!
# Line 675 | Line 735 | contains
735  
736      haveSaneForceField = .true.
737  
678 –    !! check to make sure the reaction field setting makes sense
679 –
680 –    if (FF_uses_Dipoles) then
681 –       if (electrostaticSummationMethod == REACTION_FIELD) then
682 –          dielect = getDielect()
683 –          call initialize_rf(dielect)
684 –       endif
685 –    else
686 –       if (electrostaticSummationMethod == REACTION_FIELD) then
687 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
688 –          thisStat = -1
689 –          haveSaneForceField = .false.
690 –          return
691 –       endif
692 –    endif
693 –
738      if (FF_uses_EAM) then
739         call init_EAM_FF(my_status)
740         if (my_status /= 0) then
# Line 701 | Line 745 | contains
745         end if
746      endif
747  
704 –    if (FF_uses_GayBerne) then
705 –       call check_gb_pair_FF(my_status)
706 –       if (my_status .ne. 0) then
707 –          thisStat = -1
708 –          haveSaneForceField = .false.
709 –          return
710 –       endif
711 –    endif
712 –
748      if (.not. haveNeighborList) then
749         !! Create neighbor lists
750         call expandNeighborList(nLocal, my_status)
# Line 743 | Line 778 | contains
778  
779      !! Stress Tensor
780      real( kind = dp), dimension(9) :: tau  
781 <    real ( kind = dp ) :: pot
781 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
782      logical ( kind = 2) :: do_pot_c, do_stress_c
783      logical :: do_pot
784      logical :: do_stress
785      logical :: in_switching_region
786   #ifdef IS_MPI
787 <    real( kind = DP ) :: pot_local
787 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
788      integer :: nAtomsInRow
789      integer :: nAtomsInCol
790      integer :: nprocs
# Line 764 | Line 799 | contains
799      integer :: nlist
800      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
801      real( kind = DP ) :: sw, dswdr, swderiv, mf
802 +    real( kind = DP ) :: rVal
803      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
804      real(kind=dp) :: rfpot, mu_i, virial
805 +    real(kind=dp):: rCut
806      integer :: me_i, me_j, n_in_i, n_in_j
807      logical :: is_dp_i
808      integer :: neighborListSize
# Line 774 | Line 811 | contains
811      integer :: propPack_i, propPack_j
812      integer :: loopStart, loopEnd, loop
813      integer :: iHash
814 +    integer :: i1
815    
816  
817      !! initialize local variables  
# Line 838 | Line 876 | contains
876         ! (but only on the first time through):
877         if (loop .eq. loopStart) then
878   #ifdef IS_MPI
879 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
879 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
880                 update_nlist)
881   #else
882 <          call checkNeighborList(nGroups, q_group, listSkin, &
882 >          call checkNeighborList(nGroups, q_group, skinThickness, &
883                 update_nlist)
884   #endif
885         endif
# Line 921 | Line 959 | contains
959  
960                     list(nlist) = j
961                  endif
962 +
963  
964 <                if (loop .eq. PAIR_LOOP) then
965 <                   vij = 0.0d0
927 <                   fij(1:3) = 0.0d0
928 <                endif
929 <
930 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
931 <                     in_switching_region)
932 <
933 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
934 <
935 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
964 >                
965 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
966  
967 <                   atom1 = groupListRow(ia)
968 <
969 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
970 <
971 <                      atom2 = groupListCol(jb)
972 <
973 <                      if (skipThisPair(atom1, atom2)) cycle inner
974 <
975 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
976 <                         d_atm(1:3) = d_grp(1:3)
977 <                         ratmsq = rgrpsq
978 <                      else
967 >                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
968 >                   if (loop .eq. PAIR_LOOP) then
969 >                      vij = 0.0d0
970 >                      fij(1:3) = 0.0d0
971 >                   endif
972 >                  
973 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
974 >                        group_switch, in_switching_region)
975 >                  
976 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
977 >                  
978 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
979 >                      
980 >                      atom1 = groupListRow(ia)
981 >                      
982 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
983 >                        
984 >                         atom2 = groupListCol(jb)
985 >                        
986 >                         if (skipThisPair(atom1, atom2))  cycle inner
987 >                        
988 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
989 >                            d_atm(1:3) = d_grp(1:3)
990 >                            ratmsq = rgrpsq
991 >                         else
992   #ifdef IS_MPI
993 <                         call get_interatomic_vector(q_Row(:,atom1), &
994 <                              q_Col(:,atom2), d_atm, ratmsq)
993 >                            call get_interatomic_vector(q_Row(:,atom1), &
994 >                                 q_Col(:,atom2), d_atm, ratmsq)
995   #else
996 <                         call get_interatomic_vector(q(:,atom1), &
997 <                              q(:,atom2), d_atm, ratmsq)
996 >                            call get_interatomic_vector(q(:,atom1), &
997 >                                 q(:,atom2), d_atm, ratmsq)
998   #endif
999 <                      endif
1000 <
1001 <                      if (loop .eq. PREPAIR_LOOP) then
999 >                         endif
1000 >                        
1001 >                         if (loop .eq. PREPAIR_LOOP) then
1002   #ifdef IS_MPI                      
1003 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1004 <                              rgrpsq, d_grp, do_pot, do_stress, &
1005 <                              eFrame, A, f, t, pot_local)
1003 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1004 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1005 >                                 eFrame, A, f, t, pot_local)
1006   #else
1007 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1008 <                              rgrpsq, d_grp, do_pot, do_stress, &
1009 <                              eFrame, A, f, t, pot)
1007 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1008 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1009 >                                 eFrame, A, f, t, pot)
1010   #endif                                              
1011 <                      else
1011 >                         else
1012   #ifdef IS_MPI                      
1013 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1014 <                              do_pot, &
1015 <                              eFrame, A, f, t, pot_local, vpair, fpair)
1013 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1014 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
1015 >                                 fpair, d_grp, rgrp, rCut)
1016   #else
1017 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1018 <                              do_pot,  &
1019 <                              eFrame, A, f, t, pot, vpair, fpair)
1020 < #endif
1017 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1018 >                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1019 >                                 d_grp, rgrp, rCut)
1020 > #endif
1021 >                            vij = vij + vpair
1022 >                            fij(1:3) = fij(1:3) + fpair(1:3)
1023 >                         endif
1024 >                      enddo inner
1025 >                   enddo
1026  
1027 <                         vij = vij + vpair
1028 <                         fij(1:3) = fij(1:3) + fpair(1:3)
1029 <                      endif
1030 <                   enddo inner
1031 <                enddo
1032 <
1033 <                if (loop .eq. PAIR_LOOP) then
1034 <                   if (in_switching_region) then
1035 <                      swderiv = vij*dswdr/rgrp
1036 <                      fij(1) = fij(1) + swderiv*d_grp(1)
989 <                      fij(2) = fij(2) + swderiv*d_grp(2)
990 <                      fij(3) = fij(3) + swderiv*d_grp(3)
991 <
992 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
993 <                         atom1=groupListRow(ia)
994 <                         mf = mfactRow(atom1)
1027 >                   if (loop .eq. PAIR_LOOP) then
1028 >                      if (in_switching_region) then
1029 >                         swderiv = vij*dswdr/rgrp
1030 >                         fij(1) = fij(1) + swderiv*d_grp(1)
1031 >                         fij(2) = fij(2) + swderiv*d_grp(2)
1032 >                         fij(3) = fij(3) + swderiv*d_grp(3)
1033 >                        
1034 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
1035 >                            atom1=groupListRow(ia)
1036 >                            mf = mfactRow(atom1)
1037   #ifdef IS_MPI
1038 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1039 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1040 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1038 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1039 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1040 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1041   #else
1042 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1043 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1044 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1042 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1043 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1044 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1045   #endif
1046 <                      enddo
1047 <
1048 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
1049 <                         atom2=groupListCol(jb)
1050 <                         mf = mfactCol(atom2)
1046 >                         enddo
1047 >                        
1048 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
1049 >                            atom2=groupListCol(jb)
1050 >                            mf = mfactCol(atom2)
1051   #ifdef IS_MPI
1052 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1053 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1054 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1052 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1053 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1054 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1055   #else
1056 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1057 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1058 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1056 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1057 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1058 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1059   #endif
1060 <                      enddo
1061 <                   endif
1060 >                         enddo
1061 >                      endif
1062  
1063 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1063 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1064 >                   endif
1065                  endif
1066 <             end if
1066 >             endif
1067            enddo
1068 <
1068 >          
1069         enddo outer
1070  
1071         if (update_nlist) then
# Line 1082 | Line 1125 | contains
1125  
1126      if (do_pot) then
1127         ! scatter/gather pot_row into the members of my column
1128 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1129 <
1128 >       do i = 1,LR_POT_TYPES
1129 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1130 >       end do
1131         ! scatter/gather pot_local into all other procs
1132         ! add resultant to get total pot
1133         do i = 1, nlocal
1134 <          pot_local = pot_local + pot_Temp(i)
1134 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1135 >               + pot_Temp(1:LR_POT_TYPES,i)
1136         enddo
1137  
1138         pot_Temp = 0.0_DP
1139 <
1140 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1139 >       do i = 1,LR_POT_TYPES
1140 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1141 >       end do
1142         do i = 1, nlocal
1143 <          pot_local = pot_local + pot_Temp(i)
1143 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1144 >               + pot_Temp(1:LR_POT_TYPES,i)
1145         enddo
1146  
1147      endif
1148   #endif
1149  
1150 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1150 >    if (SIM_requires_postpair_calc) then
1151 >       do i = 1, nlocal            
1152 >          
1153 >          ! we loop only over the local atoms, so we don't need row and column
1154 >          ! lookups for the types
1155 >          
1156 >          me_i = atid(i)
1157 >          
1158 >          ! is the atom electrostatic?  See if it would have an
1159 >          ! electrostatic interaction with itself
1160 >          iHash = InteractionHash(me_i,me_i)
1161  
1162 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1106 <
1162 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1163   #ifdef IS_MPI
1164 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1165 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1110 <          do i = 1,nlocal
1111 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1112 <          end do
1113 < #endif
1114 <
1115 <          do i = 1, nLocal
1116 <
1117 <             rfpot = 0.0_DP
1118 < #ifdef IS_MPI
1119 <             me_i = atid_row(i)
1164 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1165 >                  t, do_pot)
1166   #else
1167 <             me_i = atid(i)
1167 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1168 >                  t, do_pot)
1169   #endif
1170 <             iHash = InteractionHash(me_i,me_j)
1170 >          endif
1171 >  
1172 >          
1173 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1174              
1175 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1176 <
1177 <                mu_i = getDipoleMoment(me_i)
1178 <
1179 <                !! The reaction field needs to include a self contribution
1180 <                !! to the field:
1181 <                call accumulate_self_rf(i, mu_i, eFrame)
1182 <                !! Get the reaction field contribution to the
1183 <                !! potential and torques:
1184 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1175 >             ! loop over the excludes to accumulate RF stuff we've
1176 >             ! left out of the normal pair loop
1177 >            
1178 >             do i1 = 1, nSkipsForAtom(i)
1179 >                j = skipsForAtom(i, i1)
1180 >                
1181 >                ! prevent overcounting of the skips
1182 >                if (i.lt.j) then
1183 >                   call get_interatomic_vector(q(:,i), &
1184 >                        q(:,j), d_atm, ratmsq)
1185 >                   rVal = dsqrt(ratmsq)
1186 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1187 >                        in_switching_region)
1188   #ifdef IS_MPI
1189 <                pot_local = pot_local + rfpot
1189 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1190 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1191   #else
1192 <                pot = pot + rfpot
1193 <
1192 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1193 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1194   #endif
1195 <             endif
1196 <          enddo
1197 <       endif
1198 <    endif
1199 <
1200 <
1201 < #ifdef IS_MPI
1202 <
1195 >                endif
1196 >             enddo
1197 >          endif
1198 >       enddo
1199 >    endif
1200 >    
1201 > #ifdef IS_MPI
1202 >    
1203      if (do_pot) then
1204 <       pot = pot + pot_local
1205 <       !! we assume the c code will do the allreduce to get the total potential
1152 <       !! we could do it right here if we needed to...
1204 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1205 >            mpi_comm_world,mpi_err)            
1206      endif
1207 <
1207 >    
1208      if (do_stress) then
1209         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1210              mpi_comm_world,mpi_err)
1211         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1212              mpi_comm_world,mpi_err)
1213      endif
1214 <
1214 >    
1215   #else
1216 <
1216 >    
1217      if (do_stress) then
1218         tau = tau_Temp
1219         virial = virial_Temp
1220      endif
1221 <
1221 >    
1222   #endif
1223 <
1223 >    
1224    end subroutine do_force_loop
1225  
1226    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1227 <       eFrame, A, f, t, pot, vpair, fpair)
1227 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1228  
1229 <    real( kind = dp ) :: pot, vpair, sw
1229 >    real( kind = dp ) :: vpair, sw
1230 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1231      real( kind = dp ), dimension(3) :: fpair
1232      real( kind = dp ), dimension(nLocal)   :: mfact
1233      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1184 | Line 1238 | contains
1238      logical, intent(inout) :: do_pot
1239      integer, intent(in) :: i, j
1240      real ( kind = dp ), intent(inout) :: rijsq
1241 <    real ( kind = dp )                :: r
1241 >    real ( kind = dp ), intent(inout) :: r_grp
1242      real ( kind = dp ), intent(inout) :: d(3)
1243 +    real ( kind = dp ), intent(inout) :: d_grp(3)
1244 +    real ( kind = dp ), intent(inout) :: rCut
1245 +    real ( kind = dp ) :: r
1246      integer :: me_i, me_j
1247  
1248      integer :: iHash
# Line 1203 | Line 1260 | contains
1260   #endif
1261  
1262      iHash = InteractionHash(me_i, me_j)
1263 <
1263 >    
1264      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1265 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1265 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1266 >            pot(VDW_POT), f, do_pot)
1267      endif
1268 <
1268 >    
1269      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1270 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1271 <            pot, eFrame, f, t, do_pot)
1214 <
1215 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1216 <
1217 <          ! CHECK ME (RF needs to know about all electrostatic types)
1218 <          call accumulate_rf(i, j, r, eFrame, sw)
1219 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1220 <       endif
1221 <
1270 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1271 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1272      endif
1273 <
1273 >    
1274      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1275         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1276 <            pot, A, f, t, do_pot)
1276 >            pot(HB_POT), A, f, t, do_pot)
1277      endif
1278 <
1278 >    
1279      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1280         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1281 <            pot, A, f, t, do_pot)
1281 >            pot(HB_POT), A, f, t, do_pot)
1282      endif
1283 <
1283 >    
1284      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1285         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1286 <            pot, A, f, t, do_pot)
1286 >            pot(VDW_POT), A, f, t, do_pot)
1287      endif
1288      
1289      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1290 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1291 < !           pot, A, f, t, do_pot)
1290 >       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1291 >            pot(VDW_POT), A, f, t, do_pot)
1292      endif
1293 <
1293 >    
1294      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1295 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1296 <            do_pot)
1295 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1296 >            pot(METALLIC_POT), f, do_pot)
1297      endif
1298 <
1298 >    
1299      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1300         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1301 <            pot, A, f, t, do_pot)
1301 >            pot(VDW_POT), A, f, t, do_pot)
1302      endif
1303 <
1303 >    
1304      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1305         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1306 <            pot, A, f, t, do_pot)
1306 >            pot(VDW_POT), A, f, t, do_pot)
1307      endif
1308 +
1309 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1310 +       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1311 +            pot(METALLIC_POT), f, do_pot)
1312 +    endif
1313 +
1314      
1315 +    
1316    end subroutine do_pair
1317  
1318 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1318 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1319         do_pot, do_stress, eFrame, A, f, t, pot)
1320  
1321 <    real( kind = dp ) :: pot, sw
1321 >    real( kind = dp ) :: sw
1322 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1323      real( kind = dp ), dimension(9,nLocal) :: eFrame
1324      real (kind=dp), dimension(9,nLocal) :: A
1325      real (kind=dp), dimension(3,nLocal) :: f
# Line 1269 | Line 1327 | contains
1327  
1328      logical, intent(inout) :: do_pot, do_stress
1329      integer, intent(in) :: i, j
1330 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1330 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1331      real ( kind = dp )                :: r, rc
1332      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1333  
# Line 1288 | Line 1346 | contains
1346      iHash = InteractionHash(me_i, me_j)
1347  
1348      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1349 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1349 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1350      endif
1351 +
1352 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1353 +            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1354 +    endif
1355      
1356    end subroutine do_prepair
1357  
1358  
1359    subroutine do_preforce(nlocal,pot)
1360      integer :: nlocal
1361 <    real( kind = dp ) :: pot
1361 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1362  
1363      if (FF_uses_EAM .and. SIM_uses_EAM) then
1364 <       call calc_EAM_preforce_Frho(nlocal,pot)
1364 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1365      endif
1366 +    if (FF_uses_SC .and. SIM_uses_SC) then
1367 +       call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1368 +    endif
1369  
1370  
1371    end subroutine do_preforce
# Line 1385 | Line 1450 | contains
1450      pot_Col = 0.0_dp
1451      pot_Temp = 0.0_dp
1452  
1388 –    rf_Row = 0.0_dp
1389 –    rf_Col = 0.0_dp
1390 –    rf_Temp = 0.0_dp
1391 –
1453   #endif
1454  
1455      if (FF_uses_EAM .and. SIM_uses_EAM) then
1456         call clean_EAM()
1457      endif
1458  
1398 –    rf = 0.0_dp
1459      tau_Temp = 0.0_dp
1460      virial_Temp = 0.0_dp
1461    end subroutine zero_work_arrays
# Line 1489 | Line 1549 | contains
1549  
1550    function FF_RequiresPrepairCalc() result(doesit)
1551      logical :: doesit
1552 <    doesit = FF_uses_EAM
1552 >    doesit = FF_uses_EAM .or. FF_uses_SC &
1553 >         .or. FF_uses_MEAM
1554    end function FF_RequiresPrepairCalc
1555  
1495 –  function FF_RequiresPostpairCalc() result(doesit)
1496 –    logical :: doesit
1497 –    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1498 –  end function FF_RequiresPostpairCalc
1499 –
1556   #ifdef PROFILE
1557    function getforcetime() result(totalforcetime)
1558      real(kind=dp) :: totalforcetime

Diff Legend

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