ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2381 by chrisfen, Tue Oct 18 15:01:42 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.60 2005-10-18 15:01:42 chrisfen Exp $, $Date: 2005-10-18 15:01:42 $, $Name: not supported by cvs2svn $, $Revision: 1.60 $
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 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), 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 754 | Line 755 | contains
755      integer :: propPack_i, propPack_j
756      integer :: loopStart, loopEnd, loop
757      integer :: iHash
758 +    integer :: i1
759    
760  
761      !! initialize local variables  
# Line 901 | Line 903 | contains
903  
904                     list(nlist) = j
905                  endif
906 <
907 <                if (loop .eq. PAIR_LOOP) then
906 <                   vij = 0.0d0
907 <                   fij(1:3) = 0.0d0
908 <                endif
909 <
910 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
911 <                     in_switching_region)
906 >                
907 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
908  
909 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
910 <
911 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
912 <
913 <                   atom1 = groupListRow(ia)
914 <
915 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
916 <
917 <                      atom2 = groupListCol(jb)
918 <
919 <                      if (skipThisPair(atom1, atom2)) cycle inner
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
925 < #ifdef IS_MPI
926 <                         call get_interatomic_vector(q_Row(:,atom1), &
927 <                              q_Col(:,atom2), d_atm, ratmsq)
928 < #else
929 <                         call get_interatomic_vector(q(:,atom1), &
930 <                              q(:,atom2), d_atm, ratmsq)
931 < #endif
932 <                      endif
933 <
934 <                      if (loop .eq. PREPAIR_LOOP) then
935 < #ifdef IS_MPI                      
936 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
937 <                              rgrpsq, d_grp, do_pot, do_stress, &
938 <                              eFrame, A, f, t, pot_local)
939 < #else
940 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
941 <                              rgrpsq, d_grp, do_pot, do_stress, &
942 <                              eFrame, A, f, t, pot)
947 < #endif                                              
948 <                      else
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 >                      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)
939 > #endif
940 >                         endif
941 >                        
942 >                         if (loop .eq. PREPAIR_LOOP) then
943   #ifdef IS_MPI                      
944 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
945 <                              do_pot, &
946 <                              eFrame, A, f, t, pot_local, vpair, fpair)
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_pair(atom1, atom2, ratmsq, d_atm, sw, &
949 <                              do_pot,  &
950 <                              eFrame, A, f, t, pot, vpair, fpair)
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
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, 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, 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)
969 <                      fij(2) = fij(2) + swderiv*d_grp(2)
970 <                      fij(3) = fij(3) + swderiv*d_grp(3)
971 <
972 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
973 <                         atom1=groupListRow(ia)
974 <                         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)
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
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)
1005 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1006 >                   endif
1007                  endif
1008 <             end if
1008 >             endif
1009            enddo
1010 <
1010 >          
1011         enddo outer
1012  
1013         if (update_nlist) then
# Line 1084 | Line 1089 | contains
1089      endif
1090   #endif
1091  
1092 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1092 >    if (SIM_requires_postpair_calc) then
1093 >       do i = 1, nlocal            
1094 >          
1095 >          ! we loop only over the local atoms, so we don't need row and column
1096 >          ! lookups for the types
1097 >          
1098 >          me_i = atid(i)
1099 >          
1100 >          ! is the atom electrostatic?  See if it would have an
1101 >          ! electrostatic interaction with itself
1102 >          iHash = InteractionHash(me_i,me_i)
1103  
1104 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1090 <
1104 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1105   #ifdef IS_MPI
1106 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1107 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1094 <          do i = 1,nlocal
1095 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1096 <          end do
1097 < #endif
1098 <
1099 <          do i = 1, nLocal
1100 <
1101 <             rfpot = 0.0_DP
1102 < #ifdef IS_MPI
1103 <             me_i = atid_row(i)
1106 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1107 >                  t, do_pot)
1108   #else
1109 <             me_i = atid(i)
1109 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1110 >                  t, do_pot)
1111   #endif
1112 <             iHash = InteractionHash(me_i,me_j)
1112 >          endif
1113 >  
1114 >          
1115 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1116              
1117 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1118 <
1119 <                mu_i = getDipoleMoment(me_i)
1120 <
1121 <                !! The reaction field needs to include a self contribution
1122 <                !! to the field:
1123 <                call accumulate_self_rf(i, mu_i, eFrame)
1124 <                !! Get the reaction field contribution to the
1125 <                !! potential and torques:
1126 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1117 >             ! loop over the excludes to accumulate RF stuff we've
1118 >             ! left out of the normal pair loop
1119 >            
1120 >             do i1 = 1, nSkipsForAtom(i)
1121 >                j = skipsForAtom(i, i1)
1122 >                
1123 >                ! prevent overcounting of the skips
1124 >                if (i.lt.j) then
1125 >                   call get_interatomic_vector(q(:,i), &
1126 >                        q(:,j), d_atm, ratmsq)
1127 >                   rVal = dsqrt(ratmsq)
1128 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1129 >                        in_switching_region)
1130   #ifdef IS_MPI
1131 <                pot_local(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
1131 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1132 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1133   #else
1134 <                pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1135 <
1134 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1135 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1136   #endif
1137 <             endif
1138 <          enddo
1139 <       endif
1137 >                endif
1138 >             enddo
1139 >          endif
1140 >       enddo
1141      endif
1142 <
1143 <
1144 < #ifdef IS_MPI
1132 <
1142 >    
1143 > #ifdef IS_MPI
1144 >    
1145      if (do_pot) then
1146 <       pot(1:LR_POT_TYPES) = pot(1:LR_POT_TYPES) &
1147 <            + pot_local(1:LR_POT_TYPES)
1136 <       !! we assume the c code will do the allreduce to get the total potential
1137 <       !! we could do it right here if we needed to...
1146 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1147 >            mpi_comm_world,mpi_err)            
1148      endif
1149 <
1149 >    
1150      if (do_stress) then
1151         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1152              mpi_comm_world,mpi_err)
1153         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1154              mpi_comm_world,mpi_err)
1155      endif
1156 <
1156 >    
1157   #else
1158 <
1158 >    
1159      if (do_stress) then
1160         tau = tau_Temp
1161         virial = virial_Temp
1162      endif
1163 <
1163 >    
1164   #endif
1165 <
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)
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 1170 | Line 1183 | contains
1183      logical, intent(inout) :: do_pot
1184      integer, intent(in) :: i, j
1185      real ( kind = dp ), intent(inout) :: rijsq
1186 <    real ( kind = dp )                :: r
1186 >    real ( kind = dp ), intent(inout) :: r_grp
1187      real ( kind = dp ), intent(inout) :: d(3)
1188 +    real ( kind = dp ), intent(inout) :: d_grp(3)
1189 +    real ( kind = dp ) :: r
1190      integer :: me_i, me_j
1191  
1192      integer :: iHash
# Line 1189 | Line 1204 | contains
1204   #endif
1205  
1206      iHash = InteractionHash(me_i, me_j)
1207 <
1207 >    
1208      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1209         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1210              pot(VDW_POT), f, do_pot)
1211      endif
1212 <
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)
1201 <
1202 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1203 <
1204 <          ! CHECK ME (RF needs to know about all electrostatic types)
1205 <          call accumulate_rf(i, j, r, eFrame, sw)
1206 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1207 <       endif
1208 <
1217 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot, fstrs)
1218      endif
1219 <
1219 >    
1220      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1221         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1222              pot(HB_POT), A, f, t, do_pot)
1223      endif
1224 <
1224 >    
1225      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1226         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1227              pot(HB_POT), A, f, t, do_pot)
1228      endif
1229 <
1229 >    
1230      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1231         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1232              pot(VDW_POT), A, f, t, do_pot)
# Line 1227 | Line 1236 | contains
1236         call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1237              pot(VDW_POT), A, f, t, do_pot)
1238      endif
1239 <
1239 >    
1240      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1241         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1242              pot(METALLIC_POT), f, do_pot)
1243      endif
1244 <
1244 >    
1245      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1246         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1247              pot(VDW_POT), A, f, t, do_pot)
1248      endif
1249 <
1249 >    
1250      if ( iand(iHash, SHAPE_LJ).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 <    
1254 >    
1255    end subroutine do_pair
1256  
1257    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1372 | Line 1381 | contains
1381      pot_Row = 0.0_dp
1382      pot_Col = 0.0_dp
1383      pot_Temp = 0.0_dp
1375 –
1376 –    rf_Row = 0.0_dp
1377 –    rf_Col = 0.0_dp
1378 –    rf_Temp = 0.0_dp
1384  
1385   #endif
1386  
# Line 1383 | Line 1388 | contains
1388         call clean_EAM()
1389      endif
1390  
1386 –    rf = 0.0_dp
1391      tau_Temp = 0.0_dp
1392      virial_Temp = 0.0_dp
1393    end subroutine zero_work_arrays
# Line 1480 | Line 1484 | contains
1484      doesit = FF_uses_EAM
1485    end function FF_RequiresPrepairCalc
1486  
1483 –  function FF_RequiresPostpairCalc() result(doesit)
1484 –    logical :: doesit
1485 –    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1486 –  end function FF_RequiresPostpairCalc
1487 –
1487   #ifdef PROFILE
1488    function getforcetime() result(totalforcetime)
1489      real(kind=dp) :: totalforcetime

Diff Legend

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