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 3126 by gezelter, Fri Apr 6 21:53:43 2007 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.85 2007-04-06 21:53:41 gezelter Exp $, $Date: 2007-04-06 21:53:41 $, $Name: not supported by cvs2svn $, $Revision: 1.85 $
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 71 | Line 72 | module doForces
72    PRIVATE
73  
74   #define __FORTRAN90
74 – #include "UseTheForce/fSwitchingFunction.h"
75   #include "UseTheForce/fCutoffPolicy.h"
76   #include "UseTheForce/DarkSide/fInteractionMap.h"
77   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78  
79 –
79    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
80    INTEGER, PARAMETER:: PAIR_LOOP    = 2
81  
# Line 86 | Line 85 | module doForces
85    logical, save :: haveInteractionHash = .false.
86    logical, save :: haveGtypeCutoffMap = .false.
87    logical, save :: haveDefaultCutoffs = .false.
88 <  logical, save :: haveRlist = .false.
88 >  logical, save :: haveSkinThickness = .false.
89 >  logical, save :: haveElectrostaticSummationMethod = .false.
90 >  logical, save :: haveCutoffPolicy = .false.
91 >  logical, save :: VisitCutoffsAfterComputing = .false.
92 >  logical, save :: do_box_dipole = .false.
93  
94    logical, save :: FF_uses_DirectionalAtoms
95    logical, save :: FF_uses_Dipoles
96    logical, save :: FF_uses_GayBerne
97    logical, save :: FF_uses_EAM
98 +  logical, save :: FF_uses_SC
99 +  logical, save :: FF_uses_MEAM
100 +
101  
102    logical, save :: SIM_uses_DirectionalAtoms
103    logical, save :: SIM_uses_EAM
104 +  logical, save :: SIM_uses_SC
105 +  logical, save :: SIM_uses_MEAM
106    logical, save :: SIM_requires_postpair_calc
107    logical, save :: SIM_requires_prepair_calc
108    logical, save :: SIM_uses_PBC
109 +  logical, save :: SIM_uses_AtomicVirial
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 :: setBoxDipole
123 >  public :: getBoxDipole
124 >  public :: setCutoffPolicy
125 >  public :: setSkinThickness
126    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
127  
128   #ifdef PROFILE
129    public :: getforcetime
# Line 138 | Line 151 | module doForces
151    end type gtypeCutoffs
152    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
153  
154 <  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
155 <  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
143 <  real(kind=dp),save :: listSkin
144 <  
154 >  real(kind=dp), dimension(3) :: boxDipole
155 >
156   contains
157  
158 <  subroutine createInteractionHash(status)
158 >  subroutine createInteractionHash()
159      integer :: nAtypes
149 –    integer, intent(out) :: status
160      integer :: i
161      integer :: j
162      integer :: iHash
# Line 158 | Line 168 | contains
168      logical :: i_is_GB
169      logical :: i_is_EAM
170      logical :: i_is_Shape
171 +    logical :: i_is_SC
172 +    logical :: i_is_MEAM
173      logical :: j_is_LJ
174      logical :: j_is_Elect
175      logical :: j_is_Sticky
# Line 165 | Line 177 | contains
177      logical :: j_is_GB
178      logical :: j_is_EAM
179      logical :: j_is_Shape
180 +    logical :: j_is_SC
181 +    logical :: j_is_MEAM
182      real(kind=dp) :: myRcut
183  
170 –    status = 0  
171 –
184      if (.not. associated(atypes)) then
185 <       call handleError("atype", "atypes was not present before call of createInteractionHash!")
174 <       status = -1
185 >       call handleError("doForces", "atypes was not present before call of createInteractionHash!")
186         return
187      endif
188      
189      nAtypes = getSize(atypes)
190      
191      if (nAtypes == 0) then
192 <       status = -1
192 >       call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
193         return
194      end if
195  
# Line 204 | Line 215 | contains
215         call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
216         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
217         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
218 +       call getElementProperty(atypes, i, "is_SC", i_is_SC)
219 +       call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
220  
221         do j = i, nAtypes
222  
# Line 217 | Line 230 | contains
230            call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
231            call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
232            call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
233 +          call getElementProperty(atypes, j, "is_SC", j_is_SC)
234 +          call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
235  
236            if (i_is_LJ .and. j_is_LJ) then
237               iHash = ior(iHash, LJ_PAIR)            
# Line 238 | Line 253 | contains
253               iHash = ior(iHash, EAM_PAIR)
254            endif
255  
256 +          if (i_is_SC .and. j_is_SC) then
257 +             iHash = ior(iHash, SC_PAIR)
258 +          endif
259 +
260            if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
261            if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
262            if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
# Line 257 | Line 276 | contains
276      haveInteractionHash = .true.
277    end subroutine createInteractionHash
278  
279 <  subroutine createGtypeCutoffMap(stat)
279 >  subroutine createGtypeCutoffMap()
280  
262 –    integer, intent(out), optional :: stat
281      logical :: i_is_LJ
282      logical :: i_is_Elect
283      logical :: i_is_Sticky
# Line 267 | Line 285 | contains
285      logical :: i_is_GB
286      logical :: i_is_EAM
287      logical :: i_is_Shape
288 +    logical :: i_is_SC
289      logical :: GtypeFound
290  
291      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
# Line 274 | Line 293 | contains
293      integer :: nGroupsInRow
294      integer :: nGroupsInCol
295      integer :: nGroupTypesRow,nGroupTypesCol
296 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
296 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
297      real(kind=dp) :: biggestAtypeCutoff
298  
280 –    stat = 0
299      if (.not. haveInteractionHash) then
300 <       call createInteractionHash(myStatus)      
283 <       if (myStatus .ne. 0) then
284 <          write(default_error, *) 'createInteractionHash failed in doForces!'
285 <          stat = -1
286 <          return
287 <       endif
300 >       call createInteractionHash()      
301      endif
302   #ifdef IS_MPI
303      nGroupsInRow = getNgroupsInRow(plan_group_row)
# Line 302 | Line 315 | contains
315            call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
316            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
317            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
318 <          
318 >          call getElementProperty(atypes, i, "is_SC", i_is_SC)
319  
320            if (haveDefaultCutoffs) then
321               atypeMaxCutoff(i) = defaultRcut
# Line 335 | Line 348 | contains
348                  thisRcut = getShapeCut(i)
349                  if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
350               endif
351 +             if (i_is_SC) then
352 +                thisRcut = getSCCut(i)
353 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
354 +             endif
355            endif
356 <          
340 <          
356 >                    
357            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
358               biggestAtypeCutoff = atypeMaxCutoff(i)
359            endif
360  
361         endif
362      enddo
347 –  
348 –
363      
364      istart = 1
365      jstart = 1
# Line 389 | Line 403 | contains
403         allocate(groupToGtypeCol(jend))
404      end if
405  
406 <    if(.not.associated(groupToGtypeCol)) then
407 <       allocate(groupToGtypeCol(jend))
406 >    if(.not.associated(groupMaxCutoffCol)) then
407 >       allocate(groupMaxCutoffCol(jend))
408      else
409 <       deallocate(groupToGtypeCol)
410 <       allocate(groupToGtypeCol(jend))
409 >       deallocate(groupMaxCutoffCol)
410 >       allocate(groupMaxCutoffCol(jend))
411      end if
412      if(.not.associated(gtypeMaxCutoffCol)) then
413         allocate(gtypeMaxCutoffCol(jend))
# Line 414 | Line 428 | contains
428      !! largest cutoff for any atypes present in this group.  We also
429      !! create gtypes at this point.
430      
431 <    tol = 1.0d-6
431 >    tol = 1.0e-6_dp
432      nGroupTypesRow = 0
433 <
433 >    nGroupTypesCol = 0
434      do i = istart, iend      
435         n_in_i = groupStartRow(i+1) - groupStartRow(i)
436         groupMaxCutoffRow(i) = 0.0_dp
# Line 431 | Line 445 | contains
445               groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
446            endif          
447         enddo
434 –
448         if (nGroupTypesRow.eq.0) then
449            nGroupTypesRow = nGroupTypesRow + 1
450            gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
# Line 494 | Line 507 | contains
507      groupMaxCutoffCol => groupMaxCutoffRow
508   #endif
509  
497 –
498 –
499 –
500 –
510      !! allocate the gtypeCutoffMap here.
511      allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
512      !! then we do a double loop over all the group TYPES to find the cutoff
513      !! map between groups of two types
514      tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
515  
516 <    do i = 1, nGroupTypesRow
516 >    do i = 1, nGroupTypesRow      
517         do j = 1, nGroupTypesCol
518        
519            select case(cutoffPolicy)
# Line 519 | Line 528 | contains
528               return
529            end select
530            gtypeCutoffMap(i,j)%rcut = thisRcut
531 +          
532 +          if (thisRcut.gt.largestRcut) largestRcut = thisRcut
533 +
534            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
535  
536 +          if (.not.haveSkinThickness) then
537 +             skinThickness = 1.0_dp
538 +          endif
539 +
540 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
541 +
542            ! sanity check
543  
544            if (haveDefaultCutoffs) then
# Line 533 | Line 548 | contains
548            endif
549         enddo
550      enddo
551 +
552      if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
553      if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
554      if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
# Line 546 | Line 562 | contains
562      haveGtypeCutoffMap = .true.
563     end subroutine createGtypeCutoffMap
564  
565 <   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
550 <     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
551 <     integer, intent(in) :: cutPolicy
565 >   subroutine setCutoffs(defRcut, defRsw)
566  
567 +     real(kind=dp),intent(in) :: defRcut, defRsw
568 +     character(len = statusMsgSize) :: errMsg
569 +     integer :: localError
570 +
571       defaultRcut = defRcut
572       defaultRsw = defRsw
573 <     defaultRlist = defRlist
574 <     cutoffPolicy = cutPolicy
575 <
573 >    
574 >     defaultDoShift = .false.
575 >     if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
576 >        
577 >        write(errMsg, *) &
578 >             'cutoffRadius and switchingRadius are set to the same', newline &
579 >             // tab, 'value.  OOPSE will use shifted ', newline &
580 >             // tab, 'potentials instead of switching functions.'
581 >        
582 >        call handleInfo("setCutoffs", errMsg)
583 >        
584 >        defaultDoShift = .true.
585 >        
586 >     endif
587 >    
588 >     localError = 0
589 >     call setLJDefaultCutoff( defaultRcut, defaultDoShift )
590 >     call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
591 >     call setCutoffEAM( defaultRcut )
592 >     call setCutoffSC( defaultRcut )
593 >     call set_switch(defaultRsw, defaultRcut)
594 >     call setHmatDangerousRcutValue(defaultRcut)
595 >        
596       haveDefaultCutoffs = .true.
597 <   end subroutine setDefaultCutoffs
597 >     haveGtypeCutoffMap = .false.
598  
599 <   subroutine setCutoffPolicy(cutPolicy)
599 >   end subroutine setCutoffs
600  
601 +   subroutine cWasLame()
602 +    
603 +     VisitCutoffsAfterComputing = .true.
604 +     return
605 +    
606 +   end subroutine cWasLame
607 +  
608 +   subroutine setCutoffPolicy(cutPolicy)
609 +    
610       integer, intent(in) :: cutPolicy
611 +    
612       cutoffPolicy = cutPolicy
613 <     call createGtypeCutoffMap()
613 >     haveCutoffPolicy = .true.
614 >     haveGtypeCutoffMap = .false.
615 >    
616     end subroutine setCutoffPolicy
567 –    
617      
618 <  subroutine setSimVariables()
570 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
571 <    SIM_uses_EAM = SimUsesEAM()
572 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574 <    SIM_uses_PBC = SimUsesPBC()
618 >   subroutine setBoxDipole()
619  
620 <    haveSIMvariables = .true.
620 >     do_box_dipole = .true.
621 >    
622 >   end subroutine setBoxDipole
623  
624 <    return
579 <  end subroutine setSimVariables
624 >   subroutine getBoxDipole( box_dipole )
625  
626 +     real(kind=dp), intent(inout), dimension(3) :: box_dipole
627 +
628 +     box_dipole = boxDipole
629 +
630 +   end subroutine getBoxDipole
631 +
632 +   subroutine setElectrostaticMethod( thisESM )
633 +
634 +     integer, intent(in) :: thisESM
635 +
636 +     electrostaticSummationMethod = thisESM
637 +     haveElectrostaticSummationMethod = .true.
638 +    
639 +   end subroutine setElectrostaticMethod
640 +
641 +   subroutine setSkinThickness( thisSkin )
642 +    
643 +     real(kind=dp), intent(in) :: thisSkin
644 +    
645 +     skinThickness = thisSkin
646 +     haveSkinThickness = .true.    
647 +     haveGtypeCutoffMap = .false.
648 +    
649 +   end subroutine setSkinThickness
650 +      
651 +   subroutine setSimVariables()
652 +     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
653 +     SIM_uses_EAM = SimUsesEAM()
654 +     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
655 +     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
656 +     SIM_uses_PBC = SimUsesPBC()
657 +     SIM_uses_SC = SimUsesSC()
658 +     SIM_uses_AtomicVirial = SimUsesAtomicVirial()
659 +
660 +     haveSIMvariables = .true.
661 +    
662 +     return
663 +   end subroutine setSimVariables
664 +
665    subroutine doReadyCheck(error)
666      integer, intent(out) :: error
583 –
667      integer :: myStatus
668  
669      error = 0
670  
671      if (.not. haveInteractionHash) then      
672 <       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
672 >       call createInteractionHash()      
673      endif
674  
675      if (.not. haveGtypeCutoffMap) then        
676 <       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
676 >       call createGtypeCutoffMap()      
677      endif
678  
679 +    if (VisitCutoffsAfterComputing) then
680 +       call set_switch(largestRcut, largestRcut)      
681 +       call setHmatDangerousRcutValue(largestRcut)
682 +       call setCutoffEAM(largestRcut)
683 +       call setCutoffSC(largestRcut)
684 +       VisitCutoffsAfterComputing = .false.
685 +    endif
686 +
687      if (.not. haveSIMvariables) then
688         call setSimVariables()
689      endif
690  
612 –  !  if (.not. haveRlist) then
613 –  !     write(default_error, *) 'rList has not been set in doForces!'
614 –  !     error = -1
615 –  !     return
616 –  !  endif
617 –
691      if (.not. haveNeighborList) then
692         write(default_error, *) 'neighbor list has not been initialized in doForces!'
693         error = -1
694         return
695      end if
696 <
696 >    
697      if (.not. haveSaneForceField) then
698         write(default_error, *) 'Force Field is not sane in doForces!'
699         error = -1
700         return
701      end if
702 <
702 >    
703   #ifdef IS_MPI
704      if (.not. isMPISimSet()) then
705         write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
# Line 638 | Line 711 | contains
711    end subroutine doReadyCheck
712  
713  
714 <  subroutine init_FF(thisESM, thisStat)
714 >  subroutine init_FF(thisStat)
715  
643 –    integer, intent(in) :: thisESM
716      integer, intent(out) :: thisStat  
717      integer :: my_status, nMatches
718      integer, pointer :: MatchList(:) => null()
# Line 648 | Line 720 | contains
720      !! assume things are copacetic, unless they aren't
721      thisStat = 0
722  
651 –    electrostaticSummationMethod = thisESM
652 –
723      !! init_FF is called *after* all of the atom types have been
724      !! defined in atype_module using the new_atype subroutine.
725      !!
# Line 660 | Line 730 | contains
730      FF_uses_Dipoles = .false.
731      FF_uses_GayBerne = .false.
732      FF_uses_EAM = .false.
733 +    FF_uses_SC = .false.
734  
735      call getMatchingElementList(atypes, "is_Directional", .true., &
736           nMatches, MatchList)
# Line 676 | Line 747 | contains
747      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
748      if (nMatches .gt. 0) FF_uses_EAM = .true.
749  
750 +    call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
751 +    if (nMatches .gt. 0) FF_uses_SC = .true.
752  
753 +
754      haveSaneForceField = .true.
755  
756      if (FF_uses_EAM) then
# Line 741 | Line 815 | contains
815      integer :: istart, iend
816      integer :: ia, jb, atom1, atom2
817      integer :: nlist
818 <    real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
818 >    real( kind = DP ) :: ratmsq, rgrpsq, rgrp, rag, vpair, vij
819      real( kind = DP ) :: sw, dswdr, swderiv, mf
820      real( kind = DP ) :: rVal
821 <    real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
822 <    real(kind=dp) :: rfpot, mu_i, virial
821 >    real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij, fg, dag
822 >    real(kind=dp) :: rfpot, mu_i
823 >    real(kind=dp):: rCut
824      integer :: me_i, me_j, n_in_i, n_in_j
825      logical :: is_dp_i
826      integer :: neighborListSize
# Line 755 | Line 830 | contains
830      integer :: loopStart, loopEnd, loop
831      integer :: iHash
832      integer :: i1
758 –  
833  
834 +    !! the variables for the box dipole moment
835 + #ifdef IS_MPI
836 +    integer :: pChgCount_local
837 +    integer :: nChgCount_local
838 +    real(kind=dp) :: pChg_local
839 +    real(kind=dp) :: nChg_local
840 +    real(kind=dp), dimension(3) :: pChgPos_local
841 +    real(kind=dp), dimension(3) :: nChgPos_local
842 +    real(kind=dp), dimension(3) :: dipVec_local
843 + #endif
844 +    integer :: pChgCount
845 +    integer :: nChgCount
846 +    real(kind=dp) :: pChg
847 +    real(kind=dp) :: nChg
848 +    real(kind=dp) :: chg_value
849 +    real(kind=dp), dimension(3) :: pChgPos
850 +    real(kind=dp), dimension(3) :: nChgPos
851 +    real(kind=dp), dimension(3) :: dipVec
852 +    real(kind=dp), dimension(3) :: chgVec
853 +
854 +    !! initialize box dipole variables
855 +    if (do_box_dipole) then
856 + #ifdef IS_MPI
857 +       pChg_local = 0.0_dp
858 +       nChg_local = 0.0_dp
859 +       pChgCount_local = 0
860 +       nChgCount_local = 0
861 +       do i=1, 3
862 +          pChgPos_local = 0.0_dp
863 +          nChgPos_local = 0.0_dp
864 +          dipVec_local = 0.0_dp
865 +       enddo
866 + #endif
867 +       pChg = 0.0_dp
868 +       nChg = 0.0_dp
869 +       pChgCount = 0
870 +       nChgCount = 0
871 +       chg_value = 0.0_dp
872 +      
873 +       do i=1, 3
874 +          pChgPos(i) = 0.0_dp
875 +          nChgPos(i) = 0.0_dp
876 +          dipVec(i) = 0.0_dp
877 +          chgVec(i) = 0.0_dp
878 +          boxDipole(i) = 0.0_dp
879 +       enddo
880 +    endif
881 +
882      !! initialize local variables  
883  
884   #ifdef IS_MPI
# Line 819 | Line 941 | contains
941         ! (but only on the first time through):
942         if (loop .eq. loopStart) then
943   #ifdef IS_MPI
944 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
944 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
945                 update_nlist)
946   #else
947 <          call checkNeighborList(nGroups, q_group, listSkin, &
947 >          call checkNeighborList(nGroups, q_group, skinThickness, &
948                 update_nlist)
949   #endif
950         endif
# Line 902 | Line 1024 | contains
1024  
1025                     list(nlist) = j
1026                  endif
1027 <                
1027 >                
1028                  if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
1029  
1030 +                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
1031                     if (loop .eq. PAIR_LOOP) then
1032 <                      vij = 0.0d0
1033 <                      fij(1:3) = 0.0d0
1032 >                      vij = 0.0_dp
1033 >                      fij(1) = 0.0_dp
1034 >                      fij(2) = 0.0_dp
1035 >                      fij(3) = 0.0_dp
1036                     endif
1037                    
1038 <                   call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
914 <                        in_switching_region)
1038 >                   call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
1039                    
1040                     n_in_j = groupStartCol(j+1) - groupStartCol(j)
1041                    
# Line 926 | Line 1050 | contains
1050                           if (skipThisPair(atom1, atom2))  cycle inner
1051                          
1052                           if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1053 <                            d_atm(1:3) = d_grp(1:3)
1053 >                            d_atm(1) = d_grp(1)
1054 >                            d_atm(2) = d_grp(2)
1055 >                            d_atm(3) = d_grp(3)
1056                              ratmsq = rgrpsq
1057                           else
1058   #ifdef IS_MPI
# Line 941 | Line 1067 | contains
1067                           if (loop .eq. PREPAIR_LOOP) then
1068   #ifdef IS_MPI                      
1069                              call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1070 <                                 rgrpsq, d_grp, do_pot, do_stress, &
1070 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1071                                   eFrame, A, f, t, pot_local)
1072   #else
1073                              call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1074 <                                 rgrpsq, d_grp, do_pot, do_stress, &
1074 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1075                                   eFrame, A, f, t, pot)
1076   #endif                                              
1077                           else
1078   #ifdef IS_MPI                      
1079                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1080                                   do_pot, eFrame, A, f, t, pot_local, vpair, &
1081 <                                 fpair, d_grp, rgrp)
1081 >                                 fpair, d_grp, rgrp, rCut)
1082   #else
1083                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1084                                   do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1085 <                                 d_grp, rgrp)
1085 >                                 d_grp, rgrp, rCut)
1086   #endif
1087                              vij = vij + vpair
1088 <                            fij(1:3) = fij(1:3) + fpair(1:3)
1088 >                            fij(1) = fij(1) + fpair(1)
1089 >                            fij(2) = fij(2) + fpair(2)
1090 >                            fij(3) = fij(3) + fpair(3)
1091 >                            if (do_stress.and.SIM_uses_AtomicVirial) then
1092 >                               call add_stress_tensor(d_atm, fpair, tau)
1093 >                            endif
1094                           endif
1095                        enddo inner
1096                     enddo
# Line 974 | Line 1105 | contains
1105                           do ia=groupStartRow(i), groupStartRow(i+1)-1
1106                              atom1=groupListRow(ia)
1107                              mf = mfactRow(atom1)
1108 +                            ! fg is the force on atom ia due to cutoff group's
1109 +                            ! presence in switching region
1110 +                            fg = swderiv*d_grp*mf
1111   #ifdef IS_MPI
1112 <                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1113 <                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1114 <                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1112 >                            f_Row(1,atom1) = f_Row(1,atom1) + fg(1)
1113 >                            f_Row(2,atom1) = f_Row(2,atom1) + fg(2)
1114 >                            f_Row(3,atom1) = f_Row(3,atom1) + fg(3)
1115   #else
1116 <                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1117 <                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1118 <                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1116 >                            f(1,atom1) = f(1,atom1) + fg(1)
1117 >                            f(2,atom1) = f(2,atom1) + fg(2)
1118 >                            f(3,atom1) = f(3,atom1) + fg(3)
1119   #endif
1120 +                            if (do_stress.and.SIM_uses_AtomicVirial) then
1121 +                               ! find the distance between the atom and the center of
1122 +                               ! the cutoff group:
1123 + #ifdef IS_MPI
1124 +                               call get_interatomic_vector(q_Row(:,atom1), &
1125 +                                    q_group_Row(:,i), dag, rag)
1126 + #else
1127 +                               call get_interatomic_vector(q(:,atom1), &
1128 +                                    q_group(:,i), dag, rag)
1129 + #endif
1130 +                               call add_stress_tensor(dag,fg,tau)
1131 +                            endif
1132                           enddo
1133                          
1134                           do jb=groupStartCol(j), groupStartCol(j+1)-1
1135                              atom2=groupListCol(jb)
1136                              mf = mfactCol(atom2)
1137 +                            ! fg is the force on atom jb due to cutoff group's
1138 +                            ! presence in switching region
1139 +                            fg = -swderiv*d_grp*mf
1140   #ifdef IS_MPI
1141 <                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1142 <                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1143 <                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1141 >                            f_Col(1,atom2) = f_Col(1,atom2) + fg(1)
1142 >                            f_Col(2,atom2) = f_Col(2,atom2) + fg(2)
1143 >                            f_Col(3,atom2) = f_Col(3,atom2) + fg(3)
1144   #else
1145 <                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1146 <                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1147 <                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1145 >                            f(1,atom2) = f(1,atom2) + fg(1)
1146 >                            f(2,atom2) = f(2,atom2) + fg(2)
1147 >                            f(3,atom2) = f(3,atom2) + fg(3)
1148   #endif
1149 +                            if (do_stress.and.SIM_uses_AtomicVirial) then
1150 +                               ! find the distance between the atom and the center of
1151 +                               ! the cutoff group:
1152 + #ifdef IS_MPI
1153 +                               call get_interatomic_vector(q_Col(:,atom2), &
1154 +                                    q_group_Col(:,j), dag, rag)
1155 + #else
1156 +                               call get_interatomic_vector(q(:,atom2), &
1157 +                                    q_group(:,j), dag, rag)
1158 + #endif
1159 +                               call add_stress_tensor(dag,fg,tau)                              
1160 +                            endif
1161 +                            
1162                           enddo
1163                        endif
1164 <
1165 <                      if (do_stress) call add_stress_tensor(d_grp, fij)
1164 >                      if (do_stress.and.(.not.SIM_uses_AtomicVirial)) then
1165 >                         call add_stress_tensor(d_grp, fij, tau)
1166 >                      endif
1167                     endif
1168                  endif
1169               endif
# Line 1120 | Line 1283 | contains
1283                  
1284                  ! prevent overcounting of the skips
1285                  if (i.lt.j) then
1286 <                   call get_interatomic_vector(q(:,i), &
1287 <                        q(:,j), d_atm, ratmsq)
1288 <                   rVal = dsqrt(ratmsq)
1126 <                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1127 <                        in_switching_region)
1286 >                   call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1287 >                   rVal = sqrt(ratmsq)
1288 >                   call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1289   #ifdef IS_MPI
1290                     call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1291                          vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
# Line 1135 | Line 1296 | contains
1296                  endif
1297               enddo
1298            endif
1299 +
1300 +          if (do_box_dipole) then
1301 + #ifdef IS_MPI
1302 +             call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1303 +                  nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1304 +                  pChgCount_local, nChgCount_local)
1305 + #else
1306 +             call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1307 +                  pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1308 + #endif
1309 +          endif
1310         enddo
1311      endif
1312 <    
1312 >
1313   #ifdef IS_MPI
1142 –    
1314      if (do_pot) then
1315 <       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1315 > #ifdef SINGLE_PRECISION
1316 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1317              mpi_comm_world,mpi_err)            
1318 + #else
1319 +       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1320 +            mpi_sum, mpi_comm_world,mpi_err)            
1321 + #endif
1322      endif
1323 <    
1324 <    if (do_stress) then
1325 <       call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1326 <            mpi_comm_world,mpi_err)
1327 <       call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1328 <            mpi_comm_world,mpi_err)
1329 <    endif
1330 <    
1323 >        
1324 >    if (do_box_dipole) then
1325 >
1326 > #ifdef SINGLE_PRECISION
1327 >       call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1328 >            mpi_comm_world, mpi_err)
1329 >       call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1330 >            mpi_comm_world, mpi_err)
1331 >       call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1332 >            mpi_comm_world, mpi_err)
1333 >       call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1334 >            mpi_comm_world, mpi_err)
1335 >       call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1336 >            mpi_comm_world, mpi_err)
1337 >       call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1338 >            mpi_comm_world, mpi_err)
1339 >       call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1340 >            mpi_comm_world, mpi_err)
1341   #else
1342 <    
1343 <    if (do_stress) then
1344 <       tau = tau_Temp
1345 <       virial = virial_Temp
1342 >       call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1343 >            mpi_comm_world, mpi_err)
1344 >       call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1345 >            mpi_comm_world, mpi_err)
1346 >       call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1347 >            mpi_sum, mpi_comm_world, mpi_err)
1348 >       call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1349 >            mpi_sum, mpi_comm_world, mpi_err)
1350 >       call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1351 >            mpi_sum, mpi_comm_world, mpi_err)
1352 >       call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1353 >            mpi_sum, mpi_comm_world, mpi_err)
1354 >       call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1355 >            mpi_sum, mpi_comm_world, mpi_err)
1356 > #endif
1357 >
1358      endif
1359      
1360   #endif
1361 <    
1361 >
1362 >    if (do_box_dipole) then
1363 >       ! first load the accumulated dipole moment (if dipoles were present)
1364 >       boxDipole(1) = dipVec(1)
1365 >       boxDipole(2) = dipVec(2)
1366 >       boxDipole(3) = dipVec(3)
1367 >
1368 >       ! now include the dipole moment due to charges
1369 >       ! use the lesser of the positive and negative charge totals
1370 >       if (nChg .le. pChg) then
1371 >          chg_value = nChg
1372 >       else
1373 >          chg_value = pChg
1374 >       endif
1375 >      
1376 >       ! find the average positions
1377 >       if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1378 >          pChgPos = pChgPos / pChgCount
1379 >          nChgPos = nChgPos / nChgCount
1380 >       endif
1381 >
1382 >       ! dipole is from the negative to the positive (physics notation)
1383 >       chgVec(1) = pChgPos(1) - nChgPos(1)
1384 >       chgVec(2) = pChgPos(2) - nChgPos(2)
1385 >       chgVec(3) = pChgPos(3) - nChgPos(3)
1386 >
1387 >       boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1388 >       boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1389 >       boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1390 >
1391 >    endif
1392 >
1393    end subroutine do_force_loop
1394  
1395    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1396 <       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1396 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1397  
1398      real( kind = dp ) :: vpair, sw
1399      real( kind = dp ), dimension(LR_POT_TYPES) :: pot
# Line 1181 | Line 1410 | contains
1410      real ( kind = dp ), intent(inout) :: r_grp
1411      real ( kind = dp ), intent(inout) :: d(3)
1412      real ( kind = dp ), intent(inout) :: d_grp(3)
1413 +    real ( kind = dp ), intent(inout) :: rCut
1414      real ( kind = dp ) :: r
1415 +    real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1416      integer :: me_i, me_j
1417 +    integer :: k
1418  
1419      integer :: iHash
1420  
1421      r = sqrt(rijsq)
1422 <    vpair = 0.0d0
1423 <    fpair(1:3) = 0.0d0
1422 >    
1423 >    vpair = 0.0_dp
1424 >    fpair(1:3) = 0.0_dp
1425  
1426   #ifdef IS_MPI
1427      me_i = atid_row(i)
# Line 1201 | Line 1434 | contains
1434      iHash = InteractionHash(me_i, me_j)
1435      
1436      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1437 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1437 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1438              pot(VDW_POT), f, do_pot)
1439      endif
1440      
1441      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1442 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1442 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1443              pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1444      endif
1445      
# Line 1226 | Line 1459 | contains
1459      endif
1460      
1461      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1462 <       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1462 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1463              pot(VDW_POT), A, f, t, do_pot)
1464      endif
1465      
# Line 1244 | Line 1477 | contains
1477         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1478              pot(VDW_POT), A, f, t, do_pot)
1479      endif
1480 +
1481 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1482 +       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1483 +            pot(METALLIC_POT), f, do_pot)
1484 +    endif
1485      
1486    end subroutine do_pair
1487  
1488 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1488 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1489         do_pot, do_stress, eFrame, A, f, t, pot)
1490  
1491      real( kind = dp ) :: sw
# Line 1259 | Line 1497 | contains
1497  
1498      logical, intent(inout) :: do_pot, do_stress
1499      integer, intent(in) :: i, j
1500 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1500 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1501      real ( kind = dp )                :: r, rc
1502      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1503  
1504      integer :: me_i, me_j, iHash
1505  
1506      r = sqrt(rijsq)
1507 <
1507 >    
1508   #ifdef IS_MPI  
1509      me_i = atid_row(i)
1510      me_j = atid_col(j)  
# Line 1278 | Line 1516 | contains
1516      iHash = InteractionHash(me_i, me_j)
1517  
1518      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1519 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1519 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1520      endif
1521 +
1522 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1523 +            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1524 +    endif
1525      
1526    end subroutine do_prepair
1527  
# Line 1291 | Line 1533 | contains
1533      if (FF_uses_EAM .and. SIM_uses_EAM) then
1534         call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1535      endif
1536 <
1537 <
1536 >    if (FF_uses_SC .and. SIM_uses_SC) then
1537 >       call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1538 >    endif
1539    end subroutine do_preforce
1540  
1541  
# Line 1304 | Line 1547 | contains
1547      real( kind = dp ) :: d(3), scaled(3)
1548      integer i
1549  
1550 <    d(1:3) = q_j(1:3) - q_i(1:3)
1550 >    d(1) = q_j(1) - q_i(1)
1551 >    d(2) = q_j(2) - q_i(2)
1552 >    d(3) = q_j(3) - q_i(3)
1553  
1554      ! Wrap back into periodic box if necessary
1555      if ( SIM_uses_PBC ) then
1556  
1557         if( .not.boxIsOrthorhombic ) then
1558            ! calc the scaled coordinates.
1559 +          ! scaled = matmul(HmatInv, d)
1560  
1561 <          scaled = matmul(HmatInv, d)
1562 <
1561 >          scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1562 >          scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1563 >          scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1564 >          
1565            ! wrap the scaled coordinates
1566  
1567 <          scaled = scaled  - anint(scaled)
1567 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1568 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1569 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1570  
1321 –
1571            ! calc the wrapped real coordinates from the wrapped scaled
1572            ! coordinates
1573 +          ! d = matmul(Hmat,scaled)
1574 +          d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1575 +          d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1576 +          d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1577  
1325 –          d = matmul(Hmat,scaled)
1326 –
1578         else
1579            ! calc the scaled coordinates.
1580  
1581 <          do i = 1, 3
1582 <             scaled(i) = d(i) * HmatInv(i,i)
1581 >          scaled(1) = d(1) * HmatInv(1,1)
1582 >          scaled(2) = d(2) * HmatInv(2,2)
1583 >          scaled(3) = d(3) * HmatInv(3,3)
1584 >          
1585 >          ! wrap the scaled coordinates
1586 >          
1587 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1588 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1589 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1590  
1591 <             ! wrap the scaled coordinates
1591 >          ! calc the wrapped real coordinates from the wrapped scaled
1592 >          ! coordinates
1593  
1594 <             scaled(i) = scaled(i) - anint(scaled(i))
1594 >          d(1) = scaled(1)*Hmat(1,1)
1595 >          d(2) = scaled(2)*Hmat(2,2)
1596 >          d(3) = scaled(3)*Hmat(3,3)
1597  
1337 –             ! calc the wrapped real coordinates from the wrapped scaled
1338 –             ! coordinates
1339 –
1340 –             d(i) = scaled(i)*Hmat(i,i)
1341 –          enddo
1598         endif
1599  
1600      endif
1601  
1602 <    r_sq = dot_product(d,d)
1602 >    r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1603  
1604    end subroutine get_interatomic_vector
1605  
# Line 1381 | Line 1637 | contains
1637         call clean_EAM()
1638      endif
1639  
1384 –    tau_Temp = 0.0_dp
1385 –    virial_Temp = 0.0_dp
1640    end subroutine zero_work_arrays
1641  
1642    function skipThisPair(atom1, atom2) result(skip_it)
# Line 1474 | Line 1728 | contains
1728  
1729    function FF_RequiresPrepairCalc() result(doesit)
1730      logical :: doesit
1731 <    doesit = FF_uses_EAM
1731 >    doesit = FF_uses_EAM .or. FF_uses_SC &
1732 >         .or. FF_uses_MEAM
1733    end function FF_RequiresPrepairCalc
1734  
1735   #ifdef PROFILE
# Line 1486 | Line 1741 | contains
1741  
1742    !! This cleans componets of force arrays belonging only to fortran
1743  
1744 <  subroutine add_stress_tensor(dpair, fpair)
1744 >  subroutine add_stress_tensor(dpair, fpair, tau)
1745  
1746      real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1747 +    real( kind = dp ), dimension(9), intent(inout) :: tau
1748  
1749      ! because the d vector is the rj - ri vector, and
1750      ! because fx, fy, fz are the force on atom i, we need a
1751      ! negative sign here:  
1752  
1753 <    tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1754 <    tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1755 <    tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1756 <    tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1757 <    tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1758 <    tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1759 <    tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1760 <    tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1761 <    tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1753 >    tau(1) = tau(1) - dpair(1) * fpair(1)
1754 >    tau(2) = tau(2) - dpair(1) * fpair(2)
1755 >    tau(3) = tau(3) - dpair(1) * fpair(3)
1756 >    tau(4) = tau(4) - dpair(2) * fpair(1)
1757 >    tau(5) = tau(5) - dpair(2) * fpair(2)
1758 >    tau(6) = tau(6) - dpair(2) * fpair(3)
1759 >    tau(7) = tau(7) - dpair(3) * fpair(1)
1760 >    tau(8) = tau(8) - dpair(3) * fpair(2)
1761 >    tau(9) = tau(9) - dpair(3) * fpair(3)
1762  
1507 –    virial_Temp = virial_Temp + &
1508 –         (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1509 –
1763    end subroutine add_stress_tensor
1764  
1765   end module doForces

Diff Legend

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