| 45 |
|
|
| 46 |
|
!! @author Charles F. Vardeman II |
| 47 |
|
!! @author Matthew Meineke |
| 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 $ |
| 48 |
> |
!! @version $Id: doForces.F90,v 1.47 2005-09-19 23:21:46 chrisfen Exp $, $Date: 2005-09-19 23:21:46 $, $Name: not supported by cvs2svn $, $Revision: 1.47 $ |
| 49 |
|
|
| 50 |
|
|
| 51 |
|
module doForces |
| 58 |
|
use lj |
| 59 |
|
use sticky |
| 60 |
|
use electrostatic_module |
| 61 |
< |
use reaction_field |
| 61 |
> |
use reaction_field_module |
| 62 |
|
use gb_pair |
| 63 |
|
use shapes |
| 64 |
|
use vector_class |
| 73 |
|
|
| 74 |
|
#define __FORTRAN90 |
| 75 |
|
#include "UseTheForce/fSwitchingFunction.h" |
| 76 |
+ |
#include "UseTheForce/fCutoffPolicy.h" |
| 77 |
|
#include "UseTheForce/DarkSide/fInteractionMap.h" |
| 78 |
+ |
#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" |
| 79 |
|
|
| 80 |
+ |
|
| 81 |
|
INTEGER, PARAMETER:: PREPAIR_LOOP = 1 |
| 82 |
|
INTEGER, PARAMETER:: PAIR_LOOP = 2 |
| 83 |
|
|
| 81 |
– |
logical, save :: haveRlist = .false. |
| 84 |
|
logical, save :: haveNeighborList = .false. |
| 85 |
|
logical, save :: haveSIMvariables = .false. |
| 86 |
|
logical, save :: haveSaneForceField = .false. |
| 87 |
< |
logical, save :: haveInteractionMap = .false. |
| 87 |
> |
logical, save :: haveInteractionHash = .false. |
| 88 |
> |
logical, save :: haveGtypeCutoffMap = .false. |
| 89 |
> |
logical, save :: haveRlist = .false. |
| 90 |
|
|
| 91 |
|
logical, save :: FF_uses_DirectionalAtoms |
| 88 |
– |
logical, save :: FF_uses_LennardJones |
| 89 |
– |
logical, save :: FF_uses_Electrostatics |
| 90 |
– |
logical, save :: FF_uses_Charges |
| 92 |
|
logical, save :: FF_uses_Dipoles |
| 92 |
– |
logical, save :: FF_uses_Quadrupoles |
| 93 |
– |
logical, save :: FF_uses_Sticky |
| 94 |
– |
logical, save :: FF_uses_StickyPower |
| 93 |
|
logical, save :: FF_uses_GayBerne |
| 94 |
|
logical, save :: FF_uses_EAM |
| 97 |
– |
logical, save :: FF_uses_Shapes |
| 98 |
– |
logical, save :: FF_uses_FLARB |
| 99 |
– |
logical, save :: FF_uses_RF |
| 95 |
|
|
| 96 |
|
logical, save :: SIM_uses_DirectionalAtoms |
| 102 |
– |
logical, save :: SIM_uses_LennardJones |
| 103 |
– |
logical, save :: SIM_uses_Electrostatics |
| 104 |
– |
logical, save :: SIM_uses_Charges |
| 105 |
– |
logical, save :: SIM_uses_Dipoles |
| 106 |
– |
logical, save :: SIM_uses_Quadrupoles |
| 107 |
– |
logical, save :: SIM_uses_Sticky |
| 108 |
– |
logical, save :: SIM_uses_StickyPower |
| 109 |
– |
logical, save :: SIM_uses_GayBerne |
| 97 |
|
logical, save :: SIM_uses_EAM |
| 111 |
– |
logical, save :: SIM_uses_Shapes |
| 112 |
– |
logical, save :: SIM_uses_FLARB |
| 113 |
– |
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 |
| 117 |
– |
logical, save :: SIM_uses_molecular_cutoffs |
| 101 |
|
|
| 102 |
< |
!!!GO AWAY--------- |
| 120 |
< |
!!!!!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 |
| 108 |
< |
!public :: addInteraction |
| 109 |
< |
!public :: setInteractionHash |
| 110 |
< |
!public :: getInteractionHash |
| 111 |
< |
public :: createInteractionMap |
| 112 |
< |
public :: createRcuts |
| 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 |
| 118 |
|
real :: forceTimeInitial, forceTimeFinal |
| 119 |
|
integer :: nLoops |
| 120 |
|
#endif |
| 137 |
– |
|
| 138 |
– |
type, public :: Interaction |
| 139 |
– |
integer :: InteractionHash |
| 140 |
– |
real(kind=dp) :: rList = 0.0_dp |
| 141 |
– |
real(kind=dp) :: rListSq = 0.0_dp |
| 142 |
– |
end type Interaction |
| 121 |
|
|
| 122 |
< |
type(Interaction), dimension(:,:),allocatable :: InteractionMap |
| 123 |
< |
|
| 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 :: groupMaxCutoff |
| 127 |
> |
integer, dimension(:), allocatable :: groupToGtype |
| 128 |
> |
real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff |
| 129 |
> |
type ::gtypeCutoffs |
| 130 |
> |
real(kind=dp) :: rcut |
| 131 |
> |
real(kind=dp) :: rcutsq |
| 132 |
> |
real(kind=dp) :: rlistsq |
| 133 |
> |
end type gtypeCutoffs |
| 134 |
> |
type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap |
| 135 |
|
|
| 136 |
+ |
integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY |
| 137 |
+ |
real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist |
| 138 |
|
|
| 139 |
|
contains |
| 140 |
|
|
| 141 |
< |
|
| 151 |
< |
subroutine createInteractionMap(status) |
| 141 |
> |
subroutine createInteractionHash(status) |
| 142 |
|
integer :: nAtypes |
| 143 |
|
integer, intent(out) :: status |
| 144 |
|
integer :: i |
| 145 |
|
integer :: j |
| 146 |
< |
integer :: ihash |
| 157 |
< |
real(kind=dp) :: myRcut |
| 146 |
> |
integer :: iHash |
| 147 |
|
!! Test Types |
| 148 |
|
logical :: i_is_LJ |
| 149 |
|
logical :: i_is_Elect |
| 159 |
|
logical :: j_is_GB |
| 160 |
|
logical :: j_is_EAM |
| 161 |
|
logical :: j_is_Shape |
| 162 |
< |
|
| 163 |
< |
status = 0 |
| 164 |
< |
|
| 162 |
> |
real(kind=dp) :: myRcut |
| 163 |
> |
|
| 164 |
> |
status = 0 |
| 165 |
> |
|
| 166 |
|
if (.not. associated(atypes)) then |
| 167 |
< |
call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!") |
| 167 |
> |
call handleError("atype", "atypes was not present before call of createInteractionHash!") |
| 168 |
|
status = -1 |
| 169 |
|
return |
| 170 |
|
endif |
| 176 |
|
return |
| 177 |
|
end if |
| 178 |
|
|
| 179 |
< |
if (.not. allocated(InteractionMap)) then |
| 180 |
< |
allocate(InteractionMap(nAtypes,nAtypes)) |
| 179 |
> |
if (.not. allocated(InteractionHash)) then |
| 180 |
> |
allocate(InteractionHash(nAtypes,nAtypes)) |
| 181 |
|
endif |
| 182 |
+ |
|
| 183 |
+ |
if (.not. allocated(atypeMaxCutoff)) then |
| 184 |
+ |
allocate(atypeMaxCutoff(nAtypes)) |
| 185 |
+ |
endif |
| 186 |
|
|
| 187 |
|
do i = 1, nAtypes |
| 188 |
|
call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) |
| 235 |
|
if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ) |
| 236 |
|
|
| 237 |
|
|
| 238 |
< |
InteractionMap(i,j)%InteractionHash = iHash |
| 239 |
< |
InteractionMap(j,i)%InteractionHash = iHash |
| 238 |
> |
InteractionHash(i,j) = iHash |
| 239 |
> |
InteractionHash(j,i) = iHash |
| 240 |
|
|
| 241 |
|
end do |
| 242 |
|
|
| 243 |
|
end do |
| 244 |
|
|
| 245 |
< |
haveInteractionMap = .true. |
| 246 |
< |
end subroutine createInteractionMap |
| 245 |
> |
haveInteractionHash = .true. |
| 246 |
> |
end subroutine createInteractionHash |
| 247 |
|
|
| 248 |
< |
! 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 |
| 248 |
> |
subroutine createGtypeCutoffMap(stat) |
| 249 |
|
|
| 250 |
< |
stat = 0 |
| 251 |
< |
if (.not. haveInteractionMap) then |
| 250 |
> |
integer, intent(out), optional :: stat |
| 251 |
> |
logical :: i_is_LJ |
| 252 |
> |
logical :: i_is_Elect |
| 253 |
> |
logical :: i_is_Sticky |
| 254 |
> |
logical :: i_is_StickyP |
| 255 |
> |
logical :: i_is_GB |
| 256 |
> |
logical :: i_is_EAM |
| 257 |
> |
logical :: i_is_Shape |
| 258 |
> |
logical :: GtypeFound |
| 259 |
|
|
| 260 |
< |
call createInteractionMap(myStatus) |
| 260 |
> |
integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend |
| 261 |
> |
integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes |
| 262 |
> |
integer :: nGroupsInRow |
| 263 |
> |
real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin |
| 264 |
> |
real(kind=dp) :: biggestAtypeCutoff |
| 265 |
|
|
| 266 |
+ |
stat = 0 |
| 267 |
+ |
if (.not. haveInteractionHash) then |
| 268 |
+ |
call createInteractionHash(myStatus) |
| 269 |
|
if (myStatus .ne. 0) then |
| 270 |
< |
write(default_error, *) 'createInteractionMap failed in doForces!' |
| 270 |
> |
write(default_error, *) 'createInteractionHash failed in doForces!' |
| 271 |
|
stat = -1 |
| 272 |
|
return |
| 273 |
|
endif |
| 274 |
|
endif |
| 275 |
< |
|
| 276 |
< |
|
| 275 |
> |
#ifdef IS_MPI |
| 276 |
> |
nGroupsInRow = getNgroupsInRow(plan_group_row) |
| 277 |
> |
#endif |
| 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 |
| 279 |
> |
! Set all of the initial cutoffs to zero. |
| 280 |
> |
atypeMaxCutoff = 0.0_dp |
| 281 |
> |
do i = 1, nAtypes |
| 282 |
> |
if (SimHasAtype(i)) then |
| 283 |
> |
call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) |
| 284 |
> |
call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect) |
| 285 |
> |
call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky) |
| 286 |
> |
call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP) |
| 287 |
> |
call getElementProperty(atypes, i, "is_GayBerne", i_is_GB) |
| 288 |
> |
call getElementProperty(atypes, i, "is_EAM", i_is_EAM) |
| 289 |
> |
call getElementProperty(atypes, i, "is_Shape", i_is_Shape) |
| 290 |
|
|
| 291 |
< |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
| 292 |
< |
! thisRCut = getLJCutOff(map_i,map_j) |
| 293 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 291 |
> |
|
| 292 |
> |
if (i_is_LJ) then |
| 293 |
> |
thisRcut = getSigma(i) * 2.5_dp |
| 294 |
> |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 295 |
|
endif |
| 296 |
< |
|
| 297 |
< |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 298 |
< |
! thisRCut = getElectrostaticCutOff(map_i,map_j) |
| 298 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 296 |
> |
if (i_is_Elect) then |
| 297 |
> |
thisRcut = defaultRcut |
| 298 |
> |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 299 |
|
endif |
| 300 |
+ |
if (i_is_Sticky) then |
| 301 |
+ |
thisRcut = getStickyCut(i) |
| 302 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 303 |
+ |
endif |
| 304 |
+ |
if (i_is_StickyP) then |
| 305 |
+ |
thisRcut = getStickyPowerCut(i) |
| 306 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 307 |
+ |
endif |
| 308 |
+ |
if (i_is_GB) then |
| 309 |
+ |
thisRcut = getGayBerneCut(i) |
| 310 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 311 |
+ |
endif |
| 312 |
+ |
if (i_is_EAM) then |
| 313 |
+ |
thisRcut = getEAMCut(i) |
| 314 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 315 |
+ |
endif |
| 316 |
+ |
if (i_is_Shape) then |
| 317 |
+ |
thisRcut = getShapeCut(i) |
| 318 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 319 |
+ |
endif |
| 320 |
|
|
| 321 |
< |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
| 322 |
< |
! thisRCut = getStickyCutOff(map_i,map_j) |
| 323 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 324 |
< |
endif |
| 325 |
< |
|
| 326 |
< |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
| 327 |
< |
! thisRCut = getStickyPowerCutOff(map_i,map_j) |
| 328 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 329 |
< |
endif |
| 330 |
< |
|
| 331 |
< |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
| 332 |
< |
! thisRCut = getGayberneCutOff(map_i,map_j) |
| 333 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 334 |
< |
endif |
| 335 |
< |
|
| 336 |
< |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
| 337 |
< |
! thisRCut = getGaybrneLJCutOff(map_i,map_j) |
| 338 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 339 |
< |
endif |
| 340 |
< |
|
| 341 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 342 |
< |
! thisRCut = getEAMCutOff(map_i,map_j) |
| 343 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 344 |
< |
endif |
| 345 |
< |
|
| 346 |
< |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
| 347 |
< |
! thisRCut = getShapeCutOff(map_i,map_j) |
| 348 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 349 |
< |
endif |
| 350 |
< |
|
| 351 |
< |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
| 352 |
< |
! thisRCut = getShapeLJCutOff(map_i,map_j) |
| 353 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 354 |
< |
endif |
| 355 |
< |
InteractionMap(map_i, map_j)%rList = actualCutoff |
| 356 |
< |
InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff |
| 357 |
< |
end do |
| 358 |
< |
end do |
| 359 |
< |
haveRlist = .true. |
| 360 |
< |
end subroutine createRcuts |
| 361 |
< |
|
| 362 |
< |
|
| 363 |
< |
!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF |
| 344 |
< |
!!$ subroutine setRlistDF( this_rlist ) |
| 345 |
< |
!!$ |
| 346 |
< |
!!$ real(kind=dp) :: this_rlist |
| 347 |
< |
!!$ |
| 348 |
< |
!!$ rlist = this_rlist |
| 349 |
< |
!!$ rlistsq = rlist * rlist |
| 350 |
< |
!!$ |
| 351 |
< |
!!$ haveRlist = .true. |
| 352 |
< |
!!$ |
| 353 |
< |
!!$ end subroutine setRlistDF |
| 321 |
> |
if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then |
| 322 |
> |
biggestAtypeCutoff = atypeMaxCutoff(i) |
| 323 |
> |
endif |
| 324 |
> |
endif |
| 325 |
> |
enddo |
| 326 |
> |
|
| 327 |
> |
nGroupTypes = 0 |
| 328 |
> |
|
| 329 |
> |
istart = 1 |
| 330 |
> |
#ifdef IS_MPI |
| 331 |
> |
iend = nGroupsInRow |
| 332 |
> |
#else |
| 333 |
> |
iend = nGroups |
| 334 |
> |
#endif |
| 335 |
> |
|
| 336 |
> |
!! allocate the groupToGtype and gtypeMaxCutoff here. |
| 337 |
> |
if(.not.allocated(groupToGtype)) then |
| 338 |
> |
allocate(groupToGtype(iend)) |
| 339 |
> |
allocate(groupMaxCutoff(iend)) |
| 340 |
> |
allocate(gtypeMaxCutoff(iend)) |
| 341 |
> |
groupMaxCutoff = 0.0_dp |
| 342 |
> |
gtypeMaxCutoff = 0.0_dp |
| 343 |
> |
endif |
| 344 |
> |
!! first we do a single loop over the cutoff groups to find the |
| 345 |
> |
!! largest cutoff for any atypes present in this group. We also |
| 346 |
> |
!! create gtypes at this point. |
| 347 |
> |
|
| 348 |
> |
tol = 1.0d-6 |
| 349 |
> |
|
| 350 |
> |
do i = istart, iend |
| 351 |
> |
n_in_i = groupStartRow(i+1) - groupStartRow(i) |
| 352 |
> |
groupMaxCutoff(i) = 0.0_dp |
| 353 |
> |
do ia = groupStartRow(i), groupStartRow(i+1)-1 |
| 354 |
> |
atom1 = groupListRow(ia) |
| 355 |
> |
#ifdef IS_MPI |
| 356 |
> |
me_i = atid_row(atom1) |
| 357 |
> |
#else |
| 358 |
> |
me_i = atid(atom1) |
| 359 |
> |
#endif |
| 360 |
> |
if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then |
| 361 |
> |
groupMaxCutoff(i)=atypeMaxCutoff(me_i) |
| 362 |
> |
endif |
| 363 |
> |
enddo |
| 364 |
|
|
| 365 |
+ |
if (nGroupTypes.eq.0) then |
| 366 |
+ |
nGroupTypes = nGroupTypes + 1 |
| 367 |
+ |
gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) |
| 368 |
+ |
groupToGtype(i) = nGroupTypes |
| 369 |
+ |
else |
| 370 |
+ |
GtypeFound = .false. |
| 371 |
+ |
do g = 1, nGroupTypes |
| 372 |
+ |
if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then |
| 373 |
+ |
groupToGtype(i) = g |
| 374 |
+ |
GtypeFound = .true. |
| 375 |
+ |
endif |
| 376 |
+ |
enddo |
| 377 |
+ |
if (.not.GtypeFound) then |
| 378 |
+ |
nGroupTypes = nGroupTypes + 1 |
| 379 |
+ |
gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) |
| 380 |
+ |
groupToGtype(i) = nGroupTypes |
| 381 |
+ |
endif |
| 382 |
+ |
endif |
| 383 |
+ |
enddo |
| 384 |
|
|
| 385 |
+ |
!! allocate the gtypeCutoffMap here. |
| 386 |
+ |
allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes)) |
| 387 |
+ |
!! then we do a double loop over all the group TYPES to find the cutoff |
| 388 |
+ |
!! map between groups of two types |
| 389 |
+ |
|
| 390 |
+ |
do i = 1, nGroupTypes |
| 391 |
+ |
do j = 1, nGroupTypes |
| 392 |
+ |
|
| 393 |
+ |
select case(cutoffPolicy) |
| 394 |
+ |
case(TRADITIONAL_CUTOFF_POLICY) |
| 395 |
+ |
thisRcut = maxval(gtypeMaxCutoff) |
| 396 |
+ |
case(MIX_CUTOFF_POLICY) |
| 397 |
+ |
thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j)) |
| 398 |
+ |
case(MAX_CUTOFF_POLICY) |
| 399 |
+ |
thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j)) |
| 400 |
+ |
case default |
| 401 |
+ |
call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy") |
| 402 |
+ |
return |
| 403 |
+ |
end select |
| 404 |
+ |
gtypeCutoffMap(i,j)%rcut = thisRcut |
| 405 |
+ |
gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut |
| 406 |
+ |
skin = defaultRlist - defaultRcut |
| 407 |
+ |
gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2 |
| 408 |
+ |
|
| 409 |
+ |
enddo |
| 410 |
+ |
enddo |
| 411 |
+ |
|
| 412 |
+ |
haveGtypeCutoffMap = .true. |
| 413 |
+ |
end subroutine createGtypeCutoffMap |
| 414 |
+ |
|
| 415 |
+ |
subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy) |
| 416 |
+ |
real(kind=dp),intent(in) :: defRcut, defRsw, defRlist |
| 417 |
+ |
integer, intent(in) :: cutPolicy |
| 418 |
+ |
|
| 419 |
+ |
defaultRcut = defRcut |
| 420 |
+ |
defaultRsw = defRsw |
| 421 |
+ |
defaultRlist = defRlist |
| 422 |
+ |
cutoffPolicy = cutPolicy |
| 423 |
+ |
end subroutine setDefaultCutoffs |
| 424 |
+ |
|
| 425 |
+ |
subroutine setCutoffPolicy(cutPolicy) |
| 426 |
+ |
|
| 427 |
+ |
integer, intent(in) :: cutPolicy |
| 428 |
+ |
cutoffPolicy = cutPolicy |
| 429 |
+ |
call createGtypeCutoffMap() |
| 430 |
+ |
end subroutine setCutoffPolicy |
| 431 |
+ |
|
| 432 |
+ |
|
| 433 |
|
subroutine setSimVariables() |
| 434 |
|
SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms() |
| 358 |
– |
SIM_uses_LennardJones = SimUsesLennardJones() |
| 359 |
– |
SIM_uses_Electrostatics = SimUsesElectrostatics() |
| 360 |
– |
SIM_uses_Charges = SimUsesCharges() |
| 361 |
– |
SIM_uses_Dipoles = SimUsesDipoles() |
| 362 |
– |
SIM_uses_Sticky = SimUsesSticky() |
| 363 |
– |
SIM_uses_StickyPower = SimUsesStickyPower() |
| 364 |
– |
SIM_uses_GayBerne = SimUsesGayBerne() |
| 435 |
|
SIM_uses_EAM = SimUsesEAM() |
| 366 |
– |
SIM_uses_Shapes = SimUsesShapes() |
| 367 |
– |
SIM_uses_FLARB = SimUsesFLARB() |
| 368 |
– |
SIM_uses_RF = SimUsesRF() |
| 436 |
|
SIM_requires_postpair_calc = SimRequiresPostpairCalc() |
| 437 |
|
SIM_requires_prepair_calc = SimRequiresPrepairCalc() |
| 438 |
|
SIM_uses_PBC = SimUsesPBC() |
| 449 |
|
|
| 450 |
|
error = 0 |
| 451 |
|
|
| 452 |
< |
if (.not. haveInteractionMap) then |
| 452 |
> |
if (.not. haveInteractionHash) then |
| 453 |
> |
myStatus = 0 |
| 454 |
> |
call createInteractionHash(myStatus) |
| 455 |
> |
if (myStatus .ne. 0) then |
| 456 |
> |
write(default_error, *) 'createInteractionHash failed in doForces!' |
| 457 |
> |
error = -1 |
| 458 |
> |
return |
| 459 |
> |
endif |
| 460 |
> |
endif |
| 461 |
|
|
| 462 |
< |
myStatus = 0 |
| 463 |
< |
|
| 464 |
< |
call createInteractionMap(myStatus) |
| 390 |
< |
|
| 462 |
> |
if (.not. haveGtypeCutoffMap) then |
| 463 |
> |
myStatus = 0 |
| 464 |
> |
call createGtypeCutoffMap(myStatus) |
| 465 |
|
if (myStatus .ne. 0) then |
| 466 |
< |
write(default_error, *) 'createInteractionMap failed in doForces!' |
| 466 |
> |
write(default_error, *) 'createGtypeCutoffMap failed in doForces!' |
| 467 |
|
error = -1 |
| 468 |
|
return |
| 469 |
|
endif |
| 473 |
|
call setSimVariables() |
| 474 |
|
endif |
| 475 |
|
|
| 476 |
< |
if (.not. haveRlist) then |
| 477 |
< |
write(default_error, *) 'rList has not been set in doForces!' |
| 478 |
< |
error = -1 |
| 479 |
< |
return |
| 480 |
< |
endif |
| 476 |
> |
! if (.not. haveRlist) then |
| 477 |
> |
! write(default_error, *) 'rList has not been set in doForces!' |
| 478 |
> |
! error = -1 |
| 479 |
> |
! return |
| 480 |
> |
! endif |
| 481 |
|
|
| 482 |
|
if (.not. haveNeighborList) then |
| 483 |
|
write(default_error, *) 'neighbor list has not been initialized in doForces!' |
| 502 |
|
end subroutine doReadyCheck |
| 503 |
|
|
| 504 |
|
|
| 505 |
< |
subroutine init_FF(use_RF_c, thisStat) |
| 505 |
> |
subroutine init_FF(thisESM, thisStat) |
| 506 |
|
|
| 507 |
< |
logical, intent(in) :: use_RF_c |
| 434 |
< |
|
| 507 |
> |
integer, intent(in) :: thisESM |
| 508 |
|
integer, intent(out) :: thisStat |
| 509 |
|
integer :: my_status, nMatches |
| 510 |
|
integer, pointer :: MatchList(:) => null() |
| 513 |
|
!! assume things are copacetic, unless they aren't |
| 514 |
|
thisStat = 0 |
| 515 |
|
|
| 516 |
< |
!! Fortran's version of a cast: |
| 444 |
< |
FF_uses_RF = use_RF_c |
| 516 |
> |
electrostaticSummationMethod = thisESM |
| 517 |
|
|
| 518 |
|
!! init_FF is called *after* all of the atom types have been |
| 519 |
|
!! defined in atype_module using the new_atype subroutine. |
| 522 |
|
!! interactions are used by the force field. |
| 523 |
|
|
| 524 |
|
FF_uses_DirectionalAtoms = .false. |
| 453 |
– |
FF_uses_LennardJones = .false. |
| 454 |
– |
FF_uses_Electrostatics = .false. |
| 455 |
– |
FF_uses_Charges = .false. |
| 525 |
|
FF_uses_Dipoles = .false. |
| 457 |
– |
FF_uses_Sticky = .false. |
| 458 |
– |
FF_uses_StickyPower = .false. |
| 526 |
|
FF_uses_GayBerne = .false. |
| 527 |
|
FF_uses_EAM = .false. |
| 461 |
– |
FF_uses_Shapes = .false. |
| 462 |
– |
FF_uses_FLARB = .false. |
| 528 |
|
|
| 529 |
|
call getMatchingElementList(atypes, "is_Directional", .true., & |
| 530 |
|
nMatches, MatchList) |
| 531 |
|
if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true. |
| 467 |
– |
|
| 468 |
– |
call getMatchingElementList(atypes, "is_LennardJones", .true., & |
| 469 |
– |
nMatches, MatchList) |
| 470 |
– |
if (nMatches .gt. 0) FF_uses_LennardJones = .true. |
| 471 |
– |
|
| 472 |
– |
call getMatchingElementList(atypes, "is_Electrostatic", .true., & |
| 473 |
– |
nMatches, MatchList) |
| 474 |
– |
if (nMatches .gt. 0) then |
| 475 |
– |
FF_uses_Electrostatics = .true. |
| 476 |
– |
endif |
| 532 |
|
|
| 478 |
– |
call getMatchingElementList(atypes, "is_Charge", .true., & |
| 479 |
– |
nMatches, MatchList) |
| 480 |
– |
if (nMatches .gt. 0) then |
| 481 |
– |
FF_uses_Charges = .true. |
| 482 |
– |
FF_uses_Electrostatics = .true. |
| 483 |
– |
endif |
| 484 |
– |
|
| 533 |
|
call getMatchingElementList(atypes, "is_Dipole", .true., & |
| 534 |
|
nMatches, MatchList) |
| 535 |
< |
if (nMatches .gt. 0) then |
| 488 |
< |
FF_uses_Dipoles = .true. |
| 489 |
< |
FF_uses_Electrostatics = .true. |
| 490 |
< |
FF_uses_DirectionalAtoms = .true. |
| 491 |
< |
endif |
| 492 |
< |
|
| 493 |
< |
call getMatchingElementList(atypes, "is_Quadrupole", .true., & |
| 494 |
< |
nMatches, MatchList) |
| 495 |
< |
if (nMatches .gt. 0) then |
| 496 |
< |
FF_uses_Quadrupoles = .true. |
| 497 |
< |
FF_uses_Electrostatics = .true. |
| 498 |
< |
FF_uses_DirectionalAtoms = .true. |
| 499 |
< |
endif |
| 500 |
< |
|
| 501 |
< |
call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, & |
| 502 |
< |
MatchList) |
| 503 |
< |
if (nMatches .gt. 0) then |
| 504 |
< |
FF_uses_Sticky = .true. |
| 505 |
< |
FF_uses_DirectionalAtoms = .true. |
| 506 |
< |
endif |
| 507 |
< |
|
| 508 |
< |
call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, & |
| 509 |
< |
MatchList) |
| 510 |
< |
if (nMatches .gt. 0) then |
| 511 |
< |
FF_uses_StickyPower = .true. |
| 512 |
< |
FF_uses_DirectionalAtoms = .true. |
| 513 |
< |
endif |
| 535 |
> |
if (nMatches .gt. 0) FF_uses_Dipoles = .true. |
| 536 |
|
|
| 537 |
|
call getMatchingElementList(atypes, "is_GayBerne", .true., & |
| 538 |
|
nMatches, MatchList) |
| 539 |
< |
if (nMatches .gt. 0) then |
| 518 |
< |
FF_uses_GayBerne = .true. |
| 519 |
< |
FF_uses_DirectionalAtoms = .true. |
| 520 |
< |
endif |
| 539 |
> |
if (nMatches .gt. 0) FF_uses_GayBerne = .true. |
| 540 |
|
|
| 541 |
|
call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList) |
| 542 |
|
if (nMatches .gt. 0) FF_uses_EAM = .true. |
| 543 |
|
|
| 525 |
– |
call getMatchingElementList(atypes, "is_Shape", .true., & |
| 526 |
– |
nMatches, MatchList) |
| 527 |
– |
if (nMatches .gt. 0) then |
| 528 |
– |
FF_uses_Shapes = .true. |
| 529 |
– |
FF_uses_DirectionalAtoms = .true. |
| 530 |
– |
endif |
| 544 |
|
|
| 532 |
– |
call getMatchingElementList(atypes, "is_FLARB", .true., & |
| 533 |
– |
nMatches, MatchList) |
| 534 |
– |
if (nMatches .gt. 0) FF_uses_FLARB = .true. |
| 535 |
– |
|
| 536 |
– |
!! Assume sanity (for the sake of argument) |
| 545 |
|
haveSaneForceField = .true. |
| 546 |
|
|
| 547 |
< |
!! check to make sure the FF_uses_RF setting makes sense |
| 547 |
> |
!! check to make sure the reaction field setting makes sense |
| 548 |
|
|
| 549 |
< |
if (FF_uses_dipoles) then |
| 550 |
< |
if (FF_uses_RF) then |
| 549 |
> |
if (FF_uses_Dipoles) then |
| 550 |
> |
if (electrostaticSummationMethod == REACTION_FIELD) then |
| 551 |
|
dielect = getDielect() |
| 552 |
|
call initialize_rf(dielect) |
| 553 |
|
endif |
| 554 |
|
else |
| 555 |
< |
if (FF_uses_RF) then |
| 555 |
> |
if (electrostaticSummationMethod == REACTION_FIELD) then |
| 556 |
|
write(default_error,*) 'Using Reaction Field with no dipoles? Huh?' |
| 557 |
|
thisStat = -1 |
| 558 |
|
haveSaneForceField = .false. |
| 559 |
|
return |
| 560 |
|
endif |
| 561 |
|
endif |
| 554 |
– |
|
| 555 |
– |
!sticky module does not contain check_sticky_FF anymore |
| 556 |
– |
!if (FF_uses_sticky) then |
| 557 |
– |
! call check_sticky_FF(my_status) |
| 558 |
– |
! if (my_status /= 0) then |
| 559 |
– |
! thisStat = -1 |
| 560 |
– |
! haveSaneForceField = .false. |
| 561 |
– |
! return |
| 562 |
– |
! end if |
| 563 |
– |
!endif |
| 562 |
|
|
| 563 |
|
if (FF_uses_EAM) then |
| 564 |
|
call init_EAM_FF(my_status) |
| 579 |
|
endif |
| 580 |
|
endif |
| 581 |
|
|
| 584 |
– |
if (FF_uses_GayBerne .and. FF_uses_LennardJones) then |
| 585 |
– |
endif |
| 586 |
– |
|
| 582 |
|
if (.not. haveNeighborList) then |
| 583 |
|
!! Create neighbor lists |
| 584 |
|
call expandNeighborList(nLocal, my_status) |
| 642 |
|
integer :: localError |
| 643 |
|
integer :: propPack_i, propPack_j |
| 644 |
|
integer :: loopStart, loopEnd, loop |
| 645 |
< |
integer :: iMap |
| 645 |
> |
integer :: iHash |
| 646 |
|
real(kind=dp) :: listSkin = 1.0 |
| 647 |
|
|
| 648 |
|
!! initialize local variables |
| 733 |
|
iend = nGroups - 1 |
| 734 |
|
#endif |
| 735 |
|
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 |
| 736 |
|
|
| 737 |
|
if (update_nlist) point(i) = nlist + 1 |
| 738 |
|
|
| 770 |
|
q_group(:,j), d_grp, rgrpsq) |
| 771 |
|
#endif |
| 772 |
|
|
| 773 |
< |
if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then |
| 773 |
> |
if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then |
| 774 |
|
if (update_nlist) then |
| 775 |
|
nlist = nlist + 1 |
| 776 |
|
|
| 970 |
|
|
| 971 |
|
if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then |
| 972 |
|
|
| 973 |
< |
if (FF_uses_RF .and. SIM_uses_RF) then |
| 973 |
> |
if (electrostaticSummationMethod == REACTION_FIELD) then |
| 974 |
|
|
| 975 |
|
#ifdef IS_MPI |
| 976 |
|
call scatter(rf_Row,rf,plan_atom_row_3d) |
| 988 |
|
#else |
| 989 |
|
me_i = atid(i) |
| 990 |
|
#endif |
| 991 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 991 |
> |
iHash = InteractionHash(me_i,me_j) |
| 992 |
|
|
| 993 |
< |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 993 |
> |
if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then |
| 994 |
|
|
| 995 |
|
mu_i = getDipoleMoment(me_i) |
| 996 |
|
|
| 1054 |
|
real ( kind = dp ), intent(inout) :: rijsq |
| 1055 |
|
real ( kind = dp ) :: r |
| 1056 |
|
real ( kind = dp ), intent(inout) :: d(3) |
| 1068 |
– |
real ( kind = dp ) :: ebalance |
| 1057 |
|
integer :: me_i, me_j |
| 1058 |
|
|
| 1059 |
< |
integer :: iMap |
| 1059 |
> |
integer :: iHash |
| 1060 |
|
|
| 1061 |
|
r = sqrt(rijsq) |
| 1062 |
|
vpair = 0.0d0 |
| 1070 |
|
me_j = atid(j) |
| 1071 |
|
#endif |
| 1072 |
|
|
| 1073 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 1073 |
> |
iHash = InteractionHash(me_i, me_j) |
| 1074 |
|
|
| 1075 |
< |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
| 1075 |
> |
if ( iand(iHash, LJ_PAIR).ne.0 ) then |
| 1076 |
|
call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) |
| 1077 |
|
endif |
| 1078 |
|
|
| 1079 |
< |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1079 |
> |
if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1080 |
|
call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1081 |
|
pot, eFrame, f, t, do_pot) |
| 1082 |
|
|
| 1083 |
< |
if (FF_uses_RF .and. SIM_uses_RF) then |
| 1083 |
> |
if (electrostaticSummationMethod == REACTION_FIELD) then |
| 1084 |
|
|
| 1085 |
|
! CHECK ME (RF needs to know about all electrostatic types) |
| 1086 |
|
call accumulate_rf(i, j, r, eFrame, sw) |
| 1089 |
|
|
| 1090 |
|
endif |
| 1091 |
|
|
| 1092 |
< |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
| 1092 |
> |
if ( iand(iHash, STICKY_PAIR).ne.0 ) then |
| 1093 |
|
call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1094 |
|
pot, A, f, t, do_pot) |
| 1095 |
|
endif |
| 1096 |
|
|
| 1097 |
< |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
| 1097 |
> |
if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then |
| 1098 |
|
call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1099 |
|
pot, A, f, t, do_pot) |
| 1100 |
|
endif |
| 1101 |
|
|
| 1102 |
< |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
| 1102 |
> |
if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then |
| 1103 |
|
call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1104 |
|
pot, A, f, t, do_pot) |
| 1105 |
|
endif |
| 1106 |
|
|
| 1107 |
< |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
| 1107 |
> |
if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then |
| 1108 |
|
! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1109 |
|
! pot, A, f, t, do_pot) |
| 1110 |
|
endif |
| 1111 |
|
|
| 1112 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 1112 |
> |
if ( iand(iHash, EAM_PAIR).ne.0 ) then |
| 1113 |
|
call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, & |
| 1114 |
|
do_pot) |
| 1115 |
|
endif |
| 1116 |
|
|
| 1117 |
< |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
| 1117 |
> |
if ( iand(iHash, SHAPE_PAIR).ne.0 ) then |
| 1118 |
|
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1119 |
|
pot, A, f, t, do_pot) |
| 1120 |
|
endif |
| 1121 |
|
|
| 1122 |
< |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
| 1122 |
> |
if ( iand(iHash, SHAPE_LJ).ne.0 ) then |
| 1123 |
|
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1124 |
|
pot, A, f, t, do_pot) |
| 1125 |
|
endif |
| 1141 |
|
real ( kind = dp ) :: r, rc |
| 1142 |
|
real ( kind = dp ), intent(inout) :: d(3), dc(3) |
| 1143 |
|
|
| 1144 |
< |
integer :: me_i, me_j, iMap |
| 1144 |
> |
integer :: me_i, me_j, iHash |
| 1145 |
|
|
| 1146 |
+ |
r = sqrt(rijsq) |
| 1147 |
+ |
|
| 1148 |
|
#ifdef IS_MPI |
| 1149 |
|
me_i = atid_row(i) |
| 1150 |
|
me_j = atid_col(j) |
| 1153 |
|
me_j = atid(j) |
| 1154 |
|
#endif |
| 1155 |
|
|
| 1156 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 1156 |
> |
iHash = InteractionHash(me_i, me_j) |
| 1157 |
|
|
| 1158 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 1158 |
> |
if ( iand(iHash, EAM_PAIR).ne.0 ) then |
| 1159 |
|
call calc_EAM_prepair_rho(i, j, d, r, rijsq ) |
| 1160 |
|
endif |
| 1161 |
|
|
| 1352 |
|
|
| 1353 |
|
function FF_UsesDirectionalAtoms() result(doesit) |
| 1354 |
|
logical :: doesit |
| 1355 |
< |
doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. & |
| 1366 |
< |
FF_uses_Quadrupoles .or. FF_uses_Sticky .or. & |
| 1367 |
< |
FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes |
| 1355 |
> |
doesit = FF_uses_DirectionalAtoms |
| 1356 |
|
end function FF_UsesDirectionalAtoms |
| 1357 |
|
|
| 1358 |
|
function FF_RequiresPrepairCalc() result(doesit) |
| 1362 |
|
|
| 1363 |
|
function FF_RequiresPostpairCalc() result(doesit) |
| 1364 |
|
logical :: doesit |
| 1365 |
< |
doesit = FF_uses_RF |
| 1365 |
> |
if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true. |
| 1366 |
|
end function FF_RequiresPostpairCalc |
| 1367 |
|
|
| 1368 |
|
#ifdef PROFILE |
