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 2351 by chuckv, Mon Oct 10 21:34:54 2005 UTC vs.
Revision 2398 by chrisfen, Wed Oct 26 23:31:18 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.52 2005-10-10 21:34:54 chuckv Exp $, $Date: 2005-10-10 21:34:54 $, $Name: not supported by cvs2svn $, $Revision: 1.52 $
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 $
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 185 | Line 184 | contains
184  
185      if (.not. allocated(InteractionHash)) then
186         allocate(InteractionHash(nAtypes,nAtypes))
187 +    else
188 +       deallocate(InteractionHash)
189 +       allocate(InteractionHash(nAtypes,nAtypes))
190      endif
191  
192      if (.not. allocated(atypeMaxCutoff)) then
193 +       allocate(atypeMaxCutoff(nAtypes))
194 +    else
195 +       deallocate(atypeMaxCutoff)
196         allocate(atypeMaxCutoff(nAtypes))
197      endif
198          
# Line 674 | Line 679 | contains
679  
680  
681      haveSaneForceField = .true.
677
678    !! check to make sure the reaction field setting makes sense
679
680    if (FF_uses_Dipoles) then
681       if (electrostaticSummationMethod == REACTION_FIELD) then
682          dielect = getDielect()
683          call initialize_rf(dielect)
684       endif
685    else
686       if (electrostaticSummationMethod == REACTION_FIELD) then
687          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
688          thisStat = -1
689          haveSaneForceField = .false.
690          return
691       endif
692    endif
682  
683      if (FF_uses_EAM) then
684         call init_EAM_FF(my_status)
# Line 701 | Line 690 | contains
690         end if
691      endif
692  
704    if (FF_uses_GayBerne) then
705       call check_gb_pair_FF(my_status)
706       if (my_status .ne. 0) then
707          thisStat = -1
708          haveSaneForceField = .false.
709          return
710       endif
711    endif
712
693      if (.not. haveNeighborList) then
694         !! Create neighbor lists
695         call expandNeighborList(nLocal, my_status)
# Line 743 | Line 723 | contains
723  
724      !! Stress Tensor
725      real( kind = dp), dimension(9) :: tau  
726 <    real ( kind = dp ) :: pot
726 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
727      logical ( kind = 2) :: do_pot_c, do_stress_c
728      logical :: do_pot
729      logical :: do_stress
730      logical :: in_switching_region
731   #ifdef IS_MPI
732 <    real( kind = DP ) :: pot_local
732 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
733      integer :: nAtomsInRow
734      integer :: nAtomsInCol
735      integer :: nprocs
# Line 764 | Line 744 | contains
744      integer :: nlist
745      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
746      real( kind = DP ) :: sw, dswdr, swderiv, mf
747 +    real( kind = DP ) :: rVal
748      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
749      real(kind=dp) :: rfpot, mu_i, virial
750      integer :: me_i, me_j, n_in_i, n_in_j
# Line 774 | Line 755 | contains
755      integer :: propPack_i, propPack_j
756      integer :: loopStart, loopEnd, loop
757      integer :: iHash
758 +    integer :: i1
759 +    logical :: indirect_only
760    
761  
762      !! initialize local variables  
# Line 940 | Line 923 | contains
923  
924                        atom2 = groupListCol(jb)
925  
926 <                      if (skipThisPair(atom1, atom2)) cycle inner
926 >                      indirect_only = .false.
927 >    
928 >                      if (skipThisPair(atom1, atom2)) then
929 >                         if (electrostaticSummationMethod.ne.REACTION_FIELD) then
930 >                            cycle inner
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)
# Line 968 | Line 960 | contains
960                        else
961   #ifdef IS_MPI                      
962                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
963 <                              do_pot, &
964 <                              eFrame, A, f, t, pot_local, vpair, fpair)
963 >                              do_pot, eFrame, A, f, t, pot_local, vpair, &
964 >                              fpair, d_grp, rgrp, indirect_only)
965   #else
966                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
967 <                              do_pot,  &
968 <                              eFrame, A, f, t, pot, vpair, fpair)
967 >                              do_pot, eFrame, A, f, t, pot, vpair, fpair, &
968 >                              d_grp, rgrp, indirect_only)
969   #endif
970  
971                           vij = vij + vpair
# Line 1082 | Line 1074 | contains
1074  
1075      if (do_pot) then
1076         ! scatter/gather pot_row into the members of my column
1077 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1078 <
1077 >       do i = 1,LR_POT_TYPES
1078 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1079 >       end do
1080         ! scatter/gather pot_local into all other procs
1081         ! add resultant to get total pot
1082         do i = 1, nlocal
1083 <          pot_local = pot_local + pot_Temp(i)
1083 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1084 >               + pot_Temp(1:LR_POT_TYPES,i)
1085         enddo
1086  
1087         pot_Temp = 0.0_DP
1088 <
1089 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1088 >       do i = 1,LR_POT_TYPES
1089 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1090 >       end do
1091         do i = 1, nlocal
1092 <          pot_local = pot_local + pot_Temp(i)
1092 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1093 >               + pot_Temp(1:LR_POT_TYPES,i)
1094         enddo
1095  
1096      endif
1097   #endif
1098  
1099 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1099 >    if (SIM_requires_postpair_calc) then
1100 >       do i = 1, nlocal            
1101 >          
1102 >          ! we loop only over the local atoms, so we don't need row and column
1103 >          ! lookups for the types
1104 >          
1105 >          me_i = atid(i)
1106 >          
1107 >          ! is the atom electrostatic?  See if it would have an
1108 >          ! electrostatic interaction with itself
1109 >          iHash = InteractionHash(me_i,me_i)
1110  
1111 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1106 <
1111 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1112   #ifdef IS_MPI
1113 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1114 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1110 <          do i = 1,nlocal
1111 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1112 <          end do
1113 < #endif
1114 <
1115 <          do i = 1, nLocal
1116 <
1117 <             rfpot = 0.0_DP
1118 < #ifdef IS_MPI
1119 <             me_i = atid_row(i)
1113 >             call rf_self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1114 >                  t, do_pot)
1115   #else
1116 <             me_i = atid(i)
1116 >             call rf_self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1117 >                  t, do_pot)
1118   #endif
1119 <             iHash = InteractionHash(me_i,me_j)
1120 <            
1121 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1119 >          endif
1120 >  
1121 >          ! loop over the excludes to accumulate any additional RF components
1122  
1123 <                mu_i = getDipoleMoment(me_i)
1124 <
1125 <                !! The reaction field needs to include a self contribution
1126 <                !! to the field:
1127 <                call accumulate_self_rf(i, mu_i, eFrame)
1128 <                !! Get the reaction field contribution to the
1129 <                !! potential and torques:
1130 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1123 >          do i1 = 1, nSkipsForAtom(i)
1124 >             j = skipsForAtom(i, i1)
1125 >            
1126 >             ! prevent overcounting of the skips
1127 >             if (i.lt.j) then
1128 >                call get_interatomic_vector(q(:,i), &
1129 >                     q(:,j), d_atm, ratmsq)
1130 >                rVal = dsqrt(ratmsq)
1131 >                call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1132 >                     in_switching_region)
1133   #ifdef IS_MPI
1134 <                pot_local = pot_local + rfpot
1134 >                call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, vpair, &
1135 >                     pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1136   #else
1137 <                pot = pot + rfpot
1138 <
1137 >                call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, vpair, &
1138 >                     pot(ELECTROSTATIC_POT), f, t, do_pot)
1139   #endif
1140               endif
1141            enddo
1142 <       endif
1142 >       enddo      
1143      endif
1144 <
1146 <
1144 >    
1145   #ifdef IS_MPI
1146 <
1146 >    
1147      if (do_pot) then
1148 <       pot = pot + pot_local
1149 <       !! we assume the c code will do the allreduce to get the total potential
1152 <       !! we could do it right here if we needed to...
1148 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1149 >            mpi_comm_world,mpi_err)            
1150      endif
1151 <
1151 >    
1152      if (do_stress) then
1153         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1154              mpi_comm_world,mpi_err)
1155         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1156              mpi_comm_world,mpi_err)
1157      endif
1158 <
1158 >    
1159   #else
1160 <
1160 >    
1161      if (do_stress) then
1162         tau = tau_Temp
1163         virial = virial_Temp
1164      endif
1165 <
1165 >    
1166   #endif
1167 <
1167 >    
1168    end subroutine do_force_loop
1169  
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, indirect_only)
1172  
1173 <    real( kind = dp ) :: pot, vpair, sw
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(nLocal)   :: mfact
1177      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1182 | Line 1180 | contains
1180      real( kind = dp ), dimension(3,nLocal) :: t
1181  
1182      logical, intent(inout) :: do_pot
1183 +    logical, intent(inout) :: indirect_only
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 1204 | Line 1205 | contains
1205  
1206      iHash = InteractionHash(me_i, me_j)
1207  
1208 <    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1209 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1210 <    endif
1211 <
1211 <    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1212 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213 <            pot, eFrame, f, t, do_pot)
1214 <
1215 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1216 <
1217 <          ! CHECK ME (RF needs to know about all electrostatic types)
1218 <          call accumulate_rf(i, j, r, eFrame, sw)
1219 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1208 >    if (indirect_only) then
1209 >       if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1210 >          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1211 >               pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot, indirect_only)
1212         endif
1213 <
1214 <    endif
1215 <
1216 <    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1217 <       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1218 <            pot, A, f, t, do_pot)
1227 <    endif
1228 <
1229 <    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1230 <       call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1231 <            pot, A, f, t, do_pot)
1232 <    endif
1233 <
1234 <    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1235 <       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1236 <            pot, A, f, t, do_pot)
1237 <    endif
1238 <    
1239 <    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1240 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1241 < !           pot, A, f, t, do_pot)
1242 <    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, EAM_PAIR).ne.0 ) then      
1221 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1222 <            do_pot)
1223 <    endif
1224 <
1225 <    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1226 <       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1227 <            pot, A, f, t, do_pot)
1228 <    endif
1229 <
1230 <    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1231 <       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1232 <            pot, A, f, t, do_pot)
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
1259      endif
1260      
1261    end subroutine do_pair
# Line 1261 | Line 1263 | contains
1263    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1264         do_pot, do_stress, eFrame, A, f, t, pot)
1265  
1266 <    real( kind = dp ) :: pot, sw
1266 >    real( kind = dp ) :: sw
1267 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1268      real( kind = dp ), dimension(9,nLocal) :: eFrame
1269      real (kind=dp), dimension(9,nLocal) :: A
1270      real (kind=dp), dimension(3,nLocal) :: f
# Line 1296 | Line 1299 | contains
1299  
1300    subroutine do_preforce(nlocal,pot)
1301      integer :: nlocal
1302 <    real( kind = dp ) :: pot
1302 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1303  
1304      if (FF_uses_EAM .and. SIM_uses_EAM) then
1305 <       call calc_EAM_preforce_Frho(nlocal,pot)
1305 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1306      endif
1307  
1308  
# Line 1492 | Line 1495 | contains
1495      doesit = FF_uses_EAM
1496    end function FF_RequiresPrepairCalc
1497  
1495  function FF_RequiresPostpairCalc() result(doesit)
1496    logical :: doesit
1497    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1498  end function FF_RequiresPostpairCalc
1499
1498   #ifdef PROFILE
1499    function getforcetime() result(totalforcetime)
1500      real(kind=dp) :: totalforcetime

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines