| 689 |
|
d2vpdrdr = 0.0E0_DP |
| 690 |
|
|
| 691 |
|
if (rij .lt. EAM_rcut) then |
| 692 |
+ |
|
| 693 |
|
#ifdef IS_MPI |
| 693 |
– |
!!!!! FIX ME |
| 694 |
|
mytype_atom1 = atid_row(atom1) |
| 695 |
+ |
mytype_atom2 = atid_col(atom2) |
| 696 |
|
#else |
| 697 |
|
mytype_atom1 = atid(atom1) |
| 698 |
+ |
mytype_atom2 = atid(atom2) |
| 699 |
|
#endif |
| 700 |
+ |
! get cutoff for atom 1 |
| 701 |
+ |
rci = EAMList%EAMParams(mytype_atom1)%eam_rcut |
| 702 |
+ |
! get type specific cutoff for atom 2 |
| 703 |
+ |
rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut |
| 704 |
|
|
| 705 |
|
drdx = d(1)/rij |
| 706 |
|
drdy = d(2)/rij |
| 707 |
|
drdz = d(3)/rij |
| 708 |
|
|
| 709 |
< |
|
| 710 |
< |
call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, & |
| 709 |
> |
if (rij.lt.rci) then |
| 710 |
> |
call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, & |
| 711 |
|
EAMList%EAMParams(mytype_atom1)%eam_rvals, & |
| 712 |
|
EAMList%EAMParams(mytype_atom1)%eam_rho_r, & |
| 713 |
|
EAMList%EAMParams(mytype_atom1)%eam_rho_r_pp, & |
| 714 |
|
rij, rha,drha,d2rha) |
| 715 |
< |
|
| 716 |
< |
!! Calculate Phi(r) for atom1. |
| 711 |
< |
call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, & |
| 715 |
> |
!! Calculate Phi(r) for atom1. |
| 716 |
> |
call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, & |
| 717 |
|
EAMList%EAMParams(mytype_atom1)%eam_rvals, & |
| 718 |
|
EAMList%EAMParams(mytype_atom1)%eam_phi_r, & |
| 719 |
|
EAMList%EAMParams(mytype_atom1)%eam_phi_r_pp, & |
| 720 |
|
rij, pha,dpha,d2pha) |
| 721 |
+ |
endif |
| 722 |
|
|
| 723 |
< |
|
| 724 |
< |
! get cutoff for atom 1 |
| 725 |
< |
rci = EAMList%EAMParams(mytype_atom1)%eam_rcut |
| 720 |
< |
#ifdef IS_MPI |
| 721 |
< |
mytype_atom2 = atid_col(atom2) |
| 722 |
< |
#else |
| 723 |
< |
mytype_atom2 = atid(atom2) |
| 724 |
< |
#endif |
| 725 |
< |
|
| 726 |
< |
! Calculate rho,drho and d2rho for atom1 |
| 727 |
< |
call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, & |
| 723 |
> |
if (rij.lt.rcj) then |
| 724 |
> |
! Calculate rho,drho and d2rho for atom1 |
| 725 |
> |
call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, & |
| 726 |
|
EAMList%EAMParams(mytype_atom2)%eam_rvals, & |
| 727 |
|
EAMList%EAMParams(mytype_atom2)%eam_rho_r, & |
| 728 |
|
EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, & |
| 729 |
|
rij, rhb,drhb,d2rhb) |
| 730 |
< |
|
| 731 |
< |
!! Calculate Phi(r) for atom2. |
| 732 |
< |
call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, & |
| 730 |
> |
|
| 731 |
> |
!! Calculate Phi(r) for atom2. |
| 732 |
> |
call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, & |
| 733 |
|
EAMList%EAMParams(mytype_atom2)%eam_rvals, & |
| 734 |
|
EAMList%EAMParams(mytype_atom2)%eam_phi_r, & |
| 735 |
|
EAMList%EAMParams(mytype_atom2)%eam_phi_r_pp, & |
| 736 |
|
rij, phb,dphb,d2phb) |
| 737 |
+ |
endif |
| 738 |
|
|
| 740 |
– |
|
| 741 |
– |
! get type specific cutoff for atom 2 |
| 742 |
– |
rcj = EAMList%EAMParams(mytype_atom1)%eam_rcut |
| 743 |
– |
|
| 744 |
– |
|
| 745 |
– |
|
| 739 |
|
if (rij.lt.rci) then |
| 740 |
|
phab = phab + 0.5E0_DP*(rhb/rha)*pha |
| 741 |
|
dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + & |
| 746 |
|
(2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha))) |
| 747 |
|
endif |
| 748 |
|
|
| 756 |
– |
|
| 749 |
|
if (rij.lt.rcj) then |
| 750 |
|
phab = phab + 0.5E0_DP*(rha/rhb)*phb |
| 751 |
|
dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + & |
| 872 |
|
|
| 873 |
|
! check to make sure we're inside the spline range: |
| 874 |
|
if ((j.gt.nx).or.(j.lt.1)) then |
| 875 |
< |
write(errMSG,*) "EAM_splint: x is outside bounds of spline" |
| 875 |
> |
write(errMSG,*) "EAM_splint: x is outside bounds of spline: ",x,j |
| 876 |
|
call handleError(routineName,errMSG) |
| 877 |
|
endif |
| 878 |
|
! check to make sure we haven't screwed up the calculation of j: |
