module mod_atomistic
use mod_global
use mod_fsge, only : Fsge
type NeighObject_type
integer(c_int),allocatable :: neighborList(:,:)
end type NeighObject_type
type :: Atomistic
real(c_double),allocatable :: coords(:,:)
integer(c_int) :: numParticles,numContributingParticles,numSpecies,&
processdEdrFlag,processd2Edr2Flag
type(neighObject_type) :: neighObject
integer(c_int),allocatable :: particleSpecies(:),particle_contributing(:)
real(c_double) :: boxSideLengths(DIM),cutoff,energy
type(Fsge) :: fsgeType,fsgeTypeNumerical
end type Atomistic
type(Atomistic), pointer :: pA
contains
subroutine constructAtomistic(aS,numParticles)
use mod_fsge, only : initializeFsge
implicit none
integer(c_int), intent(in) :: numParticles
type(Atomistic),intent(inout) :: aS
aS%numParticles = numParticles
! Allocate memory to atomistic system
allocate(aS%coords(DIM,numParticles))
allocate(aS%particleSpecies(numParticles))
allocate(aS%particle_contributing(numParticles))
allocate(aS%neighObject%neighborList(numParticles+1, numParticles))
call initializeFsge(aS%fsgeType)
call initializeFsge(aS%fsgeTypeNumerical)
end subroutine constructAtomistic
subroutine process_dedr(data_object,deidr,r,pdx,i,j,ierr) bind(c)
use mod_matrix
implicit none
! Passed variables
type(c_ptr), intent(in),value :: data_object,pdx
real(c_double),intent(in),value :: deidr,r
integer(c_int),intent(in),value :: i,j
integer(c_int),intent(out) :: ierr
! Local variables
real(c_double), pointer :: dx(:)
real(c_double) :: R3(DIM)
call c_f_pointer(pdx,dx,[DIM])
R3 = (pA%coords(:,i) + pA%coords(:,j))/2._cd
pA%fsgeType%Tijm = pA%fsgeType%Tijm + deidr * tensor_prod(dx,dx,R3)/r
pA%fsgeType%Eijkln = pA%fsgeType%Eijkln - deidr*tensor_prod(dx,dx,dx,dx,R3)/(r**3)
pA%fsgeType%Cijkl = pA%fsgeType%Cijkl - deidr * tensor_prod(dx,dx,dx,dx)/(r**3)
pA%fsgeType%Dijmkln = pA%fsgeType%Dijmkln - deidr * tensor_prod(dx,dx,R3,dx,dx,R3)/(r**3)
pA%fsgeType%Dijmkln = pA%fsgeType%Dijmkln + deidr * tensorG_dot_tensorG(tensorG(dx,R3)) / r
ierr = 0
return
end subroutine process_dedr
subroutine process_d2edr2(data_object,d2edr2,pr,pdx,pi,pj,ierr) bind(c)
use mod_matrix
implicit none
! Passed variables
type(c_ptr), intent(in),value :: data_object,pdx,pr
real(c_double),intent(in),value :: d2edr2
type(c_ptr), intent(in),value :: pi,pj
integer(c_int),intent(out) :: ierr
! Local variables
integer(c_int), pointer :: i(:);
integer(c_int), pointer :: j(:);
real(c_double), pointer :: dx(:,:)
real(c_double), pointer :: r(:)
real(c_double) :: R1(DIM),R2(DIM),R3(DIM),R4(DIM)
call c_f_pointer(pi, i, [2])
call c_f_pointer(pj, j, [2])
call c_f_pointer(pdx,dx,[DIM,2])
call c_f_pointer(pr,r,[2])
R1 = dx(:,1)
R2 = dx(:,2)
R3 = (pA%coords(:,i(1)) + pA%coords(:,j(1)))/2._cd
R4 = (pA%coords(:,i(2)) + pA%coords(:,j(2)))/2._cd
pA%fsgeType%Eijkln = pA%fsgeType%Eijkln + d2edr2*tensor_prod(R1,R1,R2,R2,R4)/(r(1)*r(2))
pA%fsgeType%Cijkl = pA%fsgeType%Cijkl + d2edr2*tensor_prod(R1,R1,R2,R2)/(r(1)*r(2))
pA%fsgeType%Dijmkln= pA%fsgeType%Dijmkln + d2edr2*tensor_prod(R1,R1,R3,R2,R2,R4)/(r(1)*r(2))
ierr = 0
return
end subroutine process_d2edr2
function tensorG(dx,R3)
implicit none
real(c_double) ::tensorG(DIM,DIM,DIM,DIM)
real(c_double), intent(in) :: dx(DIM), R3(DIM)
integer :: p,i,j,m
tensorG = 0._cd
do m=1,DIM
do j=1,DIM
do i=1,DIM
do p=1,DIM
if (p .eq. i) tensorG(p,i,j,m) = tensorG(p,i,j,m) + dx(j)*R3(m) + dx(m)*R3(j)
if (p .eq. j) tensorG(p,i,j,m) = tensorG(p,i,j,m) + dx(i)*R3(m) + dx(m)*R3(i)
if (p .eq. m) tensorG(p,i,j,m) = tensorG(p,i,j,m) - dx(i)*R3(j) - dx(j)*R3(i)
enddo
enddo
enddo
enddo
tensorG = 0.5_cd * tensorG
end function tensorG
function tensorG_dot_tensorG(tensorG)
implicit none
real(c_double) ::tensorG_dot_tensorG(DIM,DIM,DIM,DIM,DIM,DIM)
real(c_double), intent(in) :: tensorG(DIM,DIM,DIM,DIM)
integer :: i,j,m,k,l,n,p
tensorG_dot_tensorG = 0._cd
do i = 1,DIM
do j = 1,DIM
do m = 1,DIM
do k = 1,DIM
do l = 1,DIM
do n = 1,DIM
do p = 1,DIM
tensorG_dot_tensorG(i,j,m,k,l,n) = tensorG_dot_tensorG(i,j,m,k,l,n) + &
tensorG(p,i,j,m) * tensorG(p,k,l,n)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
end function tensorG_dot_tensorG
subroutine constructNeighList(aS)
implicit none
type(Atomistic) :: aS
!call NEIGH_PURE_neighborlist(.false.,aS)
call NEIGH_PURE_neighborlist(aS)
end subroutine constructNeighList
subroutine NEIGH_PURE_neighborlist(aS)
use mod_global
implicit none
!-- Transferred variables
type(Atomistic) :: aS
!-- Local variables
integer(c_int) i, j, a
real(c_double) dx(DIM)
real(c_double) r2
real(c_double) cutoff2
cutoff2 = aS%cutoff**2
do i=1,aS%numParticles
a = 1
do j=1,aS%numParticles
dx(:) = aS%coords(:, j) - aS%coords(:, i)
r2 = dot_product(dx, dx)
if (r2.le.cutoff2) then
if (i .le. aS%numContributingParticles .and. i.ne.j) then
a = a+1
aS%neighObject%neighborList(a,i) = j
endif
endif
enddo
! part i has a-1 neighbors
aS%neighObject%neighborList(1,i) = a-1
enddo
return
end subroutine NEIGH_PURE_neighborlist
subroutine destroyAtomistic(aS)
implicit none
type(Atomistic) :: aS
if (allocated(aS%coords)) deallocate(aS%coords)
if (allocated(aS%particleSpecies)) deallocate(aS%particleSpecies)
if (allocated(aS%particle_contributing)) deallocate(aS%particle_contributing)
if (allocated(aS%neighObject%neighborList)) deallocate(aS%neighObject%neighborList)
if (allocated(aS%fsgeType%IndComponentsC)) deallocate(aS%fsgeType%IndComponentsC)
if (allocated(aS%fsgeType%IndComponentsD)) deallocate(aS%fsgeType%IndComponentsD)
if (allocated(aS%fsgeTypeNumerical%IndComponentsC)) deallocate(aS%fsgeTypeNumerical%IndComponentsC)
if (allocated(aS%fsgeTypeNumerical%IndComponentsD)) deallocate(aS%fsgeTypeNumerical%IndComponentsD)
end subroutine destroyAtomistic
end module mod_atomistic