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 2261 by gezelter, Tue Jun 28 13:58:45 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.22 2005-06-28 13:58:45 gezelter Exp $, $Date: 2005-06-28 13:58:45 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $
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
62 <  use gb_pair
61 >  use gayberne
62    use shapes
63    use vector_class
64    use eam
# Line 73 | Line 72 | module doForces
72  
73   #define __FORTRAN90
74   #include "UseTheForce/fSwitchingFunction.h"
75 + #include "UseTheForce/fCutoffPolicy.h"
76   #include "UseTheForce/DarkSide/fInteractionMap.h"
77 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78  
79 +
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82  
81 –  logical, save :: haveRlist = .false.
83    logical, save :: haveNeighborList = .false.
84    logical, save :: haveSIMvariables = .false.
85    logical, save :: haveSaneForceField = .false.
86 +  logical, save :: haveInteractionHash = .false.
87 +  logical, save :: haveGtypeCutoffMap = .false.
88 +  logical, save :: haveDefaultCutoffs = .false.
89 +  logical, save :: haveRlist = .false.
90  
91    logical, save :: FF_uses_DirectionalAtoms
87 –  logical, save :: FF_uses_LennardJones
88 –  logical, save :: FF_uses_Electrostatics
89 –  logical, save :: FF_uses_Charges
92    logical, save :: FF_uses_Dipoles
91 –  logical, save :: FF_uses_Quadrupoles
92 –  logical, save :: FF_uses_Sticky
93 –  logical, save :: FF_uses_StickyPower
93    logical, save :: FF_uses_GayBerne
94    logical, save :: FF_uses_EAM
96 –  logical, save :: FF_uses_Shapes
97 –  logical, save :: FF_uses_FLARB
98 –  logical, save :: FF_uses_RF
95  
96    logical, save :: SIM_uses_DirectionalAtoms
101 –  logical, save :: SIM_uses_LennardJones
102 –  logical, save :: SIM_uses_Electrostatics
103 –  logical, save :: SIM_uses_Charges
104 –  logical, save :: SIM_uses_Dipoles
105 –  logical, save :: SIM_uses_Quadrupoles
106 –  logical, save :: SIM_uses_Sticky
107 –  logical, save :: SIM_uses_StickyPower
108 –  logical, save :: SIM_uses_GayBerne
97    logical, save :: SIM_uses_EAM
110 –  logical, save :: SIM_uses_Shapes
111 –  logical, save :: SIM_uses_FLARB
112 –  logical, save :: SIM_uses_RF
98    logical, save :: SIM_requires_postpair_calc
99    logical, save :: SIM_requires_prepair_calc
100    logical, save :: SIM_uses_PBC
116 –  logical, save :: SIM_uses_molecular_cutoffs
101  
102 <  !!!GO AWAY---------
119 <  !!!!!real(kind=dp), save :: rlist, rlistsq
102 >  integer, save :: electrostaticSummationMethod
103  
104    public :: init_FF
105 +  public :: setDefaultCutoffs
106    public :: do_force_loop
107 < !  public :: setRlistDF
107 >  public :: createInteractionHash
108 >  public :: createGtypeCutoffMap
109 >  public :: getStickyCut
110 >  public :: getStickyPowerCut
111 >  public :: getGayBerneCut
112 >  public :: getEAMCut
113 >  public :: getShapeCut
114  
115   #ifdef PROFILE
116    public :: getforcetime
# Line 128 | Line 118 | module doForces
118    real :: forceTimeInitial, forceTimeFinal
119    integer :: nLoops
120   #endif
121 +  
122 +  !! Variables for cutoff mapping and interaction mapping
123 +  ! Bit hash to determine pair-pair interactions.
124 +  integer, dimension(:,:), allocatable :: InteractionHash
125 +  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
126 +  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
127 +  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
128  
129 <  type, public :: Interaction
130 <     integer :: InteractionHash
131 <     real(kind=dp) :: rCut
132 <  end type Interaction
129 >  integer, dimension(:), allocatable, target :: groupToGtypeRow
130 >  integer, dimension(:), pointer :: groupToGtypeCol => null()
131 >
132 >  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
133 >  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
134 >  type ::gtypeCutoffs
135 >     real(kind=dp) :: rcut
136 >     real(kind=dp) :: rcutsq
137 >     real(kind=dp) :: rlistsq
138 >  end type gtypeCutoffs
139 >  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
140 >
141 >  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
142 >  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
143 >  real(kind=dp),save :: listSkin
144    
137 –  type(Interaction), dimension(:,:),allocatable :: InteractionMap
138 –  
139 –  !public :: addInteraction
140 –  !public :: setInteractionHash
141 –  !public :: getInteractionHash
142 –  public :: createInteractionMap
143 –  
145   contains
146  
147 <
147 <  subroutine createInteractionMap(status)
147 >  subroutine createInteractionHash(status)
148      integer :: nAtypes
149 <    integer :: status
149 >    integer, intent(out) :: status
150      integer :: i
151      integer :: j
152 <    integer :: ihash
153 <    real(kind=dp) :: myRcut
154 < ! Test Types
152 >    integer :: iHash
153 >    !! Test Types
154      logical :: i_is_LJ
155      logical :: i_is_Elect
156      logical :: i_is_Sticky
# Line 166 | Line 165 | contains
165      logical :: j_is_GB
166      logical :: j_is_EAM
167      logical :: j_is_Shape
168 <    
169 <    
168 >    real(kind=dp) :: myRcut
169 >
170 >    status = 0  
171 >
172      if (.not. associated(atypes)) then
173 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
173 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
174         status = -1
175         return
176      endif
# Line 181 | Line 182 | contains
182         return
183      end if
184  
185 <    if (.not. allocated(InteractionMap)) then
186 <       allocate(InteractionMap(nAtypes,nAtypes))
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 236 | Line 247 | contains
247            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
248  
249  
250 <          InteractionMap(i,j)%InteractionHash = iHash
251 <          InteractionMap(j,i)%InteractionHash = iHash
250 >          InteractionHash(i,j) = iHash
251 >          InteractionHash(j,i) = iHash
252  
253         end do
254  
255      end do
245 –  end subroutine createInteractionMap
256  
257 +    haveInteractionHash = .true.
258 +  end subroutine createInteractionHash
259  
260 +  subroutine createGtypeCutoffMap(stat)
261  
262 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
263 < !!$  subroutine setRlistDF( this_rlist )
264 < !!$
265 < !!$   real(kind=dp) :: this_rlist
266 < !!$
267 < !!$    rlist = this_rlist
268 < !!$    rlistsq = rlist * rlist
269 < !!$
270 < !!$    haveRlist = .true.
271 < !!$
272 < !!$  end subroutine setRlistDF
262 >    integer, intent(out), optional :: stat
263 >    logical :: i_is_LJ
264 >    logical :: i_is_Elect
265 >    logical :: i_is_Sticky
266 >    logical :: i_is_StickyP
267 >    logical :: i_is_GB
268 >    logical :: i_is_EAM
269 >    logical :: i_is_Shape
270 >    logical :: GtypeFound
271 >
272 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
273 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
274 >    integer :: nGroupsInRow
275 >    integer :: nGroupsInCol
276 >    integer :: nGroupTypesRow,nGroupTypesCol
277 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
278 >    real(kind=dp) :: biggestAtypeCutoff
279 >
280 >    stat = 0
281 >    if (.not. haveInteractionHash) then
282 >       call createInteractionHash(myStatus)      
283 >       if (myStatus .ne. 0) then
284 >          write(default_error, *) 'createInteractionHash failed in doForces!'
285 >          stat = -1
286 >          return
287 >       endif
288 >    endif
289 > #ifdef IS_MPI
290 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
291 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
292 > #endif
293 >    nAtypes = getSize(atypes)
294 > ! Set all of the initial cutoffs to zero.
295 >    atypeMaxCutoff = 0.0_dp
296 >    do i = 1, nAtypes
297 >       if (SimHasAtype(i)) then    
298 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
299 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
300 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
301 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
302 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
303 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
304 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
305 >          
306 >
307 >          if (haveDefaultCutoffs) then
308 >             atypeMaxCutoff(i) = defaultRcut
309 >          else
310 >             if (i_is_LJ) then          
311 >                thisRcut = getSigma(i) * 2.5_dp
312 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313 >             endif
314 >             if (i_is_Elect) then
315 >                thisRcut = defaultRcut
316 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317 >             endif
318 >             if (i_is_Sticky) then
319 >                thisRcut = getStickyCut(i)
320 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321 >             endif
322 >             if (i_is_StickyP) then
323 >                thisRcut = getStickyPowerCut(i)
324 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325 >             endif
326 >             if (i_is_GB) then
327 >                thisRcut = getGayBerneCut(i)
328 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329 >             endif
330 >             if (i_is_EAM) then
331 >                thisRcut = getEAMCut(i)
332 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333 >             endif
334 >             if (i_is_Shape) then
335 >                thisRcut = getShapeCut(i)
336 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337 >             endif
338 >          endif
339 >          
340 >          
341 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
342 >             biggestAtypeCutoff = atypeMaxCutoff(i)
343 >          endif
344  
345 +       endif
346 +    enddo
347 +  
348  
349 +    
350 +    istart = 1
351 +    jstart = 1
352 + #ifdef IS_MPI
353 +    iend = nGroupsInRow
354 +    jend = nGroupsInCol
355 + #else
356 +    iend = nGroups
357 +    jend = nGroups
358 + #endif
359 +    
360 +    !! allocate the groupToGtype and gtypeMaxCutoff here.
361 +    if(.not.allocated(groupToGtypeRow)) then
362 +     !  allocate(groupToGtype(iend))
363 +       allocate(groupToGtypeRow(iend))
364 +    else
365 +       deallocate(groupToGtypeRow)
366 +       allocate(groupToGtypeRow(iend))
367 +    endif
368 +    if(.not.allocated(groupMaxCutoffRow)) then
369 +       allocate(groupMaxCutoffRow(iend))
370 +    else
371 +       deallocate(groupMaxCutoffRow)
372 +       allocate(groupMaxCutoffRow(iend))
373 +    end if
374 +
375 +    if(.not.allocated(gtypeMaxCutoffRow)) then
376 +       allocate(gtypeMaxCutoffRow(iend))
377 +    else
378 +       deallocate(gtypeMaxCutoffRow)
379 +       allocate(gtypeMaxCutoffRow(iend))
380 +    endif
381 +
382 +
383 + #ifdef IS_MPI
384 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
385 +    if(.not.associated(groupToGtypeCol)) then
386 +       allocate(groupToGtypeCol(jend))
387 +    else
388 +       deallocate(groupToGtypeCol)
389 +       allocate(groupToGtypeCol(jend))
390 +    end if
391 +
392 +    if(.not.associated(groupToGtypeCol)) then
393 +       allocate(groupToGtypeCol(jend))
394 +    else
395 +       deallocate(groupToGtypeCol)
396 +       allocate(groupToGtypeCol(jend))
397 +    end if
398 +    if(.not.associated(gtypeMaxCutoffCol)) then
399 +       allocate(gtypeMaxCutoffCol(jend))
400 +    else
401 +       deallocate(gtypeMaxCutoffCol)      
402 +       allocate(gtypeMaxCutoffCol(jend))
403 +    end if
404 +
405 +       groupMaxCutoffCol = 0.0_dp
406 +       gtypeMaxCutoffCol = 0.0_dp
407 +
408 + #endif
409 +       groupMaxCutoffRow = 0.0_dp
410 +       gtypeMaxCutoffRow = 0.0_dp
411 +
412 +
413 +    !! first we do a single loop over the cutoff groups to find the
414 +    !! largest cutoff for any atypes present in this group.  We also
415 +    !! create gtypes at this point.
416 +    
417 +    tol = 1.0d-6
418 +    nGroupTypesRow = 0
419 +
420 +    do i = istart, iend      
421 +       n_in_i = groupStartRow(i+1) - groupStartRow(i)
422 +       groupMaxCutoffRow(i) = 0.0_dp
423 +       do ia = groupStartRow(i), groupStartRow(i+1)-1
424 +          atom1 = groupListRow(ia)
425 + #ifdef IS_MPI
426 +          me_i = atid_row(atom1)
427 + #else
428 +          me_i = atid(atom1)
429 + #endif          
430 +          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
431 +             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
432 +          endif          
433 +       enddo
434 +
435 +       if (nGroupTypesRow.eq.0) then
436 +          nGroupTypesRow = nGroupTypesRow + 1
437 +          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
438 +          groupToGtypeRow(i) = nGroupTypesRow
439 +       else
440 +          GtypeFound = .false.
441 +          do g = 1, nGroupTypesRow
442 +             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
443 +                groupToGtypeRow(i) = g
444 +                GtypeFound = .true.
445 +             endif
446 +          enddo
447 +          if (.not.GtypeFound) then            
448 +             nGroupTypesRow = nGroupTypesRow + 1
449 +             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
450 +             groupToGtypeRow(i) = nGroupTypesRow
451 +          endif
452 +       endif
453 +    enddo    
454 +
455 + #ifdef IS_MPI
456 +    do j = jstart, jend      
457 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
458 +       groupMaxCutoffCol(j) = 0.0_dp
459 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
460 +          atom1 = groupListCol(ja)
461 +
462 +          me_j = atid_col(atom1)
463 +
464 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
465 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
466 +          endif          
467 +       enddo
468 +
469 +       if (nGroupTypesCol.eq.0) then
470 +          nGroupTypesCol = nGroupTypesCol + 1
471 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
472 +          groupToGtypeCol(j) = nGroupTypesCol
473 +       else
474 +          GtypeFound = .false.
475 +          do g = 1, nGroupTypesCol
476 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
477 +                groupToGtypeCol(j) = g
478 +                GtypeFound = .true.
479 +             endif
480 +          enddo
481 +          if (.not.GtypeFound) then            
482 +             nGroupTypesCol = nGroupTypesCol + 1
483 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
484 +             groupToGtypeCol(j) = nGroupTypesCol
485 +          endif
486 +       endif
487 +    enddo    
488 +
489 + #else
490 + ! Set pointers to information we just found
491 +    nGroupTypesCol = nGroupTypesRow
492 +    groupToGtypeCol => groupToGtypeRow
493 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
494 +    groupMaxCutoffCol => groupMaxCutoffRow
495 + #endif
496 +
497 +
498 +
499 +
500 +
501 +    !! allocate the gtypeCutoffMap here.
502 +    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
503 +    !! then we do a double loop over all the group TYPES to find the cutoff
504 +    !! map between groups of two types
505 +    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
506 +
507 +    do i = 1, nGroupTypesRow
508 +       do j = 1, nGroupTypesCol
509 +      
510 +          select case(cutoffPolicy)
511 +          case(TRADITIONAL_CUTOFF_POLICY)
512 +             thisRcut = tradRcut
513 +          case(MIX_CUTOFF_POLICY)
514 +             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
515 +          case(MAX_CUTOFF_POLICY)
516 +             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
517 +          case default
518 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
519 +             return
520 +          end select
521 +          gtypeCutoffMap(i,j)%rcut = thisRcut
522 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
523 +          skin = defaultRlist - defaultRcut
524 +          listSkin = skin ! set neighbor list skin thickness
525 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
526 +
527 +          ! sanity check
528 +
529 +          if (haveDefaultCutoffs) then
530 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
531 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
532 +             endif
533 +          endif
534 +       enddo
535 +    enddo
536 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
537 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
538 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
539 + #ifdef IS_MPI
540 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
541 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
542 + #endif
543 +    groupMaxCutoffCol => null()
544 +    gtypeMaxCutoffCol => null()
545 +    
546 +    haveGtypeCutoffMap = .true.
547 +   end subroutine createGtypeCutoffMap
548 +
549 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
550 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
551 +     integer, intent(in) :: cutPolicy
552 +
553 +     defaultRcut = defRcut
554 +     defaultRsw = defRsw
555 +     defaultRlist = defRlist
556 +     cutoffPolicy = cutPolicy
557 +
558 +     haveDefaultCutoffs = .true.
559 +   end subroutine setDefaultCutoffs
560 +
561 +   subroutine setCutoffPolicy(cutPolicy)
562 +
563 +     integer, intent(in) :: cutPolicy
564 +     cutoffPolicy = cutPolicy
565 +     call createGtypeCutoffMap()
566 +   end subroutine setCutoffPolicy
567 +    
568 +    
569    subroutine setSimVariables()
570      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
264 –    SIM_uses_LennardJones = SimUsesLennardJones()
265 –    SIM_uses_Electrostatics = SimUsesElectrostatics()
266 –    SIM_uses_Charges = SimUsesCharges()
267 –    SIM_uses_Dipoles = SimUsesDipoles()
268 –    SIM_uses_Sticky = SimUsesSticky()
269 –    SIM_uses_StickyPower = SimUsesStickyPower()
270 –    SIM_uses_GayBerne = SimUsesGayBerne()
571      SIM_uses_EAM = SimUsesEAM()
272 –    SIM_uses_Shapes = SimUsesShapes()
273 –    SIM_uses_FLARB = SimUsesFLARB()
274 –    SIM_uses_RF = SimUsesRF()
572      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574      SIM_uses_PBC = SimUsesPBC()
# Line 288 | Line 585 | contains
585  
586      error = 0
587  
588 <    if (.not. haveInteractionMap) then
588 >    if (.not. haveInteractionHash) then      
589 >       myStatus = 0      
590 >       call createInteractionHash(myStatus)      
591 >       if (myStatus .ne. 0) then
592 >          write(default_error, *) 'createInteractionHash failed in doForces!'
593 >          error = -1
594 >          return
595 >       endif
596 >    endif
597  
598 <       myStatus = 0
599 <
600 <       call createInteractionMap(myStatus)
296 <
598 >    if (.not. haveGtypeCutoffMap) then        
599 >       myStatus = 0      
600 >       call createGtypeCutoffMap(myStatus)      
601         if (myStatus .ne. 0) then
602 <          write(default_error, *) 'createInteractionMap failed in doForces!'
602 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
603            error = -1
604            return
605         endif
# Line 305 | Line 609 | contains
609         call setSimVariables()
610      endif
611  
612 <    if (.not. haveRlist) then
613 <       write(default_error, *) 'rList has not been set in doForces!'
614 <       error = -1
615 <       return
616 <    endif
612 >  !  if (.not. haveRlist) then
613 >  !     write(default_error, *) 'rList has not been set in doForces!'
614 >  !     error = -1
615 >  !     return
616 >  !  endif
617  
618      if (.not. haveNeighborList) then
619         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 334 | Line 638 | contains
638    end subroutine doReadyCheck
639  
640  
641 <  subroutine init_FF(use_RF_c, thisStat)
641 >  subroutine init_FF(thisESM, thisStat)
642  
643 <    logical, intent(in) :: use_RF_c
340 <
643 >    integer, intent(in) :: thisESM
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
646      integer, pointer :: MatchList(:) => null()
# Line 346 | Line 649 | contains
649      !! assume things are copacetic, unless they aren't
650      thisStat = 0
651  
652 <    !! Fortran's version of a cast:
350 <    FF_uses_RF = use_RF_c
652 >    electrostaticSummationMethod = thisESM
653  
654      !! init_FF is called *after* all of the atom types have been
655      !! defined in atype_module using the new_atype subroutine.
# Line 356 | Line 658 | contains
658      !! interactions are used by the force field.    
659  
660      FF_uses_DirectionalAtoms = .false.
359 –    FF_uses_LennardJones = .false.
360 –    FF_uses_Electrostatics = .false.
361 –    FF_uses_Charges = .false.    
661      FF_uses_Dipoles = .false.
363 –    FF_uses_Sticky = .false.
364 –    FF_uses_StickyPower = .false.
662      FF_uses_GayBerne = .false.
663      FF_uses_EAM = .false.
367 –    FF_uses_Shapes = .false.
368 –    FF_uses_FLARB = .false.
664  
665      call getMatchingElementList(atypes, "is_Directional", .true., &
666           nMatches, MatchList)
667      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
668  
374 –    call getMatchingElementList(atypes, "is_LennardJones", .true., &
375 –         nMatches, MatchList)
376 –    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
377 –
378 –    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
379 –         nMatches, MatchList)
380 –    if (nMatches .gt. 0) then
381 –       FF_uses_Electrostatics = .true.
382 –    endif
383 –
384 –    call getMatchingElementList(atypes, "is_Charge", .true., &
385 –         nMatches, MatchList)
386 –    if (nMatches .gt. 0) then
387 –       FF_uses_Charges = .true.  
388 –       FF_uses_Electrostatics = .true.
389 –    endif
390 –
669      call getMatchingElementList(atypes, "is_Dipole", .true., &
670           nMatches, MatchList)
671 <    if (nMatches .gt. 0) then
394 <       FF_uses_Dipoles = .true.
395 <       FF_uses_Electrostatics = .true.
396 <       FF_uses_DirectionalAtoms = .true.
397 <    endif
398 <
399 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
400 <         nMatches, MatchList)
401 <    if (nMatches .gt. 0) then
402 <       FF_uses_Quadrupoles = .true.
403 <       FF_uses_Electrostatics = .true.
404 <       FF_uses_DirectionalAtoms = .true.
405 <    endif
406 <
407 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
408 <         MatchList)
409 <    if (nMatches .gt. 0) then
410 <       FF_uses_Sticky = .true.
411 <       FF_uses_DirectionalAtoms = .true.
412 <    endif
413 <
414 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
415 <         MatchList)
416 <    if (nMatches .gt. 0) then
417 <       FF_uses_StickyPower = .true.
418 <       FF_uses_DirectionalAtoms = .true.
419 <    endif
671 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
672      
673      call getMatchingElementList(atypes, "is_GayBerne", .true., &
674           nMatches, MatchList)
675 <    if (nMatches .gt. 0) then
424 <       FF_uses_GayBerne = .true.
425 <       FF_uses_DirectionalAtoms = .true.
426 <    endif
675 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
676  
677      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
678      if (nMatches .gt. 0) FF_uses_EAM = .true.
679  
431 –    call getMatchingElementList(atypes, "is_Shape", .true., &
432 –         nMatches, MatchList)
433 –    if (nMatches .gt. 0) then
434 –       FF_uses_Shapes = .true.
435 –       FF_uses_DirectionalAtoms = .true.
436 –    endif
680  
438 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
439 –         nMatches, MatchList)
440 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
441 –
442 –    !! Assume sanity (for the sake of argument)
681      haveSaneForceField = .true.
682  
445 –    !! check to make sure the FF_uses_RF setting makes sense
446 –
447 –    if (FF_uses_dipoles) then
448 –       if (FF_uses_RF) then
449 –          dielect = getDielect()
450 –          call initialize_rf(dielect)
451 –       endif
452 –    else
453 –       if (FF_uses_RF) then          
454 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
455 –          thisStat = -1
456 –          haveSaneForceField = .false.
457 –          return
458 –       endif
459 –    endif
460 –
461 –    !sticky module does not contain check_sticky_FF anymore
462 –    !if (FF_uses_sticky) then
463 –    !   call check_sticky_FF(my_status)
464 –    !   if (my_status /= 0) then
465 –    !      thisStat = -1
466 –    !      haveSaneForceField = .false.
467 –    !      return
468 –    !   end if
469 –    !endif
470 –
683      if (FF_uses_EAM) then
684         call init_EAM_FF(my_status)
685         if (my_status /= 0) then
# Line 476 | Line 688 | contains
688            haveSaneForceField = .false.
689            return
690         end if
479 –    endif
480 –
481 –    if (FF_uses_GayBerne) then
482 –       call check_gb_pair_FF(my_status)
483 –       if (my_status .ne. 0) then
484 –          thisStat = -1
485 –          haveSaneForceField = .false.
486 –          return
487 –       endif
488 –    endif
489 –
490 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
691      endif
692  
693      if (.not. haveNeighborList) then
# Line 523 | Line 723 | contains
723  
724      !! Stress Tensor
725      real( kind = dp), dimension(9) :: tau  
726 <    real ( kind = dp ) :: pot
726 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
727      logical ( kind = 2) :: do_pot_c, do_stress_c
728      logical :: do_pot
729      logical :: do_stress
730      logical :: in_switching_region
731   #ifdef IS_MPI
732 <    real( kind = DP ) :: pot_local
732 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
733      integer :: nAtomsInRow
734      integer :: nAtomsInCol
735      integer :: nprocs
# Line 553 | Line 753 | contains
753      integer :: localError
754      integer :: propPack_i, propPack_j
755      integer :: loopStart, loopEnd, loop
756 +    integer :: iHash
757 +  
758  
557 –    real(kind=dp) :: listSkin = 1.0  
558 –
759      !! initialize local variables  
760  
761   #ifdef IS_MPI
# Line 672 | Line 872 | contains
872               endif
873  
874   #ifdef IS_MPI
875 +             me_j = atid_col(j)
876               call get_interatomic_vector(q_group_Row(:,i), &
877                    q_group_Col(:,j), d_grp, rgrpsq)
878   #else
879 +             me_j = atid(j)
880               call get_interatomic_vector(q_group(:,i), &
881                    q_group(:,j), d_grp, rgrpsq)
882 < #endif
882 > #endif      
883  
884 <             if (rgrpsq < rlistsq) then
884 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
885                  if (update_nlist) then
886                     nlist = nlist + 1
887  
# Line 800 | Line 1002 | contains
1002                  endif
1003               end if
1004            enddo
1005 +
1006         enddo outer
1007  
1008         if (update_nlist) then
# Line 859 | Line 1062 | contains
1062  
1063      if (do_pot) then
1064         ! scatter/gather pot_row into the members of my column
1065 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1066 <
1065 >       do i = 1,LR_POT_TYPES
1066 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1067 >       end do
1068         ! scatter/gather pot_local into all other procs
1069         ! add resultant to get total pot
1070         do i = 1, nlocal
1071 <          pot_local = pot_local + pot_Temp(i)
1071 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1072 >               + pot_Temp(1:LR_POT_TYPES,i)
1073         enddo
1074  
1075         pot_Temp = 0.0_DP
1076 <
1077 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1076 >       do i = 1,LR_POT_TYPES
1077 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1078 >       end do
1079         do i = 1, nlocal
1080 <          pot_local = pot_local + pot_Temp(i)
1080 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1081 >               + pot_Temp(1:LR_POT_TYPES,i)
1082         enddo
1083  
1084      endif
1085   #endif
1086  
1087 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1088 <
882 <       if (FF_uses_RF .and. SIM_uses_RF) then
883 <
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)
1101 >          me_i = atid_row(i)
1102   #else
1103 <             me_i = atid(i)
1103 >          me_i = atid(i)
1104   #endif
1105 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
1105 >          iHash = InteractionHash(me_i,me_j)
1106 >          
1107 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1108              
1109 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1110 <
1111 <                mu_i = getDipoleMoment(me_i)
1112 <
1113 <                !! The reaction field needs to include a self contribution
1114 <                !! to the field:
1115 <                call accumulate_self_rf(i, mu_i, eFrame)
1116 <                !! Get the reaction field contribution to the
910 <                !! potential and torques:
911 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
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 = pot_local + rfpot
1118 >             pot_local(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
1119   #else
1120 <                pot = pot + rfpot
1121 <
1120 >             pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1121 >            
1122   #endif
1123 <             endif
1124 <          enddo
920 <       endif
1123 >          endif
1124 >       enddo
1125      endif
1126 <
923 <
1126 >    
1127   #ifdef IS_MPI
1128  
1129      if (do_pot) then
1130 <       pot = pot + pot_local
1130 >       pot(1:LR_POT_TYPES) = pot(1:LR_POT_TYPES) &
1131 >            + pot_local(1:LR_POT_TYPES)
1132         !! we assume the c code will do the allreduce to get the total potential
1133         !! we could do it right here if we needed to...
1134      endif
# Line 950 | Line 1154 | contains
1154    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1155         eFrame, A, f, t, pot, vpair, fpair)
1156  
1157 <    real( kind = dp ) :: pot, vpair, sw
1157 >    real( kind = dp ) :: vpair, sw
1158 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1159      real( kind = dp ), dimension(3) :: fpair
1160      real( kind = dp ), dimension(nLocal)   :: mfact
1161      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 963 | Line 1168 | contains
1168      real ( kind = dp ), intent(inout) :: rijsq
1169      real ( kind = dp )                :: r
1170      real ( kind = dp ), intent(inout) :: d(3)
966 –    real ( kind = dp ) :: ebalance
1171      integer :: me_i, me_j
1172  
1173 <    integer :: iMap
1173 >    integer :: iHash
1174  
1175      r = sqrt(rijsq)
1176      vpair = 0.0d0
# Line 980 | Line 1184 | contains
1184      me_j = atid(j)
1185   #endif
1186  
1187 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1187 >    iHash = InteractionHash(me_i, me_j)
1188  
1189 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1190 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1189 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1190 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1191 >            pot(VDW_POT), f, do_pot)
1192      endif
1193  
1194 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1194 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1195         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1196 <            pot, eFrame, f, t, do_pot)
1196 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1197  
1198 <       if (FF_uses_RF .and. SIM_uses_RF) then
1198 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1199  
1200            ! CHECK ME (RF needs to know about all electrostatic types)
1201            call accumulate_rf(i, j, r, eFrame, sw)
# Line 999 | Line 1204 | contains
1204  
1205      endif
1206  
1207 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1207 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1208         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1209 <            pot, A, f, t, do_pot)
1209 >            pot(HB_POT), A, f, t, do_pot)
1210      endif
1211  
1212 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1212 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1213         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1214 <            pot, A, f, t, do_pot)
1214 >            pot(HB_POT), A, f, t, do_pot)
1215      endif
1216  
1217 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1217 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1218         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1219 <            pot, A, f, t, do_pot)
1219 >            pot(VDW_POT), A, f, t, do_pot)
1220      endif
1221      
1222 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1223 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1224 <            pot, A, f, t, do_pot)
1222 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1223 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1224 >            pot(VDW_POT), A, f, t, do_pot)
1225      endif
1226  
1227 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1228 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1229 <            do_pot)
1227 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1228 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1229 >            pot(METALLIC_POT), f, do_pot)
1230      endif
1231  
1232 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1232 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1233         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1234 <            pot, A, f, t, do_pot)
1234 >            pot(VDW_POT), A, f, t, do_pot)
1235      endif
1236  
1237 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1237 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1238         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1239 <            pot, A, f, t, do_pot)
1239 >            pot(VDW_POT), A, f, t, do_pot)
1240      endif
1241      
1242    end subroutine do_pair
# Line 1039 | Line 1244 | contains
1244    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1245         do_pot, do_stress, eFrame, A, f, t, pot)
1246  
1247 <    real( kind = dp ) :: pot, sw
1247 >    real( kind = dp ) :: sw
1248 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1249      real( kind = dp ), dimension(9,nLocal) :: eFrame
1250      real (kind=dp), dimension(9,nLocal) :: A
1251      real (kind=dp), dimension(3,nLocal) :: f
# Line 1051 | Line 1257 | contains
1257      real ( kind = dp )                :: r, rc
1258      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1259  
1260 <    integer :: me_i, me_j, iMap
1260 >    integer :: me_i, me_j, iHash
1261  
1262 +    r = sqrt(rijsq)
1263 +
1264   #ifdef IS_MPI  
1265      me_i = atid_row(i)
1266      me_j = atid_col(j)  
# Line 1061 | Line 1269 | contains
1269      me_j = atid(j)  
1270   #endif
1271  
1272 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1272 >    iHash = InteractionHash(me_i, me_j)
1273  
1274 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1274 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1275              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1276      endif
1277      
# Line 1072 | Line 1280 | contains
1280  
1281    subroutine do_preforce(nlocal,pot)
1282      integer :: nlocal
1283 <    real( kind = dp ) :: pot
1283 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1284  
1285      if (FF_uses_EAM .and. SIM_uses_EAM) then
1286 <       call calc_EAM_preforce_Frho(nlocal,pot)
1286 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1287      endif
1288  
1289  
# Line 1260 | Line 1468 | contains
1468  
1469    function FF_UsesDirectionalAtoms() result(doesit)
1470      logical :: doesit
1471 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1264 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1265 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1471 >    doesit = FF_uses_DirectionalAtoms
1472    end function FF_UsesDirectionalAtoms
1473  
1474    function FF_RequiresPrepairCalc() result(doesit)
# Line 1270 | Line 1476 | contains
1476      doesit = FF_uses_EAM
1477    end function FF_RequiresPrepairCalc
1478  
1273 –  function FF_RequiresPostpairCalc() result(doesit)
1274 –    logical :: doesit
1275 –    doesit = FF_uses_RF
1276 –  end function FF_RequiresPostpairCalc
1277 –
1479   #ifdef PROFILE
1480    function getforcetime() result(totalforcetime)
1481      real(kind=dp) :: totalforcetime

Diff Legend

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