| 45 |
|
|
| 46 |
|
!! @author Charles F. Vardeman II |
| 47 |
|
!! @author Matthew Meineke |
| 48 |
< |
!! @version $Id: doForces.F90,v 1.23 2005-07-03 20:53:43 chuckv Exp $, $Date: 2005-07-03 20:53:43 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $ |
| 48 |
> |
!! @version $Id: doForces.F90,v 1.37 2005-09-07 19:18:54 gezelter Exp $, $Date: 2005-09-07 19:18:54 $, $Name: not supported by cvs2svn $, $Revision: 1.37 $ |
| 49 |
|
|
| 50 |
|
|
| 51 |
|
module doForces |
| 73 |
|
|
| 74 |
|
#define __FORTRAN90 |
| 75 |
|
#include "UseTheForce/fSwitchingFunction.h" |
| 76 |
+ |
#include "UseTheForce/fCutoffPolicy.h" |
| 77 |
|
#include "UseTheForce/DarkSide/fInteractionMap.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 :: haveInteractionMap = .false. |
| 86 |
> |
logical, save :: haveInteractionHash = .false. |
| 87 |
> |
logical, save :: haveGtypeCutoffMap = .false. |
| 88 |
> |
logical, save :: haveRlist = .false. |
| 89 |
|
|
| 90 |
|
logical, save :: FF_uses_DirectionalAtoms |
| 88 |
– |
logical, save :: FF_uses_LennardJones |
| 89 |
– |
logical, save :: FF_uses_Electrostatics |
| 90 |
– |
logical, save :: FF_uses_Charges |
| 91 |
|
logical, save :: FF_uses_Dipoles |
| 92 |
– |
logical, save :: FF_uses_Quadrupoles |
| 93 |
– |
logical, save :: FF_uses_Sticky |
| 94 |
– |
logical, save :: FF_uses_StickyPower |
| 92 |
|
logical, save :: FF_uses_GayBerne |
| 93 |
|
logical, save :: FF_uses_EAM |
| 97 |
– |
logical, save :: FF_uses_Shapes |
| 98 |
– |
logical, save :: FF_uses_FLARB |
| 94 |
|
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 |
| 98 |
|
logical, save :: SIM_uses_RF |
| 99 |
|
logical, save :: SIM_requires_postpair_calc |
| 100 |
|
logical, save :: SIM_requires_prepair_calc |
| 101 |
|
logical, save :: SIM_uses_PBC |
| 117 |
– |
logical, save :: SIM_uses_molecular_cutoffs |
| 102 |
|
|
| 103 |
< |
!!!GO AWAY--------- |
| 120 |
< |
!!!!!real(kind=dp), save :: rlist, rlistsq |
| 103 |
> |
integer, save :: corrMethod |
| 104 |
|
|
| 105 |
|
public :: init_FF |
| 106 |
+ |
public :: setDefaultCutoffs |
| 107 |
|
public :: do_force_loop |
| 108 |
< |
! public :: setRlistDF |
| 109 |
< |
!public :: addInteraction |
| 110 |
< |
!public :: setInteractionHash |
| 111 |
< |
!public :: getInteractionHash |
| 112 |
< |
public :: createInteractionMap |
| 113 |
< |
public :: createRcuts |
| 108 |
> |
public :: createInteractionHash |
| 109 |
> |
public :: createGtypeCutoffMap |
| 110 |
> |
public :: getStickyCut |
| 111 |
> |
public :: getStickyPowerCut |
| 112 |
> |
public :: getGayBerneCut |
| 113 |
> |
public :: getEAMCut |
| 114 |
> |
public :: getShapeCut |
| 115 |
|
|
| 116 |
|
#ifdef PROFILE |
| 117 |
|
public :: getforcetime |
| 119 |
|
real :: forceTimeInitial, forceTimeFinal |
| 120 |
|
integer :: nLoops |
| 121 |
|
#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 |
| 122 |
|
|
| 123 |
< |
type(Interaction), dimension(:,:),allocatable :: InteractionMap |
| 124 |
< |
|
| 123 |
> |
!! Variables for cutoff mapping and interaction mapping |
| 124 |
> |
! Bit hash to determine pair-pair interactions. |
| 125 |
> |
integer, dimension(:,:), allocatable :: InteractionHash |
| 126 |
> |
real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff |
| 127 |
> |
real(kind=dp), dimension(:), allocatable :: groupMaxCutoff |
| 128 |
> |
integer, dimension(:), allocatable :: groupToGtype |
| 129 |
> |
real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff |
| 130 |
> |
type ::gtypeCutoffs |
| 131 |
> |
real(kind=dp) :: rcut |
| 132 |
> |
real(kind=dp) :: rcutsq |
| 133 |
> |
real(kind=dp) :: rlistsq |
| 134 |
> |
end type gtypeCutoffs |
| 135 |
> |
type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap |
| 136 |
|
|
| 137 |
+ |
integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY |
| 138 |
+ |
real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist |
| 139 |
|
|
| 140 |
|
contains |
| 141 |
|
|
| 142 |
< |
|
| 151 |
< |
subroutine createInteractionMap(status) |
| 142 |
> |
subroutine createInteractionHash(status) |
| 143 |
|
integer :: nAtypes |
| 144 |
< |
integer :: status |
| 144 |
> |
integer, intent(out) :: status |
| 145 |
|
integer :: i |
| 146 |
|
integer :: j |
| 147 |
< |
integer :: ihash |
| 148 |
< |
real(kind=dp) :: myRcut |
| 158 |
< |
! Test Types |
| 147 |
> |
integer :: iHash |
| 148 |
> |
!! Test Types |
| 149 |
|
logical :: i_is_LJ |
| 150 |
|
logical :: i_is_Elect |
| 151 |
|
logical :: i_is_Sticky |
| 160 |
|
logical :: j_is_GB |
| 161 |
|
logical :: j_is_EAM |
| 162 |
|
logical :: j_is_Shape |
| 163 |
< |
|
| 164 |
< |
|
| 163 |
> |
real(kind=dp) :: myRcut |
| 164 |
> |
|
| 165 |
> |
status = 0 |
| 166 |
> |
|
| 167 |
|
if (.not. associated(atypes)) then |
| 168 |
< |
call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!") |
| 168 |
> |
call handleError("atype", "atypes was not present before call of createInteractionHash!") |
| 169 |
|
status = -1 |
| 170 |
|
return |
| 171 |
|
endif |
| 177 |
|
return |
| 178 |
|
end if |
| 179 |
|
|
| 180 |
< |
if (.not. allocated(InteractionMap)) then |
| 181 |
< |
allocate(InteractionMap(nAtypes,nAtypes)) |
| 180 |
> |
if (.not. allocated(InteractionHash)) then |
| 181 |
> |
allocate(InteractionHash(nAtypes,nAtypes)) |
| 182 |
|
endif |
| 183 |
+ |
|
| 184 |
+ |
if (.not. allocated(atypeMaxCutoff)) then |
| 185 |
+ |
allocate(atypeMaxCutoff(nAtypes)) |
| 186 |
+ |
endif |
| 187 |
|
|
| 188 |
|
do i = 1, nAtypes |
| 189 |
|
call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) |
| 236 |
|
if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ) |
| 237 |
|
|
| 238 |
|
|
| 239 |
< |
InteractionMap(i,j)%InteractionHash = iHash |
| 240 |
< |
InteractionMap(j,i)%InteractionHash = iHash |
| 239 |
> |
InteractionHash(i,j) = iHash |
| 240 |
> |
InteractionHash(j,i) = iHash |
| 241 |
|
|
| 242 |
|
end do |
| 243 |
|
|
| 244 |
|
end do |
| 249 |
– |
end subroutine createInteractionMap |
| 245 |
|
|
| 246 |
< |
! 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. |
| 247 |
< |
subroutine createRcuts(defaultRList) |
| 253 |
< |
real(kind=dp), intent(in), optional :: defaultRList |
| 254 |
< |
integer :: iMap |
| 255 |
< |
integer :: map_i,map_j |
| 256 |
< |
real(kind=dp) :: thisRCut = 0.0_dp |
| 257 |
< |
real(kind=dp) :: actualCutoff = 0.0_dp |
| 258 |
< |
integer :: nAtypes |
| 246 |
> |
haveInteractionHash = .true. |
| 247 |
> |
end subroutine createInteractionHash |
| 248 |
|
|
| 249 |
< |
if(.not. allocated(InteractionMap)) return |
| 249 |
> |
subroutine createGtypeCutoffMap(stat) |
| 250 |
|
|
| 251 |
< |
nAtypes = getSize(atypes) |
| 252 |
< |
! If we pass a default rcut, set all atypes to that cutoff distance |
| 253 |
< |
if(present(defaultRList)) then |
| 254 |
< |
InteractionMap(:,:)%rList = defaultRList |
| 255 |
< |
InteractionMap(:,:)%rListSq = defaultRList*defaultRList |
| 256 |
< |
haveRlist = .true. |
| 257 |
< |
return |
| 258 |
< |
end if |
| 251 |
> |
integer, intent(out), optional :: stat |
| 252 |
> |
logical :: i_is_LJ |
| 253 |
> |
logical :: i_is_Elect |
| 254 |
> |
logical :: i_is_Sticky |
| 255 |
> |
logical :: i_is_StickyP |
| 256 |
> |
logical :: i_is_GB |
| 257 |
> |
logical :: i_is_EAM |
| 258 |
> |
logical :: i_is_Shape |
| 259 |
|
|
| 260 |
< |
do map_i = 1,nAtypes |
| 261 |
< |
do map_j = map_i,nAtypes |
| 262 |
< |
iMap = InteractionMap(map_i, map_j)%InteractionHash |
| 263 |
< |
|
| 264 |
< |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
| 265 |
< |
! thisRCut = getLJCutOff(map_i,map_j) |
| 266 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 267 |
< |
endif |
| 260 |
> |
integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend |
| 261 |
> |
integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes |
| 262 |
> |
real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin |
| 263 |
> |
real(kind=dp) :: biggestAtypeCutoff |
| 264 |
> |
|
| 265 |
> |
stat = 0 |
| 266 |
> |
if (.not. haveInteractionHash) then |
| 267 |
> |
call createInteractionHash(myStatus) |
| 268 |
> |
if (myStatus .ne. 0) then |
| 269 |
> |
write(default_error, *) 'createInteractionHash failed in doForces!' |
| 270 |
> |
stat = -1 |
| 271 |
> |
return |
| 272 |
> |
endif |
| 273 |
> |
endif |
| 274 |
> |
|
| 275 |
> |
nAtypes = getSize(atypes) |
| 276 |
> |
|
| 277 |
> |
do i = 1, nAtypes |
| 278 |
> |
if (SimHasAtype(i)) then |
| 279 |
> |
call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) |
| 280 |
> |
call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect) |
| 281 |
> |
call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky) |
| 282 |
> |
call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP) |
| 283 |
> |
call getElementProperty(atypes, i, "is_GayBerne", i_is_GB) |
| 284 |
> |
call getElementProperty(atypes, i, "is_EAM", i_is_EAM) |
| 285 |
> |
call getElementProperty(atypes, i, "is_Shape", i_is_Shape) |
| 286 |
|
|
| 287 |
< |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 288 |
< |
! thisRCut = getElectrostaticCutOff(map_i,map_j) |
| 289 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 287 |
> |
atypeMaxCutoff(i) = 0.0_dp |
| 288 |
> |
if (i_is_LJ) then |
| 289 |
> |
thisRcut = getSigma(i) * 2.5_dp |
| 290 |
> |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 291 |
|
endif |
| 292 |
+ |
if (i_is_Elect) then |
| 293 |
+ |
thisRcut = defaultRcut |
| 294 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 295 |
+ |
endif |
| 296 |
+ |
if (i_is_Sticky) then |
| 297 |
+ |
thisRcut = getStickyCut(i) |
| 298 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 299 |
+ |
endif |
| 300 |
+ |
if (i_is_StickyP) then |
| 301 |
+ |
thisRcut = getStickyPowerCut(i) |
| 302 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 303 |
+ |
endif |
| 304 |
+ |
if (i_is_GB) then |
| 305 |
+ |
thisRcut = getGayBerneCut(i) |
| 306 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 307 |
+ |
endif |
| 308 |
+ |
if (i_is_EAM) then |
| 309 |
+ |
thisRcut = getEAMCut(i) |
| 310 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 311 |
+ |
endif |
| 312 |
+ |
if (i_is_Shape) then |
| 313 |
+ |
thisRcut = getShapeCut(i) |
| 314 |
+ |
if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut |
| 315 |
+ |
endif |
| 316 |
|
|
| 317 |
< |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
| 318 |
< |
! thisRCut = getStickyCutOff(map_i,map_j) |
| 319 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 320 |
< |
endif |
| 321 |
< |
|
| 322 |
< |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
| 323 |
< |
! thisRCut = getStickyPowerCutOff(map_i,map_j) |
| 324 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 325 |
< |
endif |
| 326 |
< |
|
| 327 |
< |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
| 328 |
< |
! thisRCut = getGayberneCutOff(map_i,map_j) |
| 329 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 330 |
< |
endif |
| 331 |
< |
|
| 332 |
< |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
| 333 |
< |
! thisRCut = getGaybrneLJCutOff(map_i,map_j) |
| 334 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 335 |
< |
endif |
| 336 |
< |
|
| 337 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 338 |
< |
! thisRCut = getEAMCutOff(map_i,map_j) |
| 339 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 340 |
< |
endif |
| 341 |
< |
|
| 342 |
< |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
| 343 |
< |
! thisRCut = getShapeCutOff(map_i,map_j) |
| 344 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 345 |
< |
endif |
| 346 |
< |
|
| 347 |
< |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
| 348 |
< |
! thisRCut = getShapeLJCutOff(map_i,map_j) |
| 349 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 350 |
< |
endif |
| 351 |
< |
InteractionMap(map_i, map_j)%rList = actualCutoff |
| 352 |
< |
InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff |
| 353 |
< |
end do |
| 354 |
< |
end do |
| 355 |
< |
haveRlist = .true. |
| 356 |
< |
end subroutine createRcuts |
| 357 |
< |
|
| 358 |
< |
|
| 359 |
< |
!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF |
| 360 |
< |
!!$ subroutine setRlistDF( this_rlist ) |
| 361 |
< |
!!$ |
| 362 |
< |
!!$ real(kind=dp) :: this_rlist |
| 363 |
< |
!!$ |
| 364 |
< |
!!$ rlist = this_rlist |
| 365 |
< |
!!$ rlistsq = rlist * rlist |
| 366 |
< |
!!$ |
| 367 |
< |
!!$ haveRlist = .true. |
| 368 |
< |
!!$ |
| 369 |
< |
!!$ end subroutine setRlistDF |
| 317 |
> |
if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then |
| 318 |
> |
biggestAtypeCutoff = atypeMaxCutoff(i) |
| 319 |
> |
endif |
| 320 |
> |
endif |
| 321 |
> |
enddo |
| 322 |
> |
|
| 323 |
> |
nGroupTypes = 0 |
| 324 |
> |
|
| 325 |
> |
istart = 1 |
| 326 |
> |
#ifdef IS_MPI |
| 327 |
> |
iend = nGroupsInRow |
| 328 |
> |
#else |
| 329 |
> |
iend = nGroups |
| 330 |
> |
#endif |
| 331 |
> |
|
| 332 |
> |
!! allocate the groupToGtype and gtypeMaxCutoff here. |
| 333 |
> |
if(.not.allocated(groupToGtype)) then |
| 334 |
> |
allocate(groupToGtype(iend)) |
| 335 |
> |
allocate(groupMaxCutoff(iend)) |
| 336 |
> |
allocate(gtypeMaxCutoff(iend)) |
| 337 |
> |
endif |
| 338 |
> |
!! first we do a single loop over the cutoff groups to find the |
| 339 |
> |
!! largest cutoff for any atypes present in this group. We also |
| 340 |
> |
!! create gtypes at this point. |
| 341 |
> |
|
| 342 |
> |
tol = 1.0d-6 |
| 343 |
> |
|
| 344 |
> |
do i = istart, iend |
| 345 |
> |
n_in_i = groupStartRow(i+1) - groupStartRow(i) |
| 346 |
> |
groupMaxCutoff(i) = 0.0_dp |
| 347 |
> |
do ia = groupStartRow(i), groupStartRow(i+1)-1 |
| 348 |
> |
atom1 = groupListRow(ia) |
| 349 |
> |
#ifdef IS_MPI |
| 350 |
> |
me_i = atid_row(atom1) |
| 351 |
> |
#else |
| 352 |
> |
me_i = atid(atom1) |
| 353 |
> |
#endif |
| 354 |
> |
if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then |
| 355 |
> |
groupMaxCutoff(i)=atypeMaxCutoff(me_i) |
| 356 |
> |
endif |
| 357 |
> |
enddo |
| 358 |
> |
if (nGroupTypes.eq.0) then |
| 359 |
> |
nGroupTypes = nGroupTypes + 1 |
| 360 |
> |
gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) |
| 361 |
> |
groupToGtype(i) = nGroupTypes |
| 362 |
> |
else |
| 363 |
> |
do g = 1, nGroupTypes |
| 364 |
> |
if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then |
| 365 |
> |
nGroupTypes = nGroupTypes + 1 |
| 366 |
> |
gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) |
| 367 |
> |
groupToGtype(i) = nGroupTypes |
| 368 |
> |
else |
| 369 |
> |
groupToGtype(i) = g |
| 370 |
> |
endif |
| 371 |
> |
enddo |
| 372 |
> |
endif |
| 373 |
> |
enddo |
| 374 |
> |
|
| 375 |
> |
!! allocate the gtypeCutoffMap here. |
| 376 |
> |
allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes)) |
| 377 |
> |
!! then we do a double loop over all the group TYPES to find the cutoff |
| 378 |
> |
!! map between groups of two types |
| 379 |
> |
|
| 380 |
> |
do i = 1, nGroupTypes |
| 381 |
> |
do j = 1, nGroupTypes |
| 382 |
> |
|
| 383 |
> |
select case(cutoffPolicy) |
| 384 |
> |
case(TRADITIONAL_CUTOFF_POLICY) |
| 385 |
> |
thisRcut = maxval(gtypeMaxCutoff) |
| 386 |
> |
case(MIX_CUTOFF_POLICY) |
| 387 |
> |
thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j)) |
| 388 |
> |
case(MAX_CUTOFF_POLICY) |
| 389 |
> |
thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j)) |
| 390 |
> |
case default |
| 391 |
> |
call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy") |
| 392 |
> |
return |
| 393 |
> |
end select |
| 394 |
> |
gtypeCutoffMap(i,j)%rcut = thisRcut |
| 395 |
> |
gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut |
| 396 |
> |
skin = defaultRlist - defaultRcut |
| 397 |
> |
gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2 |
| 398 |
|
|
| 399 |
+ |
enddo |
| 400 |
+ |
enddo |
| 401 |
+ |
|
| 402 |
+ |
haveGtypeCutoffMap = .true. |
| 403 |
+ |
|
| 404 |
+ |
end subroutine createGtypeCutoffMap |
| 405 |
+ |
|
| 406 |
+ |
subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy) |
| 407 |
+ |
real(kind=dp),intent(in) :: defRcut, defRsw, defRlist |
| 408 |
+ |
integer, intent(in) :: cutPolicy |
| 409 |
+ |
|
| 410 |
+ |
defaultRcut = defRcut |
| 411 |
+ |
defaultRsw = defRsw |
| 412 |
+ |
defaultRlist = defRlist |
| 413 |
+ |
cutoffPolicy = cutPolicy |
| 414 |
+ |
end subroutine setDefaultCutoffs |
| 415 |
+ |
|
| 416 |
+ |
subroutine setCutoffPolicy(cutPolicy) |
| 417 |
|
|
| 418 |
+ |
integer, intent(in) :: cutPolicy |
| 419 |
+ |
cutoffPolicy = cutPolicy |
| 420 |
+ |
call createGtypeCutoffMap() |
| 421 |
+ |
|
| 422 |
+ |
end subroutine setCutoffPolicy |
| 423 |
+ |
|
| 424 |
+ |
|
| 425 |
|
subroutine setSimVariables() |
| 426 |
|
SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms() |
| 342 |
– |
SIM_uses_LennardJones = SimUsesLennardJones() |
| 343 |
– |
SIM_uses_Electrostatics = SimUsesElectrostatics() |
| 344 |
– |
SIM_uses_Charges = SimUsesCharges() |
| 345 |
– |
SIM_uses_Dipoles = SimUsesDipoles() |
| 346 |
– |
SIM_uses_Sticky = SimUsesSticky() |
| 347 |
– |
SIM_uses_StickyPower = SimUsesStickyPower() |
| 348 |
– |
SIM_uses_GayBerne = SimUsesGayBerne() |
| 427 |
|
SIM_uses_EAM = SimUsesEAM() |
| 350 |
– |
SIM_uses_Shapes = SimUsesShapes() |
| 351 |
– |
SIM_uses_FLARB = SimUsesFLARB() |
| 428 |
|
SIM_uses_RF = SimUsesRF() |
| 429 |
|
SIM_requires_postpair_calc = SimRequiresPostpairCalc() |
| 430 |
|
SIM_requires_prepair_calc = SimRequiresPrepairCalc() |
| 442 |
|
|
| 443 |
|
error = 0 |
| 444 |
|
|
| 445 |
< |
if (.not. haveInteractionMap) then |
| 445 |
> |
if (.not. haveInteractionHash) then |
| 446 |
> |
myStatus = 0 |
| 447 |
> |
call createInteractionHash(myStatus) |
| 448 |
> |
if (myStatus .ne. 0) then |
| 449 |
> |
write(default_error, *) 'createInteractionHash failed in doForces!' |
| 450 |
> |
error = -1 |
| 451 |
> |
return |
| 452 |
> |
endif |
| 453 |
> |
endif |
| 454 |
|
|
| 455 |
< |
myStatus = 0 |
| 456 |
< |
|
| 457 |
< |
call createInteractionMap(myStatus) |
| 374 |
< |
|
| 455 |
> |
if (.not. haveGtypeCutoffMap) then |
| 456 |
> |
myStatus = 0 |
| 457 |
> |
call createGtypeCutoffMap(myStatus) |
| 458 |
|
if (myStatus .ne. 0) then |
| 459 |
< |
write(default_error, *) 'createInteractionMap failed in doForces!' |
| 459 |
> |
write(default_error, *) 'createGtypeCutoffMap failed in doForces!' |
| 460 |
|
error = -1 |
| 461 |
|
return |
| 462 |
|
endif |
| 466 |
|
call setSimVariables() |
| 467 |
|
endif |
| 468 |
|
|
| 469 |
< |
if (.not. haveRlist) then |
| 470 |
< |
write(default_error, *) 'rList has not been set in doForces!' |
| 471 |
< |
error = -1 |
| 472 |
< |
return |
| 473 |
< |
endif |
| 469 |
> |
! if (.not. haveRlist) then |
| 470 |
> |
! write(default_error, *) 'rList has not been set in doForces!' |
| 471 |
> |
! error = -1 |
| 472 |
> |
! return |
| 473 |
> |
! endif |
| 474 |
|
|
| 475 |
|
if (.not. haveNeighborList) then |
| 476 |
|
write(default_error, *) 'neighbor list has not been initialized in doForces!' |
| 495 |
|
end subroutine doReadyCheck |
| 496 |
|
|
| 497 |
|
|
| 498 |
< |
subroutine init_FF(use_RF_c, thisStat) |
| 498 |
> |
subroutine init_FF(use_RF, use_UW, use_DW, thisStat) |
| 499 |
|
|
| 500 |
< |
logical, intent(in) :: use_RF_c |
| 501 |
< |
|
| 500 |
> |
logical, intent(in) :: use_RF |
| 501 |
> |
logical, intent(in) :: use_UW |
| 502 |
> |
logical, intent(in) :: use_DW |
| 503 |
|
integer, intent(out) :: thisStat |
| 504 |
|
integer :: my_status, nMatches |
| 505 |
+ |
integer :: corrMethod |
| 506 |
|
integer, pointer :: MatchList(:) => null() |
| 507 |
|
real(kind=dp) :: rcut, rrf, rt, dielect |
| 508 |
|
|
| 510 |
|
thisStat = 0 |
| 511 |
|
|
| 512 |
|
!! Fortran's version of a cast: |
| 513 |
< |
FF_uses_RF = use_RF_c |
| 513 |
> |
FF_uses_RF = use_RF |
| 514 |
|
|
| 515 |
+ |
!! set the electrostatic correction method |
| 516 |
+ |
if (use_UW) then |
| 517 |
+ |
corrMethod = 1 |
| 518 |
+ |
elseif (use_DW) then |
| 519 |
+ |
corrMethod = 2 |
| 520 |
+ |
else |
| 521 |
+ |
corrMethod = 0 |
| 522 |
+ |
endif |
| 523 |
+ |
|
| 524 |
|
!! init_FF is called *after* all of the atom types have been |
| 525 |
|
!! defined in atype_module using the new_atype subroutine. |
| 526 |
|
!! |
| 528 |
|
!! interactions are used by the force field. |
| 529 |
|
|
| 530 |
|
FF_uses_DirectionalAtoms = .false. |
| 437 |
– |
FF_uses_LennardJones = .false. |
| 438 |
– |
FF_uses_Electrostatics = .false. |
| 439 |
– |
FF_uses_Charges = .false. |
| 531 |
|
FF_uses_Dipoles = .false. |
| 441 |
– |
FF_uses_Sticky = .false. |
| 442 |
– |
FF_uses_StickyPower = .false. |
| 532 |
|
FF_uses_GayBerne = .false. |
| 533 |
|
FF_uses_EAM = .false. |
| 445 |
– |
FF_uses_Shapes = .false. |
| 446 |
– |
FF_uses_FLARB = .false. |
| 534 |
|
|
| 535 |
|
call getMatchingElementList(atypes, "is_Directional", .true., & |
| 536 |
|
nMatches, MatchList) |
| 537 |
|
if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true. |
| 538 |
|
|
| 452 |
– |
call getMatchingElementList(atypes, "is_LennardJones", .true., & |
| 453 |
– |
nMatches, MatchList) |
| 454 |
– |
if (nMatches .gt. 0) FF_uses_LennardJones = .true. |
| 455 |
– |
|
| 456 |
– |
call getMatchingElementList(atypes, "is_Electrostatic", .true., & |
| 457 |
– |
nMatches, MatchList) |
| 458 |
– |
if (nMatches .gt. 0) then |
| 459 |
– |
FF_uses_Electrostatics = .true. |
| 460 |
– |
endif |
| 461 |
– |
|
| 462 |
– |
call getMatchingElementList(atypes, "is_Charge", .true., & |
| 463 |
– |
nMatches, MatchList) |
| 464 |
– |
if (nMatches .gt. 0) then |
| 465 |
– |
FF_uses_Charges = .true. |
| 466 |
– |
FF_uses_Electrostatics = .true. |
| 467 |
– |
endif |
| 468 |
– |
|
| 539 |
|
call getMatchingElementList(atypes, "is_Dipole", .true., & |
| 540 |
|
nMatches, MatchList) |
| 541 |
< |
if (nMatches .gt. 0) then |
| 472 |
< |
FF_uses_Dipoles = .true. |
| 473 |
< |
FF_uses_Electrostatics = .true. |
| 474 |
< |
FF_uses_DirectionalAtoms = .true. |
| 475 |
< |
endif |
| 476 |
< |
|
| 477 |
< |
call getMatchingElementList(atypes, "is_Quadrupole", .true., & |
| 478 |
< |
nMatches, MatchList) |
| 479 |
< |
if (nMatches .gt. 0) then |
| 480 |
< |
FF_uses_Quadrupoles = .true. |
| 481 |
< |
FF_uses_Electrostatics = .true. |
| 482 |
< |
FF_uses_DirectionalAtoms = .true. |
| 483 |
< |
endif |
| 484 |
< |
|
| 485 |
< |
call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, & |
| 486 |
< |
MatchList) |
| 487 |
< |
if (nMatches .gt. 0) then |
| 488 |
< |
FF_uses_Sticky = .true. |
| 489 |
< |
FF_uses_DirectionalAtoms = .true. |
| 490 |
< |
endif |
| 491 |
< |
|
| 492 |
< |
call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, & |
| 493 |
< |
MatchList) |
| 494 |
< |
if (nMatches .gt. 0) then |
| 495 |
< |
FF_uses_StickyPower = .true. |
| 496 |
< |
FF_uses_DirectionalAtoms = .true. |
| 497 |
< |
endif |
| 541 |
> |
if (nMatches .gt. 0) FF_uses_Dipoles = .true. |
| 542 |
|
|
| 543 |
|
call getMatchingElementList(atypes, "is_GayBerne", .true., & |
| 544 |
|
nMatches, MatchList) |
| 545 |
< |
if (nMatches .gt. 0) then |
| 502 |
< |
FF_uses_GayBerne = .true. |
| 503 |
< |
FF_uses_DirectionalAtoms = .true. |
| 504 |
< |
endif |
| 545 |
> |
if (nMatches .gt. 0) FF_uses_GayBerne = .true. |
| 546 |
|
|
| 547 |
|
call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList) |
| 548 |
|
if (nMatches .gt. 0) FF_uses_EAM = .true. |
| 549 |
|
|
| 509 |
– |
call getMatchingElementList(atypes, "is_Shape", .true., & |
| 510 |
– |
nMatches, MatchList) |
| 511 |
– |
if (nMatches .gt. 0) then |
| 512 |
– |
FF_uses_Shapes = .true. |
| 513 |
– |
FF_uses_DirectionalAtoms = .true. |
| 514 |
– |
endif |
| 550 |
|
|
| 516 |
– |
call getMatchingElementList(atypes, "is_FLARB", .true., & |
| 517 |
– |
nMatches, MatchList) |
| 518 |
– |
if (nMatches .gt. 0) FF_uses_FLARB = .true. |
| 519 |
– |
|
| 520 |
– |
!! Assume sanity (for the sake of argument) |
| 551 |
|
haveSaneForceField = .true. |
| 552 |
|
|
| 553 |
|
!! check to make sure the FF_uses_RF setting makes sense |
| 554 |
|
|
| 555 |
< |
if (FF_uses_dipoles) then |
| 555 |
> |
if (FF_uses_Dipoles) then |
| 556 |
|
if (FF_uses_RF) then |
| 557 |
|
dielect = getDielect() |
| 558 |
|
call initialize_rf(dielect) |
| 566 |
|
endif |
| 567 |
|
endif |
| 568 |
|
|
| 539 |
– |
!sticky module does not contain check_sticky_FF anymore |
| 540 |
– |
!if (FF_uses_sticky) then |
| 541 |
– |
! call check_sticky_FF(my_status) |
| 542 |
– |
! if (my_status /= 0) then |
| 543 |
– |
! thisStat = -1 |
| 544 |
– |
! haveSaneForceField = .false. |
| 545 |
– |
! return |
| 546 |
– |
! end if |
| 547 |
– |
!endif |
| 548 |
– |
|
| 569 |
|
if (FF_uses_EAM) then |
| 570 |
|
call init_EAM_FF(my_status) |
| 571 |
|
if (my_status /= 0) then |
| 585 |
|
endif |
| 586 |
|
endif |
| 587 |
|
|
| 568 |
– |
if (FF_uses_GayBerne .and. FF_uses_LennardJones) then |
| 569 |
– |
endif |
| 570 |
– |
|
| 588 |
|
if (.not. haveNeighborList) then |
| 589 |
|
!! Create neighbor lists |
| 590 |
|
call expandNeighborList(nLocal, my_status) |
| 648 |
|
integer :: localError |
| 649 |
|
integer :: propPack_i, propPack_j |
| 650 |
|
integer :: loopStart, loopEnd, loop |
| 651 |
< |
integer :: iMap |
| 651 |
> |
integer :: iHash |
| 652 |
|
real(kind=dp) :: listSkin = 1.0 |
| 653 |
|
|
| 654 |
|
!! initialize local variables |
| 767 |
|
endif |
| 768 |
|
|
| 769 |
|
#ifdef IS_MPI |
| 770 |
+ |
me_j = atid_col(j) |
| 771 |
|
call get_interatomic_vector(q_group_Row(:,i), & |
| 772 |
|
q_group_Col(:,j), d_grp, rgrpsq) |
| 773 |
|
#else |
| 774 |
+ |
me_j = atid(j) |
| 775 |
|
call get_interatomic_vector(q_group(:,i), & |
| 776 |
|
q_group(:,j), d_grp, rgrpsq) |
| 777 |
|
#endif |
| 778 |
|
|
| 779 |
< |
if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then |
| 779 |
> |
if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then |
| 780 |
|
if (update_nlist) then |
| 781 |
|
nlist = nlist + 1 |
| 782 |
|
|
| 994 |
|
#else |
| 995 |
|
me_i = atid(i) |
| 996 |
|
#endif |
| 997 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 997 |
> |
iHash = InteractionHash(me_i,me_j) |
| 998 |
|
|
| 999 |
< |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 999 |
> |
if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1000 |
|
|
| 1001 |
|
mu_i = getDipoleMoment(me_i) |
| 1002 |
|
|
| 1063 |
|
real ( kind = dp ) :: ebalance |
| 1064 |
|
integer :: me_i, me_j |
| 1065 |
|
|
| 1066 |
< |
integer :: iMap |
| 1066 |
> |
integer :: iHash |
| 1067 |
|
|
| 1068 |
|
r = sqrt(rijsq) |
| 1069 |
|
vpair = 0.0d0 |
| 1077 |
|
me_j = atid(j) |
| 1078 |
|
#endif |
| 1079 |
|
|
| 1080 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 1080 |
> |
iHash = InteractionHash(me_i, me_j) |
| 1081 |
|
|
| 1082 |
< |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
| 1082 |
> |
if ( iand(iHash, LJ_PAIR).ne.0 ) then |
| 1083 |
|
call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) |
| 1084 |
|
endif |
| 1085 |
|
|
| 1086 |
< |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1086 |
> |
if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1087 |
|
call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1088 |
< |
pot, eFrame, f, t, do_pot) |
| 1088 |
> |
pot, eFrame, f, t, do_pot, corrMethod) |
| 1089 |
|
|
| 1090 |
|
if (FF_uses_RF .and. SIM_uses_RF) then |
| 1091 |
|
|
| 1096 |
|
|
| 1097 |
|
endif |
| 1098 |
|
|
| 1099 |
< |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
| 1099 |
> |
if ( iand(iHash, STICKY_PAIR).ne.0 ) then |
| 1100 |
|
call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1101 |
|
pot, A, f, t, do_pot) |
| 1102 |
|
endif |
| 1103 |
|
|
| 1104 |
< |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
| 1104 |
> |
if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then |
| 1105 |
|
call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1106 |
|
pot, A, f, t, do_pot) |
| 1107 |
|
endif |
| 1108 |
|
|
| 1109 |
< |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
| 1109 |
> |
if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then |
| 1110 |
|
call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1111 |
|
pot, A, f, t, do_pot) |
| 1112 |
|
endif |
| 1113 |
|
|
| 1114 |
< |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
| 1114 |
> |
if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then |
| 1115 |
|
! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1116 |
|
! pot, A, f, t, do_pot) |
| 1117 |
|
endif |
| 1118 |
|
|
| 1119 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 1119 |
> |
if ( iand(iHash, EAM_PAIR).ne.0 ) then |
| 1120 |
|
call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, & |
| 1121 |
|
do_pot) |
| 1122 |
|
endif |
| 1123 |
|
|
| 1124 |
< |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
| 1124 |
> |
if ( iand(iHash, SHAPE_PAIR).ne.0 ) then |
| 1125 |
|
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1126 |
|
pot, A, f, t, do_pot) |
| 1127 |
|
endif |
| 1128 |
|
|
| 1129 |
< |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
| 1129 |
> |
if ( iand(iHash, SHAPE_LJ).ne.0 ) then |
| 1130 |
|
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1131 |
|
pot, A, f, t, do_pot) |
| 1132 |
|
endif |
| 1148 |
|
real ( kind = dp ) :: r, rc |
| 1149 |
|
real ( kind = dp ), intent(inout) :: d(3), dc(3) |
| 1150 |
|
|
| 1151 |
< |
integer :: me_i, me_j, iMap |
| 1151 |
> |
integer :: me_i, me_j, iHash |
| 1152 |
|
|
| 1153 |
|
#ifdef IS_MPI |
| 1154 |
|
me_i = atid_row(i) |
| 1158 |
|
me_j = atid(j) |
| 1159 |
|
#endif |
| 1160 |
|
|
| 1161 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 1161 |
> |
iHash = InteractionHash(me_i, me_j) |
| 1162 |
|
|
| 1163 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 1163 |
> |
if ( iand(iHash, EAM_PAIR).ne.0 ) then |
| 1164 |
|
call calc_EAM_prepair_rho(i, j, d, r, rijsq ) |
| 1165 |
|
endif |
| 1166 |
|
|
| 1357 |
|
|
| 1358 |
|
function FF_UsesDirectionalAtoms() result(doesit) |
| 1359 |
|
logical :: doesit |
| 1360 |
< |
doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. & |
| 1342 |
< |
FF_uses_Quadrupoles .or. FF_uses_Sticky .or. & |
| 1343 |
< |
FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes |
| 1360 |
> |
doesit = FF_uses_DirectionalAtoms |
| 1361 |
|
end function FF_UsesDirectionalAtoms |
| 1362 |
|
|
| 1363 |
|
function FF_RequiresPrepairCalc() result(doesit) |
