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 2267 by tim, Fri Jul 29 17:34:06 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.25 2005-07-29 17:34:06 tim Exp $, $Date: 2005-07-29 17:34:06 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
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
157      real(kind=dp) :: myRcut
158 < ! Test Types
158 >    !! Test Types
159      logical :: i_is_LJ
160      logical :: i_is_Elect
161      logical :: i_is_Sticky
# 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 257 | Line 247 | contains
247         end do
248  
249      end do
250 +
251 +    haveInteractionMap = .true.
252    end subroutine createInteractionMap
253  
254 + ! 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.
255 +  subroutine createRcuts(defaultRList,stat)
256 +    real(kind=dp), intent(in), optional :: defaultRList
257 +    integer :: iMap
258 +    integer :: map_i,map_j
259 +    real(kind=dp) :: thisRCut = 0.0_dp
260 +    real(kind=dp) :: actualCutoff = 0.0_dp
261 +    integer, intent(out) :: stat
262 +    integer :: nAtypes
263 +    integer :: myStatus
264  
265 +    stat = 0
266 +    if (.not. haveInteractionMap) then
267  
268 +       call createInteractionMap(myStatus)
269 +
270 +       if (myStatus .ne. 0) then
271 +          write(default_error, *) 'createInteractionMap failed in doForces!'
272 +          stat = -1
273 +          return
274 +       endif
275 +    endif
276 +
277 +
278 +    nAtypes = getSize(atypes)
279 +    !! If we pass a default rcut, set all atypes to that cutoff distance
280 +    if(present(defaultRList)) then
281 +       InteractionMap(:,:)%rList = defaultRList
282 +       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
283 +       haveRlist = .true.
284 +       return
285 +    end if
286 +
287 +    do map_i = 1,nAtypes
288 +       do map_j = map_i,nAtypes
289 +          iMap = InteractionMap(map_i, map_j)%InteractionHash
290 +          
291 +          if ( iand(iMap, LJ_PAIR).ne.0 ) then
292 +             ! thisRCut = getLJCutOff(map_i,map_j)
293 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
294 +          endif
295 +          
296 +          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
297 +             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
298 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
299 +          endif
300 +          
301 +          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
302 +             ! thisRCut = getStickyCutOff(map_i,map_j)
303 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
304 +           endif
305 +          
306 +           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
307 +              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
308 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
309 +           endif
310 +          
311 +           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
312 +              ! thisRCut = getGayberneCutOff(map_i,map_j)
313 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
314 +           endif
315 +          
316 +           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
317 + !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
318 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
319 +           endif
320 +          
321 +           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
322 + !              thisRCut = getEAMCutOff(map_i,map_j)
323 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
324 +           endif
325 +          
326 +           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
327 + !              thisRCut = getShapeCutOff(map_i,map_j)
328 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
329 +           endif
330 +          
331 +           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
332 + !              thisRCut = getShapeLJCutOff(map_i,map_j)
333 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
334 +           endif
335 +           InteractionMap(map_i, map_j)%rList = actualCutoff
336 +           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
337 +        end do
338 +     end do
339 +     haveRlist = .true.
340 +  end subroutine createRcuts
341 +
342 +
343   !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
344   !!$  subroutine setRlistDF( this_rlist )
345   !!$
# Line 272 | Line 351 | contains
351   !!$    haveRlist = .true.
352   !!$
353   !!$  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
329 –
330 –       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
331 –       PropertyMap(i)%is_FLARB = thisProperty
332 –    end do
354  
334 –    havePropertyMap = .true.
355  
336 –  end subroutine createPropertyMap
337 –
356    subroutine setSimVariables()
357      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
358      SIM_uses_LennardJones = SimUsesLennardJones()
# Line 364 | Line 382 | contains
382  
383      error = 0
384  
385 <    if (.not. havePropertyMap) then
385 >    if (.not. haveInteractionMap) then
386  
387         myStatus = 0
388  
389 <       call createPropertyMap(myStatus)
389 >       call createInteractionMap(myStatus)
390  
391         if (myStatus .ne. 0) then
392 <          write(default_error, *) 'createPropertyMap failed in doForces!'
392 >          write(default_error, *) 'createInteractionMap failed in doForces!'
393            error = -1
394            return
395         endif
# Line 629 | Line 647 | contains
647      integer :: localError
648      integer :: propPack_i, propPack_j
649      integer :: loopStart, loopEnd, loop
650 <
650 >    integer :: iMap
651      real(kind=dp) :: listSkin = 1.0  
652  
653      !! initialize local variables  
# Line 720 | Line 738 | contains
738         iend = nGroups - 1
739   #endif
740         outer: do i = istart, iend
741 +
742 + #ifdef IS_MPI
743 +             me_i = atid_row(i)
744 + #else
745 +             me_i = atid(i)
746 + #endif
747  
748            if (update_nlist) point(i) = nlist + 1
749  
# Line 748 | Line 772 | contains
772               endif
773  
774   #ifdef IS_MPI
775 +             me_j = atid_col(j)
776               call get_interatomic_vector(q_group_Row(:,i), &
777                    q_group_Col(:,j), d_grp, rgrpsq)
778   #else
779 +             me_j = atid(j)
780               call get_interatomic_vector(q_group(:,i), &
781                    q_group(:,j), d_grp, rgrpsq)
782   #endif
783  
784 <             if (rgrpsq < rlistsq) then
784 >             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
785                  if (update_nlist) then
786                     nlist = nlist + 1
787  
# Line 973 | Line 999 | contains
999   #else
1000               me_i = atid(i)
1001   #endif
1002 <
1003 <             if (PropertyMap(me_i)%is_Dipole) then
1002 >             iMap = InteractionMap(me_i, me_j)%InteractionHash
1003 >            
1004 >             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1005  
1006                  mu_i = getDipoleMoment(me_i)
1007  
# Line 1065 | Line 1092 | contains
1092         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1093              pot, eFrame, f, t, do_pot)
1094  
1095 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then                
1096 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1097 <               PropertyMap(me_j)%is_Dipole) then
1098 <             if (FF_uses_RF .and. SIM_uses_RF) then
1099 <                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
1095 >       if (FF_uses_RF .and. SIM_uses_RF) then
1096 >
1097 >          ! CHECK ME (RF needs to know about all electrostatic types)
1098 >          call accumulate_rf(i, j, r, eFrame, sw)
1099 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1100         endif
1101 +
1102      endif
1103  
1104      if ( iand(iMap, STICKY_PAIR).ne.0 ) then
# Line 1092 | Line 1117 | contains
1117      endif
1118      
1119      if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1120 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1121 <            pot, A, f, t, do_pot)
1120 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1121 > !           pot, A, f, t, do_pot)
1122      endif
1123  
1124      if ( iand(iMap, EAM_PAIR).ne.0 ) then      

Diff Legend

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