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 2356 by chuckv, Wed Oct 12 19:55:26 2005 UTC vs.
Revision 2394 by chrisfen, Sun Oct 23 21:08:08 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.55 2005-10-12 19:55:26 chuckv Exp $, $Date: 2005-10-12 19:55:26 $, $Name: not supported by cvs2svn $, $Revision: 1.55 $
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 $
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 680 | Line 679 | contains
679  
680  
681      haveSaneForceField = .true.
683 –
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
682  
683      if (FF_uses_EAM) then
684         call init_EAM_FF(my_status)
# Line 707 | Line 690 | contains
690         end if
691      endif
692  
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
717 –    endif
718 –
693      if (.not. haveNeighborList) then
694         !! Create neighbor lists
695         call expandNeighborList(nLocal, my_status)
# Line 749 | Line 723 | contains
723  
724      !! Stress Tensor
725      real( kind = dp), dimension(9) :: tau  
726 <    real ( kind = dp ),dimension(POT_ARRAY_SIZE) :: 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 ), dimension(POT_ARRAY_SIZE) :: pot_local
732 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
733      integer :: nAtomsInRow
734      integer :: nAtomsInCol
735      integer :: nprocs
# Line 780 | Line 754 | contains
754      integer :: propPack_i, propPack_j
755      integer :: loopStart, loopEnd, loop
756      integer :: iHash
757 +    integer :: ig
758    
759  
760      !! initialize local variables  
# Line 974 | Line 949 | contains
949                        else
950   #ifdef IS_MPI                      
951                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
952 <                              do_pot, &
953 <                              eFrame, A, f, t, pot_local, vpair, fpair)
952 >                              do_pot, eFrame, A, f, t, pot_local, vpair, &
953 >                              fpair, d_grp, rgrp)
954   #else
955                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
956 <                              do_pot,  &
957 <                              eFrame, A, f, t, pot, vpair, fpair)
956 >                              do_pot, eFrame, A, f, t, pot, vpair, fpair, &
957 >                              d_grp, rgrp)
958   #endif
959  
960                           vij = vij + vpair
# Line 1088 | Line 1063 | contains
1063  
1064      if (do_pot) then
1065         ! scatter/gather pot_row into the members of my column
1066 <       do i = 1,POT_ARRAY_SIZE
1066 >       do i = 1,LR_POT_TYPES
1067            call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1068         end do
1069         ! scatter/gather pot_local into all other procs
1070         ! add resultant to get total pot
1071         do i = 1, nlocal
1072 <          pot_local(1:POT_ARRAY_SIZE) = pot_local(1:POT_ARRAY_SIZE &
1073 <               + pot_Temp(1:POT_ARRAY_SIZE,i)
1072 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1073 >               + pot_Temp(1:LR_POT_TYPES,i)
1074         enddo
1075  
1076         pot_Temp = 0.0_DP
1077 <       do i = 1,POT_ARRAY_SIZE
1077 >       do i = 1,LR_POT_TYPES
1078            call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1079         end do
1080         do i = 1, nlocal
1081 <          pot_local(1:POT_ARRAY_SIZE) = pot_local(1:POT_ARRAY_SIZE)&
1082 <               + pot_Temp(1:POT_ARRAY_SIZE,i)
1081 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1082 >               + pot_Temp(1:LR_POT_TYPES,i)
1083         enddo
1084  
1085      endif
1086   #endif
1087  
1088 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1088 >    if (SIM_requires_postpair_calc) then
1089 >       do i = 1, nlocal            
1090 >          
1091 >          ! we loop only over the local atoms, so we don't need row and column
1092 >          ! lookups for the types
1093  
1094 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1095 <
1094 >          me_i = atid(i)
1095 >          
1096 >          ! is the atom electrostatic?  See if it would have an
1097 >          ! electrostatic interaction with itself
1098 >          iHash = InteractionHash(me_i,me_i)
1099 >          
1100 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1101   #ifdef IS_MPI
1102 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1103 <          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)
1102 >             call rf_self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1103 >                  t, do_pot)
1104   #else
1105 <             me_i = atid(i)
1105 >             call rf_self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1106 >                  t, do_pot)
1107   #endif
1108 <             iHash = InteractionHash(me_i,me_j)
1109 <            
1135 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1136 <
1137 <                mu_i = getDipoleMoment(me_i)
1138 <
1139 <                !! The reaction field needs to include a self contribution
1140 <                !! to the field:
1141 <                call accumulate_self_rf(i, mu_i, eFrame)
1142 <                !! Get the reaction field contribution to the
1143 <                !! potential and torques:
1144 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1145 < #ifdef IS_MPI
1146 <                pot_local(RF_POT) = pot_local(RF_POT) + rfpot
1147 < #else
1148 <                pot(RF_POT) = pot(RF_POT) + rfpot
1149 <
1150 < #endif
1151 <             endif
1152 <          enddo
1153 <       endif
1108 >          endif
1109 >       enddo
1110      endif
1111 <
1156 <
1111 >    
1112   #ifdef IS_MPI
1113 <
1113 >    
1114      if (do_pot) then
1115 <       pot(1:SIZE_POT_ARRAY) = pot(1:SIZE_POT_ARRAY) &
1116 <            + pot_local(1:SIZE_POT_ARRAY)
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...
1115 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1116 >            mpi_comm_world,mpi_err)            
1117      endif
1118 <
1118 >    
1119      if (do_stress) then
1120         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1121              mpi_comm_world,mpi_err)
1122         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1123              mpi_comm_world,mpi_err)
1124      endif
1125 <
1125 >    
1126   #else
1127 <
1127 >    
1128      if (do_stress) then
1129         tau = tau_Temp
1130         virial = virial_Temp
1131      endif
1132 <
1132 >    
1133   #endif
1134 <
1134 >    
1135    end subroutine do_force_loop
1136  
1137    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1138 <       eFrame, A, f, t, pot, vpair, fpair)
1138 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1139  
1140      real( kind = dp ) :: vpair, sw
1141 <    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1141 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1142      real( kind = dp ), dimension(3) :: fpair
1143      real( kind = dp ), dimension(nLocal)   :: mfact
1144      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1196 | Line 1149 | contains
1149      logical, intent(inout) :: do_pot
1150      integer, intent(in) :: i, j
1151      real ( kind = dp ), intent(inout) :: rijsq
1152 <    real ( kind = dp )                :: r
1152 >    real ( kind = dp ), intent(inout) :: r_grp
1153      real ( kind = dp ), intent(inout) :: d(3)
1154 +    real ( kind = dp ), intent(inout) :: d_grp(3)
1155 +    real ( kind = dp ) :: r
1156      integer :: me_i, me_j
1157  
1158      integer :: iHash
# Line 1217 | Line 1172 | contains
1172      iHash = InteractionHash(me_i, me_j)
1173  
1174      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1175 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(LJ_POT), f, do_pot)
1175 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1176 >            pot(VDW_POT), f, do_pot)
1177      endif
1178  
1179      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1180         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1181              pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1226 –
1227 –       if (electrostaticSummationMethod == REACTION_FIELD) then
1228 –
1229 –          ! CHECK ME (RF needs to know about all electrostatic types)
1230 –          call accumulate_rf(i, j, r, eFrame, sw)
1231 –          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1232 –       endif
1233 –
1182      endif
1183  
1184      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1185         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1186 <            pot(STICKY_POT), A, f, t, do_pot)
1186 >            pot(HB_POT), A, f, t, do_pot)
1187      endif
1188  
1189      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1190         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1191 <            pot(STICKYPOWER_POT), A, f, t, do_pot)
1191 >            pot(HB_POT), A, f, t, do_pot)
1192      endif
1193  
1194      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1195         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1196 <            pot(GAYBERNE_POT), A, f, t, do_pot)
1196 >            pot(VDW_POT), A, f, t, do_pot)
1197      endif
1198      
1199      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1200 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1201 < !           pot(GAYBERNE_LJ_POT), A, f, t, do_pot)
1200 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1201 >            pot(VDW_POT), A, f, t, do_pot)
1202      endif
1203  
1204      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1205 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(EAM_POT), f, &
1206 <            do_pot)
1205 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1206 >            pot(METALLIC_POT), f, do_pot)
1207      endif
1208  
1209      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1210         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1211 <            pot(SHAPE_POT), A, f, t, do_pot)
1211 >            pot(VDW_POT), A, f, t, do_pot)
1212      endif
1213  
1214      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1215         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1216 <            pot(SHAPE_LJ_POT), A, f, t, do_pot)
1216 >            pot(VDW_POT), A, f, t, do_pot)
1217      endif
1218      
1219    end subroutine do_pair
# Line 1274 | Line 1222 | contains
1222         do_pot, do_stress, eFrame, A, f, t, pot)
1223  
1224      real( kind = dp ) :: sw
1225 <    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1225 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1226      real( kind = dp ), dimension(9,nLocal) :: eFrame
1227      real (kind=dp), dimension(9,nLocal) :: A
1228      real (kind=dp), dimension(3,nLocal) :: f
# Line 1309 | Line 1257 | contains
1257  
1258    subroutine do_preforce(nlocal,pot)
1259      integer :: nlocal
1260 <    real( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
1260 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1261  
1262      if (FF_uses_EAM .and. SIM_uses_EAM) then
1263 <       call calc_EAM_preforce_Frho(nlocal,pot(EAM_POT))
1263 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1264      endif
1265  
1266  
# Line 1505 | Line 1453 | contains
1453      doesit = FF_uses_EAM
1454    end function FF_RequiresPrepairCalc
1455  
1508 –  function FF_RequiresPostpairCalc() result(doesit)
1509 –    logical :: doesit
1510 –    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1511 –  end function FF_RequiresPostpairCalc
1512 –
1456   #ifdef PROFILE
1457    function getforcetime() result(totalforcetime)
1458      real(kind=dp) :: totalforcetime

Diff Legend

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