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 2355 by chuckv, Wed Oct 12 18:59:16 2005 UTC vs.
Revision 2432 by chuckv, Tue Nov 15 16:01:06 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.54 2005-10-12 18:59:16 chuckv Exp $, $Date: 2005-10-12 18:59:16 $, $Name: not supported by cvs2svn $, $Revision: 1.54 $
48 > !! @version $Id: doForces.F90,v 1.68 2005-11-15 16:01:06 chuckv Exp $, $Date: 2005-11-15 16:01:06 $, $Name: not supported by cvs2svn $, $Revision: 1.68 $
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 93 | Line 93 | module doForces
93    logical, save :: FF_uses_Dipoles
94    logical, save :: FF_uses_GayBerne
95    logical, save :: FF_uses_EAM
96 +  logical, save :: FF_uses_SC
97 +  logical, save :: FF_uses_MEAM
98 +
99  
100    logical, save :: SIM_uses_DirectionalAtoms
101    logical, save :: SIM_uses_EAM
102 +  logical, save :: SIM_uses_SC
103 +  logical, save :: SIM_uses_MEAM
104    logical, save :: SIM_requires_postpair_calc
105    logical, save :: SIM_requires_prepair_calc
106    logical, save :: SIM_uses_PBC
# 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  
180 +
181      status = 0  
182  
183      if (.not. associated(atypes)) then
# Line 205 | 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 218 | 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 237 | Line 251 | contains
251  
252            if (i_is_EAM .and. j_is_EAM) then
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)
# Line 570 | Line 588 | contains
588    subroutine setSimVariables()
589      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
590      SIM_uses_EAM = SimUsesEAM()
591 +    SIM_uses_SC  = SimUsesSC()
592      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
593      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
594      SIM_uses_PBC = SimUsesPBC()
# Line 645 | Line 664 | contains
664      integer, intent(out) :: thisStat  
665      integer :: my_status, nMatches
666      integer, pointer :: MatchList(:) => null()
648 –    real(kind=dp) :: rcut, rrf, rt, dielect
667  
668      !! assume things are copacetic, unless they aren't
669      thisStat = 0
# Line 680 | Line 698 | contains
698  
699  
700      haveSaneForceField = .true.
683 –
684 –    !! check to make sure the reaction field setting makes sense
685 –
686 –    if (FF_uses_Dipoles) then
687 –       if (electrostaticSummationMethod == REACTION_FIELD) then
688 –          dielect = getDielect()
689 –          call initialize_rf(dielect)
690 –       endif
691 –    else
692 –       if (electrostaticSummationMethod == REACTION_FIELD) then
693 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
694 –          thisStat = -1
695 –          haveSaneForceField = .false.
696 –          return
697 –       endif
698 –    endif
701  
702      if (FF_uses_EAM) then
703         call init_EAM_FF(my_status)
# Line 707 | Line 709 | contains
709         end if
710      endif
711  
710 –    if (FF_uses_GayBerne) then
711 –       call check_gb_pair_FF(my_status)
712 –       if (my_status .ne. 0) then
713 –          thisStat = -1
714 –          haveSaneForceField = .false.
715 –          return
716 –       endif
717 –    endif
718 –
712      if (.not. haveNeighborList) then
713         !! Create neighbor lists
714         call expandNeighborList(nLocal, my_status)
# Line 749 | Line 742 | contains
742  
743      !! Stress Tensor
744      real( kind = dp), dimension(9) :: tau  
745 <    real ( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
745 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
746      logical ( kind = 2) :: do_pot_c, do_stress_c
747      logical :: do_pot
748      logical :: do_stress
749      logical :: in_switching_region
750   #ifdef IS_MPI
751 <    real( kind = DP ), dimension(POT_ARRAY_SIZE) :: pot_local
751 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
752      integer :: nAtomsInRow
753      integer :: nAtomsInCol
754      integer :: nprocs
# Line 770 | Line 763 | contains
763      integer :: nlist
764      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
765      real( kind = DP ) :: sw, dswdr, swderiv, mf
766 +    real( kind = DP ) :: rVal
767      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
768      real(kind=dp) :: rfpot, mu_i, virial
769      integer :: me_i, me_j, n_in_i, n_in_j
# Line 780 | Line 774 | contains
774      integer :: propPack_i, propPack_j
775      integer :: loopStart, loopEnd, loop
776      integer :: iHash
777 +    integer :: i1
778    
779  
780      !! initialize local variables  
# Line 927 | Line 922 | contains
922  
923                     list(nlist) = j
924                  endif
925 <
926 <                if (loop .eq. PAIR_LOOP) then
932 <                   vij = 0.0d0
933 <                   fij(1:3) = 0.0d0
934 <                endif
935 <
936 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
937 <                     in_switching_region)
938 <
939 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
940 <
941 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
925 >                
926 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
927  
928 <                   atom1 = groupListRow(ia)
929 <
930 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
931 <
932 <                      atom2 = groupListCol(jb)
933 <
934 <                      if (skipThisPair(atom1, atom2)) cycle inner
935 <
936 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
937 <                         d_atm(1:3) = d_grp(1:3)
938 <                         ratmsq = rgrpsq
939 <                      else
928 >                   if (loop .eq. PAIR_LOOP) then
929 >                      vij = 0.0d0
930 >                      fij(1:3) = 0.0d0
931 >                   endif
932 >                  
933 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
934 >                        in_switching_region)
935 >                  
936 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
937 >                  
938 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
939 >                      
940 >                      atom1 = groupListRow(ia)
941 >                      
942 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
943 >                        
944 >                         atom2 = groupListCol(jb)
945 >                        
946 >                         if (skipThisPair(atom1, atom2))  cycle inner
947 >                        
948 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
949 >                            d_atm(1:3) = d_grp(1:3)
950 >                            ratmsq = rgrpsq
951 >                         else
952   #ifdef IS_MPI
953 <                         call get_interatomic_vector(q_Row(:,atom1), &
954 <                              q_Col(:,atom2), d_atm, ratmsq)
953 >                            call get_interatomic_vector(q_Row(:,atom1), &
954 >                                 q_Col(:,atom2), d_atm, ratmsq)
955   #else
956 <                         call get_interatomic_vector(q(:,atom1), &
957 <                              q(:,atom2), d_atm, ratmsq)
956 >                            call get_interatomic_vector(q(:,atom1), &
957 >                                 q(:,atom2), d_atm, ratmsq)
958   #endif
959 <                      endif
960 <
961 <                      if (loop .eq. PREPAIR_LOOP) then
959 >                         endif
960 >                        
961 >                         if (loop .eq. PREPAIR_LOOP) then
962   #ifdef IS_MPI                      
963 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
964 <                              rgrpsq, d_grp, do_pot, do_stress, &
965 <                              eFrame, A, f, t, pot_local)
963 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
964 >                                 rgrpsq, d_grp, do_pot, do_stress, &
965 >                                 eFrame, A, f, t, pot_local)
966   #else
967 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
968 <                              rgrpsq, d_grp, do_pot, do_stress, &
969 <                              eFrame, A, f, t, pot)
970 < #endif                                              
971 <                      else
967 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
968 >                                 rgrpsq, d_grp, do_pot, do_stress, &
969 >                                 eFrame, A, f, t, pot)
970 > #endif                                              
971 >                         else
972   #ifdef IS_MPI                      
973 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
974 <                              do_pot, &
975 <                              eFrame, A, f, t, pot_local, vpair, fpair)
973 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
974 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
975 >                                 fpair, d_grp, rgrp)
976   #else
977 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
978 <                              do_pot,  &
979 <                              eFrame, A, f, t, pot, vpair, fpair)
977 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
978 >                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
979 >                                 d_grp, rgrp)
980   #endif
981 +                            vij = vij + vpair
982 +                            fij(1:3) = fij(1:3) + fpair(1:3)
983 +                         endif
984 +                      enddo inner
985 +                   enddo
986  
987 <                         vij = vij + vpair
988 <                         fij(1:3) = fij(1:3) + fpair(1:3)
989 <                      endif
990 <                   enddo inner
991 <                enddo
992 <
993 <                if (loop .eq. PAIR_LOOP) then
994 <                   if (in_switching_region) then
995 <                      swderiv = vij*dswdr/rgrp
996 <                      fij(1) = fij(1) + swderiv*d_grp(1)
995 <                      fij(2) = fij(2) + swderiv*d_grp(2)
996 <                      fij(3) = fij(3) + swderiv*d_grp(3)
997 <
998 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
999 <                         atom1=groupListRow(ia)
1000 <                         mf = mfactRow(atom1)
987 >                   if (loop .eq. PAIR_LOOP) then
988 >                      if (in_switching_region) then
989 >                         swderiv = vij*dswdr/rgrp
990 >                         fij(1) = fij(1) + swderiv*d_grp(1)
991 >                         fij(2) = fij(2) + swderiv*d_grp(2)
992 >                         fij(3) = fij(3) + swderiv*d_grp(3)
993 >                        
994 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
995 >                            atom1=groupListRow(ia)
996 >                            mf = mfactRow(atom1)
997   #ifdef IS_MPI
998 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
999 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1000 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
998 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
999 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1000 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1001   #else
1002 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1003 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1004 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1002 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1003 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1004 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1005   #endif
1006 <                      enddo
1007 <
1008 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
1009 <                         atom2=groupListCol(jb)
1010 <                         mf = mfactCol(atom2)
1006 >                         enddo
1007 >                        
1008 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
1009 >                            atom2=groupListCol(jb)
1010 >                            mf = mfactCol(atom2)
1011   #ifdef IS_MPI
1012 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1013 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1014 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1012 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1013 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1014 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1015   #else
1016 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1017 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1018 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1016 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1017 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1018 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1019   #endif
1020 <                      enddo
1021 <                   endif
1020 >                         enddo
1021 >                      endif
1022  
1023 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1023 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1024 >                   endif
1025                  endif
1026 <             end if
1026 >             endif
1027            enddo
1028 <
1028 >          
1029         enddo outer
1030  
1031         if (update_nlist) then
# Line 1088 | Line 1085 | contains
1085  
1086      if (do_pot) then
1087         ! scatter/gather pot_row into the members of my column
1088 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1089 <
1088 >       do i = 1,LR_POT_TYPES
1089 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1090 >       end do
1091         ! scatter/gather pot_local into all other procs
1092         ! add resultant to get total pot
1093         do i = 1, nlocal
1094 <          pot_local = pot_local + pot_Temp(i)
1094 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1095 >               + pot_Temp(1:LR_POT_TYPES,i)
1096         enddo
1097  
1098         pot_Temp = 0.0_DP
1099 <
1100 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1099 >       do i = 1,LR_POT_TYPES
1100 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1101 >       end do
1102         do i = 1, nlocal
1103 <          pot_local = pot_local + pot_Temp(i)
1103 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1104 >               + pot_Temp(1:LR_POT_TYPES,i)
1105         enddo
1106  
1107      endif
1108   #endif
1108 –
1109 –    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1109  
1110 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1110 >    if (SIM_requires_postpair_calc) then
1111 >       do i = 1, nlocal            
1112 >          
1113 >          ! we loop only over the local atoms, so we don't need row and column
1114 >          ! lookups for the types
1115 >          
1116 >          me_i = atid(i)
1117 >          
1118 >          ! is the atom electrostatic?  See if it would have an
1119 >          ! electrostatic interaction with itself
1120 >          iHash = InteractionHash(me_i,me_i)
1121  
1122 +          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1123   #ifdef IS_MPI
1124 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1125 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1116 <          do i = 1,nlocal
1117 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1118 <          end do
1119 < #endif
1120 <
1121 <          do i = 1, nLocal
1122 <
1123 <             rfpot = 0.0_DP
1124 < #ifdef IS_MPI
1125 <             me_i = atid_row(i)
1124 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1125 >                  t, do_pot)
1126   #else
1127 <             me_i = atid(i)
1127 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1128 >                  t, do_pot)
1129   #endif
1130 <             iHash = InteractionHash(me_i,me_j)
1130 >          endif
1131 >  
1132 >          
1133 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1134              
1135 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1136 <
1137 <                mu_i = getDipoleMoment(me_i)
1138 <
1139 <                !! The reaction field needs to include a self contribution
1140 <                !! to the field:
1141 <                call accumulate_self_rf(i, mu_i, eFrame)
1142 <                !! Get the reaction field contribution to the
1143 <                !! potential and torques:
1144 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1135 >             ! loop over the excludes to accumulate RF stuff we've
1136 >             ! left out of the normal pair loop
1137 >            
1138 >             do i1 = 1, nSkipsForAtom(i)
1139 >                j = skipsForAtom(i, i1)
1140 >                
1141 >                ! prevent overcounting of the skips
1142 >                if (i.lt.j) then
1143 >                   call get_interatomic_vector(q(:,i), &
1144 >                        q(:,j), d_atm, ratmsq)
1145 >                   rVal = dsqrt(ratmsq)
1146 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1147 >                        in_switching_region)
1148   #ifdef IS_MPI
1149 <                pot_local(RF_POT) = pot_local(RF_POT) + rfpot
1149 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1150 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1151   #else
1152 <                pot(RF_POT) = pot(RF_POT) + rfpot
1153 <
1152 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1153 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1154   #endif
1155 <             endif
1156 <          enddo
1157 <       endif
1155 >                endif
1156 >             enddo
1157 >          endif
1158 >       enddo
1159      endif
1160 <
1152 <
1160 >    
1161   #ifdef IS_MPI
1162 <
1162 >    
1163      if (do_pot) then
1164 <       pot(1:SIZE_POT_ARRAY) = pot(1:SIZE_POT_ARRAY) &
1165 <            + pot_local(1:SIZE_POT_ARRAY)
1158 <       !! we assume the c code will do the allreduce to get the total potential
1159 <       !! we could do it right here if we needed to...
1164 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1165 >            mpi_comm_world,mpi_err)            
1166      endif
1167 <
1167 >    
1168      if (do_stress) then
1169         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1170              mpi_comm_world,mpi_err)
1171         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1172              mpi_comm_world,mpi_err)
1173      endif
1174 <
1174 >    
1175   #else
1176 <
1176 >    
1177      if (do_stress) then
1178         tau = tau_Temp
1179         virial = virial_Temp
1180      endif
1181 <
1181 >    
1182   #endif
1183 <
1183 >    
1184    end subroutine do_force_loop
1185  
1186    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1187 <       eFrame, A, f, t, pot, vpair, fpair)
1187 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1188  
1189      real( kind = dp ) :: vpair, sw
1190 <    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1190 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1191      real( kind = dp ), dimension(3) :: fpair
1192      real( kind = dp ), dimension(nLocal)   :: mfact
1193      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1192 | Line 1198 | contains
1198      logical, intent(inout) :: do_pot
1199      integer, intent(in) :: i, j
1200      real ( kind = dp ), intent(inout) :: rijsq
1201 <    real ( kind = dp )                :: r
1201 >    real ( kind = dp ), intent(inout) :: r_grp
1202      real ( kind = dp ), intent(inout) :: d(3)
1203 +    real ( kind = dp ), intent(inout) :: d_grp(3)
1204 +    real ( kind = dp ) :: r
1205      integer :: me_i, me_j
1206  
1207      integer :: iHash
# Line 1211 | Line 1219 | contains
1219   #endif
1220  
1221      iHash = InteractionHash(me_i, me_j)
1222 <
1222 >    
1223      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1224 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(LJ_POT), f, do_pot)
1224 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1225 >            pot(VDW_POT), f, do_pot)
1226      endif
1227 <
1227 >    
1228      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1229         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1230              pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1222 –
1223 –       if (electrostaticSummationMethod == REACTION_FIELD) then
1224 –
1225 –          ! CHECK ME (RF needs to know about all electrostatic types)
1226 –          call accumulate_rf(i, j, r, eFrame, sw)
1227 –          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1228 –       endif
1229 –
1231      endif
1232 <
1232 >    
1233      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1234         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1235 <            pot(STICKY_POT), A, f, t, do_pot)
1235 >            pot(HB_POT), A, f, t, do_pot)
1236      endif
1237 <
1237 >    
1238      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1239         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1240 <            pot(STICKYPOWER_POT), A, f, t, do_pot)
1240 >            pot(HB_POT), A, f, t, do_pot)
1241      endif
1242 <
1242 >    
1243      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1244         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1245 <            pot(GAYBERNE_POT), A, f, t, do_pot)
1245 >            pot(VDW_POT), A, f, t, do_pot)
1246      endif
1247      
1248      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1249 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1250 < !           pot(GAYBERNE_LJ_POT), A, f, t, do_pot)
1249 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1250 >            pot(VDW_POT), A, f, t, do_pot)
1251      endif
1252 <
1252 >    
1253      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1254 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(EAM_POT), f, &
1255 <            do_pot)
1254 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1255 >            pot(METALLIC_POT), f, do_pot)
1256      endif
1257 <
1257 >    
1258      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1259         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1260 <            pot(SHAPE_POT), A, f, t, do_pot)
1260 >            pot(VDW_POT), A, f, t, do_pot)
1261      endif
1262 <
1262 >    
1263      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1264         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1265 <            pot(SHAPE_LJ_POT), A, f, t, do_pot)
1265 >            pot(VDW_POT), A, f, t, do_pot)
1266      endif
1267 +
1268 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1269 +       call do_SC_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1270 +            pot(METALLIC_POT), f, do_pot)
1271 +    endif
1272 +
1273      
1274 +    
1275    end subroutine do_pair
1276  
1277    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1278         do_pot, do_stress, eFrame, A, f, t, pot)
1279  
1280      real( kind = dp ) :: sw
1281 <    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1281 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1282      real( kind = dp ), dimension(9,nLocal) :: eFrame
1283      real (kind=dp), dimension(9,nLocal) :: A
1284      real (kind=dp), dimension(3,nLocal) :: f
# Line 1299 | Line 1307 | contains
1307      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1308              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1309      endif
1310 +
1311 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1312 +            call calc_SC_prepair_rho(i, j, d, r, rijsq )
1313 +    endif
1314      
1315    end subroutine do_prepair
1316  
1317  
1318    subroutine do_preforce(nlocal,pot)
1319      integer :: nlocal
1320 <    real( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
1320 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1321  
1322      if (FF_uses_EAM .and. SIM_uses_EAM) then
1323 <       call calc_EAM_preforce_Frho(nlocal,pot(EAM_POT))
1323 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1324      endif
1325 +    if (FF_uses_SC .and. SIM_uses_SC) then
1326 +       call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1327 +    endif
1328  
1329  
1330    end subroutine do_preforce
# Line 1394 | Line 1409 | contains
1409      pot_Col = 0.0_dp
1410      pot_Temp = 0.0_dp
1411  
1397 –    rf_Row = 0.0_dp
1398 –    rf_Col = 0.0_dp
1399 –    rf_Temp = 0.0_dp
1400 –
1412   #endif
1413  
1414      if (FF_uses_EAM .and. SIM_uses_EAM) then
1415         call clean_EAM()
1416      endif
1417  
1407 –    rf = 0.0_dp
1418      tau_Temp = 0.0_dp
1419      virial_Temp = 0.0_dp
1420    end subroutine zero_work_arrays
# Line 1498 | Line 1508 | contains
1508  
1509    function FF_RequiresPrepairCalc() result(doesit)
1510      logical :: doesit
1511 <    doesit = FF_uses_EAM
1511 >    doesit = FF_uses_EAM .or. FF_uses_SC &
1512 >         .or. FF_uses_MEAM
1513    end function FF_RequiresPrepairCalc
1514  
1504 –  function FF_RequiresPostpairCalc() result(doesit)
1505 –    logical :: doesit
1506 –    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1507 –  end function FF_RequiresPostpairCalc
1508 –
1515   #ifdef PROFILE
1516    function getforcetime() result(totalforcetime)
1517      real(kind=dp) :: totalforcetime

Diff Legend

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