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 2398 by chrisfen, Wed Oct 26 23:31:18 2005 UTC vs.
Revision 2411 by chrisfen, Wed Nov 2 21:01:21 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.63 2005-10-26 23:31:18 chrisfen Exp $, $Date: 2005-10-26 23:31:18 $, $Name: not supported by cvs2svn $, $Revision: 1.63 $
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 $
49  
50  
51   module doForces
# Line 644 | Line 644 | contains
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
646      integer, pointer :: MatchList(:) => null()
647 –    real(kind=dp) :: rcut, rrf, rt, dielect
647  
648      !! assume things are copacetic, unless they aren't
649      thisStat = 0
# Line 756 | Line 755 | contains
755      integer :: loopStart, loopEnd, loop
756      integer :: iHash
757      integer :: i1
759 –    logical :: indirect_only
758    
759  
760      !! initialize local variables  
# Line 904 | Line 902 | contains
902  
903                     list(nlist) = j
904                  endif
905 <
906 <                if (loop .eq. PAIR_LOOP) then
909 <                   vij = 0.0d0
910 <                   fij(1:3) = 0.0d0
911 <                endif
912 <
913 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
914 <                     in_switching_region)
915 <
916 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
905 >                
906 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
907  
908 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
909 <
910 <                   atom1 = groupListRow(ia)
911 <
912 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
913 <
914 <                      atom2 = groupListCol(jb)
915 <
916 <                      indirect_only = .false.
917 <    
918 <                      if (skipThisPair(atom1, atom2)) then
919 <                         if (electrostaticSummationMethod.ne.REACTION_FIELD) then
920 <                            cycle inner
908 >                   if (loop .eq. PAIR_LOOP) then
909 >                      vij = 0.0d0
910 >                      fij(1:3) = 0.0d0
911 >                   endif
912 >                  
913 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
914 >                        in_switching_region)
915 >                  
916 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
917 >                  
918 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
919 >                      
920 >                      atom1 = groupListRow(ia)
921 >                      
922 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
923 >                        
924 >                         atom2 = groupListCol(jb)
925 >                        
926 >                         if (skipThisPair(atom1, atom2))  cycle inner
927 >                        
928 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
929 >                            d_atm(1:3) = d_grp(1:3)
930 >                            ratmsq = rgrpsq
931                           else
932 –                            indirect_only = .true.
933 –                         endif
934 –                      endif
935 –    
936 –
937 –                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
938 –                         d_atm(1:3) = d_grp(1:3)
939 –                         ratmsq = rgrpsq
940 –                      else
932   #ifdef IS_MPI
933 <                         call get_interatomic_vector(q_Row(:,atom1), &
934 <                              q_Col(:,atom2), d_atm, ratmsq)
933 >                            call get_interatomic_vector(q_Row(:,atom1), &
934 >                                 q_Col(:,atom2), d_atm, ratmsq)
935   #else
936 <                         call get_interatomic_vector(q(:,atom1), &
937 <                              q(:,atom2), d_atm, ratmsq)
936 >                            call get_interatomic_vector(q(:,atom1), &
937 >                                 q(:,atom2), d_atm, ratmsq)
938   #endif
939 <                      endif
940 <
941 <                      if (loop .eq. PREPAIR_LOOP) then
939 >                         endif
940 >                        
941 >                         if (loop .eq. PREPAIR_LOOP) then
942   #ifdef IS_MPI                      
943 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
944 <                              rgrpsq, d_grp, do_pot, do_stress, &
945 <                              eFrame, A, f, t, pot_local)
943 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
944 >                                 rgrpsq, d_grp, do_pot, do_stress, &
945 >                                 eFrame, A, f, t, pot_local)
946   #else
947 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
948 <                              rgrpsq, d_grp, do_pot, do_stress, &
949 <                              eFrame, A, f, t, pot)
947 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
948 >                                 rgrpsq, d_grp, do_pot, do_stress, &
949 >                                 eFrame, A, f, t, pot)
950   #endif                                              
951 <                      else
951 >                         else
952   #ifdef IS_MPI                      
953 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
954 <                              do_pot, eFrame, A, f, t, pot_local, vpair, &
955 <                              fpair, d_grp, rgrp, indirect_only)
953 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
954 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
955 >                                 fpair, d_grp, rgrp)
956   #else
957 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
958 <                              do_pot, eFrame, A, f, t, pot, vpair, fpair, &
959 <                              d_grp, rgrp, indirect_only)
957 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
958 >                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
959 >                                 d_grp, rgrp)
960   #endif
961 +                            vij = vij + vpair
962 +                            fij(1:3) = fij(1:3) + fpair(1:3)
963 +                         endif
964 +                      enddo inner
965 +                   enddo
966  
967 <                         vij = vij + vpair
968 <                         fij(1:3) = fij(1:3) + fpair(1:3)
969 <                      endif
970 <                   enddo inner
971 <                enddo
972 <
973 <                if (loop .eq. PAIR_LOOP) then
974 <                   if (in_switching_region) then
975 <                      swderiv = vij*dswdr/rgrp
976 <                      fij(1) = fij(1) + swderiv*d_grp(1)
981 <                      fij(2) = fij(2) + swderiv*d_grp(2)
982 <                      fij(3) = fij(3) + swderiv*d_grp(3)
983 <
984 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
985 <                         atom1=groupListRow(ia)
986 <                         mf = mfactRow(atom1)
967 >                   if (loop .eq. PAIR_LOOP) then
968 >                      if (in_switching_region) then
969 >                         swderiv = vij*dswdr/rgrp
970 >                         fij(1) = fij(1) + swderiv*d_grp(1)
971 >                         fij(2) = fij(2) + swderiv*d_grp(2)
972 >                         fij(3) = fij(3) + swderiv*d_grp(3)
973 >                        
974 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
975 >                            atom1=groupListRow(ia)
976 >                            mf = mfactRow(atom1)
977   #ifdef IS_MPI
978 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
979 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
980 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
978 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
979 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
980 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
981   #else
982 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
983 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
984 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
982 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
983 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
984 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
985   #endif
986 <                      enddo
987 <
988 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
989 <                         atom2=groupListCol(jb)
990 <                         mf = mfactCol(atom2)
986 >                         enddo
987 >                        
988 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
989 >                            atom2=groupListCol(jb)
990 >                            mf = mfactCol(atom2)
991   #ifdef IS_MPI
992 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
993 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
994 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
992 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
993 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
994 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
995   #else
996 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
997 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
998 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
996 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
997 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
998 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
999   #endif
1000 <                      enddo
1001 <                   endif
1002 <
1003 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1000 >                         enddo
1001 >                      endif
1002 >
1003 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1004 >                   endif
1005                  endif
1006 <             end if
1006 >             endif
1007            enddo
1008 <
1008 >          
1009         enddo outer
1010  
1011         if (update_nlist) then
# Line 1110 | Line 1101 | contains
1101  
1102            if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1103   #ifdef IS_MPI
1104 <             call rf_self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1104 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1105                    t, do_pot)
1106   #else
1107 <             call rf_self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1107 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1108                    t, do_pot)
1109   #endif
1110            endif
1111    
1112 <          ! loop over the excludes to accumulate any additional RF components
1113 <
1123 <          do i1 = 1, nSkipsForAtom(i)
1124 <             j = skipsForAtom(i, i1)
1112 >          
1113 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1114              
1115 <             ! prevent overcounting of the skips
1116 <             if (i.lt.j) then
1117 <                call get_interatomic_vector(q(:,i), &
1118 <                     q(:,j), d_atm, ratmsq)
1119 <                rVal = dsqrt(ratmsq)
1120 <                call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1121 <                     in_switching_region)
1115 >             ! loop over the excludes to accumulate RF stuff we've
1116 >             ! left out of the normal pair loop
1117 >            
1118 >             do i1 = 1, nSkipsForAtom(i)
1119 >                j = skipsForAtom(i, i1)
1120 >                
1121 >                ! prevent overcounting of the skips
1122 >                if (i.lt.j) then
1123 >                   call get_interatomic_vector(q(:,i), &
1124 >                        q(:,j), d_atm, ratmsq)
1125 >                   rVal = dsqrt(ratmsq)
1126 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1127 >                        in_switching_region)
1128   #ifdef IS_MPI
1129 <                call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, vpair, &
1130 <                     pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1129 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1130 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1131   #else
1132 <                call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, vpair, &
1133 <                     pot(ELECTROSTATIC_POT), f, t, do_pot)
1132 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1133 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1134   #endif
1135 <             endif
1136 <          enddo
1137 <       enddo      
1135 >                endif
1136 >             enddo
1137 >          endif
1138 >       enddo
1139      endif
1140      
1141   #ifdef IS_MPI
# Line 1168 | Line 1164 | contains
1164    end subroutine do_force_loop
1165  
1166    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1167 <       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, indirect_only)
1167 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1168  
1169      real( kind = dp ) :: vpair, sw
1170      real( kind = dp ), dimension(LR_POT_TYPES) :: pot
# Line 1180 | Line 1176 | contains
1176      real( kind = dp ), dimension(3,nLocal) :: t
1177  
1178      logical, intent(inout) :: do_pot
1183 –    logical, intent(inout) :: indirect_only
1179      integer, intent(in) :: i, j
1180      real ( kind = dp ), intent(inout) :: rijsq
1181      real ( kind = dp ), intent(inout) :: r_grp
# Line 1204 | Line 1199 | contains
1199   #endif
1200  
1201      iHash = InteractionHash(me_i, me_j)
1202 <
1203 <    if (indirect_only) then
1204 <       if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1205 <          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1211 <               pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot, indirect_only)
1212 <       endif
1213 <    else
1214 <          
1215 <       if ( iand(iHash, LJ_PAIR).ne.0 ) then
1216 <          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1217 <               pot(VDW_POT), f, do_pot)
1218 <       endif
1219 <
1220 <       if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1221 <          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1222 <               pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot, indirect_only)
1223 <       endif
1224 <      
1225 <       if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1226 <          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1227 <               pot(HB_POT), A, f, t, do_pot)
1228 <       endif
1229 <      
1230 <       if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1231 <          call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1232 <               pot(HB_POT), A, f, t, do_pot)
1233 <       endif
1234 <      
1235 <       if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1236 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1237 <               pot(VDW_POT), A, f, t, do_pot)
1238 <       endif
1239 <      
1240 <       if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1241 <          call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1242 <               pot(VDW_POT), A, f, t, do_pot)
1243 <       endif
1244 <      
1245 <       if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1246 <          call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1247 <               pot(METALLIC_POT), f, do_pot)
1248 <       endif
1249 <      
1250 <       if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1251 <          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1252 <               pot(VDW_POT), A, f, t, do_pot)
1253 <       endif
1254 <      
1255 <       if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1256 <          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1257 <               pot(VDW_POT), A, f, t, do_pot)
1258 <       endif
1202 >    
1203 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1204 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1205 >            pot(VDW_POT), f, do_pot)
1206      endif
1207      
1208 +    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1209 +       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1210 +            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1211 +    endif
1212 +    
1213 +    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1214 +       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1215 +            pot(HB_POT), A, f, t, do_pot)
1216 +    endif
1217 +    
1218 +    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1219 +       call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1220 +            pot(HB_POT), A, f, t, do_pot)
1221 +    endif
1222 +    
1223 +    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1224 +       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1225 +            pot(VDW_POT), A, f, t, do_pot)
1226 +    endif
1227 +    
1228 +    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1229 +       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1230 +            pot(VDW_POT), A, f, t, do_pot)
1231 +    endif
1232 +    
1233 +    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1234 +       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1235 +            pot(METALLIC_POT), f, do_pot)
1236 +    endif
1237 +    
1238 +    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1239 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1240 +            pot(VDW_POT), A, f, t, do_pot)
1241 +    endif
1242 +    
1243 +    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1244 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1245 +            pot(VDW_POT), A, f, t, do_pot)
1246 +    endif
1247 +    
1248    end subroutine do_pair
1249  
1250    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1387 | Line 1374 | contains
1374      pot_Row = 0.0_dp
1375      pot_Col = 0.0_dp
1376      pot_Temp = 0.0_dp
1390 –
1391 –    rf_Row = 0.0_dp
1392 –    rf_Col = 0.0_dp
1393 –    rf_Temp = 0.0_dp
1377  
1378   #endif
1379  
# Line 1398 | Line 1381 | contains
1381         call clean_EAM()
1382      endif
1383  
1401 –    rf = 0.0_dp
1384      tau_Temp = 0.0_dp
1385      virial_Temp = 0.0_dp
1386    end subroutine zero_work_arrays

Diff Legend

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