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 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.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.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 705 | Line 688 | contains
688            haveSaneForceField = .false.
689            return
690         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
691      endif
692  
693      if (.not. haveNeighborList) then
# Line 1109 | Line 1083 | contains
1083  
1084      endif
1085   #endif
1112
1113    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1086  
1087 <       if (electrostaticSummationMethod == REACTION_FIELD) then
1088 <
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)
1130 < #else
1131 <             me_i = atid(i)
1132 < #endif
1133 <             iHash = InteractionHash(me_i,me_j)
1134 <            
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(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
1101 >          me_i = atid_row(i)
1102   #else
1103 <                pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1149 <
1103 >          me_i = atid(i)
1104   #endif
1105 <             endif
1106 <          enddo
1107 <       endif
1105 >          iHash = InteractionHash(me_i,me_j)
1106 >          
1107 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1108 >            
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(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
1119 > #else
1120 >             pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1121 >            
1122 > #endif
1123 >          endif
1124 >       enddo
1125      endif
1126 <
1156 <
1126 >    
1127   #ifdef IS_MPI
1128  
1129      if (do_pot) then
# Line 1506 | Line 1476 | contains
1476      doesit = FF_uses_EAM
1477    end function FF_RequiresPrepairCalc
1478  
1509  function FF_RequiresPostpairCalc() result(doesit)
1510    logical :: doesit
1511    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1512  end function FF_RequiresPostpairCalc
1513
1479   #ifdef PROFILE
1480    function getforcetime() result(totalforcetime)
1481      real(kind=dp) :: totalforcetime

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines