git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12363 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
@ -1,184 +1,200 @@
|
||||
! This program calculates elastic constants
|
||||
! Must be used with in.comb-elastic script, and reads output log.lammps
|
||||
! Written by T-R Shan, MSE, UF, Feb 2010
|
||||
! This program calculates elastic constants
|
||||
! Must be used with in.comb-elastic script,
|
||||
! and reads a LAMMPS log file from standard input:
|
||||
!
|
||||
! f95 -O -o elastic.x elastic.f90
|
||||
! ./elastic.x < log.lammps
|
||||
!
|
||||
! Written by T-R Shan, MSE, UF, Feb 2010
|
||||
|
||||
program main
|
||||
implicit none
|
||||
|
||||
integer :: i,j,k,l
|
||||
real*8 :: box(6,11),force(6,11),vol(11),eps(11),strbox(6,11)
|
||||
real*8 :: sumx,sumx2,sumy,sumxy,c11,c12,c13,c14,c15,c16,c33,c44,c66,bulk,shear
|
||||
character*7 :: header1
|
||||
PROGRAM main
|
||||
IMPLICIT NONE
|
||||
|
||||
open(5,status='old',file='log.lammps')
|
||||
INTEGER, PARAMETER :: dbl=8
|
||||
INTEGER :: i,j,k,l
|
||||
REAL(kind=dbl) :: box(6,11),force(6,11),vol(11),eps(11),strbox(6,11)
|
||||
REAL(kind=dbl) :: sumx,sumx2,sumy,sumxy
|
||||
REAL(kind=dbl) :: c11,c12,c13,c14,c15,c16,c33,c44,c66
|
||||
REAL(kind=dbl) :: bulk,shear
|
||||
CHARACTER*7 :: header1
|
||||
|
||||
do i=1,500
|
||||
read(5,'(a7)') header1
|
||||
if(header1.eq.'Step Lx')then
|
||||
do j=1,11
|
||||
read(5,*) k,(box(l,j),l=1,6),(force(l,j),l=1,6)
|
||||
vol(j)=box(1,j)*box(2,j)*box(3,j)*dsqrt(1-(box(4,j)/box(1,j))*(box(4,j)/box(1,j))-(box(6,j)/box(3,j))*(box(6,j)/box(3,j)-(box(5,j)/box(3,j))*(box(5,j)/box(3,j))))
|
||||
eps(j)=(box(1,j)-box(1,1))/box(1,1)
|
||||
do l=1,6
|
||||
strbox(l,j)=force(l,j)/vol(j)/1.0d4
|
||||
enddo
|
||||
enddo
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
! open(5,status='old',file='log.lammps')
|
||||
|
||||
!! C11
|
||||
sumx = 0.0; sumx2 = 0.0; sumy = 0.0; sumxy = 0.0
|
||||
do j=2,11
|
||||
sumx = sumx + eps(j)
|
||||
sumx2= sumx2+ eps(j)*eps(j)
|
||||
sumy = sumy + strbox(1,j)
|
||||
sumxy= sumxy+ strbox(1,j)*eps(j)
|
||||
enddo
|
||||
c11 = (sumxy-sumx*sumy/10.0)/(sumx2-sumx*sumx/10.0)
|
||||
!! C12
|
||||
sumy = 0.0; sumxy = 0.0
|
||||
do j=2,11
|
||||
sumy = sumy + strbox(2,j)
|
||||
sumxy= sumxy+ strbox(2,j)*eps(j)
|
||||
enddo
|
||||
c12 = (sumxy-sumx*sumy/10.0)/(sumx2-sumx*sumx/10.0)
|
||||
!! C13
|
||||
sumy = 0.0; sumxy = 0.0
|
||||
do j=2,11
|
||||
sumy = sumy + strbox(3,j)
|
||||
sumxy= sumxy+ strbox(3,j)*eps(j)
|
||||
enddo
|
||||
c13 = (sumxy-sumx*sumy/10.0)/(sumx2-sumx*sumx/10.0)
|
||||
!! C14
|
||||
sumy = 0.0; sumxy = 0.0
|
||||
do j=2,11
|
||||
sumy = sumy + strbox(4,j)
|
||||
sumxy= sumxy+ strbox(4,j)*eps(j)
|
||||
enddo
|
||||
c14 = (sumxy-sumx*sumy/10.0)/(sumx2-sumx*sumx/10.0)
|
||||
!! C15
|
||||
sumy = 0.0; sumxy = 0.0
|
||||
do j=2,11
|
||||
sumy = sumy + strbox(5,j)
|
||||
sumxy= sumxy+ strbox(5,j)*eps(j)
|
||||
enddo
|
||||
c15 = (sumxy-sumx*sumy/10.0)/(sumx2-sumx*sumx/10.0)
|
||||
!! C16
|
||||
sumy = 0.0; sumxy = 0.0
|
||||
do j=2,11
|
||||
sumy = sumy + strbox(6,j)
|
||||
sumxy= sumxy+ strbox(6,j)*eps(j)
|
||||
enddo
|
||||
c16 = (sumxy-sumx*sumy/10.0)/(sumx2-sumx*sumx/10.0)
|
||||
|
||||
!! C33
|
||||
do i=1,500
|
||||
read(5,'(a7)') header1
|
||||
if(header1.eq.'Step Lx')then
|
||||
do j=1,11
|
||||
read(5,*) k,(box(l,j),l=1,6),(force(l,j),l=1,6)
|
||||
vol(j)=box(1,j)*box(2,j)*box(3,j)*dsqrt(1-(box(4,j)/box(1,j))*(box(4,j)/box(1,j))-(box(6,j)/box(3,j))*(box(6,j)/box(3,j)-(box(5,j)/box(3,j))*(box(5,j)/box(3,j))))
|
||||
eps(j)=(box(3,j)-box(3,1))/box(3,1)
|
||||
do l=1,6
|
||||
strbox(l,j)=force(l,j)/vol(j)/1.0d4
|
||||
enddo
|
||||
enddo
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
DO i=1,500
|
||||
READ(5,'(a7)') header1
|
||||
IF(header1.EQ.'Step Lx')THEN
|
||||
DO j=1,11
|
||||
READ(5,*) k,(box(l,j),l=1,6),(force(l,j),l=1,6)
|
||||
vol(j)=box(1,j)*box(2,j)*box(3,j) * dsqrt(1-(box(4,j)/box(1,j)) &
|
||||
*(box(4,j)/box(1,j))-(box(6,j)/box(3,j))*(box(6,j)/box(3,j) &
|
||||
-(box(5,j)/box(3,j))*(box(5,j)/box(3,j))))
|
||||
eps(j)=(box(1,j)-box(1,1))/box(1,1)
|
||||
DO l=1,6
|
||||
strbox(l,j)=force(l,j)/vol(j)/1.0d4
|
||||
ENDDO
|
||||
ENDDO
|
||||
EXIT
|
||||
ENDIF
|
||||
ENDDO
|
||||
|
||||
sumx = 0.0; sumx2 = 0.0; sumy = 0.0; sumxy = 0.0
|
||||
do j=2,11
|
||||
sumx = sumx + eps(j)
|
||||
sumx2= sumx2+ eps(j)*eps(j)
|
||||
sumy = sumy + strbox(3,j)
|
||||
sumxy= sumxy+ strbox(3,j)*eps(j)
|
||||
enddo
|
||||
c33 = (sumxy-sumx*sumy/10.0)/(sumx2-sumx*sumx/10.0)
|
||||
!! C11
|
||||
sumx = 0.0_dbl; sumx2 = 0.0_dbl; sumy = 0.0_dbl; sumxy = 0.0_dbl
|
||||
DO j=2,11
|
||||
sumx = sumx + eps(j)
|
||||
sumx2= sumx2+ eps(j)*eps(j)
|
||||
sumy = sumy + strbox(1,j)
|
||||
sumxy= sumxy+ strbox(1,j)*eps(j)
|
||||
ENDDO
|
||||
c11 = (sumxy-sumx*sumy/10.0_dbl)/(sumx2-sumx*sumx/10.0_dbl)
|
||||
!! C12
|
||||
sumy = 0.0_dbl; sumxy = 0.0_dbl
|
||||
DO j=2,11
|
||||
sumy = sumy + strbox(2,j)
|
||||
sumxy= sumxy+ strbox(2,j)*eps(j)
|
||||
ENDDO
|
||||
c12 = (sumxy-sumx*sumy/10.0_dbl)/(sumx2-sumx*sumx/10.0_dbl)
|
||||
!! C13
|
||||
sumy = 0.0_dbl; sumxy = 0.0_dbl
|
||||
DO j=2,11
|
||||
sumy = sumy + strbox(3,j)
|
||||
sumxy= sumxy+ strbox(3,j)*eps(j)
|
||||
ENDDO
|
||||
c13 = (sumxy-sumx*sumy/10.0_dbl)/(sumx2-sumx*sumx/10.0_dbl)
|
||||
!! C14
|
||||
sumy = 0.0_dbl; sumxy = 0.0_dbl
|
||||
DO j=2,11
|
||||
sumy = sumy + strbox(4,j)
|
||||
sumxy= sumxy+ strbox(4,j)*eps(j)
|
||||
ENDDO
|
||||
c14 = (sumxy-sumx*sumy/10.0_dbl)/(sumx2-sumx*sumx/10.0_dbl)
|
||||
!! C15
|
||||
sumy = 0.0_dbl; sumxy = 0.0_dbl
|
||||
DO j=2,11
|
||||
sumy = sumy + strbox(5,j)
|
||||
sumxy= sumxy+ strbox(5,j)*eps(j)
|
||||
ENDDO
|
||||
c15 = (sumxy-sumx*sumy/10.0_dbl)/(sumx2-sumx*sumx/10.0_dbl)
|
||||
!! C16
|
||||
sumy = 0.0_dbl; sumxy = 0.0_dbl
|
||||
DO j=2,11
|
||||
sumy = sumy + strbox(6,j)
|
||||
sumxy= sumxy+ strbox(6,j)*eps(j)
|
||||
ENDDO
|
||||
c16 = (sumxy-sumx*sumy/10.0_dbl)/(sumx2-sumx*sumx/10.0_dbl)
|
||||
|
||||
!! C44
|
||||
do i=1,500
|
||||
read(5,'(a7)') header1
|
||||
if(header1.eq.'Step Lx')then
|
||||
do j=1,11
|
||||
read(5,*) k,(box(l,j),l=1,6),(force(l,j),l=1,6)
|
||||
vol(j)=box(1,j)*box(2,j)*box(3,j)*dsqrt(1.0-(box(4,j)/box(1,j))*(box(4,j)/box(1,j))-(box(6,j)/box(3,j))*(box(6,j)/box(3,j)-(box(5,j)/box(3,j))*(box(5,j)/box(3,j))))
|
||||
do l=1,6
|
||||
strbox(l,j)=force(l,j)/vol(j)/1.0d4
|
||||
enddo
|
||||
enddo
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
|
||||
sumx = 0.0; sumx2 = 0.0; sumy = 0.0; sumxy = 0.0
|
||||
do j=2,11
|
||||
eps(j)=asin(box(6,j)/box(3,j))
|
||||
sumx = sumx + eps(j)
|
||||
sumx2= sumx2+ eps(j)*eps(j)
|
||||
sumy = sumy + strbox(6,j)
|
||||
sumxy= sumxy+ strbox(6,j)*eps(j)
|
||||
enddo
|
||||
c44 = (sumxy-sumx*sumy/10.0)/(sumx2-sumx*sumx/10.0)
|
||||
!! C33
|
||||
DO i=1,500
|
||||
READ(5,'(a7)') header1
|
||||
IF(header1.EQ.'Step Lx')THEN
|
||||
DO j=1,11
|
||||
READ(5,*) k,(box(l,j),l=1,6),(force(l,j),l=1,6)
|
||||
vol(j)=box(1,j)*box(2,j)*box(3,j) * dsqrt(1-(box(4,j)/box(1,j)) &
|
||||
*(box(4,j)/box(1,j))-(box(6,j)/box(3,j))*(box(6,j)/box(3,j) &
|
||||
-(box(5,j)/box(3,j))*(box(5,j)/box(3,j))))
|
||||
eps(j)=(box(3,j)-box(3,1))/box(3,1)
|
||||
DO l=1,6
|
||||
strbox(l,j)=force(l,j)/vol(j)/1.0d4
|
||||
ENDDO
|
||||
ENDDO
|
||||
EXIT
|
||||
ENDIF
|
||||
ENDDO
|
||||
|
||||
!! C66
|
||||
do i=1,500
|
||||
read(5,'(a7)') header1
|
||||
if(header1.eq.'Step Lx')then
|
||||
do j=1,11
|
||||
read(5,*) k,(box(l,j),l=1,6),(force(l,j),l=1,6)
|
||||
vol(j)=box(1,j)*box(2,j)*box(3,j)*dsqrt(1.0-(box(4,j)/box(1,j))*(box(4,j)/box(1,j))-(box(6,j)/box(3,j))*(box(6,j)/box(3,j)-(box(5,j)/box(3,j))*(box(5,j)/box(3,j))))
|
||||
do l=1,6
|
||||
strbox(l,j)=force(l,j)/vol(j)/1.0d4
|
||||
enddo
|
||||
enddo
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
sumx = 0.0_dbl; sumx2 = 0.0_dbl; sumy = 0.0_dbl; sumxy = 0.0_dbl
|
||||
DO j=2,11
|
||||
sumx = sumx + eps(j)
|
||||
sumx2= sumx2+ eps(j)*eps(j)
|
||||
sumy = sumy + strbox(3,j)
|
||||
sumxy= sumxy+ strbox(3,j)*eps(j)
|
||||
ENDDO
|
||||
c33 = (sumxy-sumx*sumy/10.0_dbl)/(sumx2-sumx*sumx/10.0_dbl)
|
||||
|
||||
sumx = 0.0; sumx2 = 0.0; sumy = 0.0; sumxy = 0.0
|
||||
do j=2,11
|
||||
eps(j)=asin(box(4,j)/box(1,j))
|
||||
sumx = sumx + eps(j)
|
||||
sumx2= sumx2+ eps(j)*eps(j)
|
||||
sumy = sumy + strbox(4,j)
|
||||
sumxy= sumxy+ strbox(4,j)*eps(j)
|
||||
enddo
|
||||
c66 = (sumxy-sumx*sumy/10.0)/(sumx2-sumx*sumx/10.0)
|
||||
!! C44
|
||||
DO i=1,500
|
||||
READ(5,'(a7)') header1
|
||||
IF(header1.EQ.'Step Lx')THEN
|
||||
DO j=1,11
|
||||
READ(5,*) k,(box(l,j),l=1,6),(force(l,j),l=1,6)
|
||||
vol(j)=box(1,j)*box(2,j)*box(3,j)*dsqrt(1.0-(box(4,j)/box(1,j)) &
|
||||
*(box(4,j)/box(1,j))-(box(6,j)/box(3,j))*(box(6,j)/box(3,j) &
|
||||
-(box(5,j)/box(3,j))*(box(5,j)/box(3,j))))
|
||||
DO l=1,6
|
||||
strbox(l,j)=force(l,j)/vol(j)/1.0d4
|
||||
ENDDO
|
||||
ENDDO
|
||||
EXIT
|
||||
ENDIF
|
||||
ENDDO
|
||||
|
||||
bulk=(c11*2.0+c33)/9.0+2.0*(c12+c13*2.0)/9.0
|
||||
shear=(2.*c11+c33-c12-2.*c13)/15.+(2.*c44+c66)/5.
|
||||
sumx = 0.0_dbl; sumx2 = 0.0_dbl; sumy = 0.0_dbl; sumxy = 0.0_dbl
|
||||
DO j=2,11
|
||||
eps(j)=ASIN(box(6,j)/box(3,j))
|
||||
sumx = sumx + eps(j)
|
||||
sumx2= sumx2+ eps(j)*eps(j)
|
||||
sumy = sumy + strbox(6,j)
|
||||
sumxy= sumxy+ strbox(6,j)*eps(j)
|
||||
enddo
|
||||
c44 = (sumxy-sumx*sumy/10.0_dbl)/(sumx2-sumx*sumx/10.0_dbl)
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'Elastic constants (GPa):'
|
||||
write(*,101) 'C11 = ',c11
|
||||
write(*,101) 'C12 = ',c12
|
||||
write(*,101) 'C13 = ',c13
|
||||
write(*,101) 'C14 = ',c14
|
||||
write(*,101) 'C15 = ',c15
|
||||
write(*,101) 'C16 = ',c16
|
||||
write(*,101) 'C33 = ',c33
|
||||
write(*,101) 'C44 = ',c44
|
||||
write(*,101) 'C66 = ',c66
|
||||
write(*,101) 'B = ',bulk
|
||||
write(*,101) 'G = ',shear
|
||||
!! C66
|
||||
DO i=1,500
|
||||
READ(5,'(a7)') header1
|
||||
IF(header1.EQ.'Step Lx')THEN
|
||||
DO j=1,11
|
||||
READ(5,*) k,(box(l,j),l=1,6),(force(l,j),l=1,6)
|
||||
vol(j)=box(1,j)*box(2,j)*box(3,j)*dsqrt(1.0-(box(4,j)/box(1,j)) &
|
||||
*(box(4,j)/box(1,j))-(box(6,j)/box(3,j))*(box(6,j)/box(3,j) &
|
||||
-(box(5,j)/box(3,j))*(box(5,j)/box(3,j))))
|
||||
DO l=1,6
|
||||
strbox(l,j)=force(l,j)/vol(j)/1.0d4
|
||||
ENDDO
|
||||
ENDDO
|
||||
EXIT
|
||||
ENDIF
|
||||
ENDDO
|
||||
|
||||
write(5,*)
|
||||
write(5,*)'Elastic constants (GPa):'
|
||||
write(5,101) 'C11 = ',c11
|
||||
write(5,101) 'C12 = ',c12
|
||||
write(5,101) 'C13 = ',c13
|
||||
write(5,101) 'C14 = ',c14
|
||||
write(5,101) 'C15 = ',c15
|
||||
write(5,101) 'C16 = ',c16
|
||||
write(5,101) 'C33 = ',c33
|
||||
write(5,101) 'C44 = ',c44
|
||||
write(5,101) 'C66 = ',c66
|
||||
write(5,101) 'B = ',bulk
|
||||
write(5,101) 'G = ',shear
|
||||
sumx = 0.0_dbl; sumx2 = 0.0_dbl; sumy = 0.0_dbl; sumxy = 0.0_dbl
|
||||
DO j=2,11
|
||||
eps(j)=ASIN(box(4,j)/box(1,j))
|
||||
sumx = sumx + eps(j)
|
||||
sumx2= sumx2+ eps(j)*eps(j)
|
||||
sumy = sumy + strbox(4,j)
|
||||
sumxy= sumxy+ strbox(4,j)*eps(j)
|
||||
ENDDO
|
||||
c66 = (sumxy-sumx*sumy/10.0_dbl)/(sumx2-sumx*sumx/10.0_dbl)
|
||||
|
||||
101 format(a6,f11.3)
|
||||
bulk=(c11*2.0_dbl+c33)/9.0_dbl+2.0_dbl*(c12+c13*2.0_dbl)/9.0_dbl
|
||||
shear=(2.0_dbl*c11+c33-c12-2.0_dbl*c13)/15.0_dbl+(2.0_dbl*c44+c66)/5.0_dbl
|
||||
|
||||
stop
|
||||
end
|
||||
WRITE(*,*)
|
||||
WRITE(*,*)'Elastic constants (GPa):'
|
||||
WRITE(*,101) 'C11 = ',c11
|
||||
WRITE(*,101) 'C12 = ',c12
|
||||
WRITE(*,101) 'C13 = ',c13
|
||||
WRITE(*,101) 'C14 = ',c14
|
||||
WRITE(*,101) 'C15 = ',c15
|
||||
WRITE(*,101) 'C16 = ',c16
|
||||
WRITE(*,101) 'C33 = ',c33
|
||||
WRITE(*,101) 'C44 = ',c44
|
||||
WRITE(*,101) 'C66 = ',c66
|
||||
WRITE(*,101) 'B = ',bulk
|
||||
WRITE(*,101) 'G = ',shear
|
||||
|
||||
! WRITE(5,*)
|
||||
! WRITE(5,*)'Elastic constants (GPa):'
|
||||
! WRITE(5,101) 'C11 = ',c11
|
||||
! WRITE(5,101) 'C12 = ',c12
|
||||
! WRITE(5,101) 'C13 = ',c13
|
||||
! WRITE(5,101) 'C14 = ',c14
|
||||
! WRITE(5,101) 'C15 = ',c15
|
||||
! WRITE(5,101) 'C16 = ',c16
|
||||
! WRITE(5,101) 'C33 = ',c33
|
||||
! WRITE(5,101) 'C44 = ',c44
|
||||
! WRITE(5,101) 'C66 = ',c66
|
||||
! WRITE(5,101) 'B = ',bulk
|
||||
! WRITE(5,101) 'G = ',shear
|
||||
|
||||
101 FORMAT(a6,f11.3)
|
||||
|
||||
STOP
|
||||
END PROGRAM main
|
||||
|
||||
Reference in New Issue
Block a user