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 2405 by chrisfen, Tue Nov 1 19:24:57 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.65 2005-11-01 19:24:54 chrisfen Exp $, $Date: 2005-11-01 19:24:54 $, $Name: not supported by cvs2svn $, $Revision: 1.65 $
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          
199      do i = 1, nAtypes
# Line 639 | Line 644 | contains
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
646      integer, pointer :: MatchList(:) => null()
642    real(kind=dp) :: rcut, rrf, rt, dielect
647  
648      !! assume things are copacetic, unless they aren't
649      thisStat = 0
# Line 674 | Line 678 | contains
678  
679  
680      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
681  
682      if (FF_uses_EAM) then
683         call init_EAM_FF(my_status)
# Line 701 | Line 689 | contains
689         end if
690      endif
691  
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
692      if (.not. haveNeighborList) then
693         !! Create neighbor lists
694         call expandNeighborList(nLocal, my_status)
# Line 743 | Line 722 | contains
722  
723      !! Stress Tensor
724      real( kind = dp), dimension(9) :: tau  
725 <    real ( kind = dp ) :: pot
725 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
726      logical ( kind = 2) :: do_pot_c, do_stress_c
727      logical :: do_pot
728      logical :: do_stress
729      logical :: in_switching_region
730   #ifdef IS_MPI
731 <    real( kind = DP ) :: pot_local
731 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
732      integer :: nAtomsInRow
733      integer :: nAtomsInCol
734      integer :: nprocs
# Line 764 | 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 774 | 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 939 | Line 920 | contains
920                     inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
921  
922                        atom2 = groupListCol(jb)
923 <
924 <                      if (skipThisPair(atom1, atom2)) cycle inner
925 <
923 >    
924 >                      if (skipThisPair(atom1, atom2))  cycle inner
925 >                      
926                        if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
927                           d_atm(1:3) = d_grp(1:3)
928                           ratmsq = rgrpsq
# Line 968 | 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 1082 | Line 1063 | contains
1063  
1064      if (do_pot) then
1065         ! scatter/gather pot_row into the members of my column
1066 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1067 <
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 = pot_local + pot_Temp(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 <
1078 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
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 = pot_local + pot_Temp(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
1089 <
1090 <       if (electrostaticSummationMethod == REACTION_FIELD) 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 >          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)
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)
1102 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1103 >                  t, do_pot)
1104   #else
1105 <             me_i = atid(i)
1105 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1106 >                  t, do_pot)
1107   #endif
1108 <             iHash = InteractionHash(me_i,me_j)
1108 >          endif
1109 >  
1110 >          
1111 > !!$          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1112              
1113 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1114 <
1115 <                mu_i = getDipoleMoment(me_i)
1116 <
1117 <                !! The reaction field needs to include a self contribution
1118 <                !! to the field:
1119 <                call accumulate_self_rf(i, mu_i, eFrame)
1120 <                !! Get the reaction field contribution to the
1121 <                !! potential and torques:
1122 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1113 >             ! loop over the excludes to accumulate RF stuff we've
1114 >             ! left out of the normal pair loop
1115 >            
1116 >             do i1 = 1, nSkipsForAtom(i)
1117 >                j = skipsForAtom(i, i1)
1118 >                
1119 >                ! prevent overcounting of the skips
1120 >                if (i.lt.j) then
1121 >                   call get_interatomic_vector(q(:,i), &
1122 >                        q(:,j), d_atm, ratmsq)
1123 >                   rVal = dsqrt(ratmsq)
1124 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1125 >                        in_switching_region)
1126   #ifdef IS_MPI
1127 <                pot_local = pot_local + rfpot
1127 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1128 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1129   #else
1130 <                pot = pot + rfpot
1131 <
1130 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1131 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1132   #endif
1133 <             endif
1134 <          enddo
1135 <       endif
1133 >                endif
1134 >             enddo
1135 > !!$          endif
1136 >       enddo
1137      endif
1138 <
1146 <
1138 >    
1139   #ifdef IS_MPI
1140 <
1140 >    
1141      if (do_pot) then
1142 <       pot = pot + pot_local
1143 <       !! 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...
1142 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1143 >            mpi_comm_world,mpi_err)            
1144      endif
1145 <
1145 >    
1146      if (do_stress) then
1147         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1148              mpi_comm_world,mpi_err)
1149         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1150              mpi_comm_world,mpi_err)
1151      endif
1152 <
1152 >    
1153   #else
1154 <
1154 >    
1155      if (do_stress) then
1156         tau = tau_Temp
1157         virial = virial_Temp
1158      endif
1159 <
1159 >    
1160   #endif
1161 <
1161 >    
1162    end subroutine do_force_loop
1163  
1164    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1165 <       eFrame, A, f, t, pot, vpair, fpair)
1165 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1166  
1167 <    real( kind = dp ) :: pot, vpair, sw
1167 >    real( kind = dp ) :: vpair, sw
1168 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1169      real( kind = dp ), dimension(3) :: fpair
1170      real( kind = dp ), dimension(nLocal)   :: mfact
1171      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1184 | Line 1176 | contains
1176      logical, intent(inout) :: do_pot
1177      integer, intent(in) :: i, j
1178      real ( kind = dp ), intent(inout) :: rijsq
1179 <    real ( kind = dp )                :: r
1179 >    real ( kind = dp ), intent(inout) :: r_grp
1180      real ( kind = dp ), intent(inout) :: d(3)
1181 +    real ( kind = dp ), intent(inout) :: d_grp(3)
1182 +    real ( kind = dp ) :: r
1183      integer :: me_i, me_j
1184  
1185      integer :: iHash
# Line 1203 | Line 1197 | contains
1197   #endif
1198  
1199      iHash = InteractionHash(me_i, me_j)
1200 <
1200 >    
1201      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1202 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1202 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1203 >            pot(VDW_POT), f, do_pot)
1204      endif
1205 <
1205 >    
1206      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1207         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1208 <            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)
1220 <       endif
1221 <
1208 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1209      endif
1210 <
1210 >    
1211      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1212         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213 <            pot, A, f, t, do_pot)
1213 >            pot(HB_POT), A, f, t, do_pot)
1214      endif
1215 <
1215 >    
1216      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1217         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1218 <            pot, A, f, t, do_pot)
1218 >            pot(HB_POT), A, f, t, do_pot)
1219      endif
1220 <
1220 >    
1221      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1222         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1223 <            pot, A, f, t, do_pot)
1223 >            pot(VDW_POT), A, f, t, do_pot)
1224      endif
1225      
1226      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1227 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1228 < !           pot, A, f, t, do_pot)
1227 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1228 >            pot(VDW_POT), A, f, t, do_pot)
1229      endif
1230 <
1230 >    
1231      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1232 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1233 <            do_pot)
1232 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1233 >            pot(METALLIC_POT), f, do_pot)
1234      endif
1235 <
1235 >    
1236      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1237         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1238 <            pot, A, f, t, do_pot)
1238 >            pot(VDW_POT), A, f, t, do_pot)
1239      endif
1240 <
1240 >    
1241      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1242         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1243 <            pot, A, f, t, do_pot)
1243 >            pot(VDW_POT), A, f, t, do_pot)
1244      endif
1245 <    
1245 >    
1246    end subroutine do_pair
1247  
1248    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1249         do_pot, do_stress, eFrame, A, f, t, pot)
1250  
1251 <    real( kind = dp ) :: pot, sw
1251 >    real( kind = dp ) :: sw
1252 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1253      real( kind = dp ), dimension(9,nLocal) :: eFrame
1254      real (kind=dp), dimension(9,nLocal) :: A
1255      real (kind=dp), dimension(3,nLocal) :: f
# Line 1296 | Line 1284 | contains
1284  
1285    subroutine do_preforce(nlocal,pot)
1286      integer :: nlocal
1287 <    real( kind = dp ) :: pot
1287 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1288  
1289      if (FF_uses_EAM .and. SIM_uses_EAM) then
1290 <       call calc_EAM_preforce_Frho(nlocal,pot)
1290 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1291      endif
1292  
1293  
# Line 1385 | Line 1373 | contains
1373      pot_Col = 0.0_dp
1374      pot_Temp = 0.0_dp
1375  
1388    rf_Row = 0.0_dp
1389    rf_Col = 0.0_dp
1390    rf_Temp = 0.0_dp
1391
1376   #endif
1377  
1378      if (FF_uses_EAM .and. SIM_uses_EAM) then
# Line 1492 | Line 1476 | contains
1476      doesit = FF_uses_EAM
1477    end function FF_RequiresPrepairCalc
1478  
1495  function FF_RequiresPostpairCalc() result(doesit)
1496    logical :: doesit
1497    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1498  end function FF_RequiresPostpairCalc
1499
1479   #ifdef PROFILE
1480    function getforcetime() result(totalforcetime)
1481      real(kind=dp) :: totalforcetime

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines