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 2355 by chuckv, Wed Oct 12 18:59:16 2005 UTC vs.
Revision 2390 by chrisfen, Wed Oct 19 19:24:40 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.54 2005-10-12 18:59:16 chuckv Exp $, $Date: 2005-10-12 18:59:16 $, $Name: not supported by cvs2svn $, $Revision: 1.54 $
48 > !! @version $Id: doForces.F90,v 1.61 2005-10-19 19:24:29 chrisfen Exp $, $Date: 2005-10-19 19:24:29 $, $Name: not supported by cvs2svn $, $Revision: 1.61 $
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 681 | Line 680 | contains
680  
681      haveSaneForceField = .true.
682  
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 –
683      if (FF_uses_EAM) then
684         call init_EAM_FF(my_status)
685         if (my_status /= 0) then
# 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 1088 | Line 1062 | contains
1062  
1063      if (do_pot) then
1064         ! scatter/gather pot_row into the members of my column
1065 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1066 <
1065 >       do i = 1,LR_POT_TYPES
1066 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1067 >       end do
1068         ! scatter/gather pot_local into all other procs
1069         ! add resultant to get total pot
1070         do i = 1, nlocal
1071 <          pot_local = pot_local + pot_Temp(i)
1071 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1072 >               + pot_Temp(1:LR_POT_TYPES,i)
1073         enddo
1074  
1075         pot_Temp = 0.0_DP
1076 <
1077 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1076 >       do i = 1,LR_POT_TYPES
1077 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1078 >       end do
1079         do i = 1, nlocal
1080 <          pot_local = pot_local + pot_Temp(i)
1080 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1081 >               + pot_Temp(1:LR_POT_TYPES,i)
1082         enddo
1083  
1084      endif
1085   #endif
1086  
1087 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1088 <
1111 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1112 <
1087 >    if (SIM_requires_postpair_calc) then
1088 >      
1089   #ifdef IS_MPI
1090 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1091 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1092 <          do i = 1,nlocal
1093 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1094 <          end do
1090 >       call scatter(rf_Row,rf,plan_atom_row_3d)
1091 >       call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1092 >       do i = 1,nlocal
1093 >          rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1094 >       end do
1095   #endif
1096  
1097 <          do i = 1, nLocal
1098 <
1099 <             rfpot = 0.0_DP
1097 >       do i = 1, nLocal
1098 >          
1099 >          rfpot = 0.0_DP
1100   #ifdef IS_MPI
1101 <             me_i = atid_row(i)
1101 >          me_i = atid_row(i)
1102   #else
1103 <             me_i = atid(i)
1103 >          me_i = atid(i)
1104   #endif
1105 <             iHash = InteractionHash(me_i,me_j)
1105 >          iHash = InteractionHash(me_i,me_j)
1106 >          
1107 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1108              
1109 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1110 <
1111 <                mu_i = getDipoleMoment(me_i)
1112 <
1113 <                !! The reaction field needs to include a self contribution
1114 <                !! to the field:
1115 <                call accumulate_self_rf(i, mu_i, eFrame)
1116 <                !! Get the reaction field contribution to the
1139 <                !! potential and torques:
1140 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1109 >             mu_i = getDipoleMoment(me_i)
1110 >            
1111 >             !! The reaction field needs to include a self contribution
1112 >             !! to the field:
1113 >             call accumulate_self_rf(i, mu_i, eFrame)
1114 >             !! Get the reaction field contribution to the
1115 >             !! potential and torques:
1116 >             call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1117   #ifdef IS_MPI
1118 <                pot_local(RF_POT) = pot_local(RF_POT) + rfpot
1118 >             pot_local(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
1119   #else
1120 <                pot(RF_POT) = pot(RF_POT) + rfpot
1121 <
1120 >             pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1121 >            
1122   #endif
1123 <             endif
1124 <          enddo
1149 <       endif
1123 >          endif
1124 >       enddo
1125      endif
1126 <
1152 <
1126 >    
1127   #ifdef IS_MPI
1128  
1129      if (do_pot) then
1130 <       pot(1:SIZE_POT_ARRAY) = pot(1:SIZE_POT_ARRAY) &
1131 <            + pot_local(1:SIZE_POT_ARRAY)
1130 >       pot(1:LR_POT_TYPES) = pot(1:LR_POT_TYPES) &
1131 >            + pot_local(1:LR_POT_TYPES)
1132         !! we assume the c code will do the allreduce to get the total potential
1133         !! we could do it right here if we needed to...
1134      endif
# Line 1181 | Line 1155 | contains
1155         eFrame, A, f, t, pot, vpair, fpair)
1156  
1157      real( kind = dp ) :: vpair, sw
1158 <    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1158 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1159      real( kind = dp ), dimension(3) :: fpair
1160      real( kind = dp ), dimension(nLocal)   :: mfact
1161      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1213 | Line 1187 | contains
1187      iHash = InteractionHash(me_i, me_j)
1188  
1189      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1190 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(LJ_POT), f, do_pot)
1190 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1191 >            pot(VDW_POT), f, do_pot)
1192      endif
1193  
1194      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
# Line 1231 | Line 1206 | contains
1206  
1207      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1208         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1209 <            pot(STICKY_POT), A, f, t, do_pot)
1209 >            pot(HB_POT), A, f, t, do_pot)
1210      endif
1211  
1212      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1213         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1214 <            pot(STICKYPOWER_POT), A, f, t, do_pot)
1214 >            pot(HB_POT), A, f, t, do_pot)
1215      endif
1216  
1217      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1218         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1219 <            pot(GAYBERNE_POT), A, f, t, do_pot)
1219 >            pot(VDW_POT), A, f, t, do_pot)
1220      endif
1221      
1222      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1223 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1224 < !           pot(GAYBERNE_LJ_POT), A, f, t, do_pot)
1223 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1224 >            pot(VDW_POT), A, f, t, do_pot)
1225      endif
1226  
1227      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1228 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(EAM_POT), f, &
1229 <            do_pot)
1228 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1229 >            pot(METALLIC_POT), f, do_pot)
1230      endif
1231  
1232      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1233         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1234 <            pot(SHAPE_POT), A, f, t, do_pot)
1234 >            pot(VDW_POT), A, f, t, do_pot)
1235      endif
1236  
1237      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1238         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1239 <            pot(SHAPE_LJ_POT), A, f, t, do_pot)
1239 >            pot(VDW_POT), A, f, t, do_pot)
1240      endif
1241      
1242    end subroutine do_pair
# Line 1270 | Line 1245 | contains
1245         do_pot, do_stress, eFrame, A, f, t, pot)
1246  
1247      real( kind = dp ) :: sw
1248 <    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1248 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1249      real( kind = dp ), dimension(9,nLocal) :: eFrame
1250      real (kind=dp), dimension(9,nLocal) :: A
1251      real (kind=dp), dimension(3,nLocal) :: f
# Line 1305 | Line 1280 | contains
1280  
1281    subroutine do_preforce(nlocal,pot)
1282      integer :: nlocal
1283 <    real( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
1283 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1284  
1285      if (FF_uses_EAM .and. SIM_uses_EAM) then
1286 <       call calc_EAM_preforce_Frho(nlocal,pot(EAM_POT))
1286 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1287      endif
1288  
1289  
# Line 1501 | Line 1476 | contains
1476      doesit = FF_uses_EAM
1477    end function FF_RequiresPrepairCalc
1478  
1504 –  function FF_RequiresPostpairCalc() result(doesit)
1505 –    logical :: doesit
1506 –    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1507 –  end function FF_RequiresPostpairCalc
1508 –
1479   #ifdef PROFILE
1480    function getforcetime() result(totalforcetime)
1481      real(kind=dp) :: totalforcetime

Diff Legend

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