v1.6.4

add a function to calculate the hardness
shahramyalameha · May 31, 2021 · 1f1581c · 1f1581c
1 parent a70619f
commit 1f1581c
Show file tree

Hide file tree

Showing 9 changed files with 569 additions and 334 deletions.
diff --git a/soc/Eatools_db.f90 b/soc/Eatools_db.f90
@@ -9,7 +9,7 @@ Subroutine databank(id,yesno)
 integer::i,j,stat
 !read(*,*)id 
      call system("clear")
- open (12,file='/home/shahram/Desktop/Cubelast/code/programMY/AAEP/soc/eatools_v1.5.2/db/All_2ID_cop.csv')
+ open (12,file='/home/shahram/Desktop/Cubelast/code/programMY/AAEP/soc/eatools_v1.6.2/db/All_2ID_cop.csv')
 do  i=1,13122
    read(12,* )id2(i)
    if (id2(i)==id )then
@@ -21,7 +21,7 @@ Subroutine databank(id,yesno)
  close (12)
 
  i=0
- open (11,file='/home/shahram/Desktop/Cubelast/code/programMY/AAEP/soc/eatools_v1.5.2/db/Cijs.binery')
+ open (11,file='/home/shahram/Desktop/Cubelast/code/programMY/AAEP/soc/eatools_v1.6.2/db/Cijs.binery')
 do  i=1,2043900
    read(11,'(Z16)',IOSTAT=stat )c1(i) ,c2(i),c3(i),c4(i),c5(i),c6(i)
  !   write(14,'(B64)')c1(j) ,c2(j),c3(j),c4(j),c5(j),c6(j)

diff --git a/soc/Eatools_hardness.f90 b/soc/Eatools_hardness.f90
@@ -0,0 +1,19 @@
+!```````````````````````````````````````````````````````````````````````````````````````````
+! Copyright (c) 2018 Shahram Yalameha <[email protected]> , <[email protected]>, `
+!               Please report bugs or suggestions to:  [email protected]                `
+!                                                                                          `
+!```````````````````````````````````````````````````````````````````````````````````````````
+! SUBROUTINE: fOR 3D MATERIAL, the core of the hardness calculation..
+! H = 0.786835118497608*G.^(0.7343).*k.^1.01993441066304; see ref:https://zenodo.org/record/3376758#.YLUzVppLjjE
+
+
+SUBROUTINE CHardness(B,E,hardvar)
+
+    DOUBLE PRECISION, PARAMETER                    ::  const1 = 0.130548175274347D0,&
+                                                       const2 = 2.2484942942017D0,  &
+                                                       const3 = -1.51675853808829D0
+    DOUBLE PRECISION                               ::  B, E, hardvar
+    hardvar = const1*(E**const2)*(B**const3)
+
+
+    end subroutine
diff --git a/soc/Eatools_main.f90 b/soc/Eatools_main.f90
@@ -4,7 +4,7 @@
 !               Please report bugs or suggestions to:  [email protected]                `
 !                                                                                          `
 !```````````````````````````````````````````````````````````````````````````````````````````
-!!  "EaTools: The code creates 3D and 2D models of an elastic property’s anisotropy. 
+!!  ElATools: A tool for analyzing anisotropic elastic properties of the 2D and 3D materials
 
 
  PROGRAM AAEP_main
@@ -29,6 +29,8 @@ PROGRAM AAEP_main
  bulkmax2d,    &
  young2dmin,   &
  young2dmax,   &
+ hard2dmax,    &
+ hard2dmin,    &
  poisson2dmax, &
  poisson2dminn,&
  poisson2dminp,&
@@ -99,6 +101,13 @@ PROGRAM AAEP_main
  G_min2_theta=0D0,          &
  G_min2_phi=0D0,            &
 
+Ha_max2 = 0D0,&
+Ha_min2 = 0.0D0,& 
+Ha_max2_theta = 0D0,& 
+Ha_min2_theta = 0D0, &
+Ha_max2_phi   = 0D0, &
+Ha_min2_phi   = 0D0,&
+
  Pugh_max=0D0,                 &
  Pugh_min=0D0,                 &
  Pugh_Ave=0D0,                 &
@@ -164,7 +173,7 @@ PROGRAM AAEP_main
  c1=0D0,                           &
  inter_theta=0D0,                  &
  phi=0,                            &
