diff --git a/soc/Eatools_bulk.f90 b/soc/Eatools_bulk.f90 new file mode 100644 index 0000000..d803185 --- /dev/null +++ b/soc/Eatools_bulk.f90 @@ -0,0 +1,22 @@ +!``````````````````````````````````````````````````````````````````````````````````````````````` +! Copyright (c) 2018-2022 Shahram Yalameha , , ` +! Please report bugs or suggestions to: yalameha93@gmail.com ` +! ` +!``````````````````````````````````````````````````````````````````````````````````````````````` +! SUBROUTINE: fOR 3D MATERIAL, the core of the bulk Modulus method II calculation.. + +SUBROUTINE bulk_method(Pratio, sheainvar, bulk_m2) + +DOUBLE PRECISION :: Pratio, sheainvar, bulk_m2 + + + + +bulk_m2 = ( 2.d0 * sheainvar * (1.d0 + Pratio) ) / (3.0d0 * (1.0d0 - 2.0d0 * Pratio)) + + + + + + +END SUBROUTINE bulk_method diff --git a/soc/Eatools_main.f90 b/soc/Eatools_main.f90 index 7188f52..5462471 100644 --- a/soc/Eatools_main.f90 +++ b/soc/Eatools_main.f90 @@ -10,14 +10,14 @@ PROGRAM AAEP_main IMPLICIT NONE ChARACTER(LEN=10) :: val,& !> for property - e1,e2 + e1,e2 INTEGER, PARAMETER :: dp = selected_real_kind(15, 307) CHARACTER(LEN=1) :: YN="Y" ,& !> for plot 3d and 2d data is "Y" yesno,yesno2,& yn_veloc,yn_ewp,& order,yn_km CHARACTER(LEN=25) :: myid !> for databank - INTEGER :: CLaS,phi_meah,theta_meah,cutmesh,npoint,ewp + INTEGER :: CLaS,phi_meah,theta_meah,cutmesh,npoint,ewp DOUBLE PRECISION, PARAMETER :: pi=3.14159265358979323846264338327950D0,ee=0.0001D0 DOUBLE PRECISION, DIMENSION(6,6) :: C1p=0D0,S1=0D0,C=0D0,S=0D0,CP=0D0,C3=0D0,CCo=0d0,Eig3d DOUBLE PRECISION, DIMENSION(80100) :: shear2dmax, Ax,& @@ -50,259 +50,267 @@ PROGRAM AAEP_main VVG_Ss_2D, & VVP_Ss_2D, & VV_Ss_PF_2D - ChARACTER(LEN=3) :: adv - DOUBLE PRECISION, DIMENSION(3,3) :: EVe=0D0, C2D,Eig2d - - DOUBLE PRECISION, DIMENSION(3) :: vec=0d0,vec_3dslic=0,& - v=0d0, & - vec1=0d0, & - EVa, & - CElastcode=0d0 - INTEGER :: loop, d2d3 ,Stable,adv_mubner !0> Stable, 1> unStable - INTEGER ::i=0, & - num=1, & - lm,MinTm, & - MaxTm, & - Ncod, & - i2, & - j=0, & - jj=0, & - ii, & - Nmesh_thata=0d0, & - Nmesh_thataf=0d0,& - Nmesh_phi=0d0, & - Nmesh_phIF=0D0, & - Nmesh_phi2=0D0, & - Nmesh_phIF2=0D0, & - kk=0D0, & - Nmesh_gamma=0D0 + ChARACTER(LEN=3) :: adv + DOUBLE PRECISION, DIMENSION(3,3) :: EVe=0D0, C2D,Eig2d + + DOUBLE PRECISION, DIMENSION(3) :: vec=0d0, & + vec_3dslic=0, & + v=0d0, & + vec1=0d0, & + EVa, & + CElastcode=0d0 + INTEGER :: loop, d2d3 ,Stable,adv_mubner !0> Stable, 1> unStable + INTEGER :: i=0, & + num=1, & + lm,MinTm, & + MaxTm, & + Ncod, & + i2, & + j=0, & + jj=0, & + ii, & + Nmesh_thata=0d0, & + Nmesh_thataf=0d0,& + Nmesh_phi=0d0, & + Nmesh_phIF=0D0, & + Nmesh_phi2=0D0, & + Nmesh_phIF2=0D0, & + kk=0D0, & + Nmesh_gamma=0D0 - DOUBLE PRECISION :: l,an, & - maxEVaLM = 0D0, & - maxEVaTM = 0D0, & - minEVaTM = 10d6, & - maxEVaLMf = 0D0, & - maxEVaTMf = 0D0, & - minEVaTMf = 10d6,test_n, new_num + DOUBLE PRECISION :: l,an, & + maxEVaLM = 0D0, & + maxEVaTM = 0D0, & + minEVaTM = 10d6, & + maxEVaLMf = 0D0, & + maxEVaTMf = 0D0, & + minEVaTMf = 10d6,test_n, new_num - DOUBLE PRECISION ::v11=0d0,& - t2d=0D0, & - v12=0d0, & - v13=0d0, & - v22=0d0, & - v23=0d0, & - v33=0d0, & - G_max=0D0, & - G_min=0D0, & - G_Ave=0D0, & - G_max2=0D0, & - G_min2=1D0, & - G_max2_theta=0D0, & - G_max2_phi=0D0, & - 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, & - Pugh_max2=0D0, & - Pugh_min2=1D0, & - Pugh_max2_theta=0D0, & - Pugh_max2_phi=0D0, & - Pugh_min2_theta=0D0, & - Pugh_min2_phi=0D0, & - Maxcomp_theta=0D0, & - Maxcomp_phi=0D0, & - Mincomp_theta=0D0, & - Mincomp_phi=0D0, & - Mincomp=10d8, & - Maxcomp=0D0, & - NPratio_max=0D0, & - Pratio=0D0, & - NPratio_min=1D0, & - NPratio_ave=0D0, & - Pratio_max=0D0, & - Pratio_min=1D0, & - Pratio_max_theta=0D0, & - Pratio_max_phi=0D0, & - Pratio_min_theta=0D0, & - Pratio_min_phi=0D0, & - v1min=0D0, & - v2min=0D0, & - v3min=0D0, & - v1max=0D0, & - v2max=0D0, & - v3max=0D0, & - v1minmin=0D0, & - v2minmin=0D0, & - v3minmin=0D0, & - v1maxmax=0D0, & - v2maxmax=0D0, & - v3maxmax=0D0, & - pmax=0.0001D0, & - pminp=0.0001D0, & - pminn=0.0001D0, & - pavep=0D0, & - paven=0.0D0, & - shmax=0.01D0, & - shminp=0.01D0, & - shminn=0.01D0, & - shavep=0D0, & - shaven=0D0, & - pughmax=0.01D0, & - pughminp=0.01D0, & - pughminn=0.01D0, & - pughavep=0D0, & - pughaven=0D0, & - comminn=0.1D0, & - comminp=0.1D0 + DOUBLE PRECISION :: v11=0d0,& + t2d=0D0, & + v12=0d0, & + v13=0d0, & + v22=0d0, & + v23=0d0, & + v33=0d0, & + G_max=0D0, & + G_min=0D0, & + G_Ave=0D0, & + G_max2=0D0, & + G_min2=1D0, & + G_max2_theta=0D0, & + G_max2_phi=0D0, & + 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, & + Pugh_max2=0D0, & + Pugh_min2=1D0, & + Pugh_max2_theta=0D0, & + Pugh_max2_phi=0D0, & + Pugh_min2_theta=0D0, & + Pugh_min2_phi=0D0, & + Maxcomp_theta=0D0, & + Maxcomp_phi=0D0, & + Mincomp_theta=0D0, & + Mincomp_phi=0D0, & + Mincomp=10d8, & + Maxcomp=0D0, & + NPratio_max=0D0, & + Pratio=0D0, & + NPratio_min=1D0, & + NPratio_ave=0D0, & + Pratio_max=0D0, & + Pratio_min=1D0, & + Pratio_max_theta=0D0, & + Pratio_max_phi=0D0, & + Pratio_min_theta=0D0, & + Pratio_min_phi=0D0, & + v1min=0D0, & + v2min=0D0, & + v3min=0D0, & + v1max=0D0, & + v2max=0D0, & + v3max=0D0, & + v1minmin=0D0, & + v2minmin=0D0, & + v3minmin=0D0, & + v1maxmax=0D0, & + v2maxmax=0D0, & + v3maxmax=0D0, & + pmax=0.0001D0, & + pminp=0.0001D0, & + pminn=0.0001D0, & + pavep=0D0, & + paven=0.0D0, & + shmax=0.01D0, & + shminp=0.01D0, & + shminn=0.01D0, & + shavep=0D0, & + shaven=0D0, & + pughmax=0.01D0, & + pughminp=0.01D0, & + pughminn=0.01D0, & + pughavep=0D0, & + pughaven=0D0, & + comminn=0.1D0, & + comminp=0.1D0 REAL(dp) ::SS=1.0_dp,SINver=0.0_dp - DOUBLE PRECISION :: inter_phi=0D0,& - BB=1, & - Maxyoung_phi=0D0, & - Maxbulk_theta=0D0, & - Maxyoung_theta=0D0, & - Maxbulk_phi=0D0, & - b_r=0D0, & - c1=0D0, & - inter_theta=0D0, & - phi=0, & - BINver=0D0,hardvar=0D0, & - a1,a2,a3 + DOUBLE PRECISION :: inter_phi=0D0, & + BB=1, & + Maxyoung_phi=0D0, & + Maxbulk_theta=0D0, & + Maxyoung_theta=0D0, & + Maxbulk_phi=0D0, & + b_r=0D0, & + c1=0D0, & + inter_theta=0D0, & + phi=0, & + BINver=0D0,hardvar=0D0, & + a1,a2,a3 - DOUBLE PRECISION :: Minyoung=10d8,& - Minbulk=0.0d0, & - Maxyoung=0D0, & - Maxbulk=0D0, & - Minyoung_theta=0D0, & - Minbulk_theta=0D0, & - Minyoung_phi=0D0, & - Minbulk_phi=0D0, & - theta=0d0, & - smkl=0D0,smkl2=0D0, & - twoDTheta=0D0, & - planetheta, & - planephi, & - mmx=1d0, & - kky=0d0, & - llz=0d0, & - mmx1=0d0, & - kky1=0d0, & - llz1=0d0 - - - DOUBLE PRECISION :: a6666=0D0,& - sheainvar,CO=0D0, & - Minsheainvar, & - Maxsheainvar, & - sheainvar_ave, & - Minpughvar, & - Maxpughvar, & - pughvar_ave, & - gamma, & - vv11=0D0, & - vv12=0D0, & - vv13=0D0, & - vv22=0D0, & - vv23=0D0, & - vv33=0D0, & - ave, & - k11, & - k12, & - k13, & - k22, & - k23, & - k33 + DOUBLE PRECISION :: Minyoung=10d8, & + Minbulk=0.0d0, & + Maxyoung=0D0, & + Maxbulk=0D0, & + Minyoung_theta=0D0, & + Minbulk_theta=0D0, & + Maxbulkm2_theta=0D0, & + Minbulkm2_theta=0D0, & + Maxbulkm2_phi=0D0, & + Minbulkm2_phi=0D0, & + Minyoung_phi=0D0, & + Minbulk_phi=0D0, & + theta=0d0, & + smkl=0D0,smkl2=0D0, & + twoDTheta=0D0, & + planetheta, & + planephi, & + mmx=1d0, & + kky=0d0, & + llz=0d0, & + mmx1=0d0, & + kky1=0d0, & + llz1=0d0 + + + DOUBLE PRECISION :: a6666=0D0, & + sheainvar,CO=0D0, & + Minsheainvar, & + Maxsheainvar, & + sheainvar_ave, & + Minpughvar, & + Maxpughvar, & + pughvar_ave, & + bulk_m2, & + Min_bulkm2=0.0d0, & + Max_bulkm2=0.0d0, & + gamma, & + vv11=0D0, & + vv12=0D0, & + vv13=0D0, & + vv22=0D0, & + vv23=0D0, & + vv33=0D0, & + ave, & + k11, & + k12, & + k13, & + k22, & + k23, & + k33 DOUBLE PRECISION, PARAMETER :: PI_C=3.14159265358979323846264338327950D0 DOUBLE PRECISION :: density,theta3,phi3,ma_avrag=1,dns,tot_at=0.0 INTEGER :: k, n, m, ac_mod, plan DOUBLE PRECISION, DIMENSION (0:9,1:3) :: poin - DOUBLE PRECISION :: VVG_P, & - VVP_P=0D0, & - VV_P_PF=0D0, & - VVG_Sf=0D0, & - VVP_Sf=0D0, & - VV_Sf_PF=0D0, & - VVG_Ss=0D0, & - VVP_Ss=0D0, & - VV_Ss_PF=0d0, & - VVG_P_max=0D0, & - VVP_P_max=0D0, & - VV_P_PF_max=0D0, & - VVG_Sf_max=0D0, & - VVP_Sf_max=0D0, & - VV_Sf_PF_max=0D0, & - VVG_Ss_max=0D0, & - VVP_Ss_max=0D0, & - VV_Ss_PF_max=0D0, & - VVG_P_min=0D0, & - VVP_P_min=0D0, & - VV_P_PF_min=0D0, & - VVG_Sf_min=0D0, & - VVP_Sf_min=0D0, & - VV_Sf_PF_min=0D0, & - VVG_Ss_min=0D0, & - VVP_Ss_min=0.0001D0, & - VV_Ss_PF_min=0d0, & - VVG_P_max_phi=0D0, & - VVP_P_max_phi=0D0, & - VV_P_PF_max_phi=0D0, & - VVG_Sf_max_phi=0D0, & - VVP_Sf_max_phi=0D0, & - VV_Sf_PF_max_phi=0D0, & - VVG_Ss_max_phi=0D0, & - VVP_Ss_max_phi=0D0, & - VVG_P_min_phi= 0D0, & - VVP_P_min_phi=0D0, & - VV_P_PF_min_phi=0D0, & - VVG_Sf_min_phi=0D0, & - VVP_Sf_min_phi=0D0, & - VV_Sf_PF_min_phi=0D0, & - VVG_Ss_min_phi=0D0, & - VVP_Ss_min_phi=0D0, & - VV_Ss_PF_min_phi=0d0, & - VV_Ss_PF_max_phi=0d0, & - VVG_P_max_theta=0D0, & - VVP_P_max_theta=0D0, & - VV_P_PF_max_theta=0D0, & - VVG_Sf_max_theta=0D0, & - VVP_Sf_max_theta=0D0, & - VV_Sf_PF_max_theta=0D0, & - VVG_Ss_max_theta=0D0, & - VVP_Ss_max_theta=0D0, & - VVG_P_min_theta=0D0, & - VVP_P_min_theta=0D0, & - VV_P_PF_min_theta=0D0, & - VVG_Sf_min_theta=0D0, & - VVP_Sf_min_theta=0D0, & - VV_Sf_PF_min_theta=0D0, & - VVG_Ss_min_theta=0D0, & - VVP_Ss_min_theta=1D0, & - VV_Ss_PF_min_theta=0d0, & - VV_Ss_PF_max_theta=0d0, & - km , & - km_max=0.d0 , & - km_min=0.0d0 , & - km_min_phi=0.0d0 , & - km_max_phi=0.0d0 , & - km_min_theta=0.0d0 , & - km_max_theta=0.0d0 + DOUBLE PRECISION ::VVG_P, & + VVP_P=0D0, & + VV_P_PF=0D0, & + VVG_Sf=0D0, & + VVP_Sf=0D0, & + VV_Sf_PF=0D0, & + VVG_Ss=0D0, & + VVP_Ss=0D0, & + VV_Ss_PF=0d0, & + VVG_P_max=0D0, & + VVP_P_max=0D0, & + VV_P_PF_max=0D0, & + VVG_Sf_max=0D0, & + VVP_Sf_max=0D0, & + VV_Sf_PF_max=0D0, & + VVG_Ss_max=0D0, & + VVP_Ss_max=0D0, & + VV_Ss_PF_max=0D0, & + VVG_P_min=0D0, & + VVP_P_min=0D0, & + VV_P_PF_min=0D0, & + VVG_Sf_min=0D0, & + VVP_Sf_min=0D0, & + VV_Sf_PF_min=0D0, & + VVG_Ss_min=0D0, & + VVP_Ss_min=0.0001D0, & + VV_Ss_PF_min=0d0, & + VVG_P_max_phi=0D0, & + VVP_P_max_phi=0D0, & + VV_P_PF_max_phi=0D0, & + VVG_Sf_max_phi=0D0, & + VVP_Sf_max_phi=0D0, & + VV_Sf_PF_max_phi=0D0, & + VVG_Ss_max_phi=0D0, & + VVP_Ss_max_phi=0D0, & + VVG_P_min_phi= 0D0, & + VVP_P_min_phi=0D0, & + VV_P_PF_min_phi=0D0, & + VVG_Sf_min_phi=0D0, & + VVP_Sf_min_phi=0D0, & + VV_Sf_PF_min_phi=0D0, & + VVG_Ss_min_phi=0D0, & + VVP_Ss_min_phi=0D0, & + VV_Ss_PF_min_phi=0d0, & + VV_Ss_PF_max_phi=0d0, & + VVG_P_max_theta=0D0, & + VVP_P_max_theta=0D0, & + VV_P_PF_max_theta=0D0, & + VVG_Sf_max_theta=0D0, & + VVP_Sf_max_theta=0D0, & + VV_Sf_PF_max_theta=0D0, & + VVG_Ss_max_theta=0D0, & + VVP_Ss_max_theta=0D0, & + VVG_P_min_theta=0D0, & + VVP_P_min_theta=0D0, & + VV_P_PF_min_theta=0D0, & + VVG_Sf_min_theta=0D0, & + VVP_Sf_min_theta=0D0, & + VV_Sf_PF_min_theta=0D0, & + VVG_Ss_min_theta=0D0, & + VVP_Ss_min_theta=1D0, & + VV_Ss_PF_min_theta=0d0, & + VV_Ss_PF_max_theta=0d0, & + km , & + km_max=0.d0 , & + km_min=0.0d0 , & + km_min_phi=0.0d0 , & + km_max_phi=0.0d0 , & + km_min_theta=0.0d0 , & + km_max_theta=0.0d0 TYPE normals DOUBLE PRECISION, DIMENSION (1:3,1:3) :: normod END TYPE normals - TYPE(normals), DIMENSION(0:8):: normp - DOUBLE PRECISION, DIMENSION(3) :: a, b - DOUBLE PRECISION :: thetar, phir, supi, supt - DOUBLE PRECISION :: cutnorm + TYPE(normals), DIMENSION(0:8) :: normp + DOUBLE PRECISION, DIMENSION(3) :: a, b + DOUBLE PRECISION :: thetar, phir, supi, supt + DOUBLE PRECISION :: cutnorm TYPE m_matrix_YLM DOUBLE PRECISION, DIMENSION (0:360,0:360) :: m_matrix END TYPE m_matrix_YLM @@ -313,22 +321,84 @@ PROGRAM AAEP_main TYPE(m_matrix_YLM), DIMENSION(3) :: normal !> Normal to the slowness surface TYPE(m_matrix_YLM) :: pow_fact !> power flow angle END TYPE mod_cal - TYPE(mod_cal), DIMENSION(1:3) :: mv - DOUBLE PRECISION :: phi2, alpha - DOUBLE PRECISION, PARAMETER :: inc = 0.00001D0 - + TYPE(mod_cal), DIMENSION(1:3) :: mv + DOUBLE PRECISION :: phi2, alpha + DOUBLE PRECISION, PARAMETER :: inc = 0.00001D0 + + + !=================v1.7.3======================= + INTEGER :: iargc + INTEGER :: argcount + CHARACTER(LEN=80) :: wq_char, & + input_m, & + input_d, & + input_ewp, & + input_mtc, & + input_opt, & + help_optin + + DOUBLE PRECISION :: input_dens + ChARACTER(len=7), Dimension(15) :: arg_mane + !=================v1.7.3======================= !>>>>>>>>>>>>>>>>>>>>>>>>END<<<<<<<<<<<<<<<<<<<<<<<<< - OPEN(99,FILE="DATA.out") !> final static data CALL SYSTEM('clear') WRITE(*,*) ' ' CALL WELCOME() - WRITE(*,*) ' ' - CALL SYSTEM('sleep 1.2') - + CALL SYSTEM('sleep 1.1') + !=================v1.7.3======================= + input_d = "N" + input_m = "N" + input_ewp = "N" + input_mtc = "N" + input_opt = "N" + input_dens = 0.0 + argcount=iargc() + DO i=1,argcount +! CALL GETARG(i,wq_char) +! arg_mane(i)=wq_char + CALL GETARG(i,wq_char) + IF ( INDEX(wq_char,'-d').NE.0 ) THEN + CALL GETARG(i+1,wq_char) + READ(wq_char,*) input_d + END IF + ! + IF ( INDEX(wq_char,'-m').NE.0 ) THEN + CALL GETARG(i+1,wq_char) + READ(wq_char,*) input_m + END IF + ! + IF ( INDEX(wq_char,'-ewp').NE.0 ) THEN + CALL GETARG(i+1,wq_char) + READ(wq_char,*) input_ewp + END IF + ! + IF ( INDEX(wq_char,'-ktc').NE.0 ) THEN + CALL GETARG(i+1,wq_char) + READ(wq_char,*) input_mtc + END IF + ! + IF ( INDEX(wq_char,'-op').NE.0 ) THEN + CALL GETARG(i+1,wq_char) + READ(wq_char,*) input_opt + END IF + IF ( INDEX(wq_char,'-h').NE.0 ) THEN + CALL GETARG(i+1,wq_char) + READ(wq_char,*) help_optin + END IF + ! + END DO + !=================v1.7.3======================= + + IF (help_optin == "help") THEN + CALL HELP_me() + STOP + ENDIF 225 CALL SYSTEM('clear') +IF (input_d == "N") THEN !!!!! FOR COMMAND LINE NEW!! ====> start + WRITE(*,*)" > Select system dimension:" !> go to type of system CALL SYSTEM('tput setaf 36;tput bold; echo " ========================";tput sgr0') !WRITE(*,*)"========================" @@ -336,98 +406,222 @@ PROGRAM AAEP_main WRITE(*,*)" 3D-Materials------ => 3" !WRITE(*,*)"========================" CALL SYSTEM('tput setaf 36;tput bold; echo " ========================";tput sgr0') - READ(*,*) d2d3 - - -IF(d2d3 == 3) then !@@@@@@@@@@@@@@@@@@@@@@@ 2D_3D system start -665 CALL SYSTEM('clear') - WRITE(*,*)" > Select the desired method (code, file and databank):" - CALL SYSTEM('tput setaf 33;tput bold; echo "============================================================";tput sgr0') - !WRITE(*,*)"============================================================" - WRITE(*,*)" IRelast-----------------------( wein2k )-=> 1" - WRITE(*,*)" Elast-------------------------( wein2k )-=> 2" - WRITE(*,*)" AELAS-------------------------( VASP )-=> 3" - WRITE(*,*)" ElaStic-----------------------( QE,Wien2k,Exciting )-=> 4" - WRITE(*,*)" Using Cij Tensor in Cij.dat---( Other codes )-=> 5" - WRITE(*,*)" Using EC Databank-------------( offline MP )-=> 6" - WRITE(*,*)" Using EC Databank-------------( online MP )-=> 7" - Write(*,*)" Back --------------------------------------------------=> 0" - !WRITE(*,*)"============================================================" - CALL SYSTEM('tput setaf 33;tput bold; echo "============================================================";tput sgr0') - READ(*,*) Ncod - IF(Ncod .EQ. 0) then; Goto 225; endif - IF (Ncod .eq. 1 .OR. Ncod .eq. 2 .OR. Ncod .eq. 3 .OR. Ncod .eq. 4 .OR. Ncod .eq. 5 .OR. Ncod .eq. 6 .OR. Ncod .eq. 7) THEN - WRITE(*,*)" > Do you want to calculate elastic wave properties? (Y/n)" !> select code for calculate of phase and group velocities - READ(*,*)yn_veloc - CALL SYSTEM('clear') -IF (yn_veloc=='Y' .or. yn_veloc=='y') then - WRITE(*,*)" > Select the desired option:" - CALL SYSTEM('tput setaf 63;tput bold; echo "================================================================";tput sgr0') - WRITE(*,"(a)")" Phase, group and PFA without Min. thermal conductivity-----=> 1" - WRITE(*,"(a)")" Phase, group and PFA with Min. thermal conductivity-----=> 2" - Write(*,"(a)")" Back ------------------------------------------------------=> 0" - CALL SYSTEM('tput setaf 63;tput bold; echo "================================================================";tput sgr0') - Read(*,*) ewp -!####################################################################### - IF (ewp== 1) then - WRITE(*,*)"Density of Compound (kg/m^3):" - WRITE(*,*)"Note: If you don't know, enter 0" - READ(*,*) density - - IF (density==0)then - CALL system("dens.x dns") - OPEN(845,file=".rhocom") - READ(845,*)density - CLOSE(845) - END iF - WRITE (99,*) " > Density of Compound =", density,"(kg/m^3)" - ENdif -!####################################################################### - IF (ewp== 2) then - yn_km='y' - WRITE(*,*)"> Density of Compound (kg/m^3):" - WRITE(*,*)"Note: If you don't know, enter 0" - READ(*,*) density - - IF (density==0)then - CALL system("dens.x dns") - OPEN(845,file=".rhocom") - READ(845,*)density - CLOSE(845) - END iF -!!!!!!!!!!!!!! - open(124,file='.mavg_com') - - WRITE(*,*)"> Mean mass of atoms in each unit cell (gr):" - WRITE(*,*)"Note: If you don't know, enter 0" - READ(*,*) ma_avrag - Write(124,*) ma_avrag - WRITE(*,*)"> Total of atoms in each unit cell:" - WRITE(*,*)"Note: If you don't know, enter 0" - READ(*,*) tot_at - WRITE(124,*) tot_at - - close(124) - IF (ma_avrag==0 .or. tot_at==0)then - CALL system("dens.x mav") - - OPEN(414,file=".mavg_com") - READ(414,*) ma_avrag - READ(414,*) tot_at - CLOSE(414) - ENDIF - WRITE (99,*) " > Mean mass of atoms =", ma_avrag, "(gr)" - - WRITE (99,*) "" + READ(*,*) d2d3 + +ELSE !v1.7.3 + + IF (input_d == "3" ) THEN + d2d3=3 + WRITE(99,*) " > System Type: 3D" + ELSE IF(input_d == "2") THEN + d2d3=2 + WRITE(99,*) " > System Type: 2D" + ELSE + WRITE(*,"(2A)") "Invalid input: -d ", input_d + STOP + ENDIF +ENDIF + !=== !v1.7.3 +IF(d2d3 == 3) THEN !@@@@@@@@@@@@@@@@@@@@@@@ 2D_3D system start + 665 CALL SYSTEM('clear') + IF (input_m == "N") THEN !=== !v1.7.3 + WRITE(*,*)" > Select the desired method (code, file and databank):" + CALL SYSTEM('tput setaf 33;tput bold; echo "============================================================";tput sgr0') + !WRITE(*,*)"============================================================" + WRITE(*,*)" IRelast-----------------------( WEIN2k )-=> 1" + WRITE(*,*)" Elast-------------------------( WEIN2k )-=> 2" + WRITE(*,*)" AELAS-------------------------( VASP )-=> 3" + WRITE(*,*)" ElaStic-----------------------( QE,WEIN2k,Exciting )-=> 4" + WRITE(*,*)" Using Cij Tensor in Cij.dat---( Other codes )-=> 5" + WRITE(*,*)" Using EC Databank-------------( Offline MP )-=> 6" + WRITE(*,*)" Using EC Databank-------------( Online MP )-=> 7" + WRITE(*,*)" Back --------------------------------------------------=> 0" + !WRITE(*,*)"============================================================" + CALL SYSTEM('tput setaf 33;tput bold; echo "============================================================";tput sgr0') + READ(*,*) Ncod + IF(Ncod .EQ. 0) THEN; Goto 225; ENDIF + ELSE !=== !v1.7.3 + IF (input_m == "IRelast" .or. input_m == "irelast" .or. input_m == "1") THEN + Ncod=1 + WRITE(99,*) " > Output code: IRelast" + ELSE IF(input_m == "Elast" .or. input_m == "elast" .or. input_m == "2") THEN + Ncod=2 + WRITE(99,*) " > Output code: Elast" + ELSE IF(input_m == "AELAS" .or. input_m == "aelas" .or. input_m == "3") THEN + Ncod=3 + WRITE(99,*) " > Output code: AELAS" + ELSE IF(input_m == "ElaStic" .or. input_m == "elastic" .or. input_m == "4") THEN + Ncod=4 + WRITE(99,*) " > Output code: ElaStic" + ELSE IF(input_m == "Cij" .or. input_m == "cij" .or. input_m == "5") THEN + Ncod=5 + WRITE(99,*) " > Output code: Cij file" + ELSE IF(input_m == "offmp" .or. input_m == "6") THEN + Ncod=6 + WRITE(99,*) " > Output code: offline Databank" + ELSE IF(input_m == "onmp" .or. input_m == "7") THEN + Ncod=7 + WRITE(99,*) " > Output code: Online Databank" + ELSE + WRITE(*,"(2A)") "Invalid input: -m ", input_m + STOP + ENDIF + ENDIF + !=== !v1.7.3 + IF(input_ewp == "n" .and. input_mtc == "n")THEN !=== !v1.7.3 + WRITE(*,"(A,A4,A,A4)") "Invalid input: -ewp ", input_ewp, "-ktc ", input_mtc + STOP ENDIF - !####################################################################### - IF (ewp== 0) then - Goto 665 - Endif + IF(input_ewp == "n" .and. input_mtc == "y")THEN !=== !v1.7.3 + WRITE(*,"(A,A4,A,A4)") "Invalid input: -ewp ", input_ewp, "-ktc ", input_mtc + STOP + ENDIF + IF(input_ewp == "N" .and. input_mtc == "y")THEN !=== !v1.7.3 + WRITE(*,"(A,A4,A,A4)") "Invalid input: -ewp ", input_ewp, "-ktc ", input_mtc + STOP + ENDIF + IF(input_ewp == "n" .and. input_mtc == "N") THEN !=== !v1.7.3 + !WRITE(*,"(A,A4,A,A4)") "Invalid input: -ewp ", input_ewp, "-ktc ", input_mtc + yn_veloc = "n" + goto 5501 + ENDIF + IF(input_ewp == "N" .or. input_mtc == "N") THEN !=== !v1.7.3 + WRITE(99,*) "Calculate elastic wave properties: On" + IF (Ncod .eq. 1 .OR. Ncod .eq. 2 .OR. Ncod .eq. 3 .OR. Ncod .eq. 4 .OR. Ncod .eq. 5 .OR. Ncod .eq. 6 .OR. Ncod .eq. 7) THEN + WRITE(*,*)" > Do you want to calculate elastic wave properties? (Y/n):" !> select code for calculate of phase and group velocities + READ(*,*) yn_veloc + CALL SYSTEM('clear') + + IF (yn_veloc == 'Y' .or. yn_veloc == 'y') THEN + WRITE(*,*)" > Select the desired option:" + CALL SYSTEM('tput setaf 63;tput bold; echo "================================================================";tput sgr0') + WRITE(*,"(a)")" Phase, group and PFA without Min. thermal conductivity-----=> 1" + WRITE(*,"(a)")" Phase, group and PFA with Min. thermal conductivity-----=> 2" + WRITE(*,"(a)")" Back ------------------------------------------------------=> 0" + CALL SYSTEM('tput setaf 63;tput bold; echo "================================================================";tput sgr0') + Read(*,*) ewp + !####################################################################### + IF (ewp == 1) THEN + WRITE(99,*) "Calculate Min. thermal conductivity: Off" + WRITE(*,*)"Density of Compound (kg/m^3):" + WRITE(*,*)"Note: IF you don't know, enter 0" + READ(*,*) density + IF (density==0)THEN + CALL system("dens.x dns") + OPEN(845,file=".rhocom") + READ(845,*)density + CLOSE(845) + END IF + WRITE (99,*) " > Density of Compound =", density,"(kg/m^3)" + ENDIF + !####################################################################### + IF (ewp== 2) THEN + WRITE(99,*) " > Calculate Min. thermal conductivity: On" + yn_km='y' + WRITE(*,*)"> Density of Compound (kg/m^3):" + WRITE(*,*)"Note: IF you don't know, enter 0" + READ(*,*) density + IF (density==0)THEN + CALL system("dens.x dns") + OPEN(845,file=".rhocom") + READ(845,*)density + CLOSE(845) + END IF + !=================================== + open(124,file='.mavg_com') + WRITE(*,*)"> Mean mass of atoms in each unit cell (gr):" + WRITE(*,*)"Note: IF you don't know, enter 0" + READ(*,*) ma_avrag + WRITE(124,*) ma_avrag + WRITE(*,*)"> Total of atoms in each unit cell:" + WRITE(*,*)"Note: IF you don't know, enter 0" + READ(*,*) tot_at + WRITE(124,*) tot_at + Close(124) + IF (ma_avrag==0 .or. tot_at==0)THEN + CALL system("dens.x mav") + OPEN(414,file=".mavg_com") + READ(414,*) ma_avrag + READ(414,*) tot_at + CLOSE(414) + ENDIF + WRITE (99,*) " > Mean mass of atoms =", ma_avrag, "(gr)" + WRITE (99,*) "" + ENDIF + !####################################################################### + IF (ewp== 0) THEN + Goto 665 + ENDIF + ENDIF + ENDIF + + ELSE !=== !v1.7.3 + + + IF (input_ewp == "y" .and. input_mtc == "n" .or. input_mtc == "N") THEN !=== !v1.7.3 + yn_veloc = "y" + ewp = 1 + WRITE(99,*) " > Calculate Min. thermal conductivity: Off" + !####################################################################### + WRITE(*,*)"Density of Compound (kg/m^3):" + WRITE(*,*)"Note: IF you don't know, enter 0" + READ(*,*) density -ENDIF - !!!!!!!!!!!!!!!!!!! - IF (Ncod .eq. 1) then + IF (density==0)THEN + CALL system("dens.x dns") + OPEN(845,file=".rhocom") + READ(845,*)density + CLOSE(845) + END IF + WRITE (99,*) " > Density of Compound =", density,"(kg/m^3)" + !####################################################################### + ENDIF + IF (input_ewp == "y" .and. input_mtc == "y" ) THEN + yn_veloc = "y" + ewp = 2 + WRITE(99,*) " > Calculate Min. thermal conductivity: On" + !####################################################################### + yn_km='y' + WRITE(*,*)"> Density of Compound (kg/m^3):" + WRITE(*,*)"Note: IF you don't know, enter 0" + READ(*,*) density + + IF (density==0)THEN + CALL system("dens.x dns") + OPEN(845,file=".rhocom") + READ(845,*)density + CLOSE(845) + END IF + !!!!!!!!!!!!!! + open(124,file='.mavg_com') + + WRITE(*,*)"> Mean mass of atoms in each unit cell (gr):" + WRITE(*,*)"Note: IF you don't know, enter 0" + READ(*,*) ma_avrag + WRITE(124,*) ma_avrag + WRITE(*,*)"> Total of atoms in each unit cell:" + WRITE(*,*)"Note: IF you don't know, enter 0" + READ(*,*) tot_at + WRITE(124,*) tot_at + + close(124) + IF(ma_avrag==0 .or. tot_at==0)THEN + CALL system("dens.x mav") + OPEN(414,file=".mavg_com") + READ(414,*) ma_avrag + READ(414,*) tot_at + CLOSE(414) + ENDIF + WRITE (99,*) " > Mean mass of atoms =", ma_avrag, "(gr)" + WRITE (99,*) "" + ENDIF + + ENDIF + + + + !=================================== +5501 IF (Ncod .eq. 1) THEN OPEN(4,FILE="INVELC-matrix",status='old',err=1367) CALL system ('cp INVELC-matrix Cij.dat') CALL C_Inv_M(6) @@ -436,15 +630,15 @@ PROGRAM AAEP_main CLOSE(4) CALL system('clear') ENDIF - IF (Ncod .eq. 2) then + IF (Ncod .eq. 2) THEN OPEN(5,FILE="elast.output",status='old',err=1366) - OPEN(6,FILE=".elsetcode") + OPEN(6,FILE=".ELSEtcode") OPEN(7,FILE="Cij.dat",status='old',err=1369) - CALL system ('cat elast.output | tail -1| cut -f2 -d=| cut -f1 -dc > .elsetcode ') !C11 - CALL system ('cat elast.output | tail -1| cut -f3 -d=| cut -f1 -dc >> .elsetcode ') !C12 - CALL system ('cat elast.output | tail -1| cut -f4 -d=| cut -f1 -dc >> .elsetcode ') !C44 + CALL system ('cat elast.output | tail -1| cut -f2 -d=| cut -f1 -dc > .ELSEtcode ') !C11 + CALL system ('cat elast.output | tail -1| cut -f3 -d=| cut -f1 -dc >> .ELSEtcode ') !C12 + CALL system ('cat elast.output | tail -1| cut -f4 -d=| cut -f1 -dc >> .ELSEtcode ') !C44 DO i2=1,3 - READ(6,*) CElastcode(i2) + READ(6,*) CElastcode(i2) ENDDO CLOSE(6) C(1,1)=CElastcode(1) ; C(1,2)=CElastcode(2);C(4,4)=CElastcode(3) @@ -465,8 +659,7 @@ PROGRAM AAEP_main CALL system('clear') ENDIF CALL system('clear') - - IF (Ncod .EQ. 3) then + IF (Ncod .EQ. 3) THEN CALL system("sed '1,2d' ELADAT > ELADAT_temp") OPEN(79,FILE="Cij.dat") OPEN(53,FILE="ELADAT_temp",status='old',err=1361) @@ -476,176 +669,162 @@ PROGRAM AAEP_main ENDDO CLOSE(53) CLOSE(79) - !call sleep(5) + !CALL sleep(5) ENDIF - + 10 format(6F11.1) IF (Ncod.EQ.4) THEN - write(*,*)"> The program does not support 3rd order of the elastic constant!" - write(*,*)"Do you continue? (Y/n)" + WRITE(*,*)"> The program does not support 3rd order of the elastic constant!" + WRITE(*,*)"Do you continue? (Y/n)" READ(*,*)order - IF(order=='y'.or.order=='Y') then - CALL system("sed '1,15d' ElaStic_2nd.out > ElaStic_2nd_temp") - OPEN(79,FILE="Cij.dat") - OPEN(53,FILE="ElaStic_2nd_temp",status='old',err=1341) - DO i=1,6 - READ(53,*) (C3(i,j),j=1,6) - WRITE(79,10) (C3(i,j),j=1,6) - ENDDO - CLOSE(53) - CLOSE(79) - else - stop - end if - end if - - IF (Ncod .EQ. 5) then - GOTO 101 - ENDIF - - IF (Ncod .EQ. 6) then - 5050 WRITE(*,*)"> ENTER ID: (EXAMPEL: mp-10 or mvc-916)" - WRITE(*,*)"============================= offline" - read(*,*)myid - CALL databank(myid,yesno) - - IF (yesno=='N')then - WRITE(*,*)"----------------------------------------------" - WRITE(*,*)"> Want to repeat again?(Y/n)" - READ(*,*)yesno2 - IF (yesno2=='Y'.or. yesno2=='y')then - CALL SYSTEM('clear') - GOTO 5050 + IF(order=='y'.or.order=='Y') THEN + CALL system("sed '1,15d' ElaStic_2nd.out > ElaStic_2nd_temp") + OPEN(79,FILE="Cij.dat") + OPEN(53,FILE="ElaStic_2nd_temp",status='old',err=1341) + DO i=1,6 + READ(53,*) (C3(i,j),j=1,6) + WRITE(79,10) (C3(i,j),j=1,6) + ENDDO + CLOSE(53) + CLOSE(79) ELSE - GOTO 1370 - END IF - END IF + stop + end IF + End IF + IF (Ncod .EQ. 5) THEN + GOTO 101 + ENDIF + + IF (Ncod .EQ. 6) THEN + 5050 WRITE(*,*)"> ENTER ID: (EXAMPEL: mp-10 or mvc-916)" + WRITE(*,*)"============================= offline" + read(*,*)myid + CALL databank(myid,yesno) - CALL sleep (2) + IF (yesno=='N')THEN + WRITE(*,*)"----------------------------------------------" + WRITE(*,*)"> Want to repeat again?(Y/n)" + READ(*,*)yesno2 + IF (yesno2=='Y'.or. yesno2=='y')THEN + CALL SYSTEM('clear') + GOTO 5050 + ELSE + GOTO 1370 + END IF + END IF + CALL sleep (2) + CALL system('clear') + CALL system('mv Cij-id.dat Cij.dat') + GOTO 101 + ENDIF + + IF (Ncod .EQ. 7) THEN + 5051 WRITE(*,*)"> ENTER ID: (EXAMPEL: mp-10 or mvc-916)" + WRITE(*,*)"============================= online" + read(*,*)myid + CALL aip_get_online(myid) + CALL system ("chmod +x aip.py; ./aip.py > Cij-id.dat") + IF (yesno=='N')THEN + WRITE(*,*)"----------------------------------------------" + WRITE(*,*)"> Want to repeat again?(Y/n)" + READ(*,*)yesno2 + IF (yesno2=='Y'.or. yesno2=='y')THEN + CALL SYSTEM('clear') + GOTO 5051 + ELSE + GOTO 1370 + END IF + END IF + + CALL sleep (2) + CALL system('clear') + CALL system('mv Cij-id.dat Cij.dat') + CALL system("rm aip.py") + GOTO 101 + ENDIF CALL system('clear') - CALL system('mv Cij-id.dat Cij.dat') - GOTO 101 - ENDIF - - - -IF (Ncod .EQ. 7) then - 5051 WRITE(*,*)"> ENTER ID: (EXAMPEL: mp-10 or mvc-916)" - WRITE(*,*)"============================= online" - read(*,*)myid - CALL aip_get_online(myid) - CALL system ("chmod +x aip.py; ./aip.py > Cij-id.dat") - IF (yesno=='N')then - WRITE(*,*)"----------------------------------------------" - WRITE(*,*)"> Want to repeat again?(Y/n)" - READ(*,*)yesno2 - IF (yesno2=='Y'.or. yesno2=='y')then - CALL SYSTEM('clear') - GOTO 5051 - ELSE - GOTO 1370 + 101 OPEN(11,FILE="Cij.dat",status='old', err=1369) ! read cij data inpout + READ(11,*) C(1,1),C(1,2),C(1,3),C(1,4),C(1,5),C(1,6) + READ(11,*) C(2,1),C(2,2),C(2,3),C(2,4),C(2,5),C(2,6) + READ(11,*) C(3,1),C(3,2),C(3,3),C(3,4),C(3,5),C(3,6) + READ(11,*) C(4,1),C(4,2),C(4,3),C(4,4),C(4,5),C(4,6) + READ(11,*) C(5,1),C(5,2),C(5,3),C(5,4),C(5,5),C(5,6) + READ(11,*) C(6,1),C(6,2),C(6,3),C(6,4),C(6,5),C(6,6) + CLOSE(11) + CALL C_Inv_M(6) + CALL system('clear') + !! for other code <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + OPEN(12,FILE="Sij.dat") ! read sij data input (generated by C_Inv_M(6) or other data and code) + READ(12,*) S1(1,1),S1(1,2),S1(1,3),S1(1,4),S1(1,5),S1(1,6) + READ(12,*) S1(2,1),S1(2,2),S1(2,3),S1(2,4),S1(2,5),S1(2,6) + READ(12,*) S1(3,1),S1(3,2),S1(3,3),S1(3,4),S1(3,5),S1(3,6) + READ(12,*) S1(4,1),S1(4,2),S1(4,3),S1(4,4),S1(4,5),S1(4,6) + READ(12,*) S1(5,1),S1(5,2),S1(5,3),S1(5,4),S1(5,5),S1(5,6) + READ(12,*) S1(6,1),S1(6,2),S1(6,3),S1(6,4),S1(6,5),S1(6,6) + CLOSE(12) + S(1,1)=S1(1,1); S(1,2)=S1(1,2); S(1,3)=S1(1,3); S(1,4)=S1(1,4); S(1,5)=S1(1,5); S(1,6)=S1(1,6) ! TESTING.... + S(2,1)=S1(2,1); S(2,2)=S1(2,2); S(2,3)=S1(2,3); S(2,4)=S1(2,4); S(2,5)=S1(2,5); S(2,6)=S1(2,6) + S(3,1)=S1(3,1); S(3,2)=S1(3,2); S(3,3)=S1(3,3); S(3,4)=S1(3,4); S(3,5)=S1(3,5); S(3,6)=S1(3,6) + S(4,1)=S1(4,1); S(4,2)=S1(4,2); S(4,3)=S1(4,3); S(4,4)=S1(4,4); S(4,5)=S1(4,5); S(4,6)=S1(4,6) + S(5,1)=S1(5,1); S(5,2)=S1(5,2); S(5,3)=S1(5,3); S(5,4)=S1(5,4); S(5,5)=S1(5,5); S(5,6)=S1(5,6) + S(6,1)=S1(6,1); S(6,2)=S1(6,2); S(6,3)=S1(6,3); S(6,4)=S1(6,4); S(6,5)=S1(6,5); S(6,6)=S1(6,6) + !> Stability + CALL stability3d( Stable, C,Eig3d ) + IF (Stable==1) THEN + CALL system ('tput setaf 9;tput bold; echo " =======================================";tput sgr0') + CALL system ('tput setaf 9;tput bold; echo " > Elastic Stability Conditions: Unstable; STOP";tput sgr0') + CALL system ('tput setaf 9;tput bold; echo " =======================================";tput sgr0') + WRITE(99,*)" ========================================" + WRITE(99,"(a)")" > Elastic Stability Conditions: Unstable; STOP" + WRITE(99,*)" ========================================" END IF - END IF - - CALL sleep (2) - CALL system('clear') - CALL system('mv Cij-id.dat Cij.dat') - CALL system("rm aip.py") - GOTO 101 -ENDIF -ELSE - WRITE(*,*)"Invalid input!"; go to 1370 -ENDIF - - -CALL system('clear') -101 OPEN(11,FILE="Cij.dat",status='old', err=1369) ! read cij data inpout -READ(11,*) C(1,1),C(1,2),C(1,3),C(1,4),C(1,5),C(1,6) -READ(11,*) C(2,1),C(2,2),C(2,3),C(2,4),C(2,5),C(2,6) -READ(11,*) C(3,1),C(3,2),C(3,3),C(3,4),C(3,5),C(3,6) -READ(11,*) C(4,1),C(4,2),C(4,3),C(4,4),C(4,5),C(4,6) -READ(11,*) C(5,1),C(5,2),C(5,3),C(5,4),C(5,5),C(5,6) -READ(11,*) C(6,1),C(6,2),C(6,3),C(6,4),C(6,5),C(6,6) -CLOSE(11) -CALL C_Inv_M(6) -CALL system('clear') - !! for other code <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< -OPEN(12,FILE="Sij.dat") ! read sij data input (generated by C_Inv_M(6) or other data and code) -READ(12,*) S1(1,1),S1(1,2),S1(1,3),S1(1,4),S1(1,5),S1(1,6) -READ(12,*) S1(2,1),S1(2,2),S1(2,3),S1(2,4),S1(2,5),S1(2,6) -READ(12,*) S1(3,1),S1(3,2),S1(3,3),S1(3,4),S1(3,5),S1(3,6) -READ(12,*) S1(4,1),S1(4,2),S1(4,3),S1(4,4),S1(4,5),S1(4,6) -READ(12,*) S1(5,1),S1(5,2),S1(5,3),S1(5,4),S1(5,5),S1(5,6) -READ(12,*) S1(6,1),S1(6,2),S1(6,3),S1(6,4),S1(6,5),S1(6,6) -CLOSE(12) -S(1,1)=S1(1,1); S(1,2)=S1(1,2); S(1,3)=S1(1,3); S(1,4)=S1(1,4); S(1,5)=S1(1,5); S(1,6)=S1(1,6) ! TESTING.... -S(2,1)=S1(2,1); S(2,2)=S1(2,2); S(2,3)=S1(2,3); S(2,4)=S1(2,4); S(2,5)=S1(2,5); S(2,6)=S1(2,6) -S(3,1)=S1(3,1); S(3,2)=S1(3,2); S(3,3)=S1(3,3); S(3,4)=S1(3,4); S(3,5)=S1(3,5); S(3,6)=S1(3,6) -S(4,1)=S1(4,1); S(4,2)=S1(4,2); S(4,3)=S1(4,3); S(4,4)=S1(4,4); S(4,5)=S1(4,5); S(4,6)=S1(4,6) -S(5,1)=S1(5,1); S(5,2)=S1(5,2); S(5,3)=S1(5,3); S(5,4)=S1(5,4); S(5,5)=S1(5,5); S(5,6)=S1(5,6) -S(6,1)=S1(6,1); S(6,2)=S1(6,2); S(6,3)=S1(6,3); S(6,4)=S1(6,4); S(6,5)=S1(6,5); S(6,6)=S1(6,6) -!> Stability -call stability3d( Stable, C,Eig3d ) -IF (Stable==1) THEN - CALL system ('tput setaf 9;tput bold; echo " =======================================";tput sgr0') - CALL system ('tput setaf 9;tput bold; echo " > Elastic Stability Conditions: Unstable; STOP";tput sgr0') - CALL system ('tput setaf 9;tput bold; echo " =======================================";tput sgr0') - Write(99,*)" ========================================" - WRITE(99,"(a)")" > Elastic Stability Conditions: Unstable; STOP" - Write(99,*)" ========================================" - STOP -END IF -IF (Stable==0) THEN - CALL system ('tput setaf 10;tput bold; echo " =======================================";tput sgr0') - CALL system ('tput setaf 10;tput bold; echo " > Elastic Stability Conditions: Stable";tput sgr0') - CALL system ('tput setaf 10;tput bold; echo " =======================================";tput sgr0') - Write(99,*)" ========================================" - WRITE(99,"(a)")" > Elastic Stability Conditions: Stable" - Write(99,*)" ========================================" - -END IF - WRITE(*,'(A,36f10.2)') " > Eigenvalues (GPa):",Eig3d(1,1),Eig3d(2,2),Eig3d(3,3) ,Eig3d(4,4) ,Eig3d(5,5) ,Eig3d(6,6) - WRITE(99,'(A,36f10.2)') " > Eigenvalues (GPa):",Eig3d(1,1),Eig3d(2,2),Eig3d(3,3) ,Eig3d(4,4) ,Eig3d(5,5) ,Eig3d(6,6) - WRITE(99,*)" " -CALL sleep(1) - !< -CALL proelast() - - -IF (yn_veloc=='Y' .or. yn_veloc=='y') WRITE (*,*) " > Density of Compound =", density,"(kg/m^3)" -IF (ewp== 2) WRITE (*,*) " > Mean mass of atoms =", ma_avrag, "(gr)" -IF (ewp== 2) CALL pro_wave(density,ma_avrag) -WRITE(*,*)"" -CALL system ('tput setaf 122;tput bold; echo " > Enter phi-meah and theta-meah between 50 and 250 (Recom.: 150*150):";tput sgr0') -!WRITE(*,*)" > Enter phi-meah and theta-meah between 50 and 250 (Recommended: 150):" -OPEN(59,file='MESH') -READ(*,*)phi_meah,theta_meah - cutmesh=(phi_meah*theta_meah)+1 - WRITE(59,*)phi_meah,theta_meah,cutmesh - close(59) - -WRITE (*,*) " " -CALL system ('tput setaf 142;tput bold; echo " > Select the (hkl) Miller indices for 2D cut:";tput sgr0') -!WRITE(*,*)"> Select the (hkl) Miller indices for 2D cut:" -READ(*,*) mmx,kky,llz -mmx1=mmx -kky1=kky -llz1=llz -CALL findpaln(mmx,kky,llz, e1,e2,yn_veloc) -!WRITE(*,*)mmx,kky,llz -!write(*,*)mmx,kky,llz -CALL sleep (1) -IF (mmx.EQ.0 .AND. kky.EQ.0 .AND. llz.EQ.0) then - mmx=1 - kky=0 - llz=0 - WRITE(*,*)'(0 0 0) is not accepted! => Chenged to (1 0 0).' - CALL sleep (1) -ENDIF -CALL system('clear; sleep 1.5') + IF (Stable==0) THEN + CALL system ('tput setaf 10;tput bold; echo " =======================================";tput sgr0') + CALL system ('tput setaf 10;tput bold; echo " > Elastic Stability Conditions: Stable";tput sgr0') + CALL system ('tput setaf 10;tput bold; echo " =======================================";tput sgr0') + WRITE(99,*)" ========================================" + WRITE(99,"(a)")" > Elastic Stability Conditions: Stable" + WRITE(99,*)" ========================================" + END IF + WRITE(*,'(A,36f10.2)') " > Eigenvalues (GPa):",Eig3d(1,1),Eig3d(2,2),Eig3d(3,3) ,Eig3d(4,4) ,Eig3d(5,5) ,Eig3d(6,6) + WRITE(99,'(A,36f10.2)') " > Eigenvalues (GPa):",Eig3d(1,1),Eig3d(2,2),Eig3d(3,3) ,Eig3d(4,4) ,Eig3d(5,5) ,Eig3d(6,6) + WRITE(99,*)" " + CALL sleep(1) + !< + CALL proelast() + + IF (yn_veloc=='Y' .or. yn_veloc=='y') WRITE (*,*) " > Density of Compound =", density,"(kg/m^3)" + IF (ewp== 2) WRITE (*,*) " > Mean mass of atoms =", ma_avrag, "(gr)" + IF (ewp== 2) CALL pro_wave(density,ma_avrag) + WRITE(*,*)"" + CALL system ('tput setaf 122;tput bold; echo " > Enter phi-meah and theta-meah between 50 and 250 (Recom.: 150*150):";tput sgr0') + !WRITE(*,*)" > Enter phi-meah and theta-meah between 50 and 250 (Recommended: 150):" + OPEN(59,file='MESH') + READ(*,*)phi_meah,theta_meah + cutmesh=(phi_meah*theta_meah)+1 + WRITE(59,*)phi_meah,theta_meah,cutmesh + close(59) + WRITE (*,*) " " + CALL system ('tput setaf 142;tput bold; echo " > Select the (hkl) Miller indices for 2D cut:";tput sgr0') + !WRITE(*,*)"> Select the (hkl) Miller indices for 2D cut:" + READ(*,*) mmx,kky,llz + mmx1=mmx + kky1=kky + llz1=llz + CALL findpaln(mmx,kky,llz, e1,e2,yn_veloc) + !WRITE(*,*)mmx,kky,llz + !WRITE(*,*)mmx,kky,llz + CALL sleep (1) + IF (mmx.EQ.0 .AND. kky.EQ.0 .AND. llz.EQ.0) THEN + mmx=1 + kky=0 + llz=0 + WRITE(*,*)'(0 0 0) is not accepted! => Chenged to (1 0 0).' + CALL sleep (1) + ENDIF + CALL system('clear; sleep 1.5') ELSE - IF(d2d3 == 2) then - + IF(d2d3 == 2) THEN OPEN(140,file="HKL") WRITE(140,*)"X" WRITE(140,*)"Y" @@ -654,34 +833,60 @@ PROGRAM AAEP_main WRITE(140,"(I2)")1 WRITE(140,*)'N' close(140) -226 call system('clear') - WRITE(*,*)" > Select the desired method (code or file):" - CALL SYSTEM('tput setaf 12;tput bold; echo " =====================================================";tput sgr0') - !WRITE(*,*)"=====================================================" - WRITE(*,*)" AELAS ( VASP ) => 1" - WRITE(*,*)" IRelast2D ( WIEN2K ) => 2" - WRITE(*,*)" Using Cij Tensor in Cij-2D.dat (other codes) => 3" - WRITE(*,*)" Back ----------------------------------------- => 0" - !WRITE(*,*)"=====================================================" - CALL SYSTEM('tput setaf 12;tput bold; echo " =====================================================";tput sgr0') - read(*,*) Ncod - IF(Ncod .EQ. 0) then; Goto 225; endif - - IF (Ncod .EQ. 1) then - call system('clear') - WRITE(*,*)" > Select the type of two-dimensional system:" - CALL SYSTEM('tput setaf 41;tput bold; echo " ====================================================";tput sgr0') - WRITE(*,*) " Default option (Hex., Squ., and Rec. systems) => 1 " - WRITE(*,*) " Advanced option ( Oblique systems ) => 2 " - WRITE(*,*) " Back ------------------------------------------=> 0 " - CALL SYSTEM('tput setaf 41;tput bold; echo " ====================================================";tput sgr0') - read(*,*)adv_mubner + 226 CALL system('clear') + IF (input_m == "N" .AND. input_d == "2")THEN + WRITE(*,*)" > Select the desired method (code or file):" + CALL SYSTEM('tput setaf 12;tput bold; echo " =====================================================";tput sgr0') + !WRITE(*,*)"=====================================================" + WRITE(*,*)" AELAS ( VASP ) => 1" + WRITE(*,*)" IRelast2D ( WIEN2K ) => 2" + WRITE(*,*)" Using Cij Tensor in Cij-2D.dat (other codes) => 3" + WRITE(*,*)" Back ----------------------------------------- => 0" + !WRITE(*,*)"=====================================================" + CALL SYSTEM('tput setaf 12;tput bold; echo " =====================================================";tput sgr0') + read(*,*) Ncod + IF(Ncod .EQ. 0) THEN; Goto 225; ENDIF + ELSE + IF (input_m == "AELAS" .or. input_m == "aelas" .or. input_m == "1") THEN + Ncod=1 + !WRITE(*,*) "Output code: Elast" + ELSE IF(input_m == "IRelast2D" .or. input_m == "irelast2D" .or. input_m == "2") THEN + Ncod=2 + !WRITE(*,*) "Output code: Elast" + ELSE IF(input_m == "Cij" .or. input_m == "cij" .or. input_m == "3") THEN + Ncod=3 + ELSE + WRITE(*,"(2A)") "Invalid input: -m ", input_m + STOP + ENDIF + ENDIF + IF (Ncod .EQ. 1) THEN + CALL system('clear') + IF(input_opt == "N" )THEN + WRITE(*,*)" > Select the type of two-dimensional system:" + CALL SYSTEM('tput setaf 41;tput bold; echo " ====================================================";tput sgr0') + WRITE(*,*) " Default option (Hex., Squ., and Rec. systems) => 1 " + WRITE(*,*) " Advanced option ( Oblique systems ) => 2 " + WRITE(*,*) " Back ------------------------------------------=> 0 " + CALL SYSTEM('tput setaf 41;tput bold; echo " ====================================================";tput sgr0') + read(*,*)adv_mubner + ELSE + IF (input_opt == "def" .or. input_opt == "1") THEN + adv_mubner=1 + !WRITE(*,*) "Output code: Elast" + ELSE IF(input_opt == "adv" .or. input_opt == "2") THEN + adv_mubner=2 + ELSE + WRITE(*,"(2A)") "Invalid input: -op ", adv_mubner + STOP + ENDIF + ENDIF IF (adv_mubner == 2 ) adv = "adv" IF (adv_mubner == 1 ) adv = "ndv" IF (adv_mubner == 0 ) goto 226 - call system("sed '1,2d' ELADAT > ELADAT_temp") + CALL system("sed '1,2d' ELADAT > ELADAT_temp") OPEN(79,FILE="Cij-2D.dat") OPEN(53,FILE="ELADAT_temp",status='old',err=1361) DO i=1,3 @@ -692,15 +897,15 @@ PROGRAM AAEP_main close(79) !CALL C_Inv_M2D(3) - !> Stability + !> Stability CALL stability2d( Stable, C2D,Eig2d ) IF (Stable == 1) THEN CALL system ('tput setaf 9;tput bold; echo " =======================================";tput sgr0') CALL system ('tput setaf 9;tput bold; echo " > Elastic Stability Conditions: Unstable; STOP";tput sgr0') CALL system ('tput setaf 9;tput bold; echo " =======================================";tput sgr0') - Write(99,*)" ========================================" - Write(99,*) " > Elastic Stability Conditions: Unstable" - Write(99,*)" ========================================" + WRITE(99,*)" ========================================" + WRITE(99,*) " > Elastic Stability Conditions: Unstable" + WRITE(99,*)" ========================================" STOP END IF @@ -708,34 +913,46 @@ PROGRAM AAEP_main CALL system ('tput setaf 10;tput bold; echo " =======================================";tput sgr0') CALL system ('tput setaf 10;tput bold; echo " > Elastic Stability Conditions: Stable";tput sgr0') CALL system ('tput setaf 10;tput bold; echo " =======================================";tput sgr0') - Write(99,*)" ========================================" - Write(99,*) " > Elastic Stability Conditions: Stable" - Write(99,*)" ========================================" + WRITE(99,*)" ========================================" + WRITE(99,*) " > Elastic Stability Conditions: Stable" + WRITE(99,*)" ========================================" END IF WRITE(*,'(A,36f10.2)')" > Eigenvalues (N/m):", Eig2d(1,1),Eig2d(2,2),Eig2d(3,3) WRITE(99,'(A,36f10.2)')" > Eigenvalues (N/m):", Eig2d(1,1),Eig2d(2,2),Eig2d(3,3) WRITE(99,*)" " - !< + !< CALL C_Inv_M2D(3) CALL proelast_2D () CALL sleep(2) - write(*,*)"" + WRITE(*,*)"" CALL system ('tput setaf 122;tput bold; echo " > Enter phi-mesh. Be divisible by 100 (e.g. 100, 200, 300, etc):";tput sgr0') !WRITE(*,*)" > Enter phi-mesh. Be divisible by 100.(e.g. 100, 200, 300, etc):" read(*,*)npoint ENDIF - IF (Ncod .eq. 2) then !> IRelast2D - - call system('clear') - WRITE(*,*)" > Select the type of two-dimensional system:" - CALL SYSTEM('tput setaf 41;tput bold; echo " ====================================================";tput sgr0') - WRITE(*,*) " Default option (Hex., Squ., and Rec. systems) => 1 " - WRITE(*,*) " Advanced option ( Oblique systems ) => 2 " - WRITE(*,*) " Back ------------------------------------------=> 0 " - CALL SYSTEM('tput setaf 41;tput bold; echo " ====================================================";tput sgr0') - read(*,*)adv_mubner + IF (Ncod .eq. 2) THEN !> IRelast2D + + CALL system('clear') + IF(input_opt == "N" )THEN + WRITE(*,*)" > Select the type of two-dimensional system:" + CALL SYSTEM('tput setaf 41;tput bold; echo " ====================================================";tput sgr0') + WRITE(*,*) " Default option (Hex., Squ., and Rec. systems) => 1 " + WRITE(*,*) " Advanced option ( Oblique systems ) => 2 " + WRITE(*,*) " Back ------------------------------------------=> 0 " + CALL SYSTEM('tput setaf 41;tput bold; echo " ====================================================";tput sgr0') + read(*,*)adv_mubner + ELSE + IF (input_opt == "def" .or. input_opt == "1") THEN + adv_mubner=1 + !WRITE(*,*) "Output code: Elast" + ELSE IF(input_opt == "adv" .or. input_opt == "2") THEN + adv_mubner=2 + ELSE + WRITE(*,"(2A)") "Invalid input: -op ", adv_mubner + STOP + ENDIF + ENDIF IF (adv_mubner == 2 ) adv = "adv" IF (adv_mubner == 1 ) adv = "ndv" @@ -758,41 +975,43 @@ PROGRAM AAEP_main ENDDO close(79) !CALL C_Inv_M2D(3) -!> Stability + !> Stability CALL stability2d( Stable, C2D,Eig2d ) IF (Stable == 1) THEN - Write(*,*)"" + WRITE(*,*)"" CALL system ('tput setaf 9;tput bold; echo " =======================================";tput sgr0') CALL system ('tput setaf 9;tput bold; echo " > Elastic Stability Conditions: Unstable; STOP";tput sgr0') CALL system ('tput setaf 9;tput bold; echo " =======================================";tput sgr0') - Write(99,*)" ========================================" - Write(99,*) " > Elastic Stability Conditions: Unstable" - Write(99,*)" ========================================" + WRITE(99,*)" ========================================" + WRITE(99,*) " > Elastic Stability Conditions: Unstable" + WRITE(99,*)" ========================================" STOP END IF IF (Stable == 0) THEN - Write(*,*)"" + WRITE(*,*)"" CALL system ('tput setaf 10;tput bold; echo " =======================================";tput sgr0') CALL system ('tput setaf 10;tput bold; echo " > Elastic Stability Conditions: Stable";tput sgr0') CALL system ('tput setaf 10;tput bold; echo " =======================================";tput sgr0') - Write(99,*)" ========================================" - Write(99,*) " > Elastic Stability Conditions: Stable" - Write(99,*)" ========================================" + WRITE(99,*)" ========================================" + WRITE(99,*) " > Elastic Stability Conditions: Stable" + WRITE(99,*)" ========================================" END IF WRITE(*,'(A,36f10.2)')" > Eigenvalues (N/m):", Eig2d(1,1),Eig2d(2,2),Eig2d(3,3) WRITE(99,'(A,36f10.2)')" > Eigenvalues (N/m):", Eig2d(1,1),Eig2d(2,2),Eig2d(3,3) WRITE(99,*)" " -!< - CALL C_Inv_M2D(3) + !> + CALL C_Inv_M2D(3) CALL proelast_2D () CALL sleep(2) - write(*,*)"" + WRITE(*,*)"" CALL system ('tput setaf 122;tput bold; echo " > Enter phi-mesh. Be divisible by 100 (e.g. 100, 200, 300, etc):";tput sgr0') !WRITE(*,*)" > Enter phi-mesh. Be divisible by 100.(e.g. 100, 200, 300, etc):" read(*,*)npoint END IF - IF (Ncod .eq. 3) then - call system('clear') + IF (Ncod .eq. 3) THEN + CALL system('clear') + IF(input_opt == "N" )THEN + WRITE(*,*)" > Select the type of two-dimensional system:" CALL SYSTEM('tput setaf 41;tput bold; echo " ====================================================";tput sgr0') WRITE(*,*) " Default option (Hex., Squ., and Rec. systems) => 1 " @@ -800,7 +1019,17 @@ PROGRAM AAEP_main WRITE(*,*) " Back ------------------------------------------=> 0 " CALL SYSTEM('tput setaf 41;tput bold; echo " ====================================================";tput sgr0') read(*,*)adv_mubner - + ELSE + IF (input_opt == "def" .or. input_opt == "1") THEN + adv_mubner=1 + !WRITE(*,*) "Output code: Elast" + ELSE IF(input_opt == "adv" .or. input_opt == "2") THEN + adv_mubner=2 + ELSE + WRITE(*,"(2A)") "Invalid input: -op ", adv_mubner + STOP + ENDIF + ENDIF IF (adv_mubner == 2 ) adv = "adv" IF (adv_mubner == 1 ) adv = "ndv" IF (adv_mubner == 0 ) goto 226 @@ -810,46 +1039,49 @@ PROGRAM AAEP_main READ(11,*) C2D(3,1),C2D(3,2),C2D(3,3) close(11) -!> Stability + !> Stability CALL stability2d( Stable, C2D,Eig2d ) IF (Stable == 1) THEN - Write(*,*)"" + WRITE(*,*)"" CALL system ('tput setaf 9;tput bold; echo " ========================================";tput sgr0') CALL system ('tput setaf 9;tput bold; echo " > Elastic Stability Conditions: Unstable; STOP";tput sgr0') CALL system ('tput setaf 9;tput bold; echo " ========================================";tput sgr0') - Write(99,*)" ========================================" - Write(99,*) " > Elastic Stability Conditions: Unstable" - Write(99,*)" ========================================" + WRITE(99,*)" ========================================" + WRITE(99,*) " > Elastic Stability Conditions: Unstable" + WRITE(99,*)" ========================================" STOP END IF IF (Stable == 0) THEN - Write(*,*)"" + WRITE(*,*)"" CALL system ('tput setaf 10;tput bold; echo " ========================================";tput sgr0') CALL system ('tput setaf 10;tput bold; echo " > Elastic Stability Conditions: Stable";tput sgr0') CALL system ('tput setaf 10;tput bold; echo " ========================================";tput sgr0') - Write(99,*)" ========================================" - Write(99,*) " > Elastic Stability Conditions: Stable" - Write(99,*)" ========================================" + WRITE(99,*)" ========================================" + WRITE(99,*) " > Elastic Stability Conditions: Stable" + WRITE(99,*)" ========================================" END IF WRITE(*,'(A,36f10.2)')" > Eigenvalues (N/m):", Eig2d(1,1),Eig2d(2,2),Eig2d(3,3) WRITE(99,'(A,36f10.2)')" > Eigenvalues (N/m):", Eig2d(1,1),Eig2d(2,2),Eig2d(3,3) WRITE(99,*)" " -!< + !< CALL C_Inv_M2D(3) CALL proelast_2D() CALL sleep(1) - write(*,*)"" + WRITE(*,*)"" CALL system ('tput setaf 122;tput bold; echo " > Enter phi-mesh. Be divisible by 100 (e.g. 100, 200, 300, etc):";tput sgr0') !WRITE(*,*)" > Enter phi-mesh. Be divisible by 100.(e.g. 100, 200, 300, etc):" read(*,*)npoint END IF - ELSE - WRITE(*,*)"Invalid Input!"; - STOP + IF(d2d3 /= 3 .and. d2d3 /= 2) THEN + WRITE(*,*)"Invalid Input!!"; + STOP + END IF END IF -END iF !@@@@@@@@@@@@@@@@@@@@@@@ 2D_3D end +END IF !@@@@@@@@@@@@@@@@@@@@@@@ 2D_3D end + + WRITE(*,*)" > Preparing data. Please wait..." -IF(d2d3 == 3) then !@@@@@@@@@@@@@@@@@@@@@@@ 2D_3D system start +IF(d2d3 == 3) THEN !@@@@@@@@@@@@@@@@@@@@@@@ 2D_3D system start !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Data clarity Nmesh_thata = theta_meah!100 Nmesh_phi = phi_meah!100 @@ -857,72 +1089,72 @@ PROGRAM AAEP_main Nmesh_thataf = 150 !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Preparing 3D data... !open(78,file='uuu') - OPEN(111, file="3d_SD.dat") - OPEN(115, file="3d_vec.dat") - OPEN(14,FILE="3d_comp.dat" ) - OPEN(28,FILE="3d_hardness.dat" ) + OPEN(111, FILE="3d_SD.dat" ) + OPEN(115, FILE="3d_vec.dat" ) + OPEN(14, FILE="3d_comp.dat" ) + OPEN(28, FILE="3d_hardness.dat" ) ! OPEN(29,FILE="CompN_Min.dat" ) - OPEN(65,FILE="3d_pugh.dat" ) + OPEN(65, FILE="3d_pugh.dat" ) ! OPEN(66,FILE="pughP_Min.dat" ) ! OPEN(67,FILE="pughN_Min.dat" ) ! OPEN(68,FILE="pughP_Ave.dat" ) - OPEN(15,FILE="3d_young.dat") + OPEN(15, FILE="3d_young.dat") - OPEN(16,FILE="3d_shear.dat" ) + OPEN(16, FILE="3d_shear.dat" ) !OPEN(17,FILE="shearP_Min.dat" ) ! OPEN(8, FILE="shearN_Min.dat" ) ! OPEN(7, FILE="shearP_Ave.dat" ) !OPEN(18,FILE="shearN_Ave.dat" ) - OPEN(20,FILE="3d_poissons.dat" ) + OPEN(20, FILE="3d_poissons.dat" ) ! OPEN(21,FILE="MinP_ratio.dat" ) ! OPEN(22,FILE="MinN_ratio.dat" ) !OPEN(23,FILE="AveP_ratio.dat" ) !OPEN(24,FILE="AveN_ratio.dat" ) - ! OPEN(30,FILE="MinP_bulk.dat" ) - OPEN(31,FILE="3d_bulk.dat" ) + OPEN(30, FILE="3d_bulk_m2.dat" ) + OPEN(31, FILE="3d_bulk.dat" ) - OPEN(32,FILE="3d_sound.dat" ) + OPEN(32, FILE="3d_sound.dat" ) !OPEN(33,FILE="MaxTm_sound.dat" ) !OPEN(34,FILE="MinTm_sound.dat" ) - If (yn_veloc=='Y' .or. yn_veloc=='y') then - OPEN(100,FILE="3d_pp.dat" ) ! 3d_pp.dat > VVP_P.dat - OPEN(101,FILE="3d_gp.dat" ) ! 3d_gp.dat > VVG_P.dat - OPEN(102,FILE="3d_pfp.dat" ) ! 3d_pfp.dat > VV_P_PF.dat - OPEN(103,FILE="3d_gf.dat" ) ! 3d_gf.dat > VVG_Sf.dat - OPEN(104,FILE="3d_pf.dat" ) ! 3d_pf.dat > VVP_Sf.dat - OPEN(105,FILE="3d_pff.dat" ) ! 3d_fpf.dat > VV_Sf_PF.dat - OPEN(106,FILE="3d_gs.dat" ) ! 3d_gs.dat > VVG_Ss.dat - OPEN(107,FILE="3d_ps.dat" ) ! 3d_ps.dat > VVP_Ss.dat - OPEN(108,FILE="3d_pfs.dat" ) ! 3d_spf.dat > VVP_Ss_PF.dat + IF (yn_veloc=='Y' .or. yn_veloc=='y') THEN + OPEN(100, FILE="3d_pp.dat" ) ! 3d_pp.dat > VVP_P.dat + OPEN(101, FILE="3d_gp.dat" ) ! 3d_gp.dat > VVG_P.dat + OPEN(102, FILE="3d_pfp.dat" ) ! 3d_pfp.dat > VV_P_PF.dat + OPEN(103, FILE="3d_gf.dat" ) ! 3d_gf.dat > VVG_Sf.dat + OPEN(104, FILE="3d_pf.dat" ) ! 3d_pf.dat > VVP_Sf.dat + OPEN(105, FILE="3d_pff.dat" ) ! 3d_fpf.dat > VV_Sf_PF.dat + OPEN(106, FILE="3d_gs.dat" ) ! 3d_gs.dat > VVG_Ss.dat + OPEN(107, FILE="3d_ps.dat" ) ! 3d_ps.dat > VVP_Ss.dat + OPEN(108, FILE="3d_pfs.dat" ) ! 3d_spf.dat > VVP_Ss_PF.dat ENDIF - If (yn_km=='Y' .or. yn_km=='y') then + IF (yn_km=='Y' .or. yn_km=='y') THEN OPEN(110,FILE="3d_km.dat" ) ENDIF !IF (MOD(Nmesh_thata,Nmesh_thataf).NE.0) THEN ! Nmesh_thata=(Nmesh_thata/Nmesh_thataf)*Nmesh_thataf - ! write(*,*)Nmesh_thata + ! WRITE(*,*)Nmesh_thata !ENDIF ! !inter_theta=Nmesh_thata/Nmesh_thataf ! IF (MOD(Nmesh_phi,Nmesh_phIF).NE.0) THEN ! Nmesh_phi=(Nmesh_phi/Nmesh_phIF)*Nmesh_phIF - ! write(*,*)Nmesh_phi + ! WRITE(*,*)Nmesh_phi !ENDIF !inter_phi=Nmesh_phi/Nmesh_phIF ! - open(59,file='.aelastpro') + open(59, FILE='.aelastpro') - If (yn_veloc=='Y' .or. yn_veloc=='y') then - open(55,file='.aelastpro2') - endif + IF (yn_veloc=='Y' .or. yn_veloc=='y') THEN + open(55, FILE='.aelastpro2') + ENDIF !! Start loop=0 for wave and elste pro.============================ - call callCij (CCO) + CALL CALLCij (CCO) i=0 j=0 theta=DBLE(i)/DBLE(Nmesh_thata)*PI @@ -963,12 +1195,12 @@ PROGRAM AAEP_main SINver=1.0_dp/SS !Young CALL CShear(G_min,G_max,G_Ave,phi,theta,v11,v12,v13,v22,v23,v33,a6666,sheainvar) CALL CHardness(BINver,SINver,sheainvar,hardvar) - !write(*,*)hardvar + !WRITE(*,*)hardvar Ha_max2 = hardvar Ha_min2 = hardvar - If (yn_veloc=='Y' .or. yn_veloc=='y') then - call wave_main_AAEP(i,j,theta,phi,vec,CCO,density, VVG_P, & + IF (yn_veloc=='Y' .or. yn_veloc=='y') THEN + CALL wave_main_AAEP(i,j,theta,phi,vec,CCO,density, VVG_P, & VVP_P, & VV_P_PF, & VVG_Sf, & @@ -992,35 +1224,35 @@ PROGRAM AAEP_main VV_Sf_PF_min = VV_Sf_PF ! VVG_Sf_max = VVG_Sf - VVP_Sf_max = VVP_Sf - VV_Sf_PF_max = VV_Sf_PF + VVP_Sf_max = VVP_Sf + VV_Sf_PF_max = VV_Sf_PF ! VVG_Ss_min = VVG_Ss VVP_Ss_min = VVP_Ss VV_Ss_PF_min = VV_Ss_PF ! VVG_Ss_max = VVG_Ss - VVP_Ss_max = VVP_Ss + VVP_Ss_max = VVP_Ss VV_Ss_PF_max = VV_Ss_PF -endif - If (yn_km=='Y' .or. yn_km=='y') then - Call Ckm_cal(density, SINver, ma_avrag, km ) +ENDIF + IF (yn_km == 'Y' .or. yn_km == 'y') THEN + CALL Ckm_cal(density, SINver, ma_avrag, km ) km_max = km km_min = km - endif + ENDIF !! Finish loop=0 for wave and elsat pro.============================ - OPEN(10,file="Cij.dat") + OPEN(10, FILE="Cij.dat") DO i=0,Nmesh_thata theta=DBLE(i)/DBLE(Nmesh_thata)*PI - Nmesh_phi2=Nmesh_phi + Nmesh_phi2 = Nmesh_phi IF (i.EQ.0 .OR. i.EQ.Nmesh_thata) Nmesh_phi2=0D0 - if ( i==Nmesh_thata/4d0 ) then ; CALL system( ' clear ' ) ;PRINT*, '|====%25====>-------------|'; ENDIF - if ( i==Nmesh_thata/3d0 ) then ; CALL system( ' clear ' ) ;PRINT*, '|======%50======>---------|'; ENDIF - if ( i==Nmesh_thata/2d0 ) then ; CALL system( ' clear ' ) ;PRINT*, '|========%75========>-----|'; ENDIF - if ( i==Nmesh_thata/1d0 ) then ; CALL system( ' clear ' ) ;PRINT*, '|==========%100==========>|'; ENDIF + IF ( i==Nmesh_thata/4d0 ) THEN ; CALL system( ' clear ' ) ;PRINT*, '|====%25====>-------------|'; ENDIF + IF ( i==Nmesh_thata/3d0 ) THEN ; CALL system( ' clear ' ) ;PRINT*, '|======%50======>---------|'; ENDIF + IF ( i==Nmesh_thata/2d0 ) THEN ; CALL system( ' clear ' ) ;PRINT*, '|========%75========>-----|'; ENDIF + IF ( i==Nmesh_thata/1d0 ) THEN ; CALL system( ' clear ' ) ;PRINT*, '|==========%100==========>|'; ENDIF DO j=0, Nmesh_phi2 phi=DBLE(j)/DBLE(Nmesh_phi)*2D0*PI @@ -1031,8 +1263,8 @@ PROGRAM AAEP_main v11=vec(1)*vec(1) ; v12=vec(1)*vec(2) v13=vec(1)*vec(3) ; v22=vec(2)*vec(2) v23=vec(2)*vec(3) ; v33=vec(3)*vec(3) - If (yn_veloc=='Y' .or. yn_veloc=='y') then - call callCij (CCO) + IF (yn_veloc == 'Y' .OR. yn_veloc == 'y') THEN + CALL CALLCij (CCO) CALL wave_main_AAEP( i,j,theta,phi,vec,CCO,density, VVG_P, & VVP_P, & @@ -1045,25 +1277,25 @@ PROGRAM AAEP_main VV_Ss_PF ) - ! write(*,*)VVP_P,VVP_Ss,VVP_Sf - IF (VVG_P .GE.VVG_P_max) THEN + ! WRITE(*,*)VVP_P,VVP_Ss,VVP_Sf + IF (VVG_P .GE. VVG_P_max) THEN VVG_P_max = VVG_P VVG_P_max_phi = (phi *180.0D0)/PI VVG_P_max_theta = (theta*180.0D0)/PI ENDIF - IF (VVG_P .LE.VVG_P_min) THEN + IF (VVG_P .LE. VVG_P_min) THEN VVG_P_min = VVG_P VVG_P_min_phi = (phi *180.0D0)/PI VVG_P_min_theta = (theta*180.0D0)/PI ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!&&&&&&&&&&&&&&&&& - IF (VVP_P .GE.VVP_P_max) THEN + IF (VVP_P .GE. VVP_P_max) THEN VVP_P_max = VVP_P VVP_P_max_phi = (phi *180.0D0)/PI VVP_P_max_theta = (theta*180.0D0)/PI ENDIF - IF (VVP_P .LE.VVP_P_min) THEN + IF (VVP_P .LE. VVP_P_min) THEN VVP_P_min = VVP_P VVP_P_min_phi = (phi *180.0D0)/PI VVP_P_min_theta = (theta*180.0D0)/PI @@ -1167,6 +1399,7 @@ PROGRAM AAEP_main ! ! Triclinic !IF(CLaS .EQ. 7) THEN +!==================================================== bulk method 1 v1.7.2 BB=( S(1,1)+S(1,2)+S(1,3) )*v11& +( S(1,6)+S(2,6)+S(3,6) )*v12& +( S(1,5)+S(2,5)+S(3,5) )*v13& @@ -1187,6 +1420,30 @@ PROGRAM AAEP_main Minbulk_phi = phi !WRITE(*,*)Minyoung ENDIF +!==================================================== bulk method 1 v1.7.2 + +!==================================================== bulk method 3 v1.7.3 + CALL bulk_method(Pratio, sheainvar*1000D0, bulk_m2) + + ! Max_bulkm2 = bulk_m2 + Min_bulkm2 = bulk_m2 + + !WRITE(*,*) bulk_m2 + IF (bulk_m2 .GE. Max_bulkm2)THEN + Max_bulkm2 = bulk_m2 + Maxbulkm2_theta = theta + Maxbulkm2_phi = phi + ! WRITE(*,*) Max_bulkm2 + ENDIF + + IF (bulk_m2 .LE. Min_bulkm2)THEN + Min_bulkm2 = bulk_m2 + Minbulkm2_theta = theta + Minbulkm2_phi = phi + !WRITE(*,*) Min_bulkm2 + ENDIF + +!==================================================== bulk method 3 v1.7.3 SS=v11*v11*S(1,1)& +2*v12*v12*S(1,2)& +2*v13*v13*S(1,3)& @@ -1220,9 +1477,7 @@ PROGRAM AAEP_main !WRITE(*,*)Minyoung ENDIF !----------------------------------------------- - If (yn_km=='Y' .OR. yn_km=='y') THEN - - + IF (yn_km=='Y' .OR. yn_km=='y') THEN CALL Ckm_cal(density, SINver, ma_avrag, km ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!&&&&&&&&&&&&&&&& > v1.7.0 @@ -1291,13 +1546,13 @@ PROGRAM AAEP_main ENDIF IF ((CO*1000D0) .LE. Mincomp) THEN - Mincomp =CO*1000D0 - Mincomp_theta =theta - Mincomp_phi =phi + Mincomp = CO*1000D0 + Mincomp_theta = theta + Mincomp_phi = phi ENDIF comminn=0.0d0 comminp=0.0d0 - ! write(41,*) Mincomp + ! WRITE(41,*) Mincomp IF (co .LE. 0.0D0) comminn = -1.0d0*CO IF (co .GE. 0.0D0) comminp = CO @@ -1351,22 +1606,28 @@ PROGRAM AAEP_main Ha_min2_theta = theta Ha_min2_phi = phi ENDIF - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - +! !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,hardvar - If (yn_veloc=='Y' .or. yn_veloc=='y') then + + 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,km - endif + ENDIF WRITE(111,*) theta, phi - write(115,"(3F30.15)")vec(1),vec(2),vec(3) + WRITE(115,"(3F30.15)")vec(1),vec(2),vec(3) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! tests +WRITE(30,"(9F30.15)") vec(1)*ABS( bulk_m2 ), vec(2)*ABS( bulk_m2 ), vec(3)*ABS( bulk_m2 ) ,& + vec(1)*ABS( Min_bulkm2 ), vec(2)*ABS( Min_bulkm2 ), vec(3)*ABS( Min_bulkm2 ),& + vec(1)*ABS( Max_bulkm2 ), vec(2)*ABS( Max_bulkm2 ), vec(3)*ABS( Max_bulkm2 ) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + 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(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 vec(1)*ABS(shminn*1000D0 ), vec(2)*ABS( shminn*1000D0), vec(3)*ABS(shminn*1000D0 ),& !>> shear @@ -1398,7 +1659,7 @@ PROGRAM AAEP_main vec(1)*ABS( maxEVaTM ), vec(2)*ABS( maxEVaTM ), vec(3)*ABS( maxEVaTM ),& !>> sound vec(1)*ABS( minEVaTM ), vec(2)*ABS( minEVaTM ), vec(3)*ABS( minEVaTM ) !>> sound - If (yn_veloc=='Y' .or. yn_veloc=='y') then + IF(yn_veloc=='Y' .or. yn_veloc=='y') THEN WRITE(100,*) vec(1)*ABS( VVP_P ), vec(2)*ABS( VVP_P ) ,vec(3)*ABS( VVP_P ) !>> VVP_P WRITE(101,*) vec(1)*ABS( VVG_P ), vec(2)*ABS( VVG_P ) ,vec(3)*ABS( VVG_P ) !>> VVG_P WRITE(102,*) vec(1)*ABS( VV_P_PF ), vec(2)*ABS( VV_P_PF ) ,vec(3)*ABS( VV_P_PF ) !>> VV_P_PF @@ -1408,57 +1669,57 @@ PROGRAM AAEP_main 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 - end if + ENDIF - If (yn_km=='Y' .or. yn_km=='y') then - WRITE(110,*) vec(1)*ABS( km ), vec(2)*ABS( km ) ,vec(3)*ABS( km ) !>> k_m + IF(yn_km=='Y' .or. yn_km=='y') THEN + WRITE(110,*) vec(1)*ABS( km ), vec(2)*ABS( km ) ,vec(3)*ABS( km ) !>> k_m ENDIF ENDDO WRITE(16,*) " " - ! WRITE(111,*) " " - !!WRITE( 8,*) " " + !WRITE(111,*) " " + !WRITE( 8,*) " " !WRITE( 7,*) " " WRITE(65,*) " " !WRITE(66,*) " " - ! WRITE(67,*) " " - ! WRITE(68,*) " " - !WRITE(18,*) " " + !WRITE(67,*) " " + !WRITE(68,*) " " + !WRITE(18,*) " " WRITE(15,*) " " WRITE(14,*) " " - WRITE(28,*) " " + WRITE(28,*) " " !WRITE(29,*) " " WRITE(20,*) " " !WRITE(21,*) " " !WRITE(22,*) " " !WRITE(23,*) " " !WRITE(24,*) " " - !WRITE(30,*) " " + !WRITE(30,*) " " WRITE(31,*) " " WRITE(32,*) " " !WRITE(33,*) " " !WRITE(34,*) " " - If (yn_veloc=='Y' .or. yn_veloc=='y') then - WRITE(100,*) " " - endif + IF (yn_veloc=='Y' .or. yn_veloc=='y') THEN + WRITE(100,*) " " + ENDIF ENDDO CALL system('sleep 0.2; clear') - CLOSE(15) - CLOSE(16) - CLOSE(111) - CLOSE(18) - CLOSE(20) + CLOSE (15) + CLOSE (16) + CLOSE (111) + CLOSE (18) + CLOSE (20) !CLOSE(21) !CLOSE(22) !CLOSE(23) !CLOSE(24) !CLOSE(7 ) !CLOSE(8 ) - !CLOSE(30) - CLOSE(31) - CLOSE(32) + CLOSE (30) + CLOSE (31) + CLOSE (32) !CLOSE(33) CLOSE (100) @@ -1470,26 +1731,26 @@ PROGRAM AAEP_main CLOSE (106) CLOSE (107) CLOSE (108) - CLOSE(110) - close(41) + CLOSE (110) + CLOSE (41) - close(59) + CLOSE (59) !CALL("touch .MaMiout2") - ! 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,& + ! 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,km_max,km_min - close(37) + CLOSE (37) ! ENDIF - open(36,file='.MaMiout') - write(36,*)Maxyoung,Minyoung,Maxcomp,Mincomp,G_max2*1000D0,G_min2*1000D0,Maxbulk,Minbulk,& + 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, Ha_max2,Ha_min2 - close(36) + CLOSE(36) !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Preparing features... WRITE(*,*)'' - call anisotropy(Maxyoung,Minyoung, Ax(1)) + CALL anisotropy(Maxyoung,Minyoung, Ax(1)) CALL SYSTEM('tput setaf 196;tput bold; echo " ==================================================> Youngs Modulus"') CALL SYSTEM('tput setaf 196;tput bold; echo " Max(GPa) Min(GPa) Anisotropy"') WRITE (*,'(3xF6.2,a,2xF6.2,a,2xF6.2)') Maxyoung," ",Minyoung," ",Ax(1) @@ -1505,7 +1766,7 @@ PROGRAM AAEP_main CALL SYSTEM('sleep 0.5') !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(*,*)'' - call anisotropy(Maxcomp,Mincomp, Ax(1)) + CALL anisotropy(Maxcomp,Mincomp, Ax(1)) CALL SYSTEM('tput setaf 202;tput bold; echo " ==================================================> Linear Compressibiliy";tput sgr0') CALL SYSTEM('tput setaf 202;tput bold; echo " Max(TPa-1) Min(TPa-1) Anisotropy"') WRITE (*,'(1xF9.3,a,F9.3,a,2xF6.3)') Maxcomp," ",Mincomp," ",Ax(1) @@ -1521,7 +1782,7 @@ PROGRAM AAEP_main CALL SYSTEM('sleep 0.5') !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(*,*)'' - call anisotropy(G_max2,G_min2, Ax(1)) + CALL anisotropy(G_max2,G_min2, Ax(1)) CALL SYSTEM('tput setaf 161;tput bold; echo " ==================================================> Shear Modulus";tput sgr0') CALL SYSTEM('tput setaf 161;tput bold; echo " Max(GPa) Min(GPa) Anisotropy"') WRITE (*,'(3xF6.2,a,2xF6.2,a,2xF6.2)') G_max2*1000D0," ",G_min2*1000D0," ",Ax(1) @@ -1537,8 +1798,8 @@ PROGRAM AAEP_main CALL SYSTEM('sleep 0.5') !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(*,*)'' - call anisotropy(Maxbulk,Minbulk, Ax(1)) - CALL SYSTEM('tput setaf 126;tput bold; echo " ==================================================> Bulk Modulus";tput sgr0') + CALL anisotropy(Maxbulk,Minbulk, Ax(1)) + CALL SYSTEM('tput setaf 126;tput bold; echo " ==================================================> Bulk Modulus method 1 (Experimental)";tput sgr0') CALL SYSTEM('tput setaf 126;tput bold; echo " Max(GPa/100) Min(GPa/100) Anisotropy"') WRITE (*,'(2xF11.2,a,1xF11.2,a,2xF6.2)') Maxbulk/100.d0," ",Minbulk/100.d0," ",Ax(1) CALL angl2cart(Maxbulk_theta,Maxbulk_phi, vec(1),vec(2),vec(3)) @@ -1552,9 +1813,25 @@ PROGRAM AAEP_main WRITE (*,*) "==================================================" CALL SYSTEM('sleep 0.5') !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - If (yn_veloc=='Y' .or. yn_veloc=='y') then + !WRITE(*,*)'' + ! CALL anisotropy(Max_bulkm2,Min_bulkm2, Ax(1)) +! CALL SYSTEM('tput setaf 126;tput bold; echo " ==================================================> Bulk Modulus method 2 (Experimental)";tput sgr0') +! CALL SYSTEM('tput setaf 126;tput bold; echo " Max(GPa/100) Min(GPa/100) Anisotropy"') +! WRITE (*,'(2xF11.2,a,1xF11.2,a,2xF6.2)') Max_bulkm2/100.d0," ",Min_bulkm2/100.d0," ",Ax(1) +! CALL angl2cart(Maxbulkm2_theta,Maxbulkm2_phi, vec(1), vec(2),vec(3)) +! CALL angl2cart(Minbulkm2_theta,Minbulkm2_phi, vec1(1),vec1(2),vec1(3)) +! WRITE(*,*)"------------------------------------------" +! WRITE (*,*) " Theta Phi"," Theta Phi" +! WRITE (*,'(1x2F6.1,11x2F6.1)') (Maxbulkm2_theta*180.0D0)/PI,(Maxbulkm2_phi*180.0D0)/PI,(Minbulkm2_theta*180.0D0)/PI,(Minbulkm2_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') + !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + IF (yn_veloc=='Y' .or. yn_veloc=='y') THEN WRITE(*,*)'' - call anisotropy(VVP_P_max,VVP_P_min, Ax(1)) + CALL anisotropy(VVP_P_max,VVP_P_min, Ax(1)) CALL SYSTEM('tput setaf 21;tput bold; echo " ==================================================> Phase P-Mode";tput sgr0') CALL SYSTEM('tput setaf 21;tput bold; echo " Max(km/a) Min(km/s) Anisotropy"') WRITE (*,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVP_P_max/1000D0," ",VVP_P_min/1000D0," ",Ax(1) @@ -1571,7 +1848,7 @@ PROGRAM AAEP_main CALL SYSTEM('sleep 0.5') !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(*,*)'' - call anisotropy(VVP_Sf_max,VVP_Sf_min, Ax(1)) + CALL anisotropy(VVP_Sf_max,VVP_Sf_min, Ax(1)) CALL SYSTEM('tput setaf 33;tput bold; echo " ==================================================> Phase Fast-Mode";tput sgr0') CALL SYSTEM('tput setaf 33;tput bold; echo " Max(km/a) Min(km/s) Anisotropy"') WRITE (*,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVP_Sf_max/1000D0," ",VVP_Sf_min/1000D0," ",Ax(1) @@ -1588,7 +1865,7 @@ PROGRAM AAEP_main CALL SYSTEM('sleep 0.5') !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(*,*)'' - call anisotropy(VVP_Ss_max,VVP_Ss_min, Ax(1)) + CALL anisotropy(VVP_Ss_max,VVP_Ss_min, Ax(1)) CALL SYSTEM('tput setaf 34;tput bold; echo " ==================================================> Phase Slow-Mode";tput sgr0') CALL SYSTEM('tput setaf 34;tput bold; echo " Max(km/a) Min(km/s) Anisotropy"') WRITE (*,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVP_Ss_max/1000D0," ",VVP_Ss_min/1000D0," ",Ax(1) @@ -1605,7 +1882,7 @@ PROGRAM AAEP_main CALL SYSTEM('sleep 0.5') !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(*,*)'' - call anisotropy(VVG_P_max,VVG_P_min, Ax(1)) + CALL anisotropy(VVG_P_max,VVG_P_min, Ax(1)) CALL SYSTEM('tput setaf 40;tput bold; echo " ==================================================> Group P-Mode";tput sgr0') CALL SYSTEM('tput setaf 40;tput bold; echo " Max(km/a) Min(km/s) Anisotropy"') WRITE (*,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVG_P_max/1000D0," ",VVG_P_min/1000D0," ",Ax(1) @@ -1622,7 +1899,7 @@ PROGRAM AAEP_main CALL SYSTEM('sleep 0.5') !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(*,*)'' - call anisotropy(VVG_Sf_max,VVG_Sf_min, Ax(1)) + CALL anisotropy(VVG_Sf_max,VVG_Sf_min, Ax(1)) CALL SYSTEM('tput setaf 226;tput bold; echo " ==================================================> Group Fast-Mode";tput sgr0') CALL SYSTEM('tput setaf 226;tput bold; echo " Max(km/a) Min(km/s) Anisotropy"') WRITE (*,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVG_Sf_max/1000D0," ",VVG_Sf_min/1000D0," ",Ax(1) @@ -1639,7 +1916,7 @@ PROGRAM AAEP_main CALL SYSTEM('sleep 0.5') !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(*,*)'' - call anisotropy(VVG_Ss_max,VVG_Ss_min, Ax(1)) + CALL anisotropy(VVG_Ss_max,VVG_Ss_min, Ax(1)) CALL SYSTEM('tput setaf 22;tput bold; echo " ==================================================> Group Slow-Mode";tput sgr0') CALL SYSTEM('tput setaf 22;tput bold; echo " Max(km/a) Min(km/s) Anisotropy"') WRITE (*,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVG_Ss_max/1000D0," ",VVG_Ss_min/1000D0," ",Ax(1) @@ -1655,9 +1932,9 @@ PROGRAM AAEP_main !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ CALL SYSTEM('sleep 0.5') !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - If (yn_km=='Y' .or. yn_km=='y') then + IF (yn_km=='Y' .or. yn_km=='y') THEN WRITE(*,*)'' - Call anisotropy(km_max,km_min, Ax(1)) + CALL anisotropy(km_max,km_min, Ax(1)) CALL SYSTEM('tput setaf 12;tput bold; echo " ==================================================> Min. thermal conductivity";tput sgr0') CALL SYSTEM('tput setaf 12;tput bold; echo " Max(W/K.m) Min(W/K.m) Anisotropy"') WRITE (*,'(2xF11.3,a,1xF11.3,a,2xF6.2)') km_max," ",km_min ," ",Ax(1) @@ -1673,10 +1950,10 @@ PROGRAM AAEP_main ENDIF !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CALL SYSTEM('sleep 0.5') - END if + END IF !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(*,*)'' - call anisotropy(Pratio_max,Pratio_min, Ax(1)) + CALL anisotropy(Pratio_max,Pratio_min, Ax(1)) CALL SYSTEM('tput setaf 201;tput bold; echo " ==================================================> Poissons Ratio";tput sgr0') CALL SYSTEM('tput setaf 201;tput bold; echo " Max Min Anisotropy"') WRITE (*,'(2xF6.3,a,2xF6.3,a,2xF6.2)') Pratio_max," ",Pratio_min," ",Ax(1) @@ -1691,7 +1968,7 @@ PROGRAM AAEP_main WRITE (*,*) "==================================================" CALL SYSTEM('sleep 0.5') WRITE(*,*)'' - call anisotropy(pugh_max2,pugh_min2, Ax(1)) + CALL anisotropy(pugh_max2,pugh_min2, Ax(1)) !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ CALL SYSTEM('tput setaf 91;tput bold; echo " ==================================================> Pugh Ratio";tput sgr0') CALL SYSTEM('tput setaf 91;tput bold; echo " Max Min Anisotropy"') @@ -1709,8 +1986,8 @@ PROGRAM AAEP_main WRITE(*,*)'' !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(*,*)'' - call anisotropy(Maxbulk,Minbulk, Ax(1)) - CALL SYSTEM('tput setaf 55;tput bold; echo " ==================================================> Hardness";tput sgr0') + CALL anisotropy(Maxbulk,Minbulk, Ax(1)) + CALL SYSTEM('tput setaf 55;tput bold; echo " ==================================================> Hardness (Experimental)";tput sgr0') CALL SYSTEM('tput setaf 55;tput bold; echo " Max Min Anisotropy"') WRITE (*,'(2xF10.3,a,1xF10.3,a,2xF6.2)') Ha_max2," ",Ha_min2," ",Ax(1) CALL angl2cart(Ha_max2_theta,Ha_max2_phi, vec(1),vec(2),vec(3)) @@ -1735,9 +2012,9 @@ PROGRAM AAEP_main !*************************************************************** 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" - IF (Ncod .EQ. 6) then + IF (Ncod .EQ. 6) THEN WRITE (99,'(2a)') " > MP ID:", myid - endif + ENDIF !WRITE (99,*) e1,e2,"plane" WRITE (99,*) '' CALL anisotropy(Maxyoung,Minyoung, Ax(1)) @@ -1756,7 +2033,7 @@ PROGRAM AAEP_main !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(99,*)' ' - call anisotropy(Maxcomp,Mincomp, Ax(1)) + CALL anisotropy(Maxcomp,Mincomp, Ax(1)) WRITE(99,*) " ==================================================> Linear Compressibiliy" WRITE(99,*) " Max(TPa-1) Min(TPa-1) Anisotropy" WRITE (99,'(1xF9.3,a,F9.3,a,2xF6.3)') Maxcomp," ",Mincomp," ",Ax(1) @@ -1772,7 +2049,7 @@ PROGRAM AAEP_main !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(99,*)' ' - call anisotropy(G_max2,G_min2, Ax(1)) + CALL anisotropy(G_max2,G_min2, Ax(1)) WRITE(99,*) " ==================================================> Shear Modulus" WRITE(99,*) " Max(GPa) Min(GPa) Anisotropy" WRITE (99,'(3xF6.2,a,2xF6.2,a,2xF6.2)') G_max2*1000D0," ",G_min2*1000D0," ",Ax(1) @@ -1788,8 +2065,8 @@ PROGRAM AAEP_main !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(99,*)' ' - call anisotropy(Maxbulk,Minbulk, Ax(1)) - WRITE(99,*) " ==================================================> Bulk Modulus" + CALL anisotropy(Maxbulk,Minbulk, Ax(1)) + WRITE(99,*) " ==================================================> Bulk Modulus (Experimental)" WRITE(99,*) " Max(GPa/100) Min(GPa/100) Anisotropy" WRITE (99,'(2xF11.2,a,1xF11.2,a,2xF6.2)') Maxbulk/100.d0," ",Minbulk/100.d0," ",Ax(1) CALL angl2cart(Maxbulk_theta,Maxbulk_phi, vec(1),vec(2),vec(3)) @@ -1803,9 +2080,9 @@ PROGRAM AAEP_main WRITE(99,*) "==================================================" !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ - If (yn_veloc=='Y' .or. yn_veloc=='y') THEN + IF (yn_veloc=='Y' .or. yn_veloc=='y') THEN WRITE(99,*)' ' - call anisotropy(VVP_P_max,VVP_P_min, Ax(1)) + CALL anisotropy(VVP_P_max,VVP_P_min, Ax(1)) WRITE(99,*) " ==================================================> Phase P-Mode" WRITE(99,*) " Max(km/s) Min(km/s) Anisotropy" WRITE (99,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVP_P_max/1000D0," ",VVP_P_min/1000D0," ",Ax(1) @@ -1821,7 +2098,7 @@ PROGRAM AAEP_main !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(99,*)' ' - call anisotropy(VVP_Sf_max,VVP_Sf_min, Ax(1)) + CALL anisotropy(VVP_Sf_max,VVP_Sf_min, Ax(1)) WRITE(99,*) " ==================================================> Phase Fast-Mode" WRITE(99,*) " Max(km/s) Min(km/s) Anisotropy" WRITE (99,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVP_Sf_max/1000D0," ",VVP_Sf_min/1000D0," ",Ax(1) @@ -1837,7 +2114,7 @@ PROGRAM AAEP_main !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(99,*)' ' - call anisotropy(VVP_Ss_max,VVP_Ss_min, Ax(1)) + CALL anisotropy(VVP_Ss_max,VVP_Ss_min, Ax(1)) WRITE(99,*) " ==================================================> Phase Slow-Mode" WRITE(99,*) " Max(km/s) Min(km/s) Anisotropy" WRITE (99,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVP_Ss_max/1000D0," ",VVP_Ss_min/1000D0," ",Ax(1) @@ -1853,7 +2130,7 @@ PROGRAM AAEP_main !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(99,*)' ' - call anisotropy(VVG_P_max,VVG_P_min, Ax(1)) + CALL anisotropy(VVG_P_max,VVG_P_min, Ax(1)) WRITE(99,*) " ==================================================> Group P-Mode" WRITE(99,*) " Max(km/s) Min(km/s) Anisotropy" WRITE (99,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVG_P_max/1000D0," ",VVG_P_min/1000D0," ",Ax(1) @@ -1869,7 +2146,7 @@ PROGRAM AAEP_main !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(99,*)' ' - call anisotropy(VVG_Sf_max,VVG_Sf_min, Ax(1)) + CALL anisotropy(VVG_Sf_max,VVG_Sf_min, Ax(1)) WRITE(99,*) " ==================================================> Group Fast-Mode" WRITE(99,*) " Max(km/s) Min(km/s) Anisotropy" WRITE (99,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVG_Sf_max/1000D0," ",VVG_Sf_min/1000D0," ",Ax(1) @@ -1885,7 +2162,7 @@ PROGRAM AAEP_main !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ WRITE(99,*)' ' - call anisotropy(VVG_Ss_max,VVG_Ss_min, Ax(1)) + CALL anisotropy(VVG_Ss_max,VVG_Ss_min, Ax(1)) WRITE(99,*) " ==================================================> Group Slow-Mode" WRITE(99,*) " Max(km/s) Min(km/s) Anisotropy" WRITE (99,'(2xF11.3,a,1xF11.3,a,2xF6.2)') VVG_Ss_max/1000D0," ",VVG_Ss_min/1000D0," ",Ax(1) @@ -1901,9 +2178,9 @@ PROGRAM AAEP_main ENDIF !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - If (yn_km=='Y' .or. yn_km=='y') then + IF (yn_km=='Y' .or. yn_km=='y') THEN WRITE(*,*)'' - Call anisotropy(km_max,km_min, Ax(1)) + CALL anisotropy(km_max,km_min, Ax(1)) WRITE(99,*) " ==================================================> Min. thermal conductivity" WRITE(99,*) " Max(W/K.m) Min(W/K.m) Anisotropy" WRITE (99,'(2xF11.3,a,1xF11.3,a,2xF6.2)') km_max," ",km_min ," ",Ax(1) @@ -1915,12 +2192,11 @@ PROGRAM AAEP_main WRITE(99,*)"------------------------------------------" WRITE (99,*) " x y z"," x y z" WRITE (99,'(1x3F6.2,3x3F6.2)') vec(1),vec(2),vec(3),vec1(1),vec1(2),vec1(3) - ENDIF !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% WRITE(99,*)'' - call anisotropy(Pratio_max,Pratio_min, Ax(1)) + CALL anisotropy(Pratio_max,Pratio_min, Ax(1)) WRITE(99,*) " ==================================================> Poisson's Ratio" WRITE(99,*) " Max Min Anisotropy" WRITE (99,'(2xF6.3,a,2xF6.3,a,2xF6.2)') Pratio_max," ",Pratio_min," ",Ax(1) @@ -1934,8 +2210,8 @@ PROGRAM AAEP_main WRITE(99,'(1x3F6.2,3x3F6.2)') vec(1),vec(2),vec(3),vec1(1),vec1(2),vec1(3) WRITE(99,*) "==================================================" WRITE(99,*)'' - call anisotropy(Maxbulk,Minbulk, Ax(1)) - WRITE(99,*) " ==================================================> Hardness" + CALL anisotropy(Maxbulk,Minbulk, Ax(1)) + WRITE(99,*) " ==================================================> Hardness (Experimental)" WRITE(99,*) " Max Min Anisotropy" WRITE (99,'(2xF10.3,a,1xF10.3,a,2xF6.2)') Ha_max2," ",Ha_min2," ",Ax(1) CALL angl2cart(Ha_max2_theta,Ha_max2_phi, vec(1),vec(2),vec(3)) @@ -1992,9 +2268,9 @@ PROGRAM AAEP_main SINver=1D0/SS IF (SINver.NE.0) THEN young2dmax(num)=SINver - ELSE + ELSE young2dmax(num)=0D0 - ENDIF + ENDIF !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! @@ -2006,23 +2282,23 @@ PROGRAM AAEP_main +( S(1,3)+S(2,3)+S(3,3) )*vv33 ! IF (CO*1000D0 .GE.0.0d0) THEN - !Maxcomp= + !Maxcomp= comMAX2d(num)=CO*1000D0 !Maxcomp ENDIF - !write(*,*)CO*1000D0 + !WRITE(*,*)CO*1000D0 IF (CO *1000D0 .LE. 0.0d0) THEN !Mincomp= comMIN2d(num)=CO*1000D0 !Mincomp ENDIF !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - BB= ( S(1,1)+S(1,2)+S(1,3) )*vv11 & - +( S(1,6)+S(2,6)+S(3,6) )*vv12 & - +( S(1,5)+S(2,5)+S(3,5) )*vv13 & - +( S(1,2)+S(2,2)+S(2,3) )*vv22 & - +( S(1,4)+S(2,4)+S(3,4) )*vv23 & - +( S(1,3)+S(2,3)+S(3,3) )*vv33 - BINver=1D0/BB + BB = ( S(1,1)+S(1,2)+S(1,3) )*vv11 & + +( S(1,6)+S(2,6)+S(3,6) )*vv12 & + +( S(1,5)+S(2,5)+S(3,5) )*vv13 & + +( S(1,2)+S(2,2)+S(2,3) )*vv22 & + +( S(1,4)+S(2,4)+S(3,4) )*vv23 & + +( S(1,3)+S(2,3)+S(3,3) )*vv33 + BINver = 1D0/BB bulkMax2d(num) = BINver ! IF (BINver.GE.0D0) THEN ! Maxbulk = BINver @@ -2030,13 +2306,13 @@ PROGRAM AAEP_main ! ENDIF - IF (BINver.LE.0D0) THEN - Minbulk = BINver - bulkMin2d(num) = Maxbulk + IF(BINver.LE.0D0) THEN + Minbulk = BINver + bulkMin2d(num) = Maxbulk ENDIF - ave=0D0 - v=0 + ave = 0D0 + v = 0 Nmesh_gamma=180 DO kk=0,Nmesh_gamma-1 gamma=DBLE(kk)/DBLE(Nmesh_gamma)*pi @@ -2071,18 +2347,14 @@ PROGRAM AAEP_main sheainvar=1D0/((4D0*a6666)) IF (kk.EQ.0.D0) THEN G_max = sheainvar; G_min = sheainvar - ELSE - IF (sheainvar.GE.G_max) G_max = sheainvar IF (sheainvar.LE.G_min) G_min = sheainvar ENDIF ave = ave+sheainvar ENDDO G_ave = ave/Nmesh_gamma - shear2dmax(num) = G_max - IF (G_max.GE.0D0) shear2dminp(num) = G_min IF (G_max.LE.0D0) shear2dminn(num) = G_min !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -2108,26 +2380,26 @@ PROGRAM AAEP_main hard2dmax(num) = 0D0 ENDIF !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - IF (yn_km=='Y' .OR. yn_km=='y') THEN - Call Ckm_cal(density, SINver, ma_avrag, km ) + IF (yn_km=='Y' .OR. yn_km=='y') THEN + CALL Ckm_cal(density, SINver, ma_avrag, km ) IF (km.NE.0) THEN km2dmax(num)=km ELSE km2dmax(num)=0D0 ENDIF - ENDIF + ENDIF !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ num=num+1 ENDDO !!!!!!!!!!!!!!#######################################!!!!!!!!!!!!!!!!!!!! - If (yn_veloc=='Y' .OR. yn_veloc=='y') THEN + IF (yn_veloc=='Y' .OR. yn_veloc=='y') THEN vec = 0D0 theta = 0.0d0 DO i=0,180 phi=0.0d0 DO j=0,360 - CALL callCij (CCO) + CALL CALLCij (CCO) CALL twoD_calc_wave(mmx,kky,llz,smkl,j,vec) CALL wave_main_AAEP( i,j,theta,phi,vec,CCo,density, VVG_P, & @@ -2139,18 +2411,18 @@ PROGRAM AAEP_main VVG_Ss, & VVP_Ss, & VV_Ss_PF ) - if (i==0)VVG_P_2D(j) = VVG_P /1000D0 - if (i==0)VVP_P_2D(j) = VVP_P /1000D0 - if (i==0)VV_P_PF_2D(j) = VV_P_PF !/1000D0 - if (i==0)VVG_Sf_2D(j) = VVG_Sf /1000D0 - if (i==0)VVP_Sf_2D(j) = VVP_Sf /1000D0 - if (i==0)VV_Sf_PF_2D(j)= VV_Sf_PF!/1000D0 - if (i==0)VVG_Ss_2D(j) = VVG_Ss /1000D0 - if (i==0)VVP_Ss_2D(j) = VVP_Ss /1000D0 - if (i==0)VV_Ss_PF_2D(j)= VV_Ss_PF!/1000D0 + IF (i==0)VVG_P_2D(j) = VVG_P /1000D0 + IF (i==0)VVP_P_2D(j) = VVP_P /1000D0 + IF (i==0)VV_P_PF_2D(j) = VV_P_PF !/1000D0 + IF (i==0)VVG_Sf_2D(j) = VVG_Sf /1000D0 + IF (i==0)VVP_Sf_2D(j) = VVP_Sf /1000D0 + IF (i==0)VV_Sf_PF_2D(j)= VV_Sf_PF!/1000D0 + IF (i==0)VVG_Ss_2D(j) = VVG_Ss /1000D0 + IF (i==0)VVP_Ss_2D(j) = VVP_Ss /1000D0 + IF (i==0)VV_Ss_PF_2D(j)= VV_Ss_PF!/1000D0 - ! WRITE(*,*)VV_Sf_PF + ! WRITE(*,*)VV_Sf_PF phi=phi+ 0.5d0 ENDDO theta=theta+0.5d0 @@ -2178,62 +2450,61 @@ PROGRAM AAEP_main ! OPEN( 73,FILE="2dpughMinn.dat" ) - If (yn_veloc=='Y' .or. yn_veloc=='y') then !----------------------------------------------------- - + IF (yn_veloc=='Y' .or. yn_veloc=='y') THEN !----------------------------------------------------- OPEN(500,FILE='2dcut_gveloc.dat') OPEN(501,FILE='2dcut_pveloc.dat') OPEN(502,FILE='2dcut_pfaveloc.dat') ENDIF !----------------------------------------------------- - If (yn_km=='Y' .or. yn_km=='y') then !!----------------------------------------------------- + IF (yn_km=='Y' .or. yn_km=='y') THEN !!----------------------------------------------------- OPEN(499,FILE='2dcut_km.dat') ENDIF !!----------------------------------------------------- - num=0 - test_n=0.0 - DO num=0,Nmesh_thata + num = 0 + test_n = 0.0 + DO num = 0,Nmesh_thata - t2d=DBLE(num)*2D0*pi/Nmesh_thata + t2d = DBLE(num)*2D0*pi/Nmesh_thata new_num = (360.0d0/Nmesh_thata)*num - call slic3D_calc(mmx,kky,llz,smkl,num, Nmesh_thata, vec_3dslic) + CALL slic3D_calc(mmx,kky,llz,smkl,num, Nmesh_thata, vec_3dslic) WRITE(50,"(7F30.20)")new_num,SQRT( (young2dmax(num)*COS(t2d) )**2D0 + (young2dmax(num)*SIN(t2d) )**2D0),& SQRT( (young2dmin(num)*COS(t2d) )**2D0 + (young2dmin(num)*SIN(t2d) )**2D0), & vec_3dslic(1)*young2dmax(num), vec_3dslic(2)*young2dmax(num) ,vec_3dslic(3)*young2dmax(num) WRITE(52,"(5F30.20)")new_num,SQRT( (shear2dmax(num)*COS(t2d) )**2D0 + (shear2dmax(num)*SIN(t2d) )**2D0),& - SQRT( (shear2dminp(num)*COS(t2d) )**2D0 + (shear2dminp(num)*SIN(t2d) )**2D0),& - SQRT( (shear2dminn(num)*COS(t2d) )**2D0 + (shear2dminn(num)*SIN(t2d) )**2D0) + SQRT( (shear2dminp(num)*COS(t2d) )**2D0 + (shear2dminp(num)*SIN(t2d) )**2D0),& + SQRT( (shear2dminn(num)*COS(t2d) )**2D0 + (shear2dminn(num)*SIN(t2d) )**2D0) - WRITE(55,"(5F25.14)")new_num,SQRT( (comMAX2d(num)*COS(t2d) )**2D0 + (comMAX2d(num)*SIN(t2d) )**2D0),& - SQRT( (comMIN2d(num)*COS(t2d) )**2D0 + (comMIN2d(num)*SIN(t2d) )**2D0) + WRITE(55,"(5F25.14)") new_num,SQRT( (comMAX2d(num)*COS(t2d) )**2D0 + (comMAX2d(num)*SIN(t2d) )**2D0),& + SQRT( (comMIN2d(num)*COS(t2d) )**2D0 + (comMIN2d(num)*SIN(t2d) )**2D0) WRITE(57,"(5F30.20)")new_num,SQRT( (poisson2dmax(num)*COS(t2d) )**2D0 + (poisson2dmax(num)*SIN(t2d) )**2D0),& - SQRT( (poisson2dminp(num)*COS(t2d))**2D0 + (poisson2dminp(num)*SIN(t2d))**2D0),& - SQRT( (poisson2dminn(num)*COS(t2d))**2D0 + (poisson2dminn(num)*SIN(t2d))**2D0) + SQRT( (poisson2dminp(num)*COS(t2d))**2D0 + (poisson2dminp(num)*SIN(t2d))**2D0),& + SQRT( (poisson2dminn(num)*COS(t2d))**2D0 + (poisson2dminn(num)*SIN(t2d))**2D0) WRITE(60,"(5F30.11)")new_num,SQRT( (bulkmin2d(num)*COS(t2d) )**2D0 + (bulkmin2d(num)*SIN(t2d) )**2D0),& - SQRT( (bulkMax2d(num)*COS(t2d) )**2D0 + (bulkmax2d(num)*SIN(t2d) )**2D0) + SQRT( (bulkMax2d(num)*COS(t2d) )**2D0 + (bulkmax2d(num)*SIN(t2d) )**2D0) WRITE(62,"(5F30.20)")new_num,SQRT( (LM2d(num)*COS(t2d) )**2D0 + (LM2d(num) *SIN(t2d) )**2D0),& - SQRT( (TMmax2d(num)*COS(t2d) )**2D0 + (TMmax2d(num)*SIN(t2d) )**2D0),& - SQRT( (TMmin2d(num)*COS(t2d) )**2D0 + (TMmin2d(num)*SIN(t2d) )**2D0) + SQRT( (TMmax2d(num)*COS(t2d) )**2D0 + (TMmax2d(num)*SIN(t2d) )**2D0),& + SQRT( (TMmin2d(num)*COS(t2d) )**2D0 + (TMmin2d(num)*SIN(t2d) )**2D0) WRITE(71,"(5F30.20)")new_num,SQRT( (pugh2dmax(num)*COS(t2d) )**2D0 + (pugh2dmax(num)*SIN(t2d) )**2D0),& - SQRT( (pugh2dminp(num)*COS(t2d) )**2D0 + (pugh2dminp(num)*SIN(t2d) )**2D0),& - SQRT( (pugh2dminn(num)*COS(t2d) )**2D0 + (pugh2dminn(num)*SIN(t2d) )**2D0) + SQRT( (pugh2dminp(num)*COS(t2d) )**2D0 + (pugh2dminp(num)*SIN(t2d) )**2D0),& + SQRT( (pugh2dminn(num)*COS(t2d) )**2D0 + (pugh2dminn(num)*SIN(t2d) )**2D0) WRITE(72,"(5F30.20)")new_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_km=='Y' .or. yn_km=='y') then + IF (yn_km == 'Y' .or. yn_km == 'y') THEN WRITE(499,"(5F30.20)")new_num,SQRT( (km2dmax(num)*COS(t2d) )**2D0 + (km2dmax(num)*SIN(t2d) )**2D0) - endif + ENDIF ENDDO - IF (yn_veloc=='Y' .or. yn_veloc=='y') then + IF (yn_veloc=='Y' .or. yn_veloc=='y') THEN DO num=0,360 t2d=DBLE(num)/360D0*2D0*PI @@ -2249,7 +2520,7 @@ PROGRAM AAEP_main WRITE(502,'(I3,3F30.15)')num,SQRT( (VV_P_PF_2D(num)*COS(t2d) )**2D0 + (VV_P_PF_2D(num)*SIN(t2d) )**2D0),& SQRT( (VV_Sf_PF_2D(num)*COS(t2d) )**2D0 + (VV_Sf_PF_2D(num)*SIN(t2d))**2D0),& SQRT( (VV_Ss_PF_2D(num)*COS(t2d) )**2D0 + (VV_Ss_PF_2D(num)*SIN(t2d))**2D0) - enddo + ENDDO ENDIF CLOSE(50) @@ -2257,23 +2528,23 @@ PROGRAM AAEP_main CLOSE(52) ! CLOSE(53) !CLOSE(55) - CLOSE(55) + CLOSE (55) ! CLOSE(56) - CLOSE(57) + CLOSE (57) !CLOSE(58) ! CLOSE(59) - CLOSE(60) + CLOSE (60) ! CLOSE(61) - CLOSE(62) + CLOSE (62) !CLOSE(63) ! CLOSE(64) - CLOSE(71) - CLOSE(72) + CLOSE (71) + CLOSE (72) !CLOSE(73) - CLOSE(500) - CLOSE(501) - CLOSE(502) - CLOSE(499) + CLOSE (500) + CLOSE (501) + CLOSE (502) + CLOSE (499) ! CLOSE(504) ! CLOSE(505) ! CLOSE(506) @@ -2283,11 +2554,11 @@ PROGRAM AAEP_main WRITE(*,*)"#==================================#" PRINT*," > Do you want to plot data? (Y/N):" READ(*,*)YN - call system('clear') + CALL system('clear') IF (YN.EQ."Y".OR. YN.EQ."y")THEN IF (YN.EQ."y") YN='Y' - CALL SYSTEM('if [ -e PicFile ]; then rm -rf PicFile; fi') - CALL SYSTEM('if [ -e DatFile ]; then rm -rf DatFile; fi') + CALL SYSTEM('IF [ -e PicFile ]; THEN rm -rf PicFile; fi') + CALL SYSTEM('IF [ -e DatFile ]; THEN rm -rf DatFile; fi') CALL SYSTEM('mkdir PicFile') CALL SYSTEM('mkdir DatFile') @@ -2295,13 +2566,13 @@ PROGRAM AAEP_main WRITE(*,*)"==================================" CALL plotpugh3d(YN) WRITE(*,*)"==================================" - if (Mincomp < 0d0) CALL plotcomp3d(YN,1) - if (Mincomp > 0d0) CALL plotcomp3d(YN,2) + IF (Mincomp < 0d0) CALL plotcomp3d(YN,1) + IF (Mincomp > 0d0) CALL plotcomp3d(YN,2) WRITE(*,*)"==================================" CALL plotbulk3d(YN) WRITE(*,*)"==================================" - if (Pratio_min < 0d0) CALL plotRatio3d(YN,1) - if (Pratio_min > 0d0) CALL plotRatio3d(YN,2) + IF (Pratio_min < 0d0) CALL plotRatio3d(YN,1) + IF (Pratio_min > 0d0) CALL plotRatio3d(YN,2) !CALL plotRatio3d(YN) WRITE(*,*)"==================================" CALL plotyoung3d(YN) @@ -2314,7 +2585,7 @@ PROGRAM AAEP_main CALL SLEEP (1) CALL SYSTEM('mv *.ps PicFile') CALL SYSTEM('mv *.png PicFile') - CALL SYSTEM('if [ -e 2DPugh.ps ]; then mv 2DPugh.ps PicFile; fi') + CALL SYSTEM('IF [ -e 2DPugh.ps ]; THEN mv 2DPugh.ps PicFile; fi') !CALL SYSTEM('cat PicFile/2Dbox.ps | epstopdf --filter > PicFile/2Dbox.pdf' ) !CALL SYSTEM('cat PicFile/2DBulk.ps | epstopdf --filter > PicFile/2DBulk.pdf' ) !CALL SYSTEM('cat PicFile/2DCompressibiliy.ps | epstopdf --filter > PicFile/2DCompressibiliy.pdf') @@ -2331,8 +2602,8 @@ PROGRAM AAEP_main CALL SYSTEM('mv .MaMiout DatFile') CALL SYSTEM('mv MESH DatFile') CALL SYSTEM('mv .MaMiout2 DatFile') - CALL SYSTEM('if [ -e .aelastpro2 ]; then mv .aelastpro2 DatFile; fi') - !CALL SYSTEM('if [ -e HKL ]; then mv HKL DatFile; fi') + CALL SYSTEM('IF [ -e .aelastpro2 ]; THEN mv .aelastpro2 DatFile; fi') + !CALL SYSTEM('IF [ -e HKL ]; THEN mv HKL DatFile; fi') ELSE !CALL zerofile() PRINT*," > Data FILEs generated:" @@ -2346,13 +2617,13 @@ PROGRAM AAEP_main CALL SYSTEM('mv .MaMiout DatFile') CALL SYSTEM('mv MESH DatFile') CALL SYSTEM('mv .MaMiout2 DatFile') - CALL SYSTEM('if [ -e .aelastpro2 ]; then mv .aelastpro2 DatFile; fi') + CALL SYSTEM('IF [ -e .aelastpro2 ]; THEN mv .aelastpro2 DatFile; fi') ENDIF !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ CALL SYSTEM('rm -rf *.gnu') ELSE - IF(d2d3 == 2) then !@@@@@@@@@@@@@@@@@@@@@@@ 2D_3D system start + IF(d2d3 == 2) THEN !@@@@@@@@@@@@@@@@@@@@@@@ 2D_3D system start CALL analiz_2D_sys(adv,npoint) @@ -2363,8 +2634,8 @@ PROGRAM AAEP_main IF (YN.EQ."Y".OR. YN.EQ."y")THEN IF (YN.EQ."y") YN='Y' - CALL SYSTEM('if [ -e PicFile ]; then rm -rf PicFile; fi') - CALL SYSTEM('if [ -e DatFile ]; then rm -rf DatFile; fi') + CALL SYSTEM('IF [ -e PicFile ]; THEN rm -rf PicFile; fi') + CALL SYSTEM('IF [ -e DatFile ]; THEN rm -rf DatFile; fi') CALL SYSTEM('mkdir PicFile') CALL SYSTEM('mkdir DatFile') @@ -2376,8 +2647,8 @@ PROGRAM AAEP_main CALL SYSTEM('cp DatFile/Cij-2D.dat . ') !CALL SYSTEM(' mv MAX1 DatFile ') CALL SYSTEM('mv .MaMiout DatFile') - !call system("dat2gnu_lapw D2") - Else + !CALL system("dat2gnu_lapw D2") + ELSE PRINT*," > Data FILEs generated:" CALL SYSTEM('mkdir DatFile') CALL SYSTEM('mv *.dat DatFile ') @@ -2385,13 +2656,13 @@ PROGRAM AAEP_main CALL SYSTEM('ls DatFile/*.dat ') !CALL SYSTEM(' mv MAX1 DatFile ') CALL SYSTEM(' mv .MaMiout DatFile ') - Endif + ENDIF END IF END IF !@@@@@@@@@@@@@@@@@@@@@@@ 2D_3D system END CLOSE(99) - call checkfolder_old(1,0,d2d3) - IF (YN.EQ."Y".OR. YN.EQ."y") then ; call checkfolder_old(0,1,d2d3); endif + CALL checkfolder_old(1,0,d2d3) + IF (YN.EQ."Y".OR. YN.EQ."y") THEN ; CALL checkfolder_old(0,1,d2d3); ENDIF CALL SYSTEM('tput setaf 55;tput bold; echo " ======================================= Visualization Tools =======================================";tput sgr0') CALL SYSTEM('tput setaf 55;tput bold; echo " - dat2gnu.x | To visualize data in polar, spherical coordinates and heat-maps diagrams: 2D sys., 2D projecton and 3D sys. ";tput sgr0') @@ -2399,19 +2670,26 @@ PROGRAM AAEP_main CALL SYSTEM('tput setaf 55;tput bold; echo " - dat2wrl.x | To visualize data in spherical coordinates : 3D sys. ";tput sgr0') CALL SYSTEM('tput setaf 55;tput bold; echo " - dat2html.x | To visualize data in spherical coordinates : 3D sys. ";tput sgr0') CALL SYSTEM('tput setaf 55;tput bold; echo " ===================================================================================================";tput sgr0') - + WRITE(*,*)"" + WRITE(*,*)"For any use of ElATools, please cite:" + WRITE(*,*)"Shahram Yalameha, Zahra Nourbakhsh, Daryoosh Vashaee," + WRITE(*,*)"ElATools: A tool for analyzing anisotropic elastic " + WRITE(*,*)"properties of the 2D and 3D materials," + WRITE(*,*)"Computer Physics Communications 271, 108195 (2022)" + WRITE(*,*)"https://doi.org/10.1016/j.cpc.2021.108195" + WRITE(*,*)"" WRITE(*,*)" ===================================================" PRINT*, " Nothing happens until something moves. " PRINT*, " >> Albert Einstion << " WRITE(*,*)" ---------------------------------------------------" - PRINT*, " > | Have A Beautiful Day | <" + PRINT*, " > | Have A BeautIFul Day | <" PRINT*, " > | Goodbye. | <" WRITE(*,*)" ===================================================" WRITE(*,*)"" - call creadfolder(1, 0,mmx1,kky1,llz1) + CALL creadfolder(1, 0,mmx1,kky1,llz1) - IF (YN.EQ."Y".OR. YN.EQ."y") then ; call creadfolder(0, 1,mmx1,kky1,llz1); endif + IF (YN.EQ."Y".OR. YN.EQ."y") THEN ; CALL creadfolder(0, 1,mmx1,kky1,llz1); ENDIF STOP 1369 WRITE(*,*) " > NOT FOUNDE Cij.dat FILE" ; STOP @@ -2421,5 +2699,4 @@ PROGRAM AAEP_main 1366 WRITE(*,*) " > NOT FOUNDE elast.output FILE " ; STOP 1361 WRITE(*,*) " > NOT FOUNDE ELADAT FILE " ; STOP 1341 WRITE(*,*) " > NOT FOUNDE ElaStic_2nd.out FILE " ; STOP -1370 END PROGRAM - +1370 END PROGRAM \ No newline at end of file diff --git a/soc/Eatools_poi.f90 b/soc/Eatools_poi.f90 index a3b4ed0..9a2ed0d 100644 --- a/soc/Eatools_poi.f90 +++ b/soc/Eatools_poi.f90 @@ -7,7 +7,7 @@ SUBROUTINE CPratio(a1111,MinPratio,MaxPratio,AvePratio,& vec1min,vec2min,vec3min,vec1max,vec2max,vec3max,& - theta,phi,k11,k12,k13,k22,k23,k33,Pratio) + theta,phi,k11,k12,k13,k22,k23,k33,Pratio) DOUBLE PRECISION, PARAMETER :: pi=3.1415926535897932384626433832795d0 DOUBLE PRECISION, DIMENSION(6,6) :: S=0D0 @@ -72,19 +72,19 @@ SUBROUTINE CPratio(a1111,MinPratio,MaxPratio,AvePratio,& + k33*v33*S(3,3) - a1122 =a1122 & + a1122 =a1122 & +(k11*v23+k23*v11)*S(1,4) & +(k11*v13+k13*v11)*S(1,5) & +(k11*v12+k12*v11)*S(1,6) & - +(k22*v23+k23*v22)*S(2,4) & + +(k22*v23+k23*v22)*S(2,4) & +(k22*v13+k13*v22)*S(2,5) & +(k22*v12+k12*v22)*S(2,6) & - +(k33*v23+k23*v33)*S(3,4) & + +(k33*v23+k23*v33)*S(3,4) & +(k33*v13+k13*v33)*S(3,5) & +(k33*v12+k12*v33)*S(3,6)+k23*v23*S(4,4) & +(k23*v13+k13*v23)*S(4,5) & +(k23*v12+k12*v23)*S(4,6)+k13*v13*S(5,5) & - +(k13*v12+k12*v13)*S(5,6)+k12*v12*S(6,6) + +(k13*v12+k12*v13)*S(5,6)+k12*v12*S(6,6) Pratio = -a1122/a1111 IF (kk.EQ.0) THEN MaxPratio=Pratio; MinPratio=Pratio; diff --git a/soc/Eatools_sh.f90 b/soc/Eatools_sh.f90 index ceebaf3..cb25e54 100644 --- a/soc/Eatools_sh.f90 +++ b/soc/Eatools_sh.f90 @@ -65,7 +65,7 @@ SUBROUTINE CShear(Minsheainvar,Maxsheainvar,sheainvar_ave,phi,theta,k11,k12,k13, !Maxsheainvar=Maxsheainvar*100D0 ENDIF IF (sheainvar.LE.Minsheainvar)THEN - Minsheainvar = (sheainvar) + Minsheainvar = sheainvar ENDIF ENDIF ave=ave+sheainvar diff --git a/soc/Eatools_welc.f90 b/soc/Eatools_welc.f90 index e7cd7e9..0c0e3f0 100644 --- a/soc/Eatools_welc.f90 +++ b/soc/Eatools_welc.f90 @@ -19,10 +19,10 @@ SUBROUTINE WELCOME() WRITE(*,'(a)')" 45 6C 41 74 6F 6F 6C 73 " WRITE(*,*)' Author : Shahram Yalameha ' WRITE(*,*)' Email : yalameha93@gmail.com ' -WRITE(*,'(a)')' v1.7.2 ' +WRITE(*,'(a)')' v1.7.3 ' WRITE(*,'(a)')' From ' WRITE(*,'(a)')' Elastic Tools Project ' -WRITE(*,'(a)')'========= (c)2018 ===========' +WRITE(*,'(a)')'========= (c)2018-2022 ==========' WRITE(*,*)' ' !!!! @@ -39,10 +39,10 @@ SUBROUTINE WELCOME() WRITE(99,'(a)')" 45 6C 41 74 6F 6F 6C 73 " WRITE(99,*)' Author : Shahram Yalameha ' WRITE(99,*)' Email : yalameha93@gmail.com ' -WRITE(99,'(a)')' v1.7.2 ' +WRITE(99,'(a)')' v1.7.3 ' WRITE(99,'(a)')' From ' WRITE(99,'(a)')' Elastic Tools Project ' -WRITE(99,'(a)')'========= (c)2018 ===========' +WRITE(99,'(a)')'========= (c)2018-2022 ==========' WRITE(99,*)' ' WRITE(99,*)' ' diff --git a/soc/Makefile b/soc/Makefile index 18e3c30..03777c2 100644 --- a/soc/Makefile +++ b/soc/Makefile @@ -13,20 +13,20 @@ #---------------- compiler ----------------- -#FC = ifort - FC = gfortran +FC = ifort +# FC = gfortran #------------- compiler options ------------ -#FOPT = -O1 -FR -mp1 -w -prec_div -pc80 -pad -ip -DINTEL_VML -traceback -assume buffered_io -I$(MKLROOT)/include # for ifort -#FGEN = -O1 -FR -mp1 -w -prec_div -pc80 -pad -ip -DINTEL_VML -traceback # for ifort + FOPT = -O1 -FR -mp1 -w -prec_div -pc80 -pad -ip -DINTEL_VML -traceback -assume buffered_io -I$(MKLROOT)/include # for ifort + FGEN = -O1 -FR -mp1 -w -prec_div -pc80 -pad -ip -DINTEL_VML -traceback # for ifort - FOPT = -ffree-form -O2 -ffree-line-length-none # for gfortran - FGEN = -ffree-form -O2 -ffree-line-length-none # for gfortran + # FOPT = -ffree-form -O2 -ffree-line-length-none # for gfortran + # FGEN = -ffree-form -O2 -ffree-line-length-none # for gfortran #-------------- loader options ------------- - LDFLAGS = $(FOPT) -lblas -llapack -lpthread # for gfortran -#LDFLAGS = $(FOPT) -L$(MKLROOT)/lib/$(MKL_TARGET_ARCH) -pthread -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -qopenmp -lpthread # for ifort + # LDFLAGS = $(FOPT) -lblas -llapack -lpthread # for gfortran +LDFLAGS = $(FOPT) -L$(MKLROOT)/lib/$(MKL_TARGET_ARCH) -pthread -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -qopenmp -lpthread # for ifort @@ -57,7 +57,7 @@ OBJS = Eatools_main.o Eatools_2Dyoung.o Eatools_2Dpoisson.o Eatools_2Dshear.o Eatools_3dPLcom.o Eatools_stability.o Eatools_db.o Eatools_3dPLsh.o Eatools_ang2car.o Eatools_dsyec3.o \ Eatools_dutester.o Eatools_proelast.o Eatools_sh.o Eatools_welc.o Eatools_pugh.o Eatools_3dPLpugh.o Eatools_hardness.o \ Eatools_km.o Eatools_wave_start.o Eatools_wave_vgvf.o Eatools_wave_rot.o Eatools_wave_cal.o Eatools_wave_main.o Eatools_2dcal_wave.o \ - Eatools_pro_wave.o Eatools_2Dadv.o Eatools_api.o Eatools_3dslic.o + Eatools_pro_wave.o Eatools_2Dadv.o Eatools_api.o Eatools_3dslic.o Eatools_bulk.o help.o OBJS1 = wrl_calylm.o wrl_colorsign.o wrl_end.o wrl_main.o \ wrl_mesh.o wrl_shape.o wrl_spherappear.o wrl_spher.o wrl_start.o wrl_setcolor.o diff --git a/soc/help.f90 b/soc/help.f90 new file mode 100644 index 0000000..836952e --- /dev/null +++ b/soc/help.f90 @@ -0,0 +1,35 @@ +SUBROUTINE help_me() + write(*,"(A)")"----------------------------------------------------------------------------------------------------------------------------------------" + write(*,"(A)")" ## Command-line Options: ##" + write(*,"(A)")"----------------------------------------------------------------------------------------------------------------------------------------" + write(*,"(A)")"Elatools.x -d [Integer] -m [Integer/String] -ewp [String] -ktc [String] -op [String]" + write(*,"(A)")" | | | | | " + write(*,"(A)")" | | | | ^ " + write(*,"(A)")" | | | | Type of 2D system " + write(*,"(A)")" | | | | 2D: def => Default option| adv => Advanced option" + write(*,"(A)")" | | | ^ " + write(*,"(A)")" | | | Calculate Min. thermal conductivity " + write(*,"(A)")" | | | 3D: y => Yes| n => No " + write(*,"(A)")" | | ^ " + write(*,"(A)")" | | Calculate elastic wave properties " + write(*,"(A)")" | | 3D: y => Yes| n => No " + write(*,"(A)")" | ^ " + write(*,"(A)")" | Output code " + write(*,"(A)")" | 3D: 1 => IRelast| 2 => Elast| 3 => AELAS| 4 => ElaStic| 5 => Cij file| 6 => offline Databank| 7 => online Databank " + write(*,"(A)")" | 2D: 1 => AELAS | 2 => IRelast2D| 3 => Cij-2D file " + write(*,"(A)")" ^" + write(*,"(A)")" System Type " + write(*,"(A)")" 3 => 3D sys. and 2 => 2D sys." + write(*,"(A)")"" + write(*,"(A)")"----------------------------------------------------------------------------------------------------------------------------------------" + + write(*,"(A)")" > Example for 3D system:" + write(*,"(A)")" Elatools.x -d 3 -m 5" + write(*,"(A)")" Elatools.x -d 3 -m 3 -ewp y" + write(*,"(A)")" Elatools.x -d 3 -m 1 -ewp y -ktc y" + write(*,"(A)")"" + write(*,"(A)")" > Example for 2D system:" + write(*,"(A)")" Elatools.x -d 2 -m 2 -op def" + write(*,"(A)")" Elatools.x -d 2 -m 3 -op adv" + +END SUBROUTINE diff --git a/soc/html_elayout.f90 b/soc/html_elayout.f90 index 042a2f1..8381ef3 100644 --- a/soc/html_elayout.f90 +++ b/soc/html_elayout.f90 @@ -90,11 +90,11 @@ SUBROUTINE buttone_polar_web(max_value) ChARACTER(len=30) :: title ChARACTER(len=10) :: namepro DOUBLE PRECISION :: max_value - if (namepro == "poi2d" ) title="Poisson\'s ratio" + if (namepro == "poi2d" .or. namepro == "2dpoi" ) title="Poisson\'s ratio" if (namepro == "pugh2d" ) title="Pugh\'s ratio" - if (namepro == "young2d") title="Young\'s modulus" + if (namepro == "young2d" .or. namepro == "2dyoung") title="Young\'s modulus" if (namepro == "bulk2d") title="Bulk modulus" - if (namepro == "shear2d") title="Shear modulus" + if (namepro == "shear2d" .or. namepro == "2dshear") title="Shear modulus" if (namepro == "com2d" ) title="Linear compressibility" if (namepro == "hard2d" ) title="Hardness" if (namepro == "pp2d" ) title="P-mode of Phase velocity" diff --git a/soc/html_gerdata_polar.f90 b/soc/html_gerdata_polar.f90 index de163c9..0c471e6 100644 --- a/soc/html_gerdata_polar.f90 +++ b/soc/html_gerdata_polar.f90 @@ -17,13 +17,44 @@ SUBROUTINE get_dataplotly_polar( namepro,rtheta, n_phif,n_thetaf,cutmesh,type_pr VVP_i_x, VVP_P_y, VVP_SF_z,VVP_SS_k,& VVG_i_x, VVG_P_y, VVG_SF_z,VVG_SS_k,& VVF_i_x, VVF_P_y, VVF_SF_z,VVF_SS_k,& - km_x, km_y + km_x, km_y ,temp_x, temp_y INTEGER, DIMENSION(190300,4) :: mesh=0 - INTEGER :: n_phif, n_thetaf,cutmesh ,num_mesh,i,ii,argl,k,& + INTEGER :: n_phif, n_thetaf,cutmesh ,num_mesh,i,ii,argl,k,io,tem_i,& start_new_reng,end_new_reng,rtheta character(len=10) :: val='',namepro ,type_pro ! type_pro : max, min, neg + +IF (namepro=="2dpoi" .OR. namepro=="2dshear" .OR. namepro=="2dyoung" .OR. namepro=="2dyou" .OR. namepro=="2dshe") THEN + tem_i=0 + OPEN(12, file="young_2d_sys.dat") + DO + tem_i=tem_i + 1 + READ(12,*, iostat = io ) temp_x, temp_y + if (io < 0) exit + ENDDO + Close(12) +!--------------------------------------------------- + IF (namepro=="2dyoung") THEN + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! loop-READer + DO ii=1,tem_i+1 + open(10, file="young_2d_sys.dat") + READ(10,*) young_x(ii),young_y(ii) + ENDDO + CLOSE(10) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end loop + If(type_pro=='max') THEN + If(rtheta==1) THEN + DO i=1,tem_i+1 + WRITE(66,"(F23.15,A)") young_x(i),"," + ENDDO + endif + Endif + ENDIF +!=========================================== + +ENDIF + !=============================================READ mesh OPEN(69, file="MESH") READ(69,*)n_phif,n_thetaf,cutmesh @@ -517,4 +548,5 @@ SUBROUTINE get_dataplotly_polar( namepro,rtheta, n_phif,n_thetaf,cutmesh,type_pr endif Endif ENDIF +!$============================================================================================ end SUBROUTINE diff --git a/soc/html_swin.f90 b/soc/html_swin.f90 index f4e4dd5..1bd14a2 100644 --- a/soc/html_swin.f90 +++ b/soc/html_swin.f90 @@ -27,12 +27,12 @@ SUBROUTINE swin_web(namepro) if (namepro=="gs" )WRITE(66,"(3a)")" Spatial dependence of Group-Slow velocity " if (namepro=="km" )WRITE(66,"(3a)")" Spatial dependence of of Min. thermal conductivity " !--------------------------------------------- - if (namepro=="young2d")WRITE(66,"(3a)")" Orientation dependence of Young's modulus " + if (namepro=="young2d" .or. namepro=="2dyoung")WRITE(66,"(3a)")" Orientation dependence of Young's modulus " if (namepro=="bulk2d" )WRITE(66,"(3a)")" Orientation dependence of Bulk modulus " if (namepro=="pugh2d" )WRITE(66,"(3a)")" Orientation dependence of Pugh's ratio " - if (namepro=="shear2d")WRITE(66,"(3a)")" Orientation dependence of Shear modulus " + if (namepro=="shear2d" .or. namepro=="2dshear")WRITE(66,"(3a)")" Orientation dependence of Shear modulus " if (namepro=="com2d" )WRITE(66,"(3a)")" Orientation dependence of linear compressibility " - if (namepro=="poi2d" )WRITE(66,"(3a)")" Orientation dependence of Poisson's ratio " + if (namepro=="poi2d" .or. namepro=="2dpoi" )WRITE(66,"(3a)")" Orientation dependence of Poisson's ratio " if (namepro=="hard2d" )WRITE(66,"(3a)")" Orientation dependence of Hardness " if (namepro=="pp2d" )WRITE(66,"(3a)")" Orientation dependence of Phase-P velocity " if (namepro=="ps2d" )WRITE(66,"(3a)")" Orientation dependence of Phase-Fast velocity " diff --git a/soc/wrl_main.f90 b/soc/wrl_main.f90 index 03e9c4d..7ea15aa 100644 --- a/soc/wrl_main.f90 +++ b/soc/wrl_main.f90 @@ -66,7 +66,7 @@ PROGRAM wrl_conv DOUBLE PRECISION, DIMENSION(1910000) :: datapoints=0d0 DOUBLE PRECISION, DIMENSION(1910000) :: G_max,shminp,shminn,shavep,SINver,CO,comminp,pugh_max,pughminp,pughminn,pughavep,& - comminn,NPratio_max,pminp,pminn,pavep,paven,& + comminn,NPratio_max,pminp,pminn,pavep,paven,& BINver,maxEVaLM1,maxEVaTM1,minEVaTM1,VVG_P,VVP_P,VV_P_PF,VVG_Sf,VVP_Sf,VV_Sf_PF,hardvar ,& VVG_Ss,VVP_Ss,VV_Ss_PF,km ChARACTER(len=7), dimension(10) :: arg_mane