| 45 |
|
|
| 46 |
|
!! @author Charles F. Vardeman II |
| 47 |
|
!! @author Matthew Meineke |
| 48 |
< |
!! @version $Id: doForces.F90,v 1.14 2005-04-21 14:12:19 chrisfen Exp $, $Date: 2005-04-21 14:12:19 $, $Name: not supported by cvs2svn $, $Revision: 1.14 $ |
| 48 |
> |
!! @version $Id: doForces.F90,v 1.25 2005-07-29 17:34:06 tim Exp $, $Date: 2005-07-29 17:34:06 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $ |
| 49 |
|
|
| 50 |
|
|
| 51 |
|
module doForces |
| 73 |
|
|
| 74 |
|
#define __FORTRAN90 |
| 75 |
|
#include "UseTheForce/fSwitchingFunction.h" |
| 76 |
+ |
#include "UseTheForce/DarkSide/fInteractionMap.h" |
| 77 |
|
|
| 78 |
|
INTEGER, PARAMETER:: PREPAIR_LOOP = 1 |
| 79 |
|
INTEGER, PARAMETER:: PAIR_LOOP = 2 |
| 81 |
|
logical, save :: haveRlist = .false. |
| 82 |
|
logical, save :: haveNeighborList = .false. |
| 83 |
|
logical, save :: haveSIMvariables = .false. |
| 83 |
– |
logical, save :: havePropertyMap = .false. |
| 84 |
|
logical, save :: haveSaneForceField = .false. |
| 85 |
+ |
logical, save :: haveInteractionMap = .false. |
| 86 |
|
|
| 87 |
|
logical, save :: FF_uses_DirectionalAtoms |
| 88 |
|
logical, save :: FF_uses_LennardJones |
| 90 |
|
logical, save :: FF_uses_Charges |
| 91 |
|
logical, save :: FF_uses_Dipoles |
| 92 |
|
logical, save :: FF_uses_Quadrupoles |
| 93 |
< |
logical, save :: FF_uses_sticky |
| 93 |
> |
logical, save :: FF_uses_Sticky |
| 94 |
> |
logical, save :: FF_uses_StickyPower |
| 95 |
|
logical, save :: FF_uses_GayBerne |
| 96 |
|
logical, save :: FF_uses_EAM |
| 97 |
|
logical, save :: FF_uses_Shapes |
| 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 |
| 110 |
|
logical, save :: SIM_uses_EAM |
| 111 |
|
logical, save :: SIM_uses_Shapes |
| 116 |
|
logical, save :: SIM_uses_PBC |
| 117 |
|
logical, save :: SIM_uses_molecular_cutoffs |
| 118 |
|
|
| 119 |
< |
real(kind=dp), save :: rlist, rlistsq |
| 119 |
> |
!!!GO AWAY--------- |
| 120 |
> |
!!!!!real(kind=dp), save :: rlist, rlistsq |
| 121 |
|
|
| 122 |
|
public :: init_FF |
| 123 |
|
public :: do_force_loop |
| 124 |
< |
public :: setRlistDF |
| 124 |
> |
! public :: setRlistDF |
| 125 |
> |
!public :: addInteraction |
| 126 |
> |
!public :: setInteractionHash |
| 127 |
> |
!public :: getInteractionHash |
| 128 |
> |
public :: createInteractionMap |
| 129 |
> |
public :: createRcuts |
| 130 |
|
|
| 131 |
|
#ifdef PROFILE |
| 132 |
|
public :: getforcetime |
| 135 |
|
integer :: nLoops |
| 136 |
|
#endif |
| 137 |
|
|
| 138 |
< |
type :: Properties |
| 139 |
< |
logical :: is_Directional = .false. |
| 140 |
< |
logical :: is_LennardJones = .false. |
| 141 |
< |
logical :: is_Electrostatic = .false. |
| 142 |
< |
logical :: is_Charge = .false. |
| 143 |
< |
logical :: is_Dipole = .false. |
| 144 |
< |
logical :: is_Quadrupole = .false. |
| 145 |
< |
logical :: is_Sticky = .false. |
| 137 |
< |
logical :: is_GayBerne = .false. |
| 138 |
< |
logical :: is_EAM = .false. |
| 139 |
< |
logical :: is_Shape = .false. |
| 140 |
< |
logical :: is_FLARB = .false. |
| 141 |
< |
end type Properties |
| 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 |
| 143 |
> |
|
| 144 |
> |
type(Interaction), dimension(:,:),allocatable :: InteractionMap |
| 145 |
> |
|
| 146 |
|
|
| 147 |
< |
type(Properties), dimension(:),allocatable :: PropertyMap |
| 144 |
< |
|
| 147 |
> |
|
| 148 |
|
contains |
| 149 |
|
|
| 150 |
< |
subroutine setRlistDF( this_rlist ) |
| 151 |
< |
|
| 149 |
< |
real(kind=dp) :: this_rlist |
| 150 |
< |
|
| 151 |
< |
rlist = this_rlist |
| 152 |
< |
rlistsq = rlist * rlist |
| 153 |
< |
|
| 154 |
< |
haveRlist = .true. |
| 155 |
< |
|
| 156 |
< |
end subroutine setRlistDF |
| 157 |
< |
|
| 158 |
< |
subroutine createPropertyMap(status) |
| 150 |
> |
|
| 151 |
> |
subroutine createInteractionMap(status) |
| 152 |
|
integer :: nAtypes |
| 153 |
< |
integer :: status |
| 153 |
> |
integer, intent(out) :: status |
| 154 |
|
integer :: i |
| 155 |
< |
logical :: thisProperty |
| 156 |
< |
real (kind=DP) :: thisDPproperty |
| 157 |
< |
|
| 155 |
> |
integer :: j |
| 156 |
> |
integer :: ihash |
| 157 |
> |
real(kind=dp) :: myRcut |
| 158 |
> |
!! Test Types |
| 159 |
> |
logical :: i_is_LJ |
| 160 |
> |
logical :: i_is_Elect |
| 161 |
> |
logical :: i_is_Sticky |
| 162 |
> |
logical :: i_is_StickyP |
| 163 |
> |
logical :: i_is_GB |
| 164 |
> |
logical :: i_is_EAM |
| 165 |
> |
logical :: i_is_Shape |
| 166 |
> |
logical :: j_is_LJ |
| 167 |
> |
logical :: j_is_Elect |
| 168 |
> |
logical :: j_is_Sticky |
| 169 |
> |
logical :: j_is_StickyP |
| 170 |
> |
logical :: j_is_GB |
| 171 |
> |
logical :: j_is_EAM |
| 172 |
> |
logical :: j_is_Shape |
| 173 |
> |
|
| 174 |
|
status = 0 |
| 175 |
< |
|
| 175 |
> |
|
| 176 |
> |
if (.not. associated(atypes)) then |
| 177 |
> |
call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!") |
| 178 |
> |
status = -1 |
| 179 |
> |
return |
| 180 |
> |
endif |
| 181 |
> |
|
| 182 |
|
nAtypes = getSize(atypes) |
| 183 |
< |
|
| 183 |
> |
|
| 184 |
|
if (nAtypes == 0) then |
| 185 |
|
status = -1 |
| 186 |
|
return |
| 187 |
|
end if |
| 188 |
|
|
| 189 |
< |
if (.not. allocated(PropertyMap)) then |
| 190 |
< |
allocate(PropertyMap(nAtypes)) |
| 189 |
> |
if (.not. allocated(InteractionMap)) then |
| 190 |
> |
allocate(InteractionMap(nAtypes,nAtypes)) |
| 191 |
|
endif |
| 192 |
< |
|
| 192 |
> |
|
| 193 |
|
do i = 1, nAtypes |
| 194 |
< |
call getElementProperty(atypes, i, "is_Directional", thisProperty) |
| 195 |
< |
PropertyMap(i)%is_Directional = thisProperty |
| 194 |
> |
call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) |
| 195 |
> |
call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect) |
| 196 |
> |
call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky) |
| 197 |
> |
call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP) |
| 198 |
> |
call getElementProperty(atypes, i, "is_GayBerne", i_is_GB) |
| 199 |
> |
call getElementProperty(atypes, i, "is_EAM", i_is_EAM) |
| 200 |
> |
call getElementProperty(atypes, i, "is_Shape", i_is_Shape) |
| 201 |
|
|
| 202 |
< |
call getElementProperty(atypes, i, "is_LennardJones", thisProperty) |
| 183 |
< |
PropertyMap(i)%is_LennardJones = thisProperty |
| 202 |
> |
do j = i, nAtypes |
| 203 |
|
|
| 204 |
< |
call getElementProperty(atypes, i, "is_Electrostatic", thisProperty) |
| 205 |
< |
PropertyMap(i)%is_Electrostatic = thisProperty |
| 204 |
> |
iHash = 0 |
| 205 |
> |
myRcut = 0.0_dp |
| 206 |
|
|
| 207 |
< |
call getElementProperty(atypes, i, "is_Charge", thisProperty) |
| 208 |
< |
PropertyMap(i)%is_Charge = thisProperty |
| 209 |
< |
|
| 210 |
< |
call getElementProperty(atypes, i, "is_Dipole", thisProperty) |
| 211 |
< |
PropertyMap(i)%is_Dipole = thisProperty |
| 207 |
> |
call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ) |
| 208 |
> |
call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect) |
| 209 |
> |
call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky) |
| 210 |
> |
call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP) |
| 211 |
> |
call getElementProperty(atypes, j, "is_GayBerne", j_is_GB) |
| 212 |
> |
call getElementProperty(atypes, j, "is_EAM", j_is_EAM) |
| 213 |
> |
call getElementProperty(atypes, j, "is_Shape", j_is_Shape) |
| 214 |
|
|
| 215 |
< |
call getElementProperty(atypes, i, "is_Quadrupole", thisProperty) |
| 216 |
< |
PropertyMap(i)%is_Quadrupole = thisProperty |
| 215 |
> |
if (i_is_LJ .and. j_is_LJ) then |
| 216 |
> |
iHash = ior(iHash, LJ_PAIR) |
| 217 |
> |
endif |
| 218 |
> |
|
| 219 |
> |
if (i_is_Elect .and. j_is_Elect) then |
| 220 |
> |
iHash = ior(iHash, ELECTROSTATIC_PAIR) |
| 221 |
> |
endif |
| 222 |
> |
|
| 223 |
> |
if (i_is_Sticky .and. j_is_Sticky) then |
| 224 |
> |
iHash = ior(iHash, STICKY_PAIR) |
| 225 |
> |
endif |
| 226 |
|
|
| 227 |
< |
call getElementProperty(atypes, i, "is_Sticky", thisProperty) |
| 228 |
< |
PropertyMap(i)%is_Sticky = thisProperty |
| 227 |
> |
if (i_is_StickyP .and. j_is_StickyP) then |
| 228 |
> |
iHash = ior(iHash, STICKYPOWER_PAIR) |
| 229 |
> |
endif |
| 230 |
|
|
| 231 |
< |
call getElementProperty(atypes, i, "is_GayBerne", thisProperty) |
| 232 |
< |
PropertyMap(i)%is_GayBerne = thisProperty |
| 231 |
> |
if (i_is_EAM .and. j_is_EAM) then |
| 232 |
> |
iHash = ior(iHash, EAM_PAIR) |
| 233 |
> |
endif |
| 234 |
|
|
| 235 |
< |
call getElementProperty(atypes, i, "is_EAM", thisProperty) |
| 236 |
< |
PropertyMap(i)%is_EAM = thisProperty |
| 235 |
> |
if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR) |
| 236 |
> |
if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ) |
| 237 |
> |
if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ) |
| 238 |
|
|
| 239 |
< |
call getElementProperty(atypes, i, "is_Shape", thisProperty) |
| 240 |
< |
PropertyMap(i)%is_Shape = thisProperty |
| 239 |
> |
if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR) |
| 240 |
> |
if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ) |
| 241 |
> |
if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ) |
| 242 |
|
|
| 243 |
< |
call getElementProperty(atypes, i, "is_FLARB", thisProperty) |
| 244 |
< |
PropertyMap(i)%is_FLARB = thisProperty |
| 243 |
> |
|
| 244 |
> |
InteractionMap(i,j)%InteractionHash = iHash |
| 245 |
> |
InteractionMap(j,i)%InteractionHash = iHash |
| 246 |
> |
|
| 247 |
> |
end do |
| 248 |
> |
|
| 249 |
|
end do |
| 250 |
|
|
| 251 |
< |
havePropertyMap = .true. |
| 251 |
> |
haveInteractionMap = .true. |
| 252 |
> |
end subroutine createInteractionMap |
| 253 |
|
|
| 254 |
< |
end subroutine createPropertyMap |
| 254 |
> |
! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff. |
| 255 |
> |
subroutine createRcuts(defaultRList,stat) |
| 256 |
> |
real(kind=dp), intent(in), optional :: defaultRList |
| 257 |
> |
integer :: iMap |
| 258 |
> |
integer :: map_i,map_j |
| 259 |
> |
real(kind=dp) :: thisRCut = 0.0_dp |
| 260 |
> |
real(kind=dp) :: actualCutoff = 0.0_dp |
| 261 |
> |
integer, intent(out) :: stat |
| 262 |
> |
integer :: nAtypes |
| 263 |
> |
integer :: myStatus |
| 264 |
> |
|
| 265 |
> |
stat = 0 |
| 266 |
> |
if (.not. haveInteractionMap) then |
| 267 |
> |
|
| 268 |
> |
call createInteractionMap(myStatus) |
| 269 |
> |
|
| 270 |
> |
if (myStatus .ne. 0) then |
| 271 |
> |
write(default_error, *) 'createInteractionMap failed in doForces!' |
| 272 |
> |
stat = -1 |
| 273 |
> |
return |
| 274 |
> |
endif |
| 275 |
> |
endif |
| 276 |
> |
|
| 277 |
> |
|
| 278 |
> |
nAtypes = getSize(atypes) |
| 279 |
> |
!! If we pass a default rcut, set all atypes to that cutoff distance |
| 280 |
> |
if(present(defaultRList)) then |
| 281 |
> |
InteractionMap(:,:)%rList = defaultRList |
| 282 |
> |
InteractionMap(:,:)%rListSq = defaultRList*defaultRList |
| 283 |
> |
haveRlist = .true. |
| 284 |
> |
return |
| 285 |
> |
end if |
| 286 |
> |
|
| 287 |
> |
do map_i = 1,nAtypes |
| 288 |
> |
do map_j = map_i,nAtypes |
| 289 |
> |
iMap = InteractionMap(map_i, map_j)%InteractionHash |
| 290 |
> |
|
| 291 |
> |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
| 292 |
> |
! thisRCut = getLJCutOff(map_i,map_j) |
| 293 |
> |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 294 |
> |
endif |
| 295 |
> |
|
| 296 |
> |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 297 |
> |
! thisRCut = getElectrostaticCutOff(map_i,map_j) |
| 298 |
> |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 299 |
> |
endif |
| 300 |
> |
|
| 301 |
> |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
| 302 |
> |
! thisRCut = getStickyCutOff(map_i,map_j) |
| 303 |
> |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 304 |
> |
endif |
| 305 |
> |
|
| 306 |
> |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
| 307 |
> |
! thisRCut = getStickyPowerCutOff(map_i,map_j) |
| 308 |
> |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 309 |
> |
endif |
| 310 |
> |
|
| 311 |
> |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
| 312 |
> |
! thisRCut = getGayberneCutOff(map_i,map_j) |
| 313 |
> |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 314 |
> |
endif |
| 315 |
> |
|
| 316 |
> |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
| 317 |
> |
! thisRCut = getGaybrneLJCutOff(map_i,map_j) |
| 318 |
> |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 319 |
> |
endif |
| 320 |
> |
|
| 321 |
> |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 322 |
> |
! thisRCut = getEAMCutOff(map_i,map_j) |
| 323 |
> |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 324 |
> |
endif |
| 325 |
> |
|
| 326 |
> |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
| 327 |
> |
! thisRCut = getShapeCutOff(map_i,map_j) |
| 328 |
> |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 329 |
> |
endif |
| 330 |
> |
|
| 331 |
> |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
| 332 |
> |
! thisRCut = getShapeLJCutOff(map_i,map_j) |
| 333 |
> |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 334 |
> |
endif |
| 335 |
> |
InteractionMap(map_i, map_j)%rList = actualCutoff |
| 336 |
> |
InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff |
| 337 |
> |
end do |
| 338 |
> |
end do |
| 339 |
> |
haveRlist = .true. |
| 340 |
> |
end subroutine createRcuts |
| 341 |
|
|
| 342 |
+ |
|
| 343 |
+ |
!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF |
| 344 |
+ |
!!$ subroutine setRlistDF( this_rlist ) |
| 345 |
+ |
!!$ |
| 346 |
+ |
!!$ real(kind=dp) :: this_rlist |
| 347 |
+ |
!!$ |
| 348 |
+ |
!!$ rlist = this_rlist |
| 349 |
+ |
!!$ rlistsq = rlist * rlist |
| 350 |
+ |
!!$ |
| 351 |
+ |
!!$ haveRlist = .true. |
| 352 |
+ |
!!$ |
| 353 |
+ |
!!$ end subroutine setRlistDF |
| 354 |
+ |
|
| 355 |
+ |
|
| 356 |
|
subroutine setSimVariables() |
| 357 |
|
SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms() |
| 358 |
|
SIM_uses_LennardJones = SimUsesLennardJones() |
| 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() |
| 365 |
|
SIM_uses_EAM = SimUsesEAM() |
| 366 |
|
SIM_uses_Shapes = SimUsesShapes() |
| 382 |
|
|
| 383 |
|
error = 0 |
| 384 |
|
|
| 385 |
< |
if (.not. havePropertyMap) then |
| 385 |
> |
if (.not. haveInteractionMap) then |
| 386 |
|
|
| 387 |
|
myStatus = 0 |
| 388 |
|
|
| 389 |
< |
call createPropertyMap(myStatus) |
| 389 |
> |
call createInteractionMap(myStatus) |
| 390 |
|
|
| 391 |
|
if (myStatus .ne. 0) then |
| 392 |
< |
write(default_error, *) 'createPropertyMap failed in doForces!' |
| 392 |
> |
write(default_error, *) 'createInteractionMap failed in doForces!' |
| 393 |
|
error = -1 |
| 394 |
|
return |
| 395 |
|
endif |
| 455 |
|
FF_uses_Charges = .false. |
| 456 |
|
FF_uses_Dipoles = .false. |
| 457 |
|
FF_uses_Sticky = .false. |
| 458 |
+ |
FF_uses_StickyPower = .false. |
| 459 |
|
FF_uses_GayBerne = .false. |
| 460 |
|
FF_uses_EAM = .false. |
| 461 |
|
FF_uses_Shapes = .false. |
| 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 |
| 514 |
+ |
|
| 515 |
|
call getMatchingElementList(atypes, "is_GayBerne", .true., & |
| 516 |
|
nMatches, MatchList) |
| 517 |
|
if (nMatches .gt. 0) then |
| 647 |
|
integer :: localError |
| 648 |
|
integer :: propPack_i, propPack_j |
| 649 |
|
integer :: loopStart, loopEnd, loop |
| 650 |
< |
|
| 650 |
> |
integer :: iMap |
| 651 |
|
real(kind=dp) :: listSkin = 1.0 |
| 652 |
|
|
| 653 |
|
!! initialize local variables |
| 738 |
|
iend = nGroups - 1 |
| 739 |
|
#endif |
| 740 |
|
outer: do i = istart, iend |
| 741 |
+ |
|
| 742 |
+ |
#ifdef IS_MPI |
| 743 |
+ |
me_i = atid_row(i) |
| 744 |
+ |
#else |
| 745 |
+ |
me_i = atid(i) |
| 746 |
+ |
#endif |
| 747 |
|
|
| 748 |
|
if (update_nlist) point(i) = nlist + 1 |
| 749 |
|
|
| 772 |
|
endif |
| 773 |
|
|
| 774 |
|
#ifdef IS_MPI |
| 775 |
+ |
me_j = atid_col(j) |
| 776 |
|
call get_interatomic_vector(q_group_Row(:,i), & |
| 777 |
|
q_group_Col(:,j), d_grp, rgrpsq) |
| 778 |
|
#else |
| 779 |
+ |
me_j = atid(j) |
| 780 |
|
call get_interatomic_vector(q_group(:,i), & |
| 781 |
|
q_group(:,j), d_grp, rgrpsq) |
| 782 |
|
#endif |
| 783 |
|
|
| 784 |
< |
if (rgrpsq < rlistsq) then |
| 784 |
> |
if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then |
| 785 |
|
if (update_nlist) then |
| 786 |
|
nlist = nlist + 1 |
| 787 |
|
|
| 999 |
|
#else |
| 1000 |
|
me_i = atid(i) |
| 1001 |
|
#endif |
| 1002 |
+ |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 1003 |
+ |
|
| 1004 |
+ |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1005 |
|
|
| 847 |
– |
if (PropertyMap(me_i)%is_Dipole) then |
| 848 |
– |
|
| 1006 |
|
mu_i = getDipoleMoment(me_i) |
| 1007 |
|
|
| 1008 |
|
!! The reaction field needs to include a self contribution |
| 1065 |
|
real ( kind = dp ), intent(inout) :: rijsq |
| 1066 |
|
real ( kind = dp ) :: r |
| 1067 |
|
real ( kind = dp ), intent(inout) :: d(3) |
| 1068 |
+ |
real ( kind = dp ) :: ebalance |
| 1069 |
|
integer :: me_i, me_j |
| 1070 |
|
|
| 1071 |
+ |
integer :: iMap |
| 1072 |
+ |
|
| 1073 |
|
r = sqrt(rijsq) |
| 1074 |
|
vpair = 0.0d0 |
| 1075 |
|
fpair(1:3) = 0.0d0 |
| 1082 |
|
me_j = atid(j) |
| 1083 |
|
#endif |
| 1084 |
|
|
| 1085 |
< |
! write(*,*) i, j, me_i, me_j |
| 1085 |
> |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 1086 |
|
|
| 1087 |
< |
if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then |
| 1088 |
< |
|
| 929 |
< |
if ( PropertyMap(me_i)%is_LennardJones .and. & |
| 930 |
< |
PropertyMap(me_j)%is_LennardJones ) then |
| 931 |
< |
call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) |
| 932 |
< |
endif |
| 933 |
< |
|
| 1087 |
> |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
| 1088 |
> |
call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) |
| 1089 |
|
endif |
| 1090 |
|
|
| 1091 |
< |
if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then |
| 1091 |
> |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1092 |
> |
call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1093 |
> |
pot, eFrame, f, t, do_pot) |
| 1094 |
|
|
| 1095 |
< |
if (PropertyMap(me_i)%is_Electrostatic .and. & |
| 939 |
< |
PropertyMap(me_j)%is_Electrostatic) then |
| 940 |
< |
call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 941 |
< |
pot, eFrame, f, t, do_pot) |
| 942 |
< |
endif |
| 1095 |
> |
if (FF_uses_RF .and. SIM_uses_RF) then |
| 1096 |
|
|
| 1097 |
< |
if (FF_uses_dipoles .and. SIM_uses_dipoles) then |
| 1098 |
< |
if ( PropertyMap(me_i)%is_Dipole .and. & |
| 1099 |
< |
PropertyMap(me_j)%is_Dipole) then |
| 947 |
< |
if (FF_uses_RF .and. SIM_uses_RF) then |
| 948 |
< |
call accumulate_rf(i, j, r, eFrame, sw) |
| 949 |
< |
call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair) |
| 950 |
< |
endif |
| 951 |
< |
endif |
| 1097 |
> |
! CHECK ME (RF needs to know about all electrostatic types) |
| 1098 |
> |
call accumulate_rf(i, j, r, eFrame, sw) |
| 1099 |
> |
call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair) |
| 1100 |
|
endif |
| 1101 |
+ |
|
| 1102 |
|
endif |
| 1103 |
|
|
| 1104 |
+ |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
| 1105 |
+ |
call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1106 |
+ |
pot, A, f, t, do_pot) |
| 1107 |
+ |
endif |
| 1108 |
|
|
| 1109 |
< |
if (FF_uses_Sticky .and. SIM_uses_sticky) then |
| 1110 |
< |
|
| 1111 |
< |
if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then |
| 959 |
< |
call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 960 |
< |
pot, A, f, t, do_pot) |
| 961 |
< |
endif |
| 962 |
< |
|
| 1109 |
> |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
| 1110 |
> |
call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1111 |
> |
pot, A, f, t, do_pot) |
| 1112 |
|
endif |
| 1113 |
|
|
| 1114 |
< |
|
| 1115 |
< |
if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then |
| 1116 |
< |
|
| 968 |
< |
if ( PropertyMap(me_i)%is_GayBerne .and. & |
| 969 |
< |
PropertyMap(me_j)%is_GayBerne) then |
| 970 |
< |
call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 971 |
< |
pot, A, f, t, do_pot) |
| 972 |
< |
endif |
| 973 |
< |
|
| 1114 |
> |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
| 1115 |
> |
call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1116 |
> |
pot, A, f, t, do_pot) |
| 1117 |
|
endif |
| 1118 |
+ |
|
| 1119 |
+ |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
| 1120 |
+ |
! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1121 |
+ |
! pot, A, f, t, do_pot) |
| 1122 |
+ |
endif |
| 1123 |
|
|
| 1124 |
< |
if (FF_uses_EAM .and. SIM_uses_EAM) then |
| 1125 |
< |
|
| 1126 |
< |
if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then |
| 979 |
< |
call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, & |
| 980 |
< |
do_pot) |
| 981 |
< |
endif |
| 982 |
< |
|
| 1124 |
> |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 1125 |
> |
call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, & |
| 1126 |
> |
do_pot) |
| 1127 |
|
endif |
| 1128 |
|
|
| 1129 |
< |
|
| 1130 |
< |
! write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape |
| 1131 |
< |
|
| 988 |
< |
if (FF_uses_Shapes .and. SIM_uses_Shapes) then |
| 989 |
< |
if ( PropertyMap(me_i)%is_Shape .and. & |
| 990 |
< |
PropertyMap(me_j)%is_Shape ) then |
| 991 |
< |
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 992 |
< |
pot, A, f, t, do_pot) |
| 993 |
< |
endif |
| 994 |
< |
if ( (PropertyMap(me_i)%is_Shape .and. & |
| 995 |
< |
PropertyMap(me_j)%is_LennardJones) .or. & |
| 996 |
< |
(PropertyMap(me_i)%is_LennardJones .and. & |
| 997 |
< |
PropertyMap(me_j)%is_Shape) ) then |
| 998 |
< |
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 999 |
< |
pot, A, f, t, do_pot) |
| 1000 |
< |
endif |
| 1129 |
> |
if ( iand(iMap, SHAPE_PAIR).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 |
| 1133 |
|
|
| 1134 |
+ |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
| 1135 |
+ |
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1136 |
+ |
pot, A, f, t, do_pot) |
| 1137 |
+ |
endif |
| 1138 |
+ |
|
| 1139 |
|
end subroutine do_pair |
| 1140 |
|
|
| 1141 |
|
subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, & |
| 1152 |
|
real ( kind = dp ), intent(inout) :: rijsq, rcijsq |
| 1153 |
|
real ( kind = dp ) :: r, rc |
| 1154 |
|
real ( kind = dp ), intent(inout) :: d(3), dc(3) |
| 1019 |
– |
|
| 1020 |
– |
logical :: is_EAM_i, is_EAM_j |
| 1021 |
– |
|
| 1022 |
– |
integer :: me_i, me_j |
| 1155 |
|
|
| 1156 |
+ |
integer :: me_i, me_j, iMap |
| 1157 |
|
|
| 1025 |
– |
r = sqrt(rijsq) |
| 1026 |
– |
if (SIM_uses_molecular_cutoffs) then |
| 1027 |
– |
rc = sqrt(rcijsq) |
| 1028 |
– |
else |
| 1029 |
– |
rc = r |
| 1030 |
– |
endif |
| 1031 |
– |
|
| 1032 |
– |
|
| 1158 |
|
#ifdef IS_MPI |
| 1159 |
|
me_i = atid_row(i) |
| 1160 |
|
me_j = atid_col(j) |
| 1163 |
|
me_j = atid(j) |
| 1164 |
|
#endif |
| 1165 |
|
|
| 1166 |
< |
if (FF_uses_EAM .and. SIM_uses_EAM) then |
| 1166 |
> |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 1167 |
|
|
| 1168 |
< |
if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) & |
| 1168 |
> |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 1169 |
|
call calc_EAM_prepair_rho(i, j, d, r, rijsq ) |
| 1045 |
– |
|
| 1170 |
|
endif |
| 1171 |
< |
|
| 1171 |
> |
|
| 1172 |
|
end subroutine do_prepair |
| 1173 |
|
|
| 1174 |
|
|
| 1364 |
|
logical :: doesit |
| 1365 |
|
doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. & |
| 1366 |
|
FF_uses_Quadrupoles .or. FF_uses_Sticky .or. & |
| 1367 |
< |
FF_uses_GayBerne .or. FF_uses_Shapes |
| 1367 |
> |
FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes |
| 1368 |
|
end function FF_UsesDirectionalAtoms |
| 1369 |
|
|
| 1370 |
|
function FF_RequiresPrepairCalc() result(doesit) |