- BINver=0D0,                       &
+ BINver=0D0,hardvar=0D0,            &
  a1,a2,a3
 
  DOUBLE PRECISION :: Minyoung=10d8,&
@@ -649,7 +658,7 @@ PROGRAM AAEP_main
 !open(78,file='uuu')
   OPEN(111, file="3d_SD.dat")
   OPEN(14,FILE="3d_comp.dat"    )
- ! OPEN(28,FILE="CompP_Min.dat"    )
+  OPEN(28,FILE="3d_harnness.dat" )
  ! OPEN(29,FILE="CompN_Min.dat"    )
 
   OPEN(65,FILE="3d_pugh.dat"    )
@@ -933,7 +942,20 @@ PROGRAM AAEP_main
       IF (G_min.ge.0D0) shminp =  G_min
       IF (G_Ave.LE.0D0) shaven = -G_Ave
       IF (G_Ave.ge.0D0) shavep =  G_Ave
-
+!!!!!!!!!!!!!!!!!!!v1.6.3
+      CALL CHardness(BINver,SINver, hardvar)
+
+      IF (hardvar.GE.Ha_max2) THEN
+        Ha_max2       = hardvar   
+        Ha_max2_theta = theta
+        Ha_max2_phi   = phi
+      ENDIF  
+      IF (hardvar.LE.Ha_min2) THEN
+        Ha_min2       = hardvar
+        Ha_min2_theta = theta
+        Ha_min2_phi   = phi
+      ENDIF
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       CALL Cpugh(pugh_min,pugh_max,pugh_Ave,phi,theta,v11,v12,v13,v22,v23,v33,a6666,sheainvar)
 
       IF (pugh_max.GE.pugh_max2) THEN
@@ -1018,13 +1040,14 @@ PROGRAM AAEP_main
       !open(41,file='o')
       WRITE(59,*) G_max*1000D0,shminp*1000D0,shminn*1000D0,shavep*1000D0,SINver,CO*1000d0,comminp*1000d0,&
       comminn*1000d0,NPratio_max,pminp,pminn,pavep,paven,BINver,maxEVaLM,maxEVaTM,minEVaTM,pugh_max*1000D0,&
-      pughminp*1000D0,pughminn*1000D0,pughavep*1000D0 
+      pughminp*1000D0,pughminn*1000D0,pughavep*1000D0,hardvar
       If (yn_veloc=='Y' .or. yn_veloc=='y') then
         WRITE(55,*) VVP_P/1000D0,VVG_P/1000D0,VVP_Sf/1000D0,VVG_Sf/1000D0,VVP_Ss/1000D0,VVG_Ss/1000D0,VV_P_PF,VV_Sf_PF,VV_Ss_PF
       endif
 
 
       WRITE(111,*) theta, phi
+      WRITE(28,"(4F30.15)")  vec(1)*ABS( hardvar ),    vec(2)*ABS( hardvar ),    vec(3)*ABS( hardvar )                     !>> hardness
       !WRITE(18,*) vec(1)*ABS( G_Ave*1000D0 ),    vec(2)*ABS( G_Ave*1000D0 ),    vec(3)*ABS( G_Ave*1000D0 )   !>> shear
       WRITE(16,"(12F30.15)") vec(1)*ABS( G_max*1000D0 ),    vec(2)*ABS( G_max*1000D0 ),    vec(3)*ABS( G_max*1000D0 ),&    !>> shear
                              vec(1)*ABS(shminp*1000D0 ),    vec(2)*ABS( shminp*1000D0),    vec(3)*ABS(shminp*1000D0 ),&    !>> shear
@@ -1066,7 +1089,7 @@ PROGRAM AAEP_main
         WRITE(105,*) vec(1)*ABS(    VV_Sf_PF),    vec(2)*ABS(   VV_Sf_PF  ) ,vec(3)*ABS(   VV_Sf_PF )   !>>  VV_Sf_PF
         WRITE(106,*) vec(1)*ABS(    VVG_Ss  ),    vec(2)*ABS(   VVG_Ss    ) ,vec(3)*ABS(   VVG_Ss   )   !>>  VVG_Ss
         WRITE(107,*) vec(1)*ABS(    VVP_Ss  ),    vec(2)*ABS(   VVP_Ss    ) ,vec(3)*ABS(   VVP_Ss   )   !>>  VVP_Ss
-        WRITE(108,*) vec(1)*ABS(  VV_Ss_PF  ),    vec(2)*ABS(   VV_Ss_PF  ) ,vec(3)*ABS(  VV_Ss_PF  )      !>>  VVP_Ss_PF
+        WRITE(108,*) vec(1)*ABS(  VV_Ss_PF  ),    vec(2)*ABS(   VV_Ss_PF  ) ,vec(3)*ABS(  VV_Ss_PF  )   !>>  VVP_Ss_PF
       end if
     ENDDO 
 
@@ -1131,13 +1154,13 @@ PROGRAM AAEP_main
  ! If (yn_veloc=='Y' .or. yn_veloc=='y') then
     open(37,file='.MaMiout2')
     write(37,*) VVP_P_max/1000D0,VVP_P_min/1000D0,VVP_Sf_max/1000D0,VVP_Sf_min/1000D0,VVP_Ss_max/1000D0,VVP_Ss_min/1000D0,VVG_P_max/1000D0,&
-             VVG_P_min/1000D0,VVG_Sf_max/1000D0,VVG_Sf_min/1000D0,VVG_Ss_max/1000D0,VVG_Ss_min/1000D0,VV_P_PF_max,VV_Sf_PF_max,VV_Ss_PF_max,&
-             VV_P_PF_min,VV_Sf_PF_min,VV_Ss_PF_min 
+                VVG_P_min/1000D0,VVG_Sf_max/1000D0,VVG_Sf_min/1000D0,VVG_Ss_max/1000D0,VVG_Ss_min/1000D0,VV_P_PF_max,VV_Sf_PF_max,VV_Ss_PF_max,&
+                VV_P_PF_min,VV_Sf_PF_min,VV_Ss_PF_min 
     close(37)
  ! ENDIF
   open(36,file='.MaMiout')
   write(36,*)Maxyoung,Minyoung,Maxcomp,Mincomp,G_max2*1000D0,G_min2*1000D0,Maxbulk,Minbulk,&
-            Pratio_max,Pratio_min,maxEVaTMf,maxEVaLM,minEVaTMf,pugh_max2*1000D0,pugh_min2*1000D0
+             Pratio_max,Pratio_min,maxEVaTMf,maxEVaLM,minEVaTMf,pugh_max2*1000D0,pugh_min2*1000D0, Ha_max2,Ha_min2
   close(36)
 
 !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@  Preparing features...
@@ -1172,7 +1195,7 @@ PROGRAM AAEP_main
 	WRITE (*,'(1x3F6.2,3x3F6.2)') vec(1),vec(2),vec(3),vec1(1),vec1(2),vec1(3)
 	WRITE (*,*) "=================================================="
 	CALL SYSTEM('sleep 0.5')
-	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+  !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   WRITE(*,*)''
   call anisotropy(G_max2,G_min2, Ax(1))
 	CALL SYSTEM('tput setaf 161;tput bold; echo " ==================================================> Shear Modulus";tput sgr0')
@@ -1188,7 +1211,7 @@ PROGRAM AAEP_main
 	WRITE (*,'(1x3F6.2,3x3F6.2)') vec(1),vec(2),vec(3),vec1(1),vec1(2),vec1(3)
 	WRITE (*,*) "=================================================="
 	CALL SYSTEM('sleep 0.5')
-	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+  !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   WRITE(*,*)''
   call anisotropy(Maxbulk,Minbulk, Ax(1))
 	CALL SYSTEM('tput setaf 126;tput bold; echo " ==================================================> Bulk Modulus";tput sgr0')
@@ -1341,15 +1364,31 @@ PROGRAM AAEP_main
 	WRITE (*,*) "=================================================="
 	CALL SYSTEM('sleep 0.5')
 	WRITE(*,*)''
-
-	CALL SYSTEM('tput setaf 56;tput bold; echo " ==================================================> Sound";tput sgr0')
+  	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+  WRITE(*,*)''
+  call anisotropy(Maxbulk,Minbulk, Ax(1))
+	CALL SYSTEM('tput setaf 55;tput bold; echo " ==================================================> Hardness";tput sgr0')
+	CALL SYSTEM('tput setaf 55;tput bold; echo "         Max                     Min  Anisotropy"')
+	WRITE (*,'(2xF10.2,a,1xF10.2,a,2xF6.2)') Ha_max2,"            ",Ha_min2," ",Ax(1)
+	CALL angl2cart(Ha_max2_theta,Ha_max2_phi, vec(1),vec(2),vec(3))
+	CALL angl2cart(Ha_min2_theta,Ha_min2_phi,vec1(1),vec1(2),vec1(3))
+	WRITE(*,*)"------------------------------------------"
+	WRITE (*,*) "  Theta   Phi","          Theta   Phi"
+	WRITE (*,'(1x2F6.1,11x2F6.1)') (Ha_max2_theta*180.0D0)/PI,(Ha_max2_phi*180.0D0)/PI,(Ha_min2_theta*180.0D0)/PI,(Ha_min2_phi*180.0D0)/PI
+	WRITE(*,*)"------------------------------------------"
+	WRITE (*,*) "   x     y     z","        x     y     z"
+	WRITE (*,'(1x3F6.2,3x3F6.2)') vec(1),vec(2),vec(3),vec1(1),vec1(2),vec1(3)
+  WRITE (*,*) "=================================================="
+  CALL SYSTEM('sleep 0.5')
+  WRITE(*,*)''
+  !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+	CALL SYSTEM('tput setaf 56;tput bold; echo " ==================================================> Sound < testing";tput sgr0')
 	CALL SYSTEM('tput setaf 56;tput bold; echo "  Transverse high  Longitudinal   Transverse low "')
 	WRITE (*,'(3xF6.2,a,2xF6.2,a,2xF6.2)') maxEVaTMf,"         ",maxEVaLM,"         ",minEVaTMf
 
 	CALL SYSTEM('tput setaf 56;tput bold; echo " ==================================================";tput sgr0')
 	WRITE(*,*)' '
 
-
 !***************************************************************
   WRITE(*,*)" > NOTE :The above information is stored in  'DATA.out' file."
   WRITE (99,'(a,F4.1,a,F4.1,a,F4.1,a)') " > Output for  (",mmx,",",kky,",",llz,") plane"
@@ -1618,6 +1657,14 @@ PROGRAM AAEP_main
         Minbulk       = BINver                       
         bulkMin2d(num) = Maxbulk 
       ENDIF
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      CALL CHardness(BINver,SINver,hardvar)  
+
+      IF (hardvar.NE.0.0) THEN
+        hard2dmax(num) = hardvar   
+      ELSE  
+        hard2dmax(num) = 0D0
+      ENDIF        
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       ave=0D0
       v=0
@@ -1740,7 +1787,7 @@ PROGRAM AAEP_main
     !OPEN(  63,FILE="2dsoundTMmax.dat" )
    ! OPEN(  64,FILE="2dsoundTMmin.dat" )
     OPEN(  71,FILE="2dcut_pugh.dat"   )
-   ! OPEN(  72,FILE="2dpughMinp.dat"  )
+    OPEN(  72,FILE="2dcut_hardness.dat")
    ! OPEN(  73,FILE="2dpughMinn.dat"  )
     If (yn_veloc=='Y' .or. yn_veloc=='y')   then !-----------------------------------------------------
 
@@ -1778,6 +1825,9 @@ PROGRAM AAEP_main
                                   SQRT( (pugh2dminp(num)*COS(t2d)  )**2D0 + (pugh2dminp(num)*SIN(t2d)    )**2D0),&
                                   SQRT( (pugh2dminn(num)*COS(t2d)  )**2D0 + (pugh2dminn(num)*SIN(t2d)    )**2D0) 
 
+      WRITE(72,"(I3,4F30.20)")num,SQRT( (hard2dmax(num)*COS(t2d)   )**2D0 + (hard2dmax(num)*SIN(t2d)   )**2D0),& 
+                                  SQRT( (hard2dmin(num)*COS(t2d)   )**2D0 + (hard2dmin(num)*SIN(t2d)   )**2D0)
+
       IF (yn_veloc=='Y' .or. yn_veloc=='y') then
         WRITE(500,'(I3,3F30.15)')num,SQRT( (VVG_P_2D(num)*COS(t2d)   )**2D0 + (VVG_P_2D(num)*SIN(t2d)      )**2D0),&
                                      SQRT( (VVG_Sf_2D(num)*COS(t2d)   )**2D0 + (VVG_Sf_2D(num)*SIN(t2d)    )**2D0),&
@@ -1808,7 +1858,7 @@ PROGRAM AAEP_main
     !CLOSE(63)
    ! CLOSE(64)
     CLOSE(71)
-    !CLOSE(72)
+    CLOSE(72)
     !CLOSE(73)
     CLOSE(500)
     CLOSE(501)