66 |
|
logical, save,:: haveMixingMap = .false. |
67 |
|
|
68 |
|
character(len = statusMsgSize) :: errMesg |
69 |
< |
integer :: eam_err |
69 |
> |
integer :: sc_err |
70 |
|
|
71 |
|
character(len = 200) :: errMsg |
72 |
|
character(len=*), parameter :: RoutineName = "Sutton-Chen MODULE" |
174 |
|
|
175 |
|
|
176 |
|
! check to see if this is the first time into |
177 |
< |
if (.not.associated(EAMList%EAMParams)) then |
177 |
> |
if (.not.associated(SCTypeList%EAMParams)) then |
178 |
|
call getMatchingElementList(atypes, "is_SuttonChen", .true., nSCtypes, MatchList) |
179 |
|
SCTypeList%nSCtypes = nSCtypes |
180 |
|
allocate(SCTypeList%SCTypes(nSCTypes)) |
215 |
|
integer :: eamID |
216 |
|
real(kind=dp) :: cutValue |
217 |
|
|
218 |
< |
eamID = EAMList%atidToEAMType(atomID) |
219 |
< |
cutValue = EAMList%EAMParams(eamID)%eam_rcut |
218 |
> |
eamID = SCTypeList%atidToEAMType(atomID) |
219 |
> |
cutValue = SCTypeList%EAMParams(eamID)%eam_rcut |
220 |
|
end function getSCCut |
221 |
|
|
222 |
|
|
432 |
|
|
433 |
|
! we don't use the derivatives, dummy variables |
434 |
|
real( kind = dp) :: drho,d2rho |
435 |
< |
integer :: eam_err |
435 |
> |
integer :: sc_err |
436 |
|
|
437 |
|
integer :: atid1,atid2 ! Global atid |
438 |
|
integer :: myid_atom1 ! EAM atid |
469 |
|
rho(atom1) = rho(atom1) + rho_j_at_i |
470 |
|
#endif |
471 |
|
|
472 |
– |
|
473 |
– |
|
472 |
|
end subroutine calc_sc_prepair_rho |
473 |
|
|
474 |
|
|
475 |
< |
!! Calculate the functional F(rho) for all local atoms |
475 |
> |
!! Calculate the rho_a for all local atoms |
476 |
|
subroutine calc_eam_preforce_Frho(nlocal,pot) |
477 |
|
integer :: nlocal |
478 |
|
real(kind=dp) :: pot |
487 |
|
cleanme = .true. |
488 |
|
!! Scatter the electron density from pre-pair calculation back to local atoms |
489 |
|
#ifdef IS_MPI |
490 |
< |
call scatter(rho_row,rho,plan_atom_row,eam_err) |
491 |
< |
if (eam_err /= 0 ) then |
490 |
> |
call scatter(rho_row,rho,plan_atom_row,sc_err) |
491 |
> |
if (sc_err /= 0 ) then |
492 |
|
write(errMsg,*) " Error scattering rho_row into rho" |
493 |
|
call handleError(RoutineName,errMesg) |
494 |
|
endif |
495 |
< |
call scatter(rho_col,rho_tmp,plan_atom_col,eam_err) |
496 |
< |
if (eam_err /= 0 ) then |
495 |
> |
call scatter(rho_col,rho_tmp,plan_atom_col,sc_err) |
496 |
> |
if (sc_err /= 0 ) then |
497 |
|
write(errMsg,*) " Error scattering rho_col into rho" |
498 |
|
call handleError(RoutineName,errMesg) |
499 |
|
|
516 |
|
|
517 |
|
#ifdef IS_MPI |
518 |
|
!! communicate f(rho) and derivatives back into row and column arrays |
519 |
< |
call gather(frho,frho_row,plan_atom_row, eam_err) |
520 |
< |
if (eam_err /= 0) then |
519 |
> |
call gather(frho,frho_row,plan_atom_row, sc_err) |
520 |
> |
if (sc_err /= 0) then |
521 |
|
call handleError("cal_eam_forces()","MPI gather frho_row failure") |
522 |
|
endif |
523 |
< |
call gather(dfrhodrho,dfrhodrho_row,plan_atom_row, eam_err) |
524 |
< |
if (eam_err /= 0) then |
523 |
> |
call gather(dfrhodrho,dfrhodrho_row,plan_atom_row, sc_err) |
524 |
> |
if (sc_err /= 0) then |
525 |
|
call handleError("cal_eam_forces()","MPI gather dfrhodrho_row failure") |
526 |
|
endif |
527 |
< |
call gather(frho,frho_col,plan_atom_col, eam_err) |
528 |
< |
if (eam_err /= 0) then |
527 |
> |
call gather(frho,frho_col,plan_atom_col, sc_err) |
528 |
> |
if (sc_err /= 0) then |
529 |
|
call handleError("cal_eam_forces()","MPI gather frho_col failure") |
530 |
|
endif |
531 |
< |
call gather(dfrhodrho,dfrhodrho_col,plan_atom_col, eam_err) |
532 |
< |
if (eam_err /= 0) then |
531 |
> |
call gather(dfrhodrho,dfrhodrho_col,plan_atom_col, sc_err) |
532 |
> |
if (sc_err /= 0) then |
533 |
|
call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure") |
534 |
|
endif |
535 |
|
|
562 |
|
real( kind = dp ) :: rha,drha,d2rha, dpha |
563 |
|
real( kind = dp ) :: rhb,drhb,d2rhb, dphb |
564 |
|
real( kind = dp ) :: dudr |
565 |
< |
real( kind = dp ) :: rci,rcj |
565 |
> |
real( kind = dp ) :: rcij |
566 |
|
real( kind = dp ) :: drhoidr,drhojdr |
567 |
|
real( kind = dp ) :: d2rhoidrdr |
568 |
|
real( kind = dp ) :: d2rhojdrdr |
581 |
|
dvpdr = 0.0E0_DP |
582 |
|
d2vpdrdr = 0.0E0_DP |
583 |
|
|
586 |
– |
if (rij .lt. EAM_rcut) then |
584 |
|
|
585 |
|
#ifdef IS_MPI |
586 |
|
atid1 = atid_row(atom1) |
590 |
|
atid2 = atid(atom2) |
591 |
|
#endif |
592 |
|
|
593 |
< |
mytype_atom1 = EAMList%atidToEAMType(atid1) |
594 |
< |
mytype_atom2 = EAMList%atidTOEAMType(atid2) |
593 |
> |
mytype_atom1 = SCTypeList%atidToSCType(atid1) |
594 |
> |
mytype_atom2 = SCTypeList%atidTOSCType(atid2) |
595 |
|
|
599 |
– |
|
600 |
– |
! get cutoff for atom 1 |
601 |
– |
rci = EAMList%EAMParams(mytype_atom1)%eam_rcut |
602 |
– |
! get type specific cutoff for atom 2 |
603 |
– |
rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut |
604 |
– |
|
596 |
|
drdx = d(1)/rij |
597 |
|
drdy = d(2)/rij |
598 |
|
drdz = d(3)/rij |
599 |
|
|
600 |
< |
if (rij.lt.rci) then |
601 |
< |
call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, & |
602 |
< |
EAMList%EAMParams(mytype_atom1)%eam_rvals, & |
603 |
< |
EAMList%EAMParams(mytype_atom1)%eam_rho_r, & |
604 |
< |
EAMList%EAMParams(mytype_atom1)%eam_rho_r_pp, & |
605 |
< |
rij, rha,drha,d2rha) |
615 |
< |
!! Calculate Phi(r) for atom1. |
616 |
< |
call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, & |
617 |
< |
EAMList%EAMParams(mytype_atom1)%eam_rvals, & |
618 |
< |
EAMList%EAMParams(mytype_atom1)%eam_phi_r, & |
619 |
< |
EAMList%EAMParams(mytype_atom1)%eam_phi_r_pp, & |
620 |
< |
rij, pha,dpha,d2pha) |
621 |
< |
endif |
622 |
< |
|
623 |
< |
if (rij.lt.rcj) then |
624 |
< |
! Calculate rho,drho and d2rho for atom1 |
625 |
< |
call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, & |
626 |
< |
EAMList%EAMParams(mytype_atom2)%eam_rvals, & |
627 |
< |
EAMList%EAMParams(mytype_atom2)%eam_rho_r, & |
628 |
< |
EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, & |
629 |
< |
rij, rhb,drhb,d2rhb) |
600 |
> |
|
601 |
> |
epsilonij = MixingMap(mytype_atom1,mytype_atom2)%epsilon |
602 |
> |
aij = MixingMap(mytype_atom1,mytype_atom2)%alpha |
603 |
> |
nij = MixingMap(mytype_atom1,mytype_atom2)%n |
604 |
> |
mij = MixingMap(mytype_atom1,mytype_atom2)%m |
605 |
> |
vcij = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot |
606 |
|
|
607 |
< |
!! Calculate Phi(r) for atom2. |
632 |
< |
call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, & |
633 |
< |
EAMList%EAMParams(mytype_atom2)%eam_rvals, & |
634 |
< |
EAMList%EAMParams(mytype_atom2)%eam_phi_r, & |
635 |
< |
EAMList%EAMParams(mytype_atom2)%eam_phi_r_pp, & |
636 |
< |
rij, phb,dphb,d2phb) |
637 |
< |
endif |
607 |
> |
vptmp = dij*((aij/r)**nij) |
608 |
|
|
639 |
– |
if (rij.lt.rci) then |
640 |
– |
phab = phab + 0.5E0_DP*(rhb/rha)*pha |
641 |
– |
dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + & |
642 |
– |
pha*((drhb/rha) - (rhb*drha/rha/rha))) |
643 |
– |
d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rhb/rha)*d2pha + & |
644 |
– |
2.0E0_DP*dpha*((drhb/rha) - (rhb*drha/rha/rha)) + & |
645 |
– |
pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + & |
646 |
– |
(2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha))) |
647 |
– |
endif |
609 |
|
|
610 |
< |
if (rij.lt.rcj) then |
611 |
< |
phab = phab + 0.5E0_DP*(rha/rhb)*phb |
612 |
< |
dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + & |
613 |
< |
phb*((drha/rhb) - (rha*drhb/rhb/rhb))) |
614 |
< |
d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rha/rhb)*d2phb + & |
615 |
< |
2.0E0_DP*dphb*((drha/rhb) - (rha*drhb/rhb/rhb)) + & |
616 |
< |
phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + & |
617 |
< |
(2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb))) |
618 |
< |
endif |
610 |
> |
dvpdr = -nij*vptmp/r |
611 |
> |
d2vpdrdr = -dvpdr*(nij+1.0d0)/r |
612 |
> |
|
613 |
> |
drhodr = -mij*((aij/r)**mij)/r |
614 |
> |
d2rhodrdr = -drhodr*(mij+1.0d0)/r |
615 |
> |
|
616 |
> |
dudr = drhodr*(dfrhodrho(i)+dfrhodrho(j)) & |
617 |
> |
+ dvpdr |
618 |
> |
|
619 |
> |
|
620 |
|
|
659 |
– |
drhoidr = drha |
660 |
– |
drhojdr = drhb |
661 |
– |
|
662 |
– |
d2rhoidrdr = d2rha |
663 |
– |
d2rhojdrdr = d2rhb |
664 |
– |
|
665 |
– |
|
621 |
|
#ifdef IS_MPI |
622 |
< |
dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) & |
622 |
> |
dudr = drhodr*(dfrhodrho_row(atom1)+dfrhodrho_col(atom2)) & |
623 |
|
+ dvpdr |
624 |
|
|
625 |
|
#else |
626 |
< |
dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) & |
626 |
> |
dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) & |
627 |
|
+ dvpdr |
673 |
– |
! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2) |
628 |
|
#endif |
629 |
|
|
630 |
+ |
|
631 |
|
fx = dudr * drdx |
632 |
|
fy = dudr * drdy |
633 |
|
fz = dudr * drdz |