| 45 |
|
|
| 46 |
|
!! @author Charles F. Vardeman II |
| 47 |
|
!! @author Matthew Meineke |
| 48 |
< |
!! @version $Id: doForces.F90,v 1.24 2005-07-28 22:12:45 chuckv Exp $, $Date: 2005-07-28 22:12:45 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $ |
| 48 |
> |
!! @version $Id: doForces.F90,v 1.64 2005-11-01 19:09:23 chrisfen Exp $, $Date: 2005-11-01 19:09:23 $, $Name: not supported by cvs2svn $, $Revision: 1.64 $ |
| 49 |
|
|
| 50 |
|
|
| 51 |
|
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 |
| 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 :: haveInteractionMap = .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 |
| 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, target :: groupMaxCutoffRow |
| 127 |
> |
real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol |
| 128 |
|
|
| 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 |
|
|
| 145 |
|
contains |
| 146 |
|
|
| 147 |
< |
|
| 151 |
< |
subroutine createInteractionMap(status) |
| 147 |
> |
subroutine createInteractionHash(status) |
| 148 |
|
integer :: nAtypes |
| 149 |
|
integer, intent(out) :: status |
| 150 |
|
integer :: i |
| 151 |
|
integer :: j |
| 152 |
< |
integer :: ihash |
| 153 |
< |
real(kind=dp) :: myRcut |
| 158 |
< |
! Test Types |
| 152 |
> |
integer :: iHash |
| 153 |
> |
!! Test Types |
| 154 |
|
logical :: i_is_LJ |
| 155 |
|
logical :: i_is_Elect |
| 156 |
|
logical :: i_is_Sticky |
| 165 |
|
logical :: j_is_GB |
| 166 |
|
logical :: j_is_EAM |
| 167 |
|
logical :: j_is_Shape |
| 168 |
< |
|
| 169 |
< |
status = 0 |
| 170 |
< |
|
| 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 |
| 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 |
| 200 |
|
call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) |
| 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 |
| 250 |
– |
end subroutine createInteractionMap |
| 256 |
|
|
| 257 |
< |
! 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. |
| 258 |
< |
subroutine createRcuts(defaultRList,stat) |
| 254 |
< |
real(kind=dp), intent(in), optional :: defaultRList |
| 255 |
< |
integer :: iMap |
| 256 |
< |
integer :: map_i,map_j |
| 257 |
< |
real(kind=dp) :: thisRCut = 0.0_dp |
| 258 |
< |
real(kind=dp) :: actualCutoff = 0.0_dp |
| 259 |
< |
integer, intent(out) :: stat |
| 260 |
< |
integer :: nAtypes |
| 261 |
< |
integer :: myStatus |
| 257 |
> |
haveInteractionHash = .true. |
| 258 |
> |
end subroutine createInteractionHash |
| 259 |
|
|
| 260 |
< |
stat = 0 |
| 264 |
< |
if (.not. haveInteractionMap) then |
| 260 |
> |
subroutine createGtypeCutoffMap(stat) |
| 261 |
|
|
| 262 |
< |
call createInteractionMap(myStatus) |
| 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, *) 'createInteractionMap failed in doForces!' |
| 284 |
> |
write(default_error, *) 'createInteractionHash failed in doForces!' |
| 285 |
|
stat = -1 |
| 286 |
|
return |
| 287 |
|
endif |
| 288 |
|
endif |
| 289 |
< |
|
| 290 |
< |
|
| 289 |
> |
#ifdef IS_MPI |
| 290 |
> |
nGroupsInRow = getNgroupsInRow(plan_group_row) |
| 291 |
> |
nGroupsInCol = getNgroupsInCol(plan_group_col) |
| 292 |
> |
#endif |
| 293 |
|
nAtypes = getSize(atypes) |
| 294 |
< |
! If we pass a default rcut, set all atypes to that cutoff distance |
| 295 |
< |
if(present(defaultRList)) then |
| 296 |
< |
InteractionMap(:,:)%rList = defaultRList |
| 297 |
< |
InteractionMap(:,:)%rListSq = defaultRList*defaultRList |
| 298 |
< |
haveRlist = .true. |
| 299 |
< |
return |
| 300 |
< |
end if |
| 301 |
< |
|
| 302 |
< |
do map_i = 1,nAtypes |
| 303 |
< |
do map_j = map_i,nAtypes |
| 304 |
< |
iMap = InteractionMap(map_i, map_j)%InteractionHash |
| 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 |
< |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
| 307 |
< |
! thisRCut = getLJCutOff(map_i,map_j) |
| 308 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 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 |
|
|
| 294 |
– |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 295 |
– |
! thisRCut = getElectrostaticCutOff(map_i,map_j) |
| 296 |
– |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 297 |
– |
endif |
| 340 |
|
|
| 341 |
< |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
| 342 |
< |
! thisRCut = getStickyCutOff(map_i,map_j) |
| 343 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 302 |
< |
endif |
| 303 |
< |
|
| 304 |
< |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
| 305 |
< |
! thisRCut = getStickyPowerCutOff(map_i,map_j) |
| 306 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 307 |
< |
endif |
| 308 |
< |
|
| 309 |
< |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
| 310 |
< |
! thisRCut = getGayberneCutOff(map_i,map_j) |
| 311 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 312 |
< |
endif |
| 313 |
< |
|
| 314 |
< |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
| 315 |
< |
! thisRCut = getGaybrneLJCutOff(map_i,map_j) |
| 316 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 317 |
< |
endif |
| 318 |
< |
|
| 319 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 320 |
< |
! thisRCut = getEAMCutOff(map_i,map_j) |
| 321 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 322 |
< |
endif |
| 323 |
< |
|
| 324 |
< |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
| 325 |
< |
! thisRCut = getShapeCutOff(map_i,map_j) |
| 326 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 327 |
< |
endif |
| 328 |
< |
|
| 329 |
< |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
| 330 |
< |
! thisRCut = getShapeLJCutOff(map_i,map_j) |
| 331 |
< |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 332 |
< |
endif |
| 333 |
< |
InteractionMap(map_i, map_j)%rList = actualCutoff |
| 334 |
< |
InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff |
| 335 |
< |
end do |
| 336 |
< |
end do |
| 337 |
< |
haveRlist = .true. |
| 338 |
< |
end subroutine createRcuts |
| 341 |
> |
if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then |
| 342 |
> |
biggestAtypeCutoff = atypeMaxCutoff(i) |
| 343 |
> |
endif |
| 344 |
|
|
| 345 |
+ |
endif |
| 346 |
+ |
enddo |
| 347 |
+ |
|
| 348 |
|
|
| 349 |
< |
!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF |
| 350 |
< |
!!$ subroutine setRlistDF( this_rlist ) |
| 351 |
< |
!!$ |
| 352 |
< |
!!$ real(kind=dp) :: this_rlist |
| 353 |
< |
!!$ |
| 354 |
< |
!!$ rlist = this_rlist |
| 355 |
< |
!!$ rlistsq = rlist * rlist |
| 356 |
< |
!!$ |
| 357 |
< |
!!$ haveRlist = .true. |
| 358 |
< |
!!$ |
| 359 |
< |
!!$ end subroutine setRlistDF |
| 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() |
| 356 |
– |
SIM_uses_LennardJones = SimUsesLennardJones() |
| 357 |
– |
SIM_uses_Electrostatics = SimUsesElectrostatics() |
| 358 |
– |
SIM_uses_Charges = SimUsesCharges() |
| 359 |
– |
SIM_uses_Dipoles = SimUsesDipoles() |
| 360 |
– |
SIM_uses_Sticky = SimUsesSticky() |
| 361 |
– |
SIM_uses_StickyPower = SimUsesStickyPower() |
| 362 |
– |
SIM_uses_GayBerne = SimUsesGayBerne() |
| 571 |
|
SIM_uses_EAM = SimUsesEAM() |
| 364 |
– |
SIM_uses_Shapes = SimUsesShapes() |
| 365 |
– |
SIM_uses_FLARB = SimUsesFLARB() |
| 366 |
– |
SIM_uses_RF = SimUsesRF() |
| 572 |
|
SIM_requires_postpair_calc = SimRequiresPostpairCalc() |
| 573 |
|
SIM_requires_prepair_calc = SimRequiresPrepairCalc() |
| 574 |
|
SIM_uses_PBC = SimUsesPBC() |
| 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) |
| 388 |
< |
|
| 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 |
| 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!' |
| 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 |
| 432 |
< |
|
| 643 |
> |
integer, intent(in) :: thisESM |
| 644 |
|
integer, intent(out) :: thisStat |
| 645 |
|
integer :: my_status, nMatches |
| 646 |
|
integer, pointer :: MatchList(:) => null() |
| 436 |
– |
real(kind=dp) :: rcut, rrf, rt, dielect |
| 647 |
|
|
| 648 |
|
!! assume things are copacetic, unless they aren't |
| 649 |
|
thisStat = 0 |
| 650 |
|
|
| 651 |
< |
!! Fortran's version of a cast: |
| 442 |
< |
FF_uses_RF = use_RF_c |
| 651 |
> |
electrostaticSummationMethod = thisESM |
| 652 |
|
|
| 653 |
|
!! init_FF is called *after* all of the atom types have been |
| 654 |
|
!! defined in atype_module using the new_atype subroutine. |
| 657 |
|
!! interactions are used by the force field. |
| 658 |
|
|
| 659 |
|
FF_uses_DirectionalAtoms = .false. |
| 451 |
– |
FF_uses_LennardJones = .false. |
| 452 |
– |
FF_uses_Electrostatics = .false. |
| 453 |
– |
FF_uses_Charges = .false. |
| 660 |
|
FF_uses_Dipoles = .false. |
| 455 |
– |
FF_uses_Sticky = .false. |
| 456 |
– |
FF_uses_StickyPower = .false. |
| 661 |
|
FF_uses_GayBerne = .false. |
| 662 |
|
FF_uses_EAM = .false. |
| 459 |
– |
FF_uses_Shapes = .false. |
| 460 |
– |
FF_uses_FLARB = .false. |
| 663 |
|
|
| 664 |
|
call getMatchingElementList(atypes, "is_Directional", .true., & |
| 665 |
|
nMatches, MatchList) |
| 666 |
|
if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true. |
| 667 |
|
|
| 466 |
– |
call getMatchingElementList(atypes, "is_LennardJones", .true., & |
| 467 |
– |
nMatches, MatchList) |
| 468 |
– |
if (nMatches .gt. 0) FF_uses_LennardJones = .true. |
| 469 |
– |
|
| 470 |
– |
call getMatchingElementList(atypes, "is_Electrostatic", .true., & |
| 471 |
– |
nMatches, MatchList) |
| 472 |
– |
if (nMatches .gt. 0) then |
| 473 |
– |
FF_uses_Electrostatics = .true. |
| 474 |
– |
endif |
| 475 |
– |
|
| 476 |
– |
call getMatchingElementList(atypes, "is_Charge", .true., & |
| 477 |
– |
nMatches, MatchList) |
| 478 |
– |
if (nMatches .gt. 0) then |
| 479 |
– |
FF_uses_Charges = .true. |
| 480 |
– |
FF_uses_Electrostatics = .true. |
| 481 |
– |
endif |
| 482 |
– |
|
| 668 |
|
call getMatchingElementList(atypes, "is_Dipole", .true., & |
| 669 |
|
nMatches, MatchList) |
| 670 |
< |
if (nMatches .gt. 0) then |
| 486 |
< |
FF_uses_Dipoles = .true. |
| 487 |
< |
FF_uses_Electrostatics = .true. |
| 488 |
< |
FF_uses_DirectionalAtoms = .true. |
| 489 |
< |
endif |
| 490 |
< |
|
| 491 |
< |
call getMatchingElementList(atypes, "is_Quadrupole", .true., & |
| 492 |
< |
nMatches, MatchList) |
| 493 |
< |
if (nMatches .gt. 0) then |
| 494 |
< |
FF_uses_Quadrupoles = .true. |
| 495 |
< |
FF_uses_Electrostatics = .true. |
| 496 |
< |
FF_uses_DirectionalAtoms = .true. |
| 497 |
< |
endif |
| 498 |
< |
|
| 499 |
< |
call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, & |
| 500 |
< |
MatchList) |
| 501 |
< |
if (nMatches .gt. 0) then |
| 502 |
< |
FF_uses_Sticky = .true. |
| 503 |
< |
FF_uses_DirectionalAtoms = .true. |
| 504 |
< |
endif |
| 505 |
< |
|
| 506 |
< |
call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, & |
| 507 |
< |
MatchList) |
| 508 |
< |
if (nMatches .gt. 0) then |
| 509 |
< |
FF_uses_StickyPower = .true. |
| 510 |
< |
FF_uses_DirectionalAtoms = .true. |
| 511 |
< |
endif |
| 670 |
> |
if (nMatches .gt. 0) FF_uses_Dipoles = .true. |
| 671 |
|
|
| 672 |
|
call getMatchingElementList(atypes, "is_GayBerne", .true., & |
| 673 |
|
nMatches, MatchList) |
| 674 |
< |
if (nMatches .gt. 0) then |
| 516 |
< |
FF_uses_GayBerne = .true. |
| 517 |
< |
FF_uses_DirectionalAtoms = .true. |
| 518 |
< |
endif |
| 674 |
> |
if (nMatches .gt. 0) FF_uses_GayBerne = .true. |
| 675 |
|
|
| 676 |
|
call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList) |
| 677 |
|
if (nMatches .gt. 0) FF_uses_EAM = .true. |
| 678 |
|
|
| 523 |
– |
call getMatchingElementList(atypes, "is_Shape", .true., & |
| 524 |
– |
nMatches, MatchList) |
| 525 |
– |
if (nMatches .gt. 0) then |
| 526 |
– |
FF_uses_Shapes = .true. |
| 527 |
– |
FF_uses_DirectionalAtoms = .true. |
| 528 |
– |
endif |
| 529 |
– |
|
| 530 |
– |
call getMatchingElementList(atypes, "is_FLARB", .true., & |
| 531 |
– |
nMatches, MatchList) |
| 532 |
– |
if (nMatches .gt. 0) FF_uses_FLARB = .true. |
| 679 |
|
|
| 534 |
– |
!! Assume sanity (for the sake of argument) |
| 680 |
|
haveSaneForceField = .true. |
| 681 |
|
|
| 537 |
– |
!! check to make sure the FF_uses_RF setting makes sense |
| 538 |
– |
|
| 539 |
– |
if (FF_uses_dipoles) then |
| 540 |
– |
if (FF_uses_RF) then |
| 541 |
– |
dielect = getDielect() |
| 542 |
– |
call initialize_rf(dielect) |
| 543 |
– |
endif |
| 544 |
– |
else |
| 545 |
– |
if (FF_uses_RF) then |
| 546 |
– |
write(default_error,*) 'Using Reaction Field with no dipoles? Huh?' |
| 547 |
– |
thisStat = -1 |
| 548 |
– |
haveSaneForceField = .false. |
| 549 |
– |
return |
| 550 |
– |
endif |
| 551 |
– |
endif |
| 552 |
– |
|
| 553 |
– |
!sticky module does not contain check_sticky_FF anymore |
| 554 |
– |
!if (FF_uses_sticky) then |
| 555 |
– |
! call check_sticky_FF(my_status) |
| 556 |
– |
! if (my_status /= 0) then |
| 557 |
– |
! thisStat = -1 |
| 558 |
– |
! haveSaneForceField = .false. |
| 559 |
– |
! return |
| 560 |
– |
! end if |
| 561 |
– |
!endif |
| 562 |
– |
|
| 682 |
|
if (FF_uses_EAM) then |
| 683 |
|
call init_EAM_FF(my_status) |
| 684 |
|
if (my_status /= 0) then |
| 689 |
|
end if |
| 690 |
|
endif |
| 691 |
|
|
| 573 |
– |
if (FF_uses_GayBerne) then |
| 574 |
– |
call check_gb_pair_FF(my_status) |
| 575 |
– |
if (my_status .ne. 0) then |
| 576 |
– |
thisStat = -1 |
| 577 |
– |
haveSaneForceField = .false. |
| 578 |
– |
return |
| 579 |
– |
endif |
| 580 |
– |
endif |
| 581 |
– |
|
| 582 |
– |
if (FF_uses_GayBerne .and. FF_uses_LennardJones) then |
| 583 |
– |
endif |
| 584 |
– |
|
| 692 |
|
if (.not. haveNeighborList) then |
| 693 |
|
!! Create neighbor lists |
| 694 |
|
call expandNeighborList(nLocal, my_status) |
| 722 |
|
|
| 723 |
|
!! Stress Tensor |
| 724 |
|
real( kind = dp), dimension(9) :: tau |
| 725 |
< |
real ( kind = dp ) :: pot |
| 725 |
> |
real ( kind = dp ),dimension(LR_POT_TYPES) :: pot |
| 726 |
|
logical ( kind = 2) :: do_pot_c, do_stress_c |
| 727 |
|
logical :: do_pot |
| 728 |
|
logical :: do_stress |
| 729 |
|
logical :: in_switching_region |
| 730 |
|
#ifdef IS_MPI |
| 731 |
< |
real( kind = DP ) :: pot_local |
| 731 |
> |
real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local |
| 732 |
|
integer :: nAtomsInRow |
| 733 |
|
integer :: nAtomsInCol |
| 734 |
|
integer :: nprocs |
| 743 |
|
integer :: nlist |
| 744 |
|
real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij |
| 745 |
|
real( kind = DP ) :: sw, dswdr, swderiv, mf |
| 746 |
+ |
real( kind = DP ) :: rVal |
| 747 |
|
real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij |
| 748 |
|
real(kind=dp) :: rfpot, mu_i, virial |
| 749 |
|
integer :: me_i, me_j, n_in_i, n_in_j |
| 753 |
|
integer :: localError |
| 754 |
|
integer :: propPack_i, propPack_j |
| 755 |
|
integer :: loopStart, loopEnd, loop |
| 756 |
< |
integer :: iMap |
| 757 |
< |
real(kind=dp) :: listSkin = 1.0 |
| 756 |
> |
integer :: iHash |
| 757 |
> |
integer :: i1 |
| 758 |
> |
|
| 759 |
|
|
| 760 |
|
!! initialize local variables |
| 761 |
|
|
| 852 |
|
|
| 853 |
|
if (update_nlist) then |
| 854 |
|
#ifdef IS_MPI |
| 746 |
– |
me_i = atid_row(i) |
| 855 |
|
jstart = 1 |
| 856 |
|
jend = nGroupsInCol |
| 857 |
|
#else |
| 750 |
– |
me_i = atid(i) |
| 858 |
|
jstart = i+1 |
| 859 |
|
jend = nGroups |
| 860 |
|
#endif |
| 880 |
|
me_j = atid(j) |
| 881 |
|
call get_interatomic_vector(q_group(:,i), & |
| 882 |
|
q_group(:,j), d_grp, rgrpsq) |
| 883 |
< |
#endif |
| 883 |
> |
#endif |
| 884 |
|
|
| 885 |
< |
if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then |
| 885 |
> |
if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then |
| 886 |
|
if (update_nlist) then |
| 887 |
|
nlist = nlist + 1 |
| 888 |
|
|
| 920 |
|
inner: do jb = groupStartCol(j), groupStartCol(j+1)-1 |
| 921 |
|
|
| 922 |
|
atom2 = groupListCol(jb) |
| 923 |
< |
|
| 924 |
< |
if (skipThisPair(atom1, atom2)) cycle inner |
| 925 |
< |
|
| 923 |
> |
|
| 924 |
> |
if (skipThisPair(atom1, atom2)) cycle inner |
| 925 |
> |
|
| 926 |
|
if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then |
| 927 |
|
d_atm(1:3) = d_grp(1:3) |
| 928 |
|
ratmsq = rgrpsq |
| 949 |
|
else |
| 950 |
|
#ifdef IS_MPI |
| 951 |
|
call do_pair(atom1, atom2, ratmsq, d_atm, sw, & |
| 952 |
< |
do_pot, & |
| 953 |
< |
eFrame, A, f, t, pot_local, vpair, fpair) |
| 952 |
> |
do_pot, eFrame, A, f, t, pot_local, vpair, & |
| 953 |
> |
fpair, d_grp, rgrp) |
| 954 |
|
#else |
| 955 |
|
call do_pair(atom1, atom2, ratmsq, d_atm, sw, & |
| 956 |
< |
do_pot, & |
| 957 |
< |
eFrame, A, f, t, pot, vpair, fpair) |
| 956 |
> |
do_pot, eFrame, A, f, t, pot, vpair, fpair, & |
| 957 |
> |
d_grp, rgrp) |
| 958 |
|
#endif |
| 959 |
|
|
| 960 |
|
vij = vij + vpair |
| 1003 |
|
endif |
| 1004 |
|
end if |
| 1005 |
|
enddo |
| 1006 |
+ |
|
| 1007 |
|
enddo outer |
| 1008 |
|
|
| 1009 |
|
if (update_nlist) then |
| 1063 |
|
|
| 1064 |
|
if (do_pot) then |
| 1065 |
|
! scatter/gather pot_row into the members of my column |
| 1066 |
< |
call scatter(pot_Row, pot_Temp, plan_atom_row) |
| 1067 |
< |
|
| 1066 |
> |
do i = 1,LR_POT_TYPES |
| 1067 |
> |
call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row) |
| 1068 |
> |
end do |
| 1069 |
|
! scatter/gather pot_local into all other procs |
| 1070 |
|
! add resultant to get total pot |
| 1071 |
|
do i = 1, nlocal |
| 1072 |
< |
pot_local = pot_local + pot_Temp(i) |
| 1072 |
> |
pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) & |
| 1073 |
> |
+ pot_Temp(1:LR_POT_TYPES,i) |
| 1074 |
|
enddo |
| 1075 |
|
|
| 1076 |
|
pot_Temp = 0.0_DP |
| 1077 |
< |
|
| 1078 |
< |
call scatter(pot_Col, pot_Temp, plan_atom_col) |
| 1077 |
> |
do i = 1,LR_POT_TYPES |
| 1078 |
> |
call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col) |
| 1079 |
> |
end do |
| 1080 |
|
do i = 1, nlocal |
| 1081 |
< |
pot_local = pot_local + pot_Temp(i) |
| 1081 |
> |
pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)& |
| 1082 |
> |
+ pot_Temp(1:LR_POT_TYPES,i) |
| 1083 |
|
enddo |
| 1084 |
|
|
| 1085 |
|
endif |
| 1086 |
|
#endif |
| 1087 |
|
|
| 1088 |
< |
if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then |
| 1088 |
> |
if (SIM_requires_postpair_calc) then |
| 1089 |
> |
do i = 1, nlocal |
| 1090 |
> |
|
| 1091 |
> |
! we loop only over the local atoms, so we don't need row and column |
| 1092 |
> |
! lookups for the types |
| 1093 |
> |
|
| 1094 |
> |
me_i = atid(i) |
| 1095 |
> |
|
| 1096 |
> |
! is the atom electrostatic? See if it would have an |
| 1097 |
> |
! electrostatic interaction with itself |
| 1098 |
> |
iHash = InteractionHash(me_i,me_i) |
| 1099 |
|
|
| 1100 |
< |
if (FF_uses_RF .and. SIM_uses_RF) then |
| 979 |
< |
|
| 1100 |
> |
if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1101 |
|
#ifdef IS_MPI |
| 1102 |
< |
call scatter(rf_Row,rf,plan_atom_row_3d) |
| 1103 |
< |
call scatter(rf_Col,rf_Temp,plan_atom_col_3d) |
| 983 |
< |
do i = 1,nlocal |
| 984 |
< |
rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i) |
| 985 |
< |
end do |
| 986 |
< |
#endif |
| 987 |
< |
|
| 988 |
< |
do i = 1, nLocal |
| 989 |
< |
|
| 990 |
< |
rfpot = 0.0_DP |
| 991 |
< |
#ifdef IS_MPI |
| 992 |
< |
me_i = atid_row(i) |
| 1102 |
> |
call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), & |
| 1103 |
> |
t, do_pot) |
| 1104 |
|
#else |
| 1105 |
< |
me_i = atid(i) |
| 1105 |
> |
call self_self(i, eFrame, pot(ELECTROSTATIC_POT), & |
| 1106 |
> |
t, do_pot) |
| 1107 |
|
#endif |
| 1108 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 1108 |
> |
endif |
| 1109 |
> |
|
| 1110 |
> |
|
| 1111 |
> |
!!$ if (electrostaticSummationMethod.eq.REACTION_FIELD) then |
| 1112 |
|
|
| 1113 |
< |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1114 |
< |
|
| 1115 |
< |
mu_i = getDipoleMoment(me_i) |
| 1116 |
< |
|
| 1117 |
< |
!! The reaction field needs to include a self contribution |
| 1118 |
< |
!! to the field: |
| 1119 |
< |
call accumulate_self_rf(i, mu_i, eFrame) |
| 1120 |
< |
!! Get the reaction field contribution to the |
| 1121 |
< |
!! potential and torques: |
| 1122 |
< |
call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot) |
| 1113 |
> |
! loop over the excludes to accumulate RF stuff we've |
| 1114 |
> |
! left out of the normal pair loop |
| 1115 |
> |
|
| 1116 |
> |
do i1 = 1, nSkipsForAtom(i) |
| 1117 |
> |
j = skipsForAtom(i, i1) |
| 1118 |
> |
|
| 1119 |
> |
! prevent overcounting of the skips |
| 1120 |
> |
if (i.lt.j) then |
| 1121 |
> |
call get_interatomic_vector(q(:,i), & |
| 1122 |
> |
q(:,j), d_atm, ratmsq) |
| 1123 |
> |
rVal = dsqrt(ratmsq) |
| 1124 |
> |
call get_switch(ratmsq, sw, dswdr, rVal, group_switch, & |
| 1125 |
> |
in_switching_region) |
| 1126 |
|
#ifdef IS_MPI |
| 1127 |
< |
pot_local = pot_local + rfpot |
| 1127 |
> |
call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, & |
| 1128 |
> |
vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot) |
| 1129 |
|
#else |
| 1130 |
< |
pot = pot + rfpot |
| 1131 |
< |
|
| 1130 |
> |
call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, & |
| 1131 |
> |
vpair, pot(ELECTROSTATIC_POT), f, t, do_pot) |
| 1132 |
|
#endif |
| 1133 |
< |
endif |
| 1134 |
< |
enddo |
| 1135 |
< |
endif |
| 1133 |
> |
endif |
| 1134 |
> |
enddo |
| 1135 |
> |
!!$ endif |
| 1136 |
> |
enddo |
| 1137 |
|
endif |
| 1138 |
< |
|
| 1019 |
< |
|
| 1138 |
> |
|
| 1139 |
|
#ifdef IS_MPI |
| 1140 |
< |
|
| 1140 |
> |
|
| 1141 |
|
if (do_pot) then |
| 1142 |
< |
pot = pot + pot_local |
| 1143 |
< |
!! we assume the c code will do the allreduce to get the total potential |
| 1025 |
< |
!! we could do it right here if we needed to... |
| 1142 |
> |
call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, & |
| 1143 |
> |
mpi_comm_world,mpi_err) |
| 1144 |
|
endif |
| 1145 |
< |
|
| 1145 |
> |
|
| 1146 |
|
if (do_stress) then |
| 1147 |
|
call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, & |
| 1148 |
|
mpi_comm_world,mpi_err) |
| 1149 |
|
call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, & |
| 1150 |
|
mpi_comm_world,mpi_err) |
| 1151 |
|
endif |
| 1152 |
< |
|
| 1152 |
> |
|
| 1153 |
|
#else |
| 1154 |
< |
|
| 1154 |
> |
|
| 1155 |
|
if (do_stress) then |
| 1156 |
|
tau = tau_Temp |
| 1157 |
|
virial = virial_Temp |
| 1158 |
|
endif |
| 1159 |
< |
|
| 1159 |
> |
|
| 1160 |
|
#endif |
| 1161 |
< |
|
| 1161 |
> |
|
| 1162 |
|
end subroutine do_force_loop |
| 1163 |
|
|
| 1164 |
|
subroutine do_pair(i, j, rijsq, d, sw, do_pot, & |
| 1165 |
< |
eFrame, A, f, t, pot, vpair, fpair) |
| 1165 |
> |
eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp) |
| 1166 |
> |
!!$ subroutine do_pair(i, j, rijsq, d, sw, do_pot, & |
| 1167 |
> |
!!$ eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, felec) |
| 1168 |
|
|
| 1169 |
< |
real( kind = dp ) :: pot, vpair, sw |
| 1169 |
> |
real( kind = dp ) :: vpair, sw |
| 1170 |
> |
real( kind = dp ), dimension(LR_POT_TYPES) :: pot |
| 1171 |
|
real( kind = dp ), dimension(3) :: fpair |
| 1172 |
+ |
real( kind = dp ), dimension(3) :: felec |
| 1173 |
|
real( kind = dp ), dimension(nLocal) :: mfact |
| 1174 |
|
real( kind = dp ), dimension(9,nLocal) :: eFrame |
| 1175 |
|
real( kind = dp ), dimension(9,nLocal) :: A |
| 1179 |
|
logical, intent(inout) :: do_pot |
| 1180 |
|
integer, intent(in) :: i, j |
| 1181 |
|
real ( kind = dp ), intent(inout) :: rijsq |
| 1182 |
< |
real ( kind = dp ) :: r |
| 1182 |
> |
real ( kind = dp ), intent(inout) :: r_grp |
| 1183 |
|
real ( kind = dp ), intent(inout) :: d(3) |
| 1184 |
< |
real ( kind = dp ) :: ebalance |
| 1184 |
> |
real ( kind = dp ), intent(inout) :: d_grp(3) |
| 1185 |
> |
real ( kind = dp ) :: r |
| 1186 |
|
integer :: me_i, me_j |
| 1187 |
|
|
| 1188 |
< |
integer :: iMap |
| 1188 |
> |
integer :: iHash |
| 1189 |
|
|
| 1190 |
|
r = sqrt(rijsq) |
| 1191 |
|
vpair = 0.0d0 |
| 1199 |
|
me_j = atid(j) |
| 1200 |
|
#endif |
| 1201 |
|
|
| 1202 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 1203 |
< |
|
| 1204 |
< |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
| 1205 |
< |
call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) |
| 1202 |
> |
iHash = InteractionHash(me_i, me_j) |
| 1203 |
> |
|
| 1204 |
> |
if ( iand(iHash, LJ_PAIR).ne.0 ) then |
| 1205 |
> |
call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1206 |
> |
pot(VDW_POT), f, do_pot) |
| 1207 |
|
endif |
| 1208 |
< |
|
| 1209 |
< |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1208 |
> |
|
| 1209 |
> |
if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then |
| 1210 |
|
call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1211 |
< |
pot, eFrame, f, t, do_pot) |
| 1212 |
< |
|
| 1213 |
< |
if (FF_uses_RF .and. SIM_uses_RF) then |
| 1214 |
< |
|
| 1215 |
< |
! CHECK ME (RF needs to know about all electrostatic types) |
| 1216 |
< |
call accumulate_rf(i, j, r, eFrame, sw) |
| 1093 |
< |
call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair) |
| 1094 |
< |
endif |
| 1095 |
< |
|
| 1096 |
< |
endif |
| 1097 |
< |
|
| 1098 |
< |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
| 1211 |
> |
pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot) |
| 1212 |
> |
!!$ call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1213 |
> |
!!$ pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot, felec) |
| 1214 |
> |
endif |
| 1215 |
> |
|
| 1216 |
> |
if ( iand(iHash, STICKY_PAIR).ne.0 ) then |
| 1217 |
|
call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1218 |
< |
pot, A, f, t, do_pot) |
| 1218 |
> |
pot(HB_POT), A, f, t, do_pot) |
| 1219 |
|
endif |
| 1220 |
< |
|
| 1221 |
< |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
| 1220 |
> |
|
| 1221 |
> |
if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then |
| 1222 |
|
call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1223 |
< |
pot, A, f, t, do_pot) |
| 1223 |
> |
pot(HB_POT), A, f, t, do_pot) |
| 1224 |
|
endif |
| 1225 |
< |
|
| 1226 |
< |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
| 1225 |
> |
|
| 1226 |
> |
if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then |
| 1227 |
|
call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1228 |
< |
pot, A, f, t, do_pot) |
| 1228 |
> |
pot(VDW_POT), A, f, t, do_pot) |
| 1229 |
|
endif |
| 1230 |
|
|
| 1231 |
< |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
| 1232 |
< |
! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1233 |
< |
! pot, A, f, t, do_pot) |
| 1231 |
> |
if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then |
| 1232 |
> |
call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1233 |
> |
pot(VDW_POT), A, f, t, do_pot) |
| 1234 |
|
endif |
| 1235 |
< |
|
| 1236 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 1237 |
< |
call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, & |
| 1238 |
< |
do_pot) |
| 1235 |
> |
|
| 1236 |
> |
if ( iand(iHash, EAM_PAIR).ne.0 ) then |
| 1237 |
> |
call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1238 |
> |
pot(METALLIC_POT), f, do_pot) |
| 1239 |
|
endif |
| 1240 |
< |
|
| 1241 |
< |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
| 1240 |
> |
|
| 1241 |
> |
if ( iand(iHash, SHAPE_PAIR).ne.0 ) then |
| 1242 |
|
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1243 |
< |
pot, A, f, t, do_pot) |
| 1243 |
> |
pot(VDW_POT), A, f, t, do_pot) |
| 1244 |
|
endif |
| 1245 |
< |
|
| 1246 |
< |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
| 1245 |
> |
|
| 1246 |
> |
if ( iand(iHash, SHAPE_LJ).ne.0 ) then |
| 1247 |
|
call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1248 |
< |
pot, A, f, t, do_pot) |
| 1248 |
> |
pot(VDW_POT), A, f, t, do_pot) |
| 1249 |
|
endif |
| 1250 |
< |
|
| 1250 |
> |
|
| 1251 |
|
end subroutine do_pair |
| 1252 |
|
|
| 1253 |
|
subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, & |
| 1254 |
|
do_pot, do_stress, eFrame, A, f, t, pot) |
| 1255 |
|
|
| 1256 |
< |
real( kind = dp ) :: pot, sw |
| 1256 |
> |
real( kind = dp ) :: sw |
| 1257 |
> |
real( kind = dp ), dimension(LR_POT_TYPES) :: pot |
| 1258 |
|
real( kind = dp ), dimension(9,nLocal) :: eFrame |
| 1259 |
|
real (kind=dp), dimension(9,nLocal) :: A |
| 1260 |
|
real (kind=dp), dimension(3,nLocal) :: f |
| 1266 |
|
real ( kind = dp ) :: r, rc |
| 1267 |
|
real ( kind = dp ), intent(inout) :: d(3), dc(3) |
| 1268 |
|
|
| 1269 |
< |
integer :: me_i, me_j, iMap |
| 1269 |
> |
integer :: me_i, me_j, iHash |
| 1270 |
|
|
| 1271 |
+ |
r = sqrt(rijsq) |
| 1272 |
+ |
|
| 1273 |
|
#ifdef IS_MPI |
| 1274 |
|
me_i = atid_row(i) |
| 1275 |
|
me_j = atid_col(j) |
| 1278 |
|
me_j = atid(j) |
| 1279 |
|
#endif |
| 1280 |
|
|
| 1281 |
< |
iMap = InteractionMap(me_i, me_j)%InteractionHash |
| 1281 |
> |
iHash = InteractionHash(me_i, me_j) |
| 1282 |
|
|
| 1283 |
< |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 1283 |
> |
if ( iand(iHash, EAM_PAIR).ne.0 ) then |
| 1284 |
|
call calc_EAM_prepair_rho(i, j, d, r, rijsq ) |
| 1285 |
|
endif |
| 1286 |
|
|
| 1289 |
|
|
| 1290 |
|
subroutine do_preforce(nlocal,pot) |
| 1291 |
|
integer :: nlocal |
| 1292 |
< |
real( kind = dp ) :: pot |
| 1292 |
> |
real( kind = dp ),dimension(LR_POT_TYPES) :: pot |
| 1293 |
|
|
| 1294 |
|
if (FF_uses_EAM .and. SIM_uses_EAM) then |
| 1295 |
< |
call calc_EAM_preforce_Frho(nlocal,pot) |
| 1295 |
> |
call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT)) |
| 1296 |
|
endif |
| 1297 |
|
|
| 1298 |
|
|
| 1378 |
|
pot_Col = 0.0_dp |
| 1379 |
|
pot_Temp = 0.0_dp |
| 1380 |
|
|
| 1260 |
– |
rf_Row = 0.0_dp |
| 1261 |
– |
rf_Col = 0.0_dp |
| 1262 |
– |
rf_Temp = 0.0_dp |
| 1263 |
– |
|
| 1381 |
|
#endif |
| 1382 |
|
|
| 1383 |
|
if (FF_uses_EAM .and. SIM_uses_EAM) then |
| 1473 |
|
|
| 1474 |
|
function FF_UsesDirectionalAtoms() result(doesit) |
| 1475 |
|
logical :: doesit |
| 1476 |
< |
doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. & |
| 1360 |
< |
FF_uses_Quadrupoles .or. FF_uses_Sticky .or. & |
| 1361 |
< |
FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes |
| 1476 |
> |
doesit = FF_uses_DirectionalAtoms |
| 1477 |
|
end function FF_UsesDirectionalAtoms |
| 1478 |
|
|
| 1479 |
|
function FF_RequiresPrepairCalc() result(doesit) |
| 1481 |
|
doesit = FF_uses_EAM |
| 1482 |
|
end function FF_RequiresPrepairCalc |
| 1483 |
|
|
| 1369 |
– |
function FF_RequiresPostpairCalc() result(doesit) |
| 1370 |
– |
logical :: doesit |
| 1371 |
– |
doesit = FF_uses_RF |
| 1372 |
– |
end function FF_RequiresPostpairCalc |
| 1373 |
– |
|
| 1484 |
|
#ifdef PROFILE |
| 1485 |
|
function getforcetime() result(totalforcetime) |
| 1486 |
|
real(kind=dp) :: totalforcetime |
