| 4 |
|
|
| 5 |
|
!! @author Charles F. Vardeman II |
| 6 |
|
!! @author Matthew Meineke |
| 7 |
< |
!! @version $Id: do_Forces.F90,v 1.29 2003-08-13 21:20:20 chuckv Exp $, $Date: 2003-08-13 21:20:20 $, $Name: not supported by cvs2svn $, $Revision: 1.29 $ |
| 7 |
> |
!! @version $Id: do_Forces.F90,v 1.43 2004-01-05 22:18:52 chuckv Exp $, $Date: 2004-01-05 22:18:52 $, $Name: not supported by cvs2svn $, $Revision: 1.43 $ |
| 8 |
|
|
| 9 |
|
module do_Forces |
| 10 |
|
use force_globals |
| 46 |
|
public :: setRlistDF |
| 47 |
|
|
| 48 |
|
#ifdef PROFILE |
| 49 |
< |
real(kind = dp) :: forceTime |
| 50 |
< |
real(kind = dp) :: forceTimeInitial, forceTimeFinal |
| 51 |
< |
real(kind = dp) :: globalForceTime |
| 52 |
< |
real(kind = dp) :: maxForceTime |
| 53 |
< |
integer, save :: nloops = 0 |
| 49 |
> |
public :: getforcetime |
| 50 |
> |
real, save :: forceTime = 0 |
| 51 |
> |
real :: forceTimeInitial, forceTimeFinal |
| 52 |
> |
integer :: nLoops |
| 53 |
|
#endif |
| 54 |
|
|
| 55 |
+ |
logical, allocatable :: propertyMapI(:,:) |
| 56 |
+ |
logical, allocatable :: propertyMapJ(:,:) |
| 57 |
+ |
|
| 58 |
|
contains |
| 59 |
|
|
| 60 |
|
subroutine setRlistDF( this_rlist ) |
| 156 |
|
|
| 157 |
|
|
| 158 |
|
if (FF_uses_EAM) then |
| 159 |
< |
call init_EAM_FF(my_status) |
| 159 |
> |
call init_EAM_FF(my_status) |
| 160 |
|
if (my_status /= 0) then |
| 161 |
+ |
write(*,*) "init_EAM_FF returned a bad status" |
| 162 |
|
thisStat = -1 |
| 163 |
|
return |
| 164 |
|
end if |
| 207 |
|
real ( kind = dp ), dimension(3,getNlocal()) :: f |
| 208 |
|
!! Torsion array provided by C, dimensioned by getNlocal |
| 209 |
|
real( kind = dp ), dimension(3,getNlocal()) :: t |
| 210 |
+ |
|
| 211 |
|
!! Stress Tensor |
| 212 |
|
real( kind = dp), dimension(9) :: tau |
| 213 |
|
real ( kind = dp ) :: pot |
| 228 |
|
real( kind = DP ) :: rijsq |
| 229 |
|
real(kind=dp),dimension(3) :: d |
| 230 |
|
real(kind=dp) :: rfpot, mu_i, virial |
| 231 |
< |
integer :: me_i |
| 231 |
> |
integer :: me_i, me_j |
| 232 |
|
logical :: is_dp_i |
| 233 |
|
integer :: neighborListSize |
| 234 |
|
integer :: listerror, error |
| 235 |
|
integer :: localError |
| 236 |
+ |
integer :: propPack_i, propPack_j |
| 237 |
|
|
| 238 |
< |
real(kind=dp) :: listSkin = 1.0 |
| 234 |
< |
|
| 238 |
> |
real(kind=dp) :: listSkin = 1.0 |
| 239 |
|
|
| 240 |
|
!! initialize local variables |
| 241 |
|
|
| 260 |
|
do_pot = do_pot_c |
| 261 |
|
do_stress = do_stress_c |
| 262 |
|
|
| 263 |
< |
|
| 263 |
> |
|
| 264 |
> |
#ifdef IS_MPI |
| 265 |
> |
if (.not.allocated(propertyMapI)) then |
| 266 |
> |
allocate(propertyMapI(5,nrow)) |
| 267 |
> |
endif |
| 268 |
> |
|
| 269 |
> |
do i = 1, nrow |
| 270 |
> |
me_i = atid_row(i) |
| 271 |
> |
#else |
| 272 |
> |
if (.not.allocated(propertyMapI)) then |
| 273 |
> |
allocate(propertyMapI(5,nlocal)) |
| 274 |
> |
endif |
| 275 |
> |
|
| 276 |
> |
do i = 1, natoms |
| 277 |
> |
me_i = atid(i) |
| 278 |
> |
#endif |
| 279 |
> |
|
| 280 |
> |
propertyMapI(1:5,i) = .false. |
| 281 |
> |
|
| 282 |
> |
call getElementProperty(atypes, me_i, "propertyPack", propPack_i) |
| 283 |
> |
|
| 284 |
> |
! unpack the properties |
| 285 |
> |
|
| 286 |
> |
if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) & |
| 287 |
> |
propertyMapI(1, i) = .true. |
| 288 |
> |
if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) & |
| 289 |
> |
propertyMapI(2, i) = .true. |
| 290 |
> |
if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) & |
| 291 |
> |
propertyMapI(3, i) = .true. |
| 292 |
> |
if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) & |
| 293 |
> |
propertyMapI(4, i) = .true. |
| 294 |
> |
if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) & |
| 295 |
> |
propertyMapI(5, i) = .true. |
| 296 |
> |
|
| 297 |
> |
end do |
| 298 |
> |
|
| 299 |
> |
#ifdef IS_MPI |
| 300 |
> |
if (.not.allocated(propertyMapJ)) then |
| 301 |
> |
allocate(propertyMapJ(5,ncol)) |
| 302 |
> |
endif |
| 303 |
> |
|
| 304 |
> |
do j = 1, ncol |
| 305 |
> |
me_j = atid_col(j) |
| 306 |
> |
#else |
| 307 |
> |
if (.not.allocated(propertyMapJ)) then |
| 308 |
> |
allocate(propertyMapJ(5,nlocal)) |
| 309 |
> |
endif |
| 310 |
> |
|
| 311 |
> |
do j = 1, natoms |
| 312 |
> |
me_j = atid(j) |
| 313 |
> |
#endif |
| 314 |
> |
|
| 315 |
> |
propertyMapJ(1:5,j) = .false. |
| 316 |
> |
|
| 317 |
> |
call getElementProperty(atypes, me_j, "propertyPack", propPack_j) |
| 318 |
> |
|
| 319 |
> |
! unpack the properties |
| 320 |
> |
|
| 321 |
> |
if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) & |
| 322 |
> |
propertyMapJ(1, j) = .true. |
| 323 |
> |
if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) & |
| 324 |
> |
propertyMapJ(2, j) = .true. |
| 325 |
> |
if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) & |
| 326 |
> |
propertyMapJ(3, j) = .true. |
| 327 |
> |
if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) & |
| 328 |
> |
propertyMapJ(4, j) = .true. |
| 329 |
> |
if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) & |
| 330 |
> |
propertyMapJ(5, j) = .true. |
| 331 |
> |
|
| 332 |
> |
end do |
| 333 |
> |
|
| 334 |
|
! Gather all information needed by all force loops: |
| 335 |
|
|
| 336 |
|
#ifdef IS_MPI |
| 518 |
|
#ifdef IS_MPI |
| 519 |
|
|
| 520 |
|
if (update_nlist) then |
| 447 |
– |
|
| 521 |
|
!! save current configuration, construct neighbor list, |
| 522 |
|
!! and calculate forces |
| 523 |
|
call saveNeighborList(nlocal, q) |
| 526 |
|
nlist = 0 |
| 527 |
|
|
| 528 |
|
do i = 1, nrow |
| 529 |
+ |
|
| 530 |
|
point(i) = nlist + 1 |
| 531 |
|
|
| 532 |
|
inner: do j = 1, ncol |
| 584 |
|
#else |
| 585 |
|
|
| 586 |
|
if (update_nlist) then |
| 587 |
< |
|
| 587 |
> |
|
| 588 |
|
! save current configuration, contruct neighbor list, |
| 589 |
|
! and calculate forces |
| 590 |
|
call saveNeighborList(natoms, q) |
| 771 |
|
tau = tau_Temp |
| 772 |
|
virial = virial_Temp |
| 773 |
|
endif |
| 774 |
< |
|
| 774 |
> |
|
| 775 |
|
#endif |
| 776 |
< |
|
| 777 |
< |
#ifdef PROFILE |
| 778 |
< |
if (do_pot) then |
| 705 |
< |
|
| 706 |
< |
#ifdef IS_MPI |
| 707 |
< |
|
| 708 |
< |
|
| 709 |
< |
call printCommTime() |
| 710 |
< |
|
| 711 |
< |
call mpi_allreduce(forceTime,globalForceTime,1,MPI_DOUBLE_PRECISION, & |
| 712 |
< |
mpi_sum,mpi_comm_world,mpi_err) |
| 713 |
< |
|
| 714 |
< |
call mpi_allreduce(forceTime,maxForceTime,1,MPI_DOUBLE_PRECISION, & |
| 715 |
< |
MPI_MAX,mpi_comm_world,mpi_err) |
| 716 |
< |
|
| 717 |
< |
call mpi_comm_size( MPI_COMM_WORLD, nprocs,mpi_err) |
| 718 |
< |
|
| 719 |
< |
if (getMyNode() == 0) then |
| 720 |
< |
write(*,*) "Total processor time spent in force calculations is: ", globalForceTime |
| 721 |
< |
write(*,*) "Total Time spent in force loop per processor is: ", globalforceTime/nprocs |
| 722 |
< |
write(*,*) "Maximum force time on any processor is: ", maxForceTime |
| 723 |
< |
end if |
| 724 |
< |
#else |
| 725 |
< |
write(*,*) "Time spent in force loop is: ", forceTime |
| 726 |
< |
#endif |
| 727 |
< |
|
| 728 |
< |
|
| 729 |
< |
endif |
| 730 |
< |
|
| 731 |
< |
#endif |
| 732 |
< |
|
| 733 |
< |
|
| 734 |
< |
|
| 776 |
> |
|
| 777 |
> |
|
| 778 |
> |
|
| 779 |
|
end subroutine do_force_loop |
| 780 |
|
|
| 781 |
|
subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot) |
| 797 |
|
logical :: is_EAM_i,is_EAM_j |
| 798 |
|
logical :: is_Sticky_i, is_Sticky_j |
| 799 |
|
integer :: me_i, me_j |
| 800 |
< |
|
| 800 |
> |
integer :: propPack_i |
| 801 |
> |
integer :: propPack_j |
| 802 |
|
r = sqrt(rijsq) |
| 803 |
|
|
| 804 |
|
#ifdef IS_MPI |
| 815 |
|
me_j = atid(j) |
| 816 |
|
|
| 817 |
|
#endif |
| 818 |
< |
|
| 818 |
> |
|
| 819 |
|
if (FF_uses_LJ .and. SimUsesLJ()) then |
| 775 |
– |
call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i) |
| 776 |
– |
call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j) |
| 820 |
|
|
| 821 |
< |
if ( is_LJ_i .and. is_LJ_j ) & |
| 821 |
> |
if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) & |
| 822 |
|
call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) |
| 823 |
+ |
|
| 824 |
|
endif |
| 825 |
|
|
| 826 |
|
if (FF_uses_dipoles .and. SimUsesDipoles()) then |
| 783 |
– |
call getElementProperty(atypes, me_i, "is_DP", is_DP_i) |
| 784 |
– |
call getElementProperty(atypes, me_j, "is_DP", is_DP_j) |
| 827 |
|
|
| 828 |
< |
if ( is_DP_i .and. is_DP_j ) then |
| 787 |
< |
|
| 828 |
> |
if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then |
| 829 |
|
call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, & |
| 830 |
|
do_pot, do_stress) |
| 831 |
|
if (FF_uses_RF .and. SimUsesRF()) then |
| 838 |
|
|
| 839 |
|
if (FF_uses_Sticky .and. SimUsesSticky()) then |
| 840 |
|
|
| 841 |
< |
call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i) |
| 801 |
< |
call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j) |
| 802 |
< |
|
| 803 |
< |
if ( is_Sticky_i .and. is_Sticky_j ) then |
| 841 |
> |
if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then |
| 842 |
|
call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, & |
| 843 |
|
do_pot, do_stress) |
| 844 |
|
endif |
| 846 |
|
|
| 847 |
|
|
| 848 |
|
if (FF_uses_GB .and. SimUsesGB()) then |
| 811 |
– |
|
| 812 |
– |
call getElementProperty(atypes, me_i, "is_GB", is_GB_i) |
| 813 |
– |
call getElementProperty(atypes, me_j, "is_GB", is_GB_j) |
| 849 |
|
|
| 850 |
< |
if ( is_GB_i .and. is_GB_j ) then |
| 850 |
> |
if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then |
| 851 |
|
call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, & |
| 852 |
|
do_pot, do_stress) |
| 853 |
|
endif |
| 854 |
+ |
|
| 855 |
|
endif |
| 856 |
|
|
| 857 |
|
|
| 858 |
|
|
| 859 |
|
if (FF_uses_EAM .and. SimUsesEAM()) then |
| 824 |
– |
call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i) |
| 825 |
– |
call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j) |
| 860 |
|
|
| 861 |
< |
if ( is_EAM_i .and. is_EAM_j ) & |
| 862 |
< |
call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) |
| 863 |
< |
endif |
| 830 |
< |
|
| 831 |
< |
|
| 861 |
> |
if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then |
| 862 |
> |
call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) |
| 863 |
> |
endif |
| 864 |
|
|
| 865 |
+ |
endif |
| 866 |
|
|
| 867 |
|
end subroutine do_pair |
| 868 |
|
|
| 1154 |
|
doesit = FF_uses_RF |
| 1155 |
|
end function FF_RequiresPostpairCalc |
| 1156 |
|
|
| 1157 |
+ |
#ifdef PROFILE |
| 1158 |
+ |
function getforcetime() result(totalforcetime) |
| 1159 |
+ |
real(kind=dp) :: totalforcetime |
| 1160 |
+ |
totalforcetime = forcetime |
| 1161 |
+ |
end function getforcetime |
| 1162 |
+ |
#endif |
| 1163 |
+ |
|
| 1164 |
|
!! This cleans componets of force arrays belonging only to fortran |
| 1165 |
|
|
| 1166 |
|
end module do_Forces |
