| 15 |
|
logical, save :: EAM_FF_initialized = .false. |
| 16 |
|
integer, save :: EAM_Mixing_Policy |
| 17 |
|
real(kind = dp), save :: EAM_rcut |
| 18 |
< |
real(kind = dp), save :: EAM_rcut_orig |
| 18 |
> |
logical, save :: haveRcut = .false. |
| 19 |
|
|
| 20 |
|
character(len = statusMsgSize) :: errMesg |
| 21 |
|
integer :: eam_err |
| 51 |
|
|
| 52 |
|
|
| 53 |
|
!! Arrays for derivatives used in force calculation |
| 54 |
< |
real( kind = dp), dimension(:), allocatable :: frho |
| 55 |
< |
real( kind = dp), dimension(:), allocatable :: rho |
| 54 |
> |
real( kind = dp),save, dimension(:), allocatable :: frho |
| 55 |
> |
real( kind = dp),save, dimension(:), allocatable :: rho |
| 56 |
|
|
| 57 |
< |
real( kind = dp), dimension(:), allocatable :: dfrhodrho |
| 58 |
< |
real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho |
| 57 |
> |
real( kind = dp),save, dimension(:), allocatable :: dfrhodrho |
| 58 |
> |
real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho |
| 59 |
|
|
| 60 |
|
|
| 61 |
|
!! Arrays for MPI storage |
| 62 |
|
#ifdef IS_MPI |
| 63 |
< |
real( kind = dp), dimension(:), allocatable :: dfrhodrho_col |
| 64 |
< |
real( kind = dp), dimension(:), allocatable :: dfrhodrho_row |
| 65 |
< |
real( kind = dp), dimension(:), allocatable :: frho_row |
| 66 |
< |
real( kind = dp), dimension(:), allocatable :: frho_col |
| 67 |
< |
real( kind = dp), dimension(:), allocatable :: rho_row |
| 68 |
< |
real( kind = dp), dimension(:), allocatable :: rho_col |
| 69 |
< |
real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho_col |
| 70 |
< |
real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho_row |
| 63 |
> |
real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col |
| 64 |
> |
real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_row |
| 65 |
> |
real( kind = dp),save, dimension(:), allocatable :: frho_row |
| 66 |
> |
real( kind = dp),save, dimension(:), allocatable :: frho_col |
| 67 |
> |
real( kind = dp),save, dimension(:), allocatable :: rho_row |
| 68 |
> |
real( kind = dp),save, dimension(:), allocatable :: rho_col |
| 69 |
> |
real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_col |
| 70 |
> |
real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row |
| 71 |
|
#endif |
| 72 |
|
|
| 73 |
|
type, private :: EAMTypeList |
| 117 |
|
|
| 118 |
|
status = 0 |
| 119 |
|
|
| 120 |
< |
write(*,*) "Adding new eamtype: ",eam_ident |
| 120 |
> |
|
| 121 |
|
!! Assume that atypes has already been set and get the total number of types in atypes |
| 122 |
|
!! Also assume that every member of atypes is a EAM model. |
| 123 |
|
|
| 222 |
|
enddo |
| 223 |
|
|
| 224 |
|
|
| 225 |
< |
current_rcut_max = EAMList%EAMParams(1)%eam_rcut |
| 226 |
< |
!! find the smallest rcut for any eam atype |
| 227 |
< |
do i = 2, EAMList%currentAddition |
| 228 |
< |
current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut) |
| 229 |
< |
end do |
| 225 |
> |
! current_rcut_max = EAMList%EAMParams(1)%eam_rcut |
| 226 |
> |
!! find the smallest rcut for any eam atype |
| 227 |
> |
! do i = 2, EAMList%currentAddition |
| 228 |
> |
! current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut) |
| 229 |
> |
! end do |
| 230 |
|
|
| 231 |
< |
EAM_rcut = current_rcut_max |
| 232 |
< |
EAM_rcut_orig = current_rcut_max |
| 231 |
> |
! EAM_rcut = current_rcut_max |
| 232 |
> |
! EAM_rcut_orig = current_rcut_max |
| 233 |
|
! do i = 1, EAMList%currentAddition |
| 234 |
|
! EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i |
| 235 |
|
! end do |
| 350 |
|
|
| 351 |
|
end subroutine allocateEAM |
| 352 |
|
|
| 353 |
+ |
!! C sets rcut to be the largest cutoff of any atype |
| 354 |
+ |
!! present in this simulation. Doesn't include all atypes |
| 355 |
+ |
!! sim knows about, just those in the simulation. |
| 356 |
|
subroutine setCutoffEAM(rcut, status) |
| 357 |
|
real(kind=dp) :: rcut |
| 358 |
|
integer :: status |
| 359 |
+ |
status = 0 |
| 360 |
|
|
| 361 |
< |
if (rcut < EAM_rcut) then |
| 358 |
< |
EAM_rcut = rcut |
| 359 |
< |
endif |
| 361 |
> |
EAM_rcut = rcut |
| 362 |
|
|
| 361 |
– |
|
| 363 |
|
end subroutine setCutoffEAM |
| 364 |
|
|
| 365 |
|
|
| 366 |
|
|
| 367 |
|
subroutine clean_EAM() |
| 368 |
< |
|
| 368 |
< |
! clean non-IS_MPI first |
| 368 |
> |
! clean non-IS_MPI first |
| 369 |
|
frho = 0.0_dp |
| 370 |
|
rho = 0.0_dp |
| 371 |
|
dfrhodrho = 0.0_dp |
| 484 |
|
integer :: myid_atom2 |
| 485 |
|
|
| 486 |
|
! check to see if we need to be cleaned at the start of a force loop |
| 487 |
– |
if (cleanme) call clean_EAM |
| 488 |
– |
cleanme = .false. |
| 487 |
|
|
| 488 |
+ |
if (cleanme) then |
| 489 |
+ |
call clean_EAM |
| 490 |
+ |
cleanme = .false. |
| 491 |
+ |
end if |
| 492 |
+ |
|
| 493 |
|
|
| 494 |
+ |
|
| 495 |
+ |
|
| 496 |
|
#ifdef IS_MPI |
| 497 |
|
myid_atom1 = atid_Row(atom1) |
| 498 |
|
myid_atom2 = atid_Col(atom2) |
| 537 |
|
#endif |
| 538 |
|
endif |
| 539 |
|
|
| 540 |
+ |
|
| 541 |
+ |
|
| 542 |
|
end subroutine calc_eam_prepair_rho |
| 543 |
|
|
| 544 |
|
|
| 556 |
|
integer :: n_rho_points |
| 557 |
|
! reset clean forces to be true at top of calc rho. |
| 558 |
|
cleanme = .true. |
| 559 |
< |
|
| 559 |
> |
|
| 560 |
|
!! Scatter the electron density from pre-pair calculation back to local atoms |
| 561 |
|
#ifdef IS_MPI |
| 562 |
|
call scatter(rho_row,rho,plan_row,eam_err) |
| 601 |
|
pot = pot + u |
| 602 |
|
enddo |
| 603 |
|
|
| 604 |
+ |
|
| 605 |
|
|
| 598 |
– |
|
| 606 |
|
#ifdef IS_MPI |
| 607 |
|
!! communicate f(rho) and derivatives back into row and column arrays |
| 608 |
|
call gather(frho,frho_row,plan_row, eam_err) |
| 668 |
|
!Local Variables |
| 669 |
|
|
| 670 |
|
|
| 671 |
< |
|
| 671 |
> |
|
| 672 |
|
phab = 0.0E0_DP |
| 673 |
|
dvpdr = 0.0E0_DP |
| 674 |
|
d2vpdrdr = 0.0E0_DP |
| 675 |
< |
|
| 669 |
< |
|
| 675 |
> |
|
| 676 |
|
if (rij .lt. EAM_rcut) then |
| 677 |
|
#ifdef IS_MPI |
| 678 |
|
!!!!! FIX ME |
| 759 |
|
#ifdef IS_MPI |
| 760 |
|
dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) & |
| 761 |
|
+ dvpdr |
| 762 |
< |
|
| 762 |
> |
|
| 763 |
|
#else |
| 764 |
|
dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) & |
| 765 |
|
+ dvpdr |
| 766 |
< |
|
| 766 |
> |
! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2) |
| 767 |
|
#endif |
| 768 |
|
|
| 769 |
|
fx = dudr * drdx |
| 785 |
|
f_Col(2,atom2) = f_Col(2,atom2) - fy |
| 786 |
|
f_Col(3,atom2) = f_Col(3,atom2) - fz |
| 787 |
|
#else |
| 788 |
< |
if(do_pot) pot = pot + phab |
| 789 |
< |
|
| 788 |
> |
|
| 789 |
> |
if(do_pot) then |
| 790 |
> |
pot = pot + phab |
| 791 |
> |
end if |
| 792 |
> |
|
| 793 |
|
f(1,atom1) = f(1,atom1) + fx |
| 794 |
|
f(2,atom1) = f(2,atom1) + fy |
| 795 |
|
f(3,atom1) = f(3,atom1) + fz |