| 45 |
|
|
| 46 |
|
!! @author Charles F. Vardeman II |
| 47 |
|
!! @author Matthew Meineke |
| 48 |
< |
!! @version $Id: doForces.F90,v 1.22 2005-06-28 13:58:45 gezelter Exp $, $Date: 2005-06-28 13:58:45 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $ |
| 48 |
> |
!! @version $Id: doForces.F90,v 1.26 2005-07-29 19:38:27 gezelter Exp $, $Date: 2005-07-29 19:38:27 $, $Name: not supported by cvs2svn $, $Revision: 1.26 $ |
| 49 |
|
|
| 50 |
|
|
| 51 |
|
module doForces |
| 82 |
|
logical, save :: haveNeighborList = .false. |
| 83 |
|
logical, save :: haveSIMvariables = .false. |
| 84 |
|
logical, save :: haveSaneForceField = .false. |
| 85 |
+ |
logical, save :: haveInteractionMap = .false. |
| 86 |
|
|
| 87 |
|
logical, save :: FF_uses_DirectionalAtoms |
| 88 |
|
logical, save :: FF_uses_LennardJones |
| 122 |
|
public :: init_FF |
| 123 |
|
public :: do_force_loop |
| 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 |
| 137 |
|
|
| 138 |
|
type, public :: Interaction |
| 139 |
|
integer :: InteractionHash |
| 140 |
< |
real(kind=dp) :: rCut |
| 140 |
> |
real(kind=dp) :: rCut = 0.0_dp |
| 141 |
> |
real(kind=dp) :: rCutSq = 0.0_dp |
| 142 |
> |
real(kind=dp) :: rListSq = 0.0_dp |
| 143 |
|
end type Interaction |
| 144 |
|
|
| 145 |
|
type(Interaction), dimension(:,:),allocatable :: InteractionMap |
| 146 |
|
|
| 147 |
< |
!public :: addInteraction |
| 140 |
< |
!public :: setInteractionHash |
| 141 |
< |
!public :: getInteractionHash |
| 142 |
< |
public :: createInteractionMap |
| 147 |
> |
|
| 148 |
|
|
| 149 |
|
contains |
| 150 |
|
|
| 151 |
|
|
| 152 |
|
subroutine createInteractionMap(status) |
| 153 |
|
integer :: nAtypes |
| 154 |
< |
integer :: status |
| 154 |
> |
integer, intent(out) :: status |
| 155 |
|
integer :: i |
| 156 |
|
integer :: j |
| 157 |
|
integer :: ihash |
| 158 |
|
real(kind=dp) :: myRcut |
| 159 |
< |
! Test Types |
| 159 |
> |
!! Test Types |
| 160 |
|
logical :: i_is_LJ |
| 161 |
|
logical :: i_is_Elect |
| 162 |
|
logical :: i_is_Sticky |
| 172 |
|
logical :: j_is_EAM |
| 173 |
|
logical :: j_is_Shape |
| 174 |
|
|
| 175 |
< |
|
| 175 |
> |
status = 0 |
| 176 |
> |
|
| 177 |
|
if (.not. associated(atypes)) then |
| 178 |
|
call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!") |
| 179 |
|
status = -1 |
| 248 |
|
end do |
| 249 |
|
|
| 250 |
|
end do |
| 251 |
+ |
|
| 252 |
+ |
haveInteractionMap = .true. |
| 253 |
|
end subroutine createInteractionMap |
| 254 |
+ |
|
| 255 |
+ |
! Query each potential and return the cutoff for that potential. We |
| 256 |
+ |
! build the neighbor list based on the largest cutoff value for that |
| 257 |
+ |
! atype. Each potential can decide whether to calculate the force for |
| 258 |
+ |
! that atype based upon it's own cutoff. |
| 259 |
+ |
|
| 260 |
+ |
subroutine createRcuts(defaultRcut, defaultSkinThickness, stat) |
| 261 |
+ |
|
| 262 |
+ |
real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness |
| 263 |
+ |
integer :: iMap |
| 264 |
+ |
integer :: map_i,map_j |
| 265 |
+ |
real(kind=dp) :: thisRCut = 0.0_dp |
| 266 |
+ |
real(kind=dp) :: actualCutoff = 0.0_dp |
| 267 |
+ |
integer, intent(out) :: stat |
| 268 |
+ |
integer :: nAtypes |
| 269 |
+ |
integer :: myStatus |
| 270 |
+ |
|
| 271 |
+ |
stat = 0 |
| 272 |
+ |
if (.not. haveInteractionMap) then |
| 273 |
+ |
|
| 274 |
+ |
call createInteractionMap(myStatus) |
| 275 |
+ |
|
| 276 |
+ |
if (myStatus .ne. 0) then |
| 277 |
+ |
write(default_error, *) 'createInteractionMap failed in doForces!' |
| 278 |
+ |
stat = -1 |
| 279 |
+ |
return |
| 280 |
+ |
endif |
| 281 |
+ |
endif |
| 282 |
+ |
|
| 283 |
+ |
nAtypes = getSize(atypes) |
| 284 |
+ |
!! If we pass a default rcut, set all atypes to that cutoff distance |
| 285 |
+ |
if(present(defaultRList)) then |
| 286 |
+ |
InteractionMap(:,:)%rCut = defaultRCut |
| 287 |
+ |
InteractionMap(:,:)%rCutSq = defaultRCut*defaultRCut |
| 288 |
+ |
InteractionMap(:,:)%rListSq = (defaultRCut+defaultSkinThickness)**2 |
| 289 |
+ |
haveRlist = .true. |
| 290 |
+ |
return |
| 291 |
+ |
end if |
| 292 |
+ |
|
| 293 |
+ |
do map_i = 1,nAtypes |
| 294 |
+ |
do map_j = map_i,nAtypes |
| 295 |
+ |
iMap = InteractionMap(map_i, map_j)%InteractionHash |
| 296 |
+ |
|
| 297 |
+ |
if ( iand(iMap, LJ_PAIR).ne.0 ) then |
| 298 |
+ |
! thisRCut = getLJCutOff(map_i,map_j) |
| 299 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 300 |
+ |
endif |
| 301 |
+ |
|
| 302 |
+ |
if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then |
| 303 |
+ |
! thisRCut = getElectrostaticCutOff(map_i,map_j) |
| 304 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 305 |
+ |
endif |
| 306 |
+ |
|
| 307 |
+ |
if ( iand(iMap, STICKY_PAIR).ne.0 ) then |
| 308 |
+ |
! thisRCut = getStickyCutOff(map_i,map_j) |
| 309 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 310 |
+ |
endif |
| 311 |
+ |
|
| 312 |
+ |
if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then |
| 313 |
+ |
! thisRCut = getStickyPowerCutOff(map_i,map_j) |
| 314 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 315 |
+ |
endif |
| 316 |
+ |
|
| 317 |
+ |
if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then |
| 318 |
+ |
! thisRCut = getGayberneCutOff(map_i,map_j) |
| 319 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 320 |
+ |
endif |
| 321 |
+ |
|
| 322 |
+ |
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
| 323 |
+ |
! thisRCut = getGaybrneLJCutOff(map_i,map_j) |
| 324 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 325 |
+ |
endif |
| 326 |
+ |
|
| 327 |
+ |
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
| 328 |
+ |
! thisRCut = getEAMCutOff(map_i,map_j) |
| 329 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 330 |
+ |
endif |
| 331 |
+ |
|
| 332 |
+ |
if ( iand(iMap, SHAPE_PAIR).ne.0 ) then |
| 333 |
+ |
! thisRCut = getShapeCutOff(map_i,map_j) |
| 334 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 335 |
+ |
endif |
| 336 |
+ |
|
| 337 |
+ |
if ( iand(iMap, SHAPE_LJ).ne.0 ) then |
| 338 |
+ |
! thisRCut = getShapeLJCutOff(map_i,map_j) |
| 339 |
+ |
if (thisRcut > actualCutoff) actualCutoff = thisRcut |
| 340 |
+ |
endif |
| 341 |
+ |
InteractionMap(map_i, map_j)%rCut = actualCutoff |
| 342 |
+ |
InteractionMap(map_i, map_j)%rCutSq = actualCutoff * actualCutoff |
| 343 |
+ |
InteractionMap(map_i, map_j)%rListSq = (actualCutoff + skinThickness)**2 |
| 344 |
|
|
| 345 |
+ |
InteractionMap(map_j, map_i)%rCut = InteractionMap(map_i, map_j)%rCut |
| 346 |
+ |
InteractionMap(map_j, map_i)%rCutSq = InteractionMap(map_i, map_j)%rCutSq |
| 347 |
+ |
InteractionMap(map_j, map_i)%rListSq = InteractionMap(map_i, map_j)%rListSq |
| 348 |
+ |
end do |
| 349 |
+ |
end do |
| 350 |
+ |
! now the groups |
| 351 |
|
|
| 352 |
|
|
| 353 |
+ |
|
| 354 |
+ |
haveRlist = .true. |
| 355 |
+ |
end subroutine createRcuts |
| 356 |
+ |
|
| 357 |
+ |
|
| 358 |
|
!!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF |
| 359 |
|
!!$ subroutine setRlistDF( this_rlist ) |
| 360 |
|
!!$ |
| 398 |
|
error = 0 |
| 399 |
|
|
| 400 |
|
if (.not. haveInteractionMap) then |
| 401 |
< |
|
| 402 |
< |
myStatus = 0 |
| 294 |
< |
|
| 401 |
> |
|
| 402 |
> |
myStatus = 0 |
| 403 |
|
call createInteractionMap(myStatus) |
| 404 |
< |
|
| 404 |
> |
|
| 405 |
|
if (myStatus .ne. 0) then |
| 406 |
|
write(default_error, *) 'createInteractionMap failed in doForces!' |
| 407 |
|
error = -1 |
| 661 |
|
integer :: localError |
| 662 |
|
integer :: propPack_i, propPack_j |
| 663 |
|
integer :: loopStart, loopEnd, loop |
| 664 |
< |
|
| 664 |
> |
integer :: iMap |
| 665 |
|
real(kind=dp) :: listSkin = 1.0 |
| 666 |
|
|
| 667 |
|
!! initialize local variables |
| 753 |
|
#endif |
| 754 |
|
outer: do i = istart, iend |
| 755 |
|
|
| 756 |
+ |
#ifdef IS_MPI |
| 757 |
+ |
me_i = atid_row(i) |
| 758 |
+ |
#else |
| 759 |
+ |
me_i = atid(i) |
| 760 |
+ |
#endif |
| 761 |
+ |
|
| 762 |
|
if (update_nlist) point(i) = nlist + 1 |
| 763 |
|
|
| 764 |
|
n_in_i = groupStartRow(i+1) - groupStartRow(i) |
| 786 |
|
endif |
| 787 |
|
|
| 788 |
|
#ifdef IS_MPI |
| 789 |
+ |
me_j = atid_col(j) |
| 790 |
|
call get_interatomic_vector(q_group_Row(:,i), & |
| 791 |
|
q_group_Col(:,j), d_grp, rgrpsq) |
| 792 |
|
#else |
| 793 |
+ |
me_j = atid(j) |
| 794 |
|
call get_interatomic_vector(q_group(:,i), & |
| 795 |
|
q_group(:,j), d_grp, rgrpsq) |
| 796 |
|
#endif |
| 797 |
|
|
| 798 |
< |
if (rgrpsq < rlistsq) then |
| 798 |
> |
if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then |
| 799 |
|
if (update_nlist) then |
| 800 |
|
nlist = nlist + 1 |
| 801 |
|
|
| 1131 |
|
endif |
| 1132 |
|
|
| 1133 |
|
if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then |
| 1134 |
< |
call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1135 |
< |
pot, A, f, t, do_pot) |
| 1134 |
> |
! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & |
| 1135 |
> |
! pot, A, f, t, do_pot) |
| 1136 |
|
endif |
| 1137 |
|
|
| 1138 |
|
if ( iand(iMap, EAM_PAIR).ne.0 ) then |
