ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2260 by chuckv, Mon Jun 27 22:21:37 2005 UTC vs.
Revision 2262 by chuckv, Sun Jul 3 20:53:43 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.21 2005-06-27 22:21:37 chuckv Exp $, $Date: 2005-06-27 22:21:37 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $
48 > !! @version $Id: doForces.F90,v 1.23 2005-07-03 20:53:43 chuckv Exp $, $Date: 2005-07-03 20:53:43 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $
49  
50  
51   module doForces
# Line 81 | Line 81 | module doForces
81    logical, save :: haveRlist = .false.
82    logical, save :: haveNeighborList = .false.
83    logical, save :: haveSIMvariables = .false.
84 –  logical, save :: havePropertyMap = .false.
84    logical, save :: haveSaneForceField = .false.
85 +  logical, save :: haveInteractionMap = .false.
86  
87    logical, save :: FF_uses_DirectionalAtoms
88    logical, save :: FF_uses_LennardJones
# Line 122 | Line 122 | module doForces
122    public :: init_FF
123    public :: do_force_loop
124   !  public :: setRlistDF
125 +  !public :: addInteraction
126 +  !public :: setInteractionHash
127 +  !public :: getInteractionHash
128 +  public :: createInteractionMap
129 +  public :: createRcuts
130  
131   #ifdef PROFILE
132    public :: getforcetime
# Line 130 | Line 135 | module doForces
135    integer :: nLoops
136   #endif
137  
133 –  type :: Properties
134 –     logical :: is_Directional   = .false.
135 –     logical :: is_LennardJones  = .false.
136 –     logical :: is_Electrostatic = .false.
137 –     logical :: is_Charge        = .false.
138 –     logical :: is_Dipole        = .false.
139 –     logical :: is_Quadrupole    = .false.
140 –     logical :: is_Sticky        = .false.
141 –     logical :: is_StickyPower   = .false.
142 –     logical :: is_GayBerne      = .false.
143 –     logical :: is_EAM           = .false.
144 –     logical :: is_Shape         = .false.
145 –     logical :: is_FLARB         = .false.
146 –  end type Properties
147 –
148 –  type(Properties), dimension(:),allocatable :: PropertyMap
149 –
150 –
151 –  
138    type, public :: Interaction
139       integer :: InteractionHash
140 <     real(kind=dp) :: rCut
140 >     real(kind=dp) :: rList = 0.0_dp
141 >     real(kind=dp) :: rListSq = 0.0_dp
142    end type Interaction
143    
144 <  type(Interaction), public, dimension(:,:), allocatable :: InteractionMap
144 >  type(Interaction), dimension(:,:),allocatable :: InteractionMap
145    
159 –  !public :: addInteraction
160 –  !public :: setInteractionHash
161 –  !public :: getInteractionHash
162 –  public :: createInteractionMap
146  
147 +  
148   contains
149  
150  
# Line 228 | Line 212 | contains
212            call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
213  
214            if (i_is_LJ .and. j_is_LJ) then
215 <             iHash = ior(iHash, LJ_PAIR)
232 <            
233 <
234 <
215 >             iHash = ior(iHash, LJ_PAIR)            
216            endif
217 +          
218 +          if (i_is_Elect .and. j_is_Elect) then
219 +             iHash = ior(iHash, ELECTROSTATIC_PAIR)
220 +          endif
221 +          
222 +          if (i_is_Sticky .and. j_is_Sticky) then
223 +             iHash = ior(iHash, STICKY_PAIR)
224 +          endif
225  
226 <
227 <
228 <          if (i_is_Elect .and. j_is_Elect) iHash = ior(iHash, ELECTROSTATIC_PAIR)
240 <          if (i_is_Sticky .and. j_is_Sticky) iHash = ior(iHash, STICKY_PAIR)
241 <          if (i_is_StickyP .and. j_is_StickyP) iHash = ior(iHash, STICKYPOWER_PAIR)
226 >          if (i_is_StickyP .and. j_is_StickyP) then
227 >             iHash = ior(iHash, STICKYPOWER_PAIR)
228 >          endif
229  
230 <          if (i_is_EAM .and. j_is_EAM) iHash = ior(iHash, EAM_PAIR)
230 >          if (i_is_EAM .and. j_is_EAM) then
231 >             iHash = ior(iHash, EAM_PAIR)
232 >          endif
233  
234            if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
235            if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
# Line 259 | Line 248 | contains
248      end do
249    end subroutine createInteractionMap
250  
251 + ! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff.
252 +  subroutine createRcuts(defaultRList)
253 +    real(kind=dp), intent(in), optional :: defaultRList
254 +    integer :: iMap
255 +    integer :: map_i,map_j
256 +    real(kind=dp) :: thisRCut = 0.0_dp
257 +    real(kind=dp) :: actualCutoff = 0.0_dp
258 +    integer :: nAtypes
259  
260 +    if(.not. allocated(InteractionMap)) return
261  
262 +    nAtypes = getSize(atypes)
263 + ! If we pass a default rcut, set all atypes to that cutoff distance
264 +    if(present(defaultRList)) then
265 +       InteractionMap(:,:)%rList = defaultRList
266 +       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
267 +       haveRlist = .true.
268 +       return
269 +    end if
270 +
271 +    do map_i = 1,nAtypes
272 +       do map_j = map_i,nAtypes
273 +          iMap = InteractionMap(map_i, map_j)%InteractionHash
274 +          
275 +          if ( iand(iMap, LJ_PAIR).ne.0 ) then
276 + !            thisRCut = getLJCutOff(map_i,map_j)
277 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
278 +          endif
279 +          
280 +          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
281 + !            thisRCut = getElectrostaticCutOff(map_i,map_j)
282 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
283 +          endif
284 +          
285 +          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
286 + !             thisRCut = getStickyCutOff(map_i,map_j)
287 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
288 +           endif
289 +          
290 +           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
291 + !              thisRCut = getStickyPowerCutOff(map_i,map_j)
292 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
293 +           endif
294 +          
295 +           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
296 + !              thisRCut = getGayberneCutOff(map_i,map_j)
297 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
298 +           endif
299 +          
300 +           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
301 + !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
302 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
303 +           endif
304 +          
305 +           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
306 + !              thisRCut = getEAMCutOff(map_i,map_j)
307 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
308 +           endif
309 +          
310 +           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
311 + !              thisRCut = getShapeCutOff(map_i,map_j)
312 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
313 +           endif
314 +          
315 +           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
316 + !              thisRCut = getShapeLJCutOff(map_i,map_j)
317 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
318 +           endif
319 +           InteractionMap(map_i, map_j)%rList = actualCutoff
320 +           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
321 +        end do
322 +     end do
323 +          haveRlist = .true.
324 +  end subroutine createRcuts
325 +
326 +
327   !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
328   !!$  subroutine setRlistDF( this_rlist )
329   !!$
# Line 272 | Line 335 | contains
335   !!$    haveRlist = .true.
336   !!$
337   !!$  end subroutine setRlistDF
275 –
276 –  subroutine createPropertyMap(status)
277 –    integer :: nAtypes
278 –    integer :: status
279 –    integer :: i
280 –    logical :: thisProperty
281 –    real (kind=DP) :: thisDPproperty
282 –
283 –    status = 0
284 –
285 –    nAtypes = getSize(atypes)
286 –
287 –    if (nAtypes == 0) then
288 –       status = -1
289 –       return
290 –    end if
291 –
292 –    if (.not. allocated(PropertyMap)) then
293 –       allocate(PropertyMap(nAtypes))
294 –    endif
295 –
296 –    do i = 1, nAtypes
297 –       call getElementProperty(atypes, i, "is_Directional", thisProperty)
298 –       PropertyMap(i)%is_Directional = thisProperty
299 –
300 –       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
301 –       PropertyMap(i)%is_LennardJones = thisProperty
302 –
303 –       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
304 –       PropertyMap(i)%is_Electrostatic = thisProperty
305 –
306 –       call getElementProperty(atypes, i, "is_Charge", thisProperty)
307 –       PropertyMap(i)%is_Charge = thisProperty
308 –
309 –       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
310 –       PropertyMap(i)%is_Dipole = thisProperty
311 –
312 –       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
313 –       PropertyMap(i)%is_Quadrupole = thisProperty
314 –
315 –       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
316 –       PropertyMap(i)%is_Sticky = thisProperty
317 –      
318 –       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
319 –       PropertyMap(i)%is_StickyPower = thisProperty
320 –
321 –       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
322 –       PropertyMap(i)%is_GayBerne = thisProperty
323 –
324 –       call getElementProperty(atypes, i, "is_EAM", thisProperty)
325 –       PropertyMap(i)%is_EAM = thisProperty
326 –
327 –       call getElementProperty(atypes, i, "is_Shape", thisProperty)
328 –       PropertyMap(i)%is_Shape = thisProperty
338  
330 –       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
331 –       PropertyMap(i)%is_FLARB = thisProperty
332 –    end do
339  
334 –    havePropertyMap = .true.
335 –
336 –  end subroutine createPropertyMap
337 –
340    subroutine setSimVariables()
341      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
342      SIM_uses_LennardJones = SimUsesLennardJones()
# Line 364 | Line 366 | contains
366  
367      error = 0
368  
369 <    if (.not. havePropertyMap) then
369 >    if (.not. haveInteractionMap) then
370  
371         myStatus = 0
372  
373 <       call createPropertyMap(myStatus)
373 >       call createInteractionMap(myStatus)
374  
375         if (myStatus .ne. 0) then
376 <          write(default_error, *) 'createPropertyMap failed in doForces!'
376 >          write(default_error, *) 'createInteractionMap failed in doForces!'
377            error = -1
378            return
379         endif
# Line 629 | Line 631 | contains
631      integer :: localError
632      integer :: propPack_i, propPack_j
633      integer :: loopStart, loopEnd, loop
634 <
634 >    integer :: iMap
635      real(kind=dp) :: listSkin = 1.0  
636  
637      !! initialize local variables  
# Line 755 | Line 757 | contains
757                    q_group(:,j), d_grp, rgrpsq)
758   #endif
759  
760 <             if (rgrpsq < rlistsq) then
760 >             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
761                  if (update_nlist) then
762                     nlist = nlist + 1
763  
# Line 973 | Line 975 | contains
975   #else
976               me_i = atid(i)
977   #endif
978 <
979 <             if (PropertyMap(me_i)%is_Dipole) then
978 >             iMap = InteractionMap(me_i, me_j)%InteractionHash
979 >            
980 >             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
981  
982                  mu_i = getDipoleMoment(me_i)
983  
# Line 1065 | Line 1068 | contains
1068         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1069              pot, eFrame, f, t, do_pot)
1070  
1071 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then                
1072 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1073 <               PropertyMap(me_j)%is_Dipole) then
1074 <             if (FF_uses_RF .and. SIM_uses_RF) then
1075 <                call accumulate_rf(i, j, r, eFrame, sw)
1073 <                call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1074 <             endif
1075 <          endif
1071 >       if (FF_uses_RF .and. SIM_uses_RF) then
1072 >
1073 >          ! CHECK ME (RF needs to know about all electrostatic types)
1074 >          call accumulate_rf(i, j, r, eFrame, sw)
1075 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1076         endif
1077 +
1078      endif
1079  
1080      if ( iand(iMap, STICKY_PAIR).ne.0 ) then
# Line 1092 | Line 1093 | contains
1093      endif
1094      
1095      if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1096 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1097 <            pot, A, f, t, do_pot)
1096 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1097 > !           pot, A, f, t, do_pot)
1098      endif
1099  
1100      if ( iand(iMap, EAM_PAIR).ne.0 ) then      

Diff Legend

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