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 2260 by chuckv, Mon Jun 27 22:21:37 2005 UTC vs.
Revision 2266 by chuckv, Thu Jul 28 22:12:45 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.24 2005-07-28 22:12:45 chuckv Exp $, $Date: 2005-07-28 22:12:45 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
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  
151    subroutine createInteractionMap(status)
152      integer :: nAtypes
153 <    integer :: status
153 >    integer, intent(out) :: status
154      integer :: i
155      integer :: j
156      integer :: ihash
# Line 187 | Line 171 | contains
171      logical :: j_is_EAM
172      logical :: j_is_Shape
173      
174 +    status = 0
175      
176      if (.not. associated(atypes)) then
177         call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
# Line 228 | Line 213 | contains
213            call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
214  
215            if (i_is_LJ .and. j_is_LJ) then
216 <             iHash = ior(iHash, LJ_PAIR)
217 <            
216 >             iHash = ior(iHash, LJ_PAIR)            
217 >          endif
218 >          
219 >          if (i_is_Elect .and. j_is_Elect) then
220 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
221 >          endif
222 >          
223 >          if (i_is_Sticky .and. j_is_Sticky) then
224 >             iHash = ior(iHash, STICKY_PAIR)
225 >          endif
226  
227 <
227 >          if (i_is_StickyP .and. j_is_StickyP) then
228 >             iHash = ior(iHash, STICKYPOWER_PAIR)
229            endif
230  
231 <
232 <
233 <          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)
242 <
243 <          if (i_is_EAM .and. j_is_EAM) iHash = ior(iHash, EAM_PAIR)
231 >          if (i_is_EAM .and. j_is_EAM) then
232 >             iHash = ior(iHash, EAM_PAIR)
233 >          endif
234  
235            if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
236            if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
# Line 259 | Line 249 | contains
249      end do
250    end subroutine createInteractionMap
251  
252 + ! 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.
253 +  subroutine createRcuts(defaultRList,stat)
254 +    real(kind=dp), intent(in), optional :: defaultRList
255 +    integer :: iMap
256 +    integer :: map_i,map_j
257 +    real(kind=dp) :: thisRCut = 0.0_dp
258 +    real(kind=dp) :: actualCutoff = 0.0_dp
259 +    integer, intent(out) :: stat
260 +    integer :: nAtypes
261 +    integer :: myStatus
262  
263 +    stat = 0
264 +    if (.not. haveInteractionMap) then
265  
266 +       call createInteractionMap(myStatus)
267 +
268 +       if (myStatus .ne. 0) then
269 +          write(default_error, *) 'createInteractionMap failed in doForces!'
270 +          stat = -1
271 +          return
272 +       endif
273 +    endif
274 +
275 +
276 +    nAtypes = getSize(atypes)
277 + ! If we pass a default rcut, set all atypes to that cutoff distance
278 +    if(present(defaultRList)) then
279 +       InteractionMap(:,:)%rList = defaultRList
280 +       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
281 +       haveRlist = .true.
282 +       return
283 +    end if
284 +
285 +    do map_i = 1,nAtypes
286 +       do map_j = map_i,nAtypes
287 +          iMap = InteractionMap(map_i, map_j)%InteractionHash
288 +          
289 +          if ( iand(iMap, LJ_PAIR).ne.0 ) then
290 + !            thisRCut = getLJCutOff(map_i,map_j)
291 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
292 +          endif
293 +          
294 +          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
295 + !            thisRCut = getElectrostaticCutOff(map_i,map_j)
296 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
297 +          endif
298 +          
299 +          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
300 + !             thisRCut = getStickyCutOff(map_i,map_j)
301 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
302 +           endif
303 +          
304 +           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
305 + !              thisRCut = getStickyPowerCutOff(map_i,map_j)
306 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
307 +           endif
308 +          
309 +           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
310 + !              thisRCut = getGayberneCutOff(map_i,map_j)
311 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
312 +           endif
313 +          
314 +           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
315 + !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
316 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
317 +           endif
318 +          
319 +           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
320 + !              thisRCut = getEAMCutOff(map_i,map_j)
321 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
322 +           endif
323 +          
324 +           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
325 + !              thisRCut = getShapeCutOff(map_i,map_j)
326 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
327 +           endif
328 +          
329 +           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
330 + !              thisRCut = getShapeLJCutOff(map_i,map_j)
331 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
332 +           endif
333 +           InteractionMap(map_i, map_j)%rList = actualCutoff
334 +           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
335 +        end do
336 +     end do
337 +          haveRlist = .true.
338 +  end subroutine createRcuts
339 +
340 +
341   !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
342   !!$  subroutine setRlistDF( this_rlist )
343   !!$
# Line 272 | Line 349 | contains
349   !!$    haveRlist = .true.
350   !!$
351   !!$  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
352  
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
353  
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
329 –
330 –       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
331 –       PropertyMap(i)%is_FLARB = thisProperty
332 –    end do
333 –
334 –    havePropertyMap = .true.
335 –
336 –  end subroutine createPropertyMap
337 –
354    subroutine setSimVariables()
355      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
356      SIM_uses_LennardJones = SimUsesLennardJones()
# Line 364 | Line 380 | contains
380  
381      error = 0
382  
383 <    if (.not. havePropertyMap) then
383 >    if (.not. haveInteractionMap) then
384  
385         myStatus = 0
386  
387 <       call createPropertyMap(myStatus)
387 >       call createInteractionMap(myStatus)
388  
389         if (myStatus .ne. 0) then
390 <          write(default_error, *) 'createPropertyMap failed in doForces!'
390 >          write(default_error, *) 'createInteractionMap failed in doForces!'
391            error = -1
392            return
393         endif
# Line 629 | Line 645 | contains
645      integer :: localError
646      integer :: propPack_i, propPack_j
647      integer :: loopStart, loopEnd, loop
648 <
648 >    integer :: iMap
649      real(kind=dp) :: listSkin = 1.0  
650  
651      !! initialize local variables  
# Line 727 | Line 743 | contains
743  
744            if (update_nlist) then
745   #ifdef IS_MPI
746 +             me_i = atid_row(i)
747               jstart = 1
748               jend = nGroupsInCol
749   #else
750 +             me_i = atid(i)
751               jstart = i+1
752               jend = nGroups
753   #endif
# Line 748 | Line 766 | contains
766               endif
767  
768   #ifdef IS_MPI
769 +             me_j = atid_col(j)
770               call get_interatomic_vector(q_group_Row(:,i), &
771                    q_group_Col(:,j), d_grp, rgrpsq)
772   #else
773 +             me_j = atid(j)
774               call get_interatomic_vector(q_group(:,i), &
775                    q_group(:,j), d_grp, rgrpsq)
776   #endif
777  
778 <             if (rgrpsq < rlistsq) then
778 >             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
779                  if (update_nlist) then
780                     nlist = nlist + 1
781  
# Line 973 | Line 993 | contains
993   #else
994               me_i = atid(i)
995   #endif
996 <
997 <             if (PropertyMap(me_i)%is_Dipole) then
996 >             iMap = InteractionMap(me_i, me_j)%InteractionHash
997 >            
998 >             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
999  
1000                  mu_i = getDipoleMoment(me_i)
1001  
# Line 1065 | Line 1086 | contains
1086         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1087              pot, eFrame, f, t, do_pot)
1088  
1089 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then                
1090 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1091 <               PropertyMap(me_j)%is_Dipole) then
1092 <             if (FF_uses_RF .and. SIM_uses_RF) then
1093 <                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
1089 >       if (FF_uses_RF .and. SIM_uses_RF) then
1090 >
1091 >          ! CHECK ME (RF needs to know about all electrostatic types)
1092 >          call accumulate_rf(i, j, r, eFrame, sw)
1093 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1094         endif
1095 +
1096      endif
1097  
1098      if ( iand(iMap, STICKY_PAIR).ne.0 ) then
# Line 1092 | Line 1111 | contains
1111      endif
1112      
1113      if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1114 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1115 <            pot, A, f, t, do_pot)
1114 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1115 > !           pot, A, f, t, do_pot)
1116      endif
1117  
1118      if ( iand(iMap, EAM_PAIR).ne.0 ) then      

Diff Legend

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