diff --git a/lib/mesont/CNTPot.f90 b/lib/mesont/CNTPot.f90 index 74fa0dd9ef..c74b485274 100644 --- a/lib/mesont/CNTPot.f90 +++ b/lib/mesont/CNTPot.f90 @@ -134,7 +134,7 @@ implicit none real(c_double) :: CNTSTRUl, CNTSTRUf ! ! Axial buckling - hysteresis approach - real(c_double) :: CNTSTREc = -0.0142d+00 ! The minimum buckling strain + real(c_double) :: CNTSTREc = -0.0142d+00 ! The minimum buckling strain real(c_double) :: CNTSTREc1 = -0.04d+00 ! Critical axial buckling strain real(c_double) :: CNTSTREc2 = -0.45d+00 ! Maximum buckling strain @@ -249,7 +249,7 @@ contains !********************************************************************** integer(c_int) function CNTSTRH1BCalc ( U, dUdL, L, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Young's modulus depends on R, see [4]. - ! Axial buckling without hysteresis. + ! Axial buckling without hysteresis. !------------------------------------------------------------------------------------------- real(c_double), intent(out) :: U, dUdL real(c_double), intent(in) :: L, R0, L0 @@ -275,7 +275,7 @@ contains !********************************************************************** U = 0.5d+00 * L0 * (d+E-CNTSTREc2) * dUdL + dUbcl - Ud CNTSTRH1BCalc = CNTPOT_STRETCHING end if - end function CNTSTRH1BCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + end function CNTSTRH1BCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Stretching without fracture, harmonic potential, with axial buckling with hysteresis @@ -289,49 +289,49 @@ contains !********************************************************************** integer(c_int), intent(in) :: ABF !------------------------------------------------------------------------------------------- real(c_double) :: E, K, dUbcl, Ebcl, Kbcl, Edu - real(c_double) :: C, DE, t + real(c_double) :: C, DE, t !------------------------------------------------------------------------------------------- E = ( L - L0 ) / L0 - K = 86.64d+00 + 100.56d+00 * R0 - Kbcl = -10.98d+00 * L0 - if ( E .gt. CNTSTREc ) then ! Harmonic potential - no buckling - dUdL = K * E + K = 86.64d+00 + 100.56d+00 * R0 + Kbcl = -10.98d+00 * L0 + if ( E .gt. CNTSTREc ) then ! Harmonic potential - no buckling + dUdL = K * E U = 0.5d+00 * L0 * E * dUdL CNTSTRH1BHCalc = CNTPOT_STRETCHING Ebuc = 0.0d+00 - else if ( E .gt. CNTSTREc1 ) then ! Above minimal buckling strain, but not at critical strain - if ( ABF .eq. 0 ) then ! Not buckled. Continue harmonic potential - dUdL = K * E - U = 0.5d+00 * L0 * E * dUdL - CNTSTRH1BHCalc = CNTPOT_STRETCHING - Ebuc = 0.0d+00 - else ! Relaxing from buckled state. Use buckling potential - dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc - U = Kbcl * E + dUbcl - dUdL = Kbcl / L0 - CNTSTRH1BHCalc = CNTPOT_SBUCKLING - Ebuc = 0.0d+00 - end if - else if( E .gt. CNTSTREc2 ) then ! Axial buckling strain region - if ( ABF .eq. 0 ) then ! Newly buckled - dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc - U = Kbcl * E + dUbcl - dUdL = Kbcl / L0 - CNTSTRH1BHCalc = CNTPOT_SBUCKLING - Ebuc = 0.5d+00 * L0 * K * CNTSTREc1 * CNTSTREc1 - Kbcl * CNTSTREc1 - dUbcl - else ! Already buckled - dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc - U = Kbcl * E + dUbcl - dUdL = Kbcl / L0 - CNTSTRH1BHCalc = CNTPOT_SBUCKLING - Ebuc = 0.0d+00 - end if - else ! Maximum strain and return to harmonic potential - dUdL = K * E - U = 0.5d+00 * L0 * E * dUdL - CNTSTRH1BHCalc = CNTPOT_STRETCHING - Ebuc = 0.0d+00 - end if + else if ( E .gt. CNTSTREc1 ) then ! Above minimal buckling strain, but not at critical strain + if ( ABF .eq. 0 ) then ! Not buckled. Continue harmonic potential + dUdL = K * E + U = 0.5d+00 * L0 * E * dUdL + CNTSTRH1BHCalc = CNTPOT_STRETCHING + Ebuc = 0.0d+00 + else ! Relaxing from buckled state. Use buckling potential + dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc + U = Kbcl * E + dUbcl + dUdL = Kbcl / L0 + CNTSTRH1BHCalc = CNTPOT_SBUCKLING + Ebuc = 0.0d+00 + end if + else if( E .gt. CNTSTREc2 ) then ! Axial buckling strain region + if ( ABF .eq. 0 ) then ! Newly buckled + dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc + U = Kbcl * E + dUbcl + dUdL = Kbcl / L0 + CNTSTRH1BHCalc = CNTPOT_SBUCKLING + Ebuc = 0.5d+00 * L0 * K * CNTSTREc1 * CNTSTREc1 - Kbcl * CNTSTREc1 - dUbcl + else ! Already buckled + dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc + U = Kbcl * E + dUbcl + dUdL = Kbcl / L0 + CNTSTRH1BHCalc = CNTPOT_SBUCKLING + Ebuc = 0.0d+00 + end if + else ! Maximum strain and return to harmonic potential + dUdL = K * E + U = 0.5d+00 * L0 * E * dUdL + CNTSTRH1BHCalc = CNTPOT_STRETCHING + Ebuc = 0.0d+00 + end if end function CNTSTRH1BHCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! @@ -449,7 +449,7 @@ contains !********************************************************************** CNTSTRUl = 0.5d+00 * CNTSTRFl * CNTSTREl CNTSTRUf = CNTSTRUl + ( CNTSTRFl + CNTSTRAA ) * ( CNTSTREf - CNTSTREl ) - CNTSTRAA * dlog ( C ) / CNTSTRBB end subroutine CNTSTRNH1Init !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - + ! ! General ! @@ -457,9 +457,9 @@ contains !********************************************************************** integer(c_int) function CNTSTRCalc ( U, dUdL, L, R0, L0 , ABF, Ebuc ) !!!!!!!!!!!!!!!!!!!!!!! real(c_double), intent(out) :: U, dUdL, Ebuc real(c_double), intent(in) :: L, R0, L0 - integer(c_int), intent(in) :: ABF + integer(c_int), intent(in) :: ABF !------------------------------------------------------------------------------------------- - Ebuc = 0.0d+00 + Ebuc = 0.0d+00 select case ( CNTSTRModel ) case ( CNTSTRMODEL_H0 ) CNTSTRCalc = CNTSTRH0Calc ( U, dUdL, L, R0, L0 ) @@ -471,9 +471,9 @@ contains !********************************************************************** CNTSTRCalc = CNTSTRNH1Calc ( U, dUdL, L, R0, L0 ) case ( CNTSTRMODEL_NH1F ) CNTSTRCalc = CNTSTRNH1FCalc ( U, dUdL, L, R0, L0 ) - case ( CNTSTRMODEL_H1B ) + case ( CNTSTRMODEL_H1B ) CNTSTRCalc = CNTSTRH1BCalc ( U, dUdL, L, R0, L0 ) - case ( CNTSTRMODEL_H1BH ) + case ( CNTSTRMODEL_H1BH ) CNTSTRCalc = CNTSTRH1BHCalc ( U, dUdL, L, R0, L0, ABF, Ebuc ) end select end function CNTSTRCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -601,7 +601,7 @@ contains !********************************************************************** CNTBNDHBFCalc = CNTPOT_BENDING end if end function CNTBNDHBFCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - + integer(c_int) function CNTBNDHBHCalc ( U, dUdC, C, R0, L0, BBF, Ebuc ) !!!!!!!!!!!!!!!!!!!! ! Bending model of type 1: Harmonic bending potential with buckling with hysteresis approach. !------------------------------------------------------------------------------------------- @@ -610,9 +610,9 @@ contains !********************************************************************** integer(c_int), intent(in) :: BBF real(c_double) :: E1, E2, C2, Kbnd, Kbcl,Theta,DUbcl, Ubcl, Cmin,Rmax !------------------------------------------------------------------------------------------- - Rmax = 340.0d+00 - Cmin = 1.0/(Rmax*Rmax) - E1 = 1.0d+00 - C + Rmax = 340.0d+00 + Cmin = 1.0/(Rmax*Rmax) + E1 = 1.0d+00 - C E2 = 1.0d+00 + C ! Calculate the square of curvature C2 = 4.0d+00 * E2 / ( L0 * L0 * E1 ) @@ -625,21 +625,21 @@ contains !********************************************************************** Ebuc = 0.0 else if ( C2 .ge. Cmin .and. C2 .lt. CNTBNDC2 ) then ! Potential depends on buckling flag of a node if ( BBF .eq. 0 ) then ! Not buckled yet. Continue harmonic bending - Kbnd = 2.0d+00 * ( 63.8d+00 * R0**2.93d+00 ) / L0 - U = Kbnd * E2 / E1 - dUdC = 2.0d+00 * Kbnd / ( E1 * E1 ) - CNTBNDHBHCalc = CNTPOT_BENDING - Ebuc = 0.0d+00 + Kbnd = 2.0d+00 * ( 63.8d+00 * R0**2.93d+00 ) / L0 + U = Kbnd * E2 / E1 + dUdC = 2.0d+00 * Kbnd / ( E1 * E1 ) + CNTBNDHBHCalc = CNTPOT_BENDING + Ebuc = 0.0d+00 else ! Already has been buckled or is buckled. Use buckling potential until Cmin. - Theta= M_PI - acos ( C ) - Kbnd = 63.8d+00 * R0**2.93d+00 - Kbcl = CNTBNDB * Kbnd / CNTBNDR - DUbcl= 2.0d+00*Kbnd * & + Theta= M_PI - acos ( C ) + Kbnd = 63.8d+00 * R0**2.93d+00 + Kbcl = CNTBNDB * Kbnd / CNTBNDR + DUbcl= 2.0d+00*Kbnd * & (1.0d+00+cos(l0/Rmax+M_PI))/(1.0d+00-cos(l0/Rmax+M_PI))/L0-Kbcl*abs(l0/Rmax)**CNTBNDN - U = Kbcl * abs( Theta )**CNTBNDN + DUbcl - dUdC = Kbcl * CNTBNDN * abs( Theta )**CNTBNDN1 / sqrt ( 1.0d+00 - C * C ) - Ebuc = 0.0d+00 - CNTBNDHBHCalc = CNTPOT_BBUCKLING + U = Kbcl * abs( Theta )**CNTBNDN + DUbcl + dUdC = Kbcl * CNTBNDN * abs( Theta )**CNTBNDN1 / sqrt ( 1.0d+00 - C * C ) + Ebuc = 0.0d+00 + CNTBNDHBHCalc = CNTPOT_BBUCKLING end if else ! Greater than buckling critical point if ( BBF .eq. 1 ) then ! Already buckled @@ -663,8 +663,8 @@ contains !********************************************************************** Ebuc = 2.0d+00*Kbnd * (1.0d+00+cos(l0/CNTBNDR+M_PI)) / (1.0d+00-cos(l0/CNTBNDR+M_PI))/L0 & - Kbcl * abs ( l0 / CNTBNDR ) ** CNTBNDN - dUbcl CNTBNDHBHCalc = CNTPOT_BBUCKLING - end if - end if + end if + end if end function CNTBNDHBHCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! diff --git a/lib/mesont/Makefile.gfortran b/lib/mesont/Makefile.gfortran index f8ffc38e6d..55d6f59b20 100644 --- a/lib/mesont/Makefile.gfortran +++ b/lib/mesont/Makefile.gfortran @@ -19,7 +19,7 @@ OBJ = $(SRC:.f90=.o) # ------ SETTINGS ------ F90 = gfortran -F90FLAGS = -O3 -fPIC -ftree-vectorize -g -Wall +F90FLAGS = -O3 -fPIC -ftree-vectorize -g ARCHIVE = ar ARCHFLAG = -rc USRLIB = diff --git a/lib/mesont/TPMM0.f90 b/lib/mesont/TPMM0.f90 index 0ec9ce6248..f5b55e487a 100644 --- a/lib/mesont/TPMM0.f90 +++ b/lib/mesont/TPMM0.f90 @@ -41,7 +41,7 @@ contains !********************************************************************** !------------------------------------------------------------------------------------------- real(c_double) :: Qa, Ua, Fd, L2 real(c_double), dimension(0:2) :: F1_1a, F1_2a, F2_1a, F2_2a, R2_3, R2, Laxis2, F - integer(c_int) :: IntSign + integer(c_int) :: IntSign !------------------------------------------------------------------------------------------- R2 = 0.5d+00 * ( R2_1 + R2_2 ) Laxis2 = R2_2 - R2_1 @@ -51,34 +51,34 @@ contains !********************************************************************** TPMInteractionFSS = TPMInteractionF ( Q, U, F1_1, F1_2, F2_1, F2_2, Fd, R1_1, R1_2, R2_1, R2_2, 1 ) R2_3 = R2_2 + R2_2 - R2_1 IntSign = TPMInteractionF ( Qa, Ua, F1_1a, F1_2a, F2_1a, F2_2a, Fd, R1_1, R1_2, R2_2, R2_3, 1 ) - if ( IntSign > 0 ) then - TPMInteractionFSS = 1 - call TPMSegmentForces ( F2_1a, F2_2a, F1_1a, F1_2a, R1_1, R1_2, R2, Laxis2, L2 ) - F = ( Fd - S_V3xV3 ( F2_2a, Laxis2 ) ) * Laxis2 - F2_2a = F2_2a + F - F2_1a = F2_1a - F - end if + if ( IntSign > 0 ) then + TPMInteractionFSS = 1 + call TPMSegmentForces ( F2_1a, F2_2a, F1_1a, F1_2a, R1_1, R1_2, R2, Laxis2, L2 ) + F = ( Fd - S_V3xV3 ( F2_2a, Laxis2 ) ) * Laxis2 + F2_2a = F2_2a + F + F2_1a = F2_1a - F + end if else TPMInteractionFSS = TPMInteractionF ( Q, U, F1_1, F1_2, F2_1, F2_2, Fd, R1_1, R1_2, R2_1, R2_2, 2 ) R2_3 = R2_1 + R2_1 - R2_2 IntSign = TPMInteractionF ( Qa, Ua, F1_1a, F1_2a, F2_1a, F2_2a, Fd, R1_1, R1_2, R2_1, R2_3, 1 ) - if ( IntSign > 0 ) then - TPMInteractionFSS = 1 - call TPMSegmentForces ( F2_1a, F2_2a, F1_1a, F1_2a, R1_1, R1_2, R2, Laxis2, L2 ) - F = ( - Fd - S_V3xV3 ( F2_1a, Laxis2 ) ) * Laxis2 - F2_1a = F2_1a + F - F2_2a = F2_2a - F - end if + if ( IntSign > 0 ) then + TPMInteractionFSS = 1 + call TPMSegmentForces ( F2_1a, F2_2a, F1_1a, F1_2a, R1_1, R1_2, R2, Laxis2, L2 ) + F = ( - Fd - S_V3xV3 ( F2_1a, Laxis2 ) ) * Laxis2 + F2_1a = F2_1a + F + F2_2a = F2_2a - F + end if end if - if ( IntSign > 0 ) then - Q = Q - Qa + if ( IntSign > 0 ) then + Q = Q - Qa if ( Q < 0.0d+00 ) Q = 0.0d+00 - U = U - Ua - F2_1 = F2_1 - F2_1a - F2_2 = F2_2 - F2_2a - F1_1 = F1_1 - F1_1a - F1_2 = F1_2 - F1_2a - end if + U = U - Ua + F2_1 = F2_1 - F2_1a + F2_2 = F2_2 - F2_2a + F1_1 = F1_1 - F1_1a + F1_2 = F1_2 - F1_2a + end if end function TPMInteractionFSS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! integer(c_int) function TPMInteractionFW0 ( QQ, U, U1, U2, UU, F1, F2, F, G1, G2, R1, R2, N, NMAX, R ) @@ -174,20 +174,20 @@ contains !********************************************************************** if ( TPMInteractionFSS ( QQ(i), Ua, F1_1a, F1_2a, F2_1a, F2_2a, R1, R2, R(0:2,i), R(0:2,i+1), & EType ) > 0 ) then - TPMInteractionFW0 = 1 - U = U + Ua - Ua = 0.25d+00 * Ua - U1 = U1 + Ua - U2 = U2 + Ua - UU(i) = UU(i) + Ua - UU(i+1) = UU(i+1) + Ua - F1 = F1 + F1_1a - F2 = F2 + F1_2a - F(0:2,i) = F(0:2,i) + F2_1a - F(0:2,i+1) = F(0:2,i+1) + F2_2a - G2(0:2,i) = F2_1a - G1(0:2,i+1) = F2_2a - end if + TPMInteractionFW0 = 1 + U = U + Ua + Ua = 0.25d+00 * Ua + U1 = U1 + Ua + U2 = U2 + Ua + UU(i) = UU(i) + Ua + UU(i+1) = UU(i+1) + Ua + F1 = F1 + F1_1a + F2 = F2 + F1_2a + F(0:2,i) = F(0:2,i) + F2_1a + F(0:2,i+1) = F(0:2,i+1) + F2_2a + G2(0:2,i) = F2_1a + G1(0:2,i+1) = F2_2a + end if end do end function TPMInteractionFW0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/lib/mesont/TubePotMono.f90 b/lib/mesont/TubePotMono.f90 index b6ff0f6d15..c97131ab11 100644 --- a/lib/mesont/TubePotMono.f90 +++ b/lib/mesont/TubePotMono.f90 @@ -1677,9 +1677,9 @@ contains !********************************************************************** call fdate( PDate ) write ( unit = TPMUnitID, fmt = '(a,a)' ) 'DATE ', PDate write ( unit = TPMUnitID, fmt = '(a,i3,a,i3,a)' ) & - 'Tabulated data of the tubular potential for (', ChiIndM, ',', ChiIndN, ') CNTs' + 'Tabulated data of the tubular potential for (', ChiIndM, ',', ChiIndN, ') CNTs' write ( unit = TPMUnitID, fmt = '(a)' ) & - 'A. N. Volkov, L. V. Zhigilei, J. Phys. Chem. C 114, 5513-5531, 2010. doi: 10.1021/jp906142h' + 'A. N. Volkov, L. V. Zhigilei, J. Phys. Chem. C 114, 5513-5531, 2010. doi: 10.1021/jp906142h' end if call TPMSSTPInit ()