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 2405 by chrisfen, Tue Nov 1 19:24:57 2005 UTC vs.
Revision 2407 by chrisfen, Wed Nov 2 20:35:34 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.65 2005-11-01 19:24:54 chrisfen Exp $, $Date: 2005-11-01 19:24:54 $, $Name: not supported by cvs2svn $, $Revision: 1.65 $
48 > !! @version $Id: doForces.F90,v 1.66 2005-11-02 20:35:34 chrisfen Exp $, $Date: 2005-11-02 20:35:34 $, $Name: not supported by cvs2svn $, $Revision: 1.66 $
49  
50  
51   module doForces
# Line 745 | Line 745 | contains
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), dimension(3) :: fstrs, f2strs
749      real(kind=dp) :: rfpot, mu_i, virial
750      integer :: me_i, me_j, n_in_i, n_in_j
751      logical :: is_dp_i
# Line 902 | Line 903 | contains
903  
904                     list(nlist) = j
905                  endif
906 <
907 <                if (loop .eq. PAIR_LOOP) then
907 <                   vij = 0.0d0
908 <                   fij(1:3) = 0.0d0
909 <                endif
910 <
911 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
912 <                     in_switching_region)
913 <
914 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
906 >                
907 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
908  
909 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
910 <
911 <                   atom1 = groupListRow(ia)
912 <
913 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
914 <
915 <                      atom2 = groupListCol(jb)
916 <    
917 <                      if (skipThisPair(atom1, atom2))  cycle inner
909 >                   if (loop .eq. PAIR_LOOP) then
910 >                      vij = 0.0d0
911 >                      fij(1:3) = 0.0d0
912 >                   endif
913 >                  
914 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
915 >                        in_switching_region)
916 >                  
917 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
918 >                  
919 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
920                        
921 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
922 <                         d_atm(1:3) = d_grp(1:3)
923 <                         ratmsq = rgrpsq
924 <                      else
921 >                      atom1 = groupListRow(ia)
922 >                      
923 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
924 >                        
925 >                         atom2 = groupListCol(jb)
926 >                        
927 >                         if (skipThisPair(atom1, atom2))  cycle inner
928 >                        
929 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
930 >                            d_atm(1:3) = d_grp(1:3)
931 >                            ratmsq = rgrpsq
932 >                         else
933   #ifdef IS_MPI
934 <                         call get_interatomic_vector(q_Row(:,atom1), &
935 <                              q_Col(:,atom2), d_atm, ratmsq)
936 < #else
937 <                         call get_interatomic_vector(q(:,atom1), &
938 <                              q(:,atom2), d_atm, ratmsq)
934 >                            call get_interatomic_vector(q_Row(:,atom1), &
935 >                                 q_Col(:,atom2), d_atm, ratmsq)
936 > #else
937 >                            call get_interatomic_vector(q(:,atom1), &
938 >                                 q(:,atom2), d_atm, ratmsq)
939   #endif
940 <                      endif
941 <
942 <                      if (loop .eq. PREPAIR_LOOP) then
940 >                         endif
941 >                        
942 >                         if (loop .eq. PREPAIR_LOOP) then
943   #ifdef IS_MPI                      
944 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
945 <                              rgrpsq, d_grp, do_pot, do_stress, &
946 <                              eFrame, A, f, t, pot_local)
944 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
945 >                                 rgrpsq, d_grp, do_pot, do_stress, &
946 >                                 eFrame, A, f, t, pot_local)
947   #else
948 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
949 <                              rgrpsq, d_grp, do_pot, do_stress, &
950 <                              eFrame, A, f, t, pot)
948 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
949 >                                 rgrpsq, d_grp, do_pot, do_stress, &
950 >                                 eFrame, A, f, t, pot)
951   #endif                                              
952 <                      else
952 >                         else
953   #ifdef IS_MPI                      
954 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
955 <                              do_pot, eFrame, A, f, t, pot_local, vpair, &
956 <                              fpair, d_grp, rgrp)
954 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
955 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
956 >                                 fpair, d_grp, rgrp, fstrs)
957   #else
958 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
959 <                              do_pot, eFrame, A, f, t, pot, vpair, fpair, &
960 <                              d_grp, rgrp)
958 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
959 >                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
960 >                                 d_grp, rgrp, fstrs)
961   #endif
962 +                            f2strs(1:3) = f2strs(1:3) + fstrs(1:3)
963 +                            vij = vij + vpair
964 +                            fij(1:3) = fij(1:3) + fpair(1:3)
965 +                         endif
966 +                      enddo inner
967 +                   enddo
968  
969 <                         vij = vij + vpair
970 <                         fij(1:3) = fij(1:3) + fpair(1:3)
971 <                      endif
972 <                   enddo inner
973 <                enddo
974 <
975 <                if (loop .eq. PAIR_LOOP) then
976 <                   if (in_switching_region) then
977 <                      swderiv = vij*dswdr/rgrp
978 <                      fij(1) = fij(1) + swderiv*d_grp(1)
970 <                      fij(2) = fij(2) + swderiv*d_grp(2)
971 <                      fij(3) = fij(3) + swderiv*d_grp(3)
972 <
973 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
974 <                         atom1=groupListRow(ia)
975 <                         mf = mfactRow(atom1)
969 >                   if (loop .eq. PAIR_LOOP) then
970 >                      if (in_switching_region) then
971 >                         swderiv = vij*dswdr/rgrp
972 >                         fij(1) = fij(1) + swderiv*d_grp(1)
973 >                         fij(2) = fij(2) + swderiv*d_grp(2)
974 >                         fij(3) = fij(3) + swderiv*d_grp(3)
975 >                        
976 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
977 >                            atom1=groupListRow(ia)
978 >                            mf = mfactRow(atom1)
979   #ifdef IS_MPI
980 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
981 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
982 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
980 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
981 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
982 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
983   #else
984 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
985 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
986 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
984 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
985 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
986 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
987   #endif
988 <                      enddo
989 <
990 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
991 <                         atom2=groupListCol(jb)
992 <                         mf = mfactCol(atom2)
993 < #ifdef IS_MPI
994 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
995 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
996 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
988 >                         enddo
989 >                        
990 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
991 >                            atom2=groupListCol(jb)
992 >                            mf = mfactCol(atom2)
993 > #ifdef IS_MPI
994 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
995 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
996 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
997   #else
998 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
999 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1000 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
998 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
999 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1000 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1001   #endif
1002 <                      enddo
1003 <                   endif
1002 >                         enddo
1003 >                      endif
1004  
1005 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1006 <                endif
1007 <             end if
1005 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1006 >                   endif
1007 >                endif
1008 >             endif
1009            enddo
1010 <
1010 >          
1011         enddo outer
1012  
1013         if (update_nlist) then
# Line 1108 | Line 1112 | contains
1112            endif
1113    
1114            
1115 < !!$          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1115 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1116              
1117               ! loop over the excludes to accumulate RF stuff we've
1118               ! left out of the normal pair loop
# Line 1132 | Line 1136 | contains
1136   #endif
1137                  endif
1138               enddo
1139 < !!$          endif
1139 >          endif
1140         enddo
1141      endif
1142      
# Line 1161 | Line 1165 | contains
1165      
1166    end subroutine do_force_loop
1167  
1168 + !!$  subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1169 + !!$       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1170    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1171 <       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1171 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, fstrs)
1172  
1173      real( kind = dp ) :: vpair, sw
1174      real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1175      real( kind = dp ), dimension(3) :: fpair
1176 +    real( kind = dp ), dimension(3) :: fstrs
1177      real( kind = dp ), dimension(nLocal)   :: mfact
1178      real( kind = dp ), dimension(9,nLocal) :: eFrame
1179      real( kind = dp ), dimension(9,nLocal) :: A
# Line 1204 | Line 1211 | contains
1211      endif
1212      
1213      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1214 + !!$       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1215 + !!$            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1216         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1217 <            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1217 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot, fstrs)
1218      endif
1219      
1220      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
# Line 1379 | Line 1388 | contains
1388         call clean_EAM()
1389      endif
1390  
1382 –    rf = 0.0_dp
1391      tau_Temp = 0.0_dp
1392      virial_Temp = 0.0_dp
1393    end subroutine zero_work_arrays

Diff Legend

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