44 |
|
use force_globals |
45 |
|
use status |
46 |
|
use atype_module |
47 |
< |
use Vector_class |
47 |
> |
use vector_class |
48 |
|
#ifdef IS_MPI |
49 |
|
use mpiSimulation |
50 |
|
#endif |
51 |
|
implicit none |
52 |
|
PRIVATE |
53 |
+ |
#define __FORTRAN90 |
54 |
+ |
#include "UseTheForce/DarkSide/fInteractionMap.h" |
55 |
|
|
56 |
|
INTEGER, PARAMETER :: DP = selected_real_kind(15) |
57 |
|
|
119 |
|
integer :: currentAddition = 0 |
120 |
|
|
121 |
|
type (EAMtype), pointer :: EAMParams(:) => null() |
122 |
+ |
integer, pointer :: atidToEAMType(:) => null() |
123 |
|
end type EAMTypeList |
124 |
|
|
125 |
|
|
136 |
|
public :: calc_eam_preforce_Frho |
137 |
|
public :: clean_EAM |
138 |
|
public :: destroyEAMTypes |
139 |
+ |
public :: getEAMCut |
140 |
|
|
141 |
|
contains |
142 |
|
|
143 |
|
|
144 |
|
subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,& |
145 |
|
eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,& |
146 |
< |
eam_ident,status) |
146 |
> |
c_ident,status) |
147 |
|
real (kind = dp ) :: lattice_constant |
148 |
|
integer :: eam_nrho |
149 |
|
real (kind = dp ) :: eam_drho |
153 |
|
real (kind = dp ), dimension(eam_nr) :: eam_Z_r |
154 |
|
real (kind = dp ), dimension(eam_nr) :: eam_rho_r |
155 |
|
real (kind = dp ), dimension(eam_nrho) :: eam_F_rho |
156 |
< |
integer :: eam_ident |
156 |
> |
integer :: c_ident |
157 |
|
integer :: status |
158 |
|
|
159 |
< |
integer :: nAtypes |
159 |
> |
integer :: nAtypes,nEAMTypes,myATID |
160 |
|
integer :: maxVals |
161 |
|
integer :: alloc_stat |
162 |
|
integer :: current |
171 |
|
|
172 |
|
! check to see if this is the first time into |
173 |
|
if (.not.associated(EAMList%EAMParams)) then |
174 |
< |
call getMatchingElementList(atypes, "is_EAM", .true., nAtypes, MatchList) |
175 |
< |
EAMList%n_eam_types = nAtypes |
176 |
< |
allocate(EAMList%EAMParams(nAtypes)) |
174 |
> |
call getMatchingElementList(atypes, "is_EAM", .true., nEAMtypes, MatchList) |
175 |
> |
EAMList%n_eam_types = nEAMtypes |
176 |
> |
allocate(EAMList%EAMParams(nEAMTypes)) |
177 |
> |
nAtypes = getSize(atypes) |
178 |
> |
allocate(EAMList%atidToEAMType(nAtypes)) |
179 |
|
end if |
180 |
|
|
181 |
|
EAMList%currentAddition = EAMList%currentAddition + 1 |
182 |
|
current = EAMList%currentAddition |
183 |
|
|
184 |
+ |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
185 |
+ |
EAMList%atidToEAMType(myATID) = current |
186 |
|
|
187 |
|
call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat) |
188 |
|
if (alloc_stat /= 0) then |
190 |
|
return |
191 |
|
end if |
192 |
|
|
193 |
< |
! this is a possible bug, we assume a correspondence between the vector atypes and |
194 |
< |
! EAMAtypes |
187 |
< |
|
188 |
< |
EAMList%EAMParams(current)%eam_atype = eam_ident |
193 |
> |
|
194 |
> |
EAMList%EAMParams(current)%eam_atype = c_ident |
195 |
|
EAMList%EAMParams(current)%eam_lattice = lattice_constant |
196 |
|
EAMList%EAMParams(current)%eam_nrho = eam_nrho |
197 |
|
EAMList%EAMParams(current)%eam_drho = eam_drho |
221 |
|
eamList%currentAddition = 0 |
222 |
|
|
223 |
|
end subroutine destroyEAMtypes |
224 |
+ |
|
225 |
+ |
function getEAMCut(atomID) result(cutValue) |
226 |
+ |
integer, intent(in) :: atomID |
227 |
+ |
integer :: eamID |
228 |
+ |
real(kind=dp) :: cutValue |
229 |
+ |
|
230 |
+ |
eamID = EAMList%atidToEAMType(atomID) |
231 |
+ |
cutValue = EAMList%EAMParams(eamID)%eam_rcut |
232 |
+ |
end function getEAMCut |
233 |
|
|
234 |
|
subroutine init_EAM_FF(status) |
235 |
|
integer :: status |
561 |
|
! we don't use the derivatives, dummy variables |
562 |
|
real( kind = dp) :: drho,d2rho |
563 |
|
integer :: eam_err |
564 |
+ |
|
565 |
+ |
integer :: atid1,atid2 ! Global atid |
566 |
+ |
integer :: myid_atom1 ! EAM atid |
567 |
+ |
integer :: myid_atom2 |
568 |
|
|
550 |
– |
integer :: myid_atom1 |
551 |
– |
integer :: myid_atom2 |
569 |
|
|
570 |
|
! check to see if we need to be cleaned at the start of a force loop |
571 |
|
|
573 |
|
|
574 |
|
|
575 |
|
#ifdef IS_MPI |
576 |
< |
myid_atom1 = atid_Row(atom1) |
577 |
< |
myid_atom2 = atid_Col(atom2) |
576 |
> |
Atid1 = Atid_row(Atom1) |
577 |
> |
Atid2 = Atid_col(Atom2) |
578 |
|
#else |
579 |
< |
myid_atom1 = atid(atom1) |
580 |
< |
myid_atom2 = atid(atom2) |
579 |
> |
Atid1 = Atid(Atom1) |
580 |
> |
Atid2 = Atid(Atom2) |
581 |
|
#endif |
582 |
|
|
583 |
+ |
Myid_atom1 = Eamlist%atidtoeamtype(Atid1) |
584 |
+ |
Myid_atom2 = Eamlist%atidtoeamtype(Atid2) |
585 |
+ |
|
586 |
|
if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then |
587 |
|
|
588 |
|
|
600 |
|
#else |
601 |
|
rho(atom2) = rho(atom2) + rho_i_at_j |
602 |
|
#endif |
603 |
< |
! write(*,*) atom1,atom2,r,rho_i_at_j |
603 |
> |
! write(*,*) atom1,atom2,r,rho_i_at_j |
604 |
|
endif |
605 |
|
|
606 |
|
if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then |
638 |
|
integer :: atom |
639 |
|
real(kind=dp) :: U,U1,U2 |
640 |
|
integer :: atype1 |
641 |
< |
integer :: me |
641 |
> |
integer :: me,atid1 |
642 |
|
integer :: n_rho_points |
643 |
|
|
644 |
|
|
663 |
|
|
664 |
|
!! Calculate F(rho) and derivative |
665 |
|
do atom = 1, nlocal |
666 |
< |
me = atid(atom) |
666 |
> |
atid1 = atid(atom) |
667 |
> |
me = eamList%atidToEAMtype(atid1) |
668 |
|
n_rho_points = EAMList%EAMParams(me)%eam_nrho |
669 |
|
! Check to see that the density is not greater than the larges rho we have calculated |
670 |
|
if (rho(atom) < EAMList%EAMParams(me)%eam_rhovals(n_rho_points)) then |
758 |
|
integer :: id1,id2 |
759 |
|
integer :: mytype_atom1 |
760 |
|
integer :: mytype_atom2 |
761 |
< |
|
761 |
> |
integer :: atid1,atid2 |
762 |
|
!Local Variables |
763 |
|
|
764 |
|
! write(*,*) "Frho: ", Frho(atom1) |
770 |
|
if (rij .lt. EAM_rcut) then |
771 |
|
|
772 |
|
#ifdef IS_MPI |
773 |
< |
mytype_atom1 = atid_row(atom1) |
774 |
< |
mytype_atom2 = atid_col(atom2) |
773 |
> |
atid1 = atid_row(atom1) |
774 |
> |
atid2 = atid_col(atom2) |
775 |
|
#else |
776 |
< |
mytype_atom1 = atid(atom1) |
777 |
< |
mytype_atom2 = atid(atom2) |
776 |
> |
atid1 = atid(atom1) |
777 |
> |
atid2 = atid(atom2) |
778 |
|
#endif |
779 |
+ |
|
780 |
+ |
mytype_atom1 = EAMList%atidToEAMType(atid1) |
781 |
+ |
mytype_atom2 = EAMList%atidTOEAMType(atid2) |
782 |
+ |
|
783 |
+ |
|
784 |
|
! get cutoff for atom 1 |
785 |
|
rci = EAMList%EAMParams(mytype_atom1)%eam_rcut |
786 |
|
! get type specific cutoff for atom 2 |
864 |
|
|
865 |
|
#ifdef IS_MPI |
866 |
|
if (do_pot) then |
867 |
< |
pot_Row(atom1) = pot_Row(atom1) + phab*0.5 |
868 |
< |
pot_Col(atom2) = pot_Col(atom2) + phab*0.5 |
867 |
> |
pot_Row(METALLIC_POT,atom1) = pot_Row(METALLIC_POT,atom1) + phab*0.5 |
868 |
> |
pot_Col(METALLIC_POT,atom2) = pot_Col(METALLIC_POT,atom2) + phab*0.5 |
869 |
|
end if |
870 |
|
|
871 |
|
f_Row(1,atom1) = f_Row(1,atom1) + fx |