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: |