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 2411 by chrisfen, Wed Nov 2 21:01:21 2005 UTC vs.
Revision 2532 by tim, Fri Dec 30 21:25:56 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.67 2005-11-02 21:01:18 chrisfen Exp $, $Date: 2005-11-02 21:01:18 $, $Name: not supported by cvs2svn $, $Revision: 1.67 $
48 > !! @version $Id: doForces.F90,v 1.73 2005-12-30 21:25:56 tim Exp $, $Date: 2005-12-30 21:25:56 $, $Name: not supported by cvs2svn $, $Revision: 1.73 $
49  
50  
51   module doForces
# Line 62 | Line 62 | module doForces
62    use shapes
63    use vector_class
64    use eam
65 +  use suttonchen
66    use status
67   #ifdef IS_MPI
68    use mpiSimulation
# Line 86 | 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
107  public :: createInteractionHash
108  public :: createGtypeCutoffMap
109  public :: getStickyCut
110  public :: getStickyPowerCut
111  public :: getGayBerneCut
112  public :: getEAMCut
113  public :: getShapeCut
125  
126   #ifdef PROFILE
127    public :: getforcetime
# Line 138 | Line 149 | module doForces
149    end type gtypeCutoffs
150    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
151  
141  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
142  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
143  real(kind=dp),save :: listSkin
144  
152   contains
153  
154 <  subroutine createInteractionHash(status)
154 >  subroutine createInteractionHash()
155      integer :: nAtypes
149    integer, intent(out) :: status
156      integer :: i
157      integer :: j
158      integer :: iHash
# Line 158 | 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 165 | 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  
170    status = 0  
171
180      if (.not. associated(atypes)) then
181 <       call handleError("atype", "atypes was not present before call of createInteractionHash!")
174 <       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  
# Line 204 | 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 217 | 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 238 | 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 257 | Line 272 | contains
272      haveInteractionHash = .true.
273    end subroutine createInteractionHash
274  
275 <  subroutine createGtypeCutoffMap(stat)
275 >  subroutine createGtypeCutoffMap()
276  
262    integer, intent(out), optional :: stat
277      logical :: i_is_LJ
278      logical :: i_is_Elect
279      logical :: i_is_Sticky
# Line 267 | 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
# Line 274 | Line 289 | contains
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  
280    stat = 0
295      if (.not. haveInteractionHash) then
296 <       call createInteractionHash(myStatus)      
283 <       if (myStatus .ne. 0) then
284 <          write(default_error, *) 'createInteractionHash failed in doForces!'
285 <          stat = -1
286 <          return
287 <       endif
296 >       call createInteractionHash()      
297      endif
298   #ifdef IS_MPI
299      nGroupsInRow = getNgroupsInRow(plan_group_row)
# Line 302 | 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 335 | 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 <          
340 <          
352 >                    
353            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
354               biggestAtypeCutoff = atypeMaxCutoff(i)
355            endif
356  
357         endif
358      enddo
347  
348
359      
360      istart = 1
361      jstart = 1
# Line 389 | Line 399 | contains
399         allocate(groupToGtypeCol(jend))
400      end if
401  
402 <    if(.not.associated(groupToGtypeCol)) then
403 <       allocate(groupToGtypeCol(jend))
402 >    if(.not.associated(groupMaxCutoffCol)) then
403 >       allocate(groupMaxCutoffCol(jend))
404      else
405 <       deallocate(groupToGtypeCol)
406 <       allocate(groupToGtypeCol(jend))
405 >       deallocate(groupMaxCutoffCol)
406 >       allocate(groupMaxCutoffCol(jend))
407      end if
408      if(.not.associated(gtypeMaxCutoffCol)) then
409         allocate(gtypeMaxCutoffCol(jend))
# Line 416 | Line 426 | contains
426      
427      tol = 1.0d-6
428      nGroupTypesRow = 0
429 <
429 >    nGroupTypesCol = 0
430      do i = istart, iend      
431         n_in_i = groupStartRow(i+1) - groupStartRow(i)
432         groupMaxCutoffRow(i) = 0.0_dp
# Line 431 | Line 441 | contains
441               groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
442            endif          
443         enddo
434
444         if (nGroupTypesRow.eq.0) then
445            nGroupTypesRow = nGroupTypesRow + 1
446            gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
# Line 494 | Line 503 | contains
503      groupMaxCutoffCol => groupMaxCutoffRow
504   #endif
505  
497
498
499
500
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 519 | 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
523          skin = defaultRlist - defaultRcut
524          listSkin = skin ! set neighbor list skin thickness
525          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 533 | 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 546 | Line 558 | contains
558      haveGtypeCutoffMap = .true.
559     end subroutine createGtypeCutoffMap
560  
561 <   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
550 <     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
551 <     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 <     haveDefaultCutoffs = .true.
585 <   end subroutine setDefaultCutoffs
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 <   subroutine setCutoffPolicy(cutPolicy)
594 >     haveDefaultCutoffs = .true.
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 +    
607       integer, intent(in) :: cutPolicy
608 +    
609       cutoffPolicy = cutPolicy
610 <     call createGtypeCutoffMap()
611 <   end subroutine setCutoffPolicy
567 <    
610 >     haveCutoffPolicy = .true.
611 >     haveGtypeCutoffMap = .false.
612      
613 <  subroutine setSimVariables()
614 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
615 <    SIM_uses_EAM = SimUsesEAM()
572 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574 <    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
648  
# Line 586 | Line 651 | contains
651      error = 0
652  
653      if (.not. haveInteractionHash) then      
654 <       myStatus = 0      
590 <       call createInteractionHash(myStatus)      
591 <       if (myStatus .ne. 0) then
592 <          write(default_error, *) 'createInteractionHash failed in doForces!'
593 <          error = -1
594 <          return
595 <       endif
654 >       call createInteractionHash()      
655      endif
656  
657      if (.not. haveGtypeCutoffMap) then        
658 <       myStatus = 0      
600 <       call createGtypeCutoffMap(myStatus)      
601 <       if (myStatus .ne. 0) then
602 <          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
603 <          error = -1
604 <          return
605 <       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 638 | Line 697 | contains
697    end subroutine doReadyCheck
698  
699  
700 <  subroutine init_FF(thisESM, thisStat)
700 >  subroutine init_FF(thisStat)
701  
643    integer, intent(in) :: thisESM
702      integer, intent(out) :: thisStat  
703      integer :: my_status, nMatches
704      integer, pointer :: MatchList(:) => null()
# Line 648 | Line 706 | contains
706      !! assume things are copacetic, unless they aren't
707      thisStat = 0
708  
651    electrostaticSummationMethod = thisESM
652
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 746 | Line 802 | contains
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 819 | 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 902 | Line 959 | contains
959  
960                     list(nlist) = j
961                  endif
962 +
963 +
964                  
965                  if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
966  
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, group_switch, &
974 <                        in_switching_region)
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                    
# Line 941 | Line 1001 | contains
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, &
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, &
1008 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1009                                   eFrame, A, f, t, pot)
1010   #endif                                              
1011                           else
1012   #ifdef IS_MPI                      
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)
1015 >                                 fpair, d_grp, rgrp, rCut)
1016   #else
1017                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1018                                   do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1019 <                                 d_grp, rgrp)
1019 >                                 d_grp, rgrp, rCut)
1020   #endif
1021                              vij = vij + vpair
1022                              fij(1:3) = fij(1:3) + fpair(1:3)
# Line 1164 | Line 1224 | contains
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, d_grp, r_grp)
1227 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1228  
1229      real( kind = dp ) :: vpair, sw
1230      real( kind = dp ), dimension(LR_POT_TYPES) :: pot
# Line 1181 | Line 1241 | contains
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  
# Line 1201 | Line 1262 | contains
1262      iHash = InteractionHash(me_i, me_j)
1263      
1264      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1265 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1265 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1266              pot(VDW_POT), f, do_pot)
1267      endif
1268      
1269      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1270 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1270 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1271              pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1272      endif
1273      
# Line 1226 | Line 1287 | contains
1287      endif
1288      
1289      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1290 <       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
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      
# Line 1244 | Line 1305 | contains
1305         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
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 ) :: sw
# Line 1259 | 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 1278 | 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
# Line 1291 | Line 1363 | contains
1363      if (FF_uses_EAM .and. SIM_uses_EAM) then
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 1474 | 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  
1556   #ifdef PROFILE

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines