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 2367 by kdaily, Fri Oct 14 15:44:37 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.58 2005-10-14 15:44:37 kdaily Exp $, $Date: 2005-10-14 15:44:37 $, $Name: not supported by cvs2svn $, $Revision: 1.58 $
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 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
# Line 645 | Line 644 | contains
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
646      integer, pointer :: MatchList(:) => null()
648    real(kind=dp) :: rcut, rrf, rt, dielect
647  
648      !! assume things are copacetic, unless they aren't
649      thisStat = 0
# Line 681 | Line 679 | contains
679  
680      haveSaneForceField = .true.
681  
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
699
682      if (FF_uses_EAM) then
683         call init_EAM_FF(my_status)
684         if (my_status /= 0) then
# Line 705 | Line 687 | contains
687            haveSaneForceField = .false.
688            return
689         end if
708    endif
709
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
690      endif
691  
692      if (.not. haveNeighborList) then
# Line 770 | 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 780 | Line 754 | contains
754      integer :: propPack_i, propPack_j
755      integer :: loopStart, loopEnd, loop
756      integer :: iHash
757 +    integer :: i1
758    
759  
760      !! initialize local variables  
# Line 927 | Line 902 | contains
902  
903                     list(nlist) = j
904                  endif
905 +                
906 +                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
907  
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
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)
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, &
955 <                              eFrame, A, f, t, pot_local, vpair, fpair)
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,  &
959 <                              eFrame, A, f, t, pot, vpair, fpair)
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)
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)
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
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 1110 | Line 1087 | contains
1087      endif
1088   #endif
1089  
1090 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1090 >    if (SIM_requires_postpair_calc) then
1091 >       do i = 1, nlocal            
1092 >          
1093 >          ! we loop only over the local atoms, so we don't need row and column
1094 >          ! lookups for the types
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  
1102 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1116 <
1102 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1103   #ifdef IS_MPI
1104 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1105 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1120 <          do i = 1,nlocal
1121 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1122 <          end do
1123 < #endif
1124 <
1125 <          do i = 1, nLocal
1126 <
1127 <             rfpot = 0.0_DP
1128 < #ifdef IS_MPI
1129 <             me_i = atid_row(i)
1104 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1105 >                  t, do_pot)
1106   #else
1107 <             me_i = atid(i)
1107 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1108 >                  t, do_pot)
1109   #endif
1110 <             iHash = InteractionHash(me_i,me_j)
1110 >          endif
1111 >  
1112 >          
1113 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1114              
1115 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1116 <
1117 <                mu_i = getDipoleMoment(me_i)
1118 <
1119 <                !! The reaction field needs to include a self contribution
1120 <                !! to the field:
1121 <                call accumulate_self_rf(i, mu_i, eFrame)
1122 <                !! Get the reaction field contribution to the
1123 <                !! potential and torques:
1124 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
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 <                pot_local(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
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 <                pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1133 <
1134 < #endif
1135 <             endif
1136 <          enddo
1137 <       endif
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
1140 <
1156 <
1140 >    
1141   #ifdef IS_MPI
1142 <
1142 >    
1143      if (do_pot) then
1144 <       pot(1:LR_POT_TYPES) = pot(1:LR_POT_TYPES) &
1145 <            + pot_local(1:LR_POT_TYPES)
1162 <       !! we assume the c code will do the allreduce to get the total potential
1163 <       !! we could do it right here if we needed to...
1144 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1145 >            mpi_comm_world,mpi_err)            
1146      endif
1147 <
1147 >    
1148      if (do_stress) then
1149         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1150              mpi_comm_world,mpi_err)
1151         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1152              mpi_comm_world,mpi_err)
1153      endif
1154 <
1154 >    
1155   #else
1156 <
1156 >    
1157      if (do_stress) then
1158         tau = tau_Temp
1159         virial = virial_Temp
1160      endif
1161 <
1161 >    
1162   #endif
1163 <
1163 >    
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)
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 1196 | Line 1178 | contains
1178      logical, intent(inout) :: do_pot
1179      integer, intent(in) :: i, j
1180      real ( kind = dp ), intent(inout) :: rijsq
1181 <    real ( kind = dp )                :: r
1181 >    real ( kind = dp ), intent(inout) :: r_grp
1182      real ( kind = dp ), intent(inout) :: d(3)
1183 +    real ( kind = dp ), intent(inout) :: d_grp(3)
1184 +    real ( kind = dp ) :: r
1185      integer :: me_i, me_j
1186  
1187      integer :: iHash
# Line 1215 | 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)
1227
1228       if (electrostaticSummationMethod == REACTION_FIELD) then
1229
1230          ! CHECK ME (RF needs to know about all electrostatic types)
1231          call accumulate_rf(i, j, r, eFrame, sw)
1232          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1233       endif
1234
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 1253 | 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 1398 | Line 1374 | contains
1374      pot_Row = 0.0_dp
1375      pot_Col = 0.0_dp
1376      pot_Temp = 0.0_dp
1401
1402    rf_Row = 0.0_dp
1403    rf_Col = 0.0_dp
1404    rf_Temp = 0.0_dp
1377  
1378   #endif
1379  
# Line 1409 | Line 1381 | contains
1381         call clean_EAM()
1382      endif
1383  
1412    rf = 0.0_dp
1384      tau_Temp = 0.0_dp
1385      virial_Temp = 0.0_dp
1386    end subroutine zero_work_arrays
# Line 1506 | Line 1477 | contains
1477      doesit = FF_uses_EAM
1478    end function FF_RequiresPrepairCalc
1479  
1509  function FF_RequiresPostpairCalc() result(doesit)
1510    logical :: doesit
1511    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1512  end function FF_RequiresPostpairCalc
1513
1480   #ifdef PROFILE
1481    function getforcetime() result(totalforcetime)
1482      real(kind=dp) :: totalforcetime

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines