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 2394 by chrisfen, Sun Oct 23 21:08:08 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.62 2005-10-23 21:08:02 chrisfen Exp $, $Date: 2005-10-23 21:08:02 $, $Name: not supported by cvs2svn $, $Revision: 1.62 $
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 744 | Line 743 | contains
743      integer :: nlist
744      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
745      real( kind = DP ) :: sw, dswdr, swderiv, mf
746 +    real( kind = DP ) :: rVal
747      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
748      real(kind=dp) :: rfpot, mu_i, virial
749      integer :: me_i, me_j, n_in_i, n_in_j
# Line 754 | Line 754 | contains
754      integer :: propPack_i, propPack_j
755      integer :: loopStart, loopEnd, loop
756      integer :: iHash
757 <    integer :: ig
757 >    integer :: i1
758    
759  
760      !! initialize local variables  
# Line 901 | Line 901 | contains
901                     endif
902  
903                     list(nlist) = j
904 –                endif
905 –
906 –                if (loop .eq. PAIR_LOOP) then
907 –                   vij = 0.0d0
908 –                   fij(1:3) = 0.0d0
904                  endif
905 <
906 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
912 <                     in_switching_region)
913 <
914 <                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 <                      if (skipThisPair(atom1, atom2)) cycle inner
917 <
918 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
919 <                         d_atm(1:3) = d_grp(1:3)
920 <                         ratmsq = rgrpsq
921 <                      else
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   #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)
938 < #endif
939 <                      endif
940 <
941 <                      if (loop .eq. PREPAIR_LOOP) then
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
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)
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)
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)
977 <                      fij(2) = fij(2) + swderiv*d_grp(2)
978 <                      fij(3) = fij(3) + swderiv*d_grp(3)
979 <
980 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
974 <                         atom1=groupListRow(ia)
975 <                         mf = mfactRow(atom1)
976 < #ifdef IS_MPI
977 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
978 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
979 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
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
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
1000 >                         enddo
1001 >                      endif
1002  
1003 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
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 1090 | Line 1092 | contains
1092            
1093            ! we loop only over the local atoms, so we don't need row and column
1094            ! lookups for the types
1095 <
1095 >          
1096            me_i = atid(i)
1097            
1098            ! is the atom electrostatic?  See if it would have an
1099            ! electrostatic interaction with itself
1100            iHash = InteractionHash(me_i,me_i)
1101 <          
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 +          
1113 +          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1114 +            
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, &
1130 +                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1131 + #else
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            endif
1138         enddo
1139      endif
# Line 1170 | Line 1199 | contains
1199   #endif
1200  
1201      iHash = InteractionHash(me_i, me_j)
1202 <
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 <
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 <
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 <
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 <
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)
# Line 1200 | Line 1229 | contains
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 <
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 <
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 <
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 <    
1247 >    
1248    end subroutine do_pair
1249  
1250    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1345 | Line 1374 | contains
1374      pot_Row = 0.0_dp
1375      pot_Col = 0.0_dp
1376      pot_Temp = 0.0_dp
1348 –
1349 –    rf_Row = 0.0_dp
1350 –    rf_Col = 0.0_dp
1351 –    rf_Temp = 0.0_dp
1377  
1378   #endif
1379  
# Line 1356 | Line 1381 | contains
1381         call clean_EAM()
1382      endif
1383  
1359 –    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