ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

apply whitespace checking/fixing also to free-format Fortran files

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-08-25 23:46:48 -04:00
parent 45e599cb33
commit 61855c5058
14 changed files with 564 additions and 563 deletions
Download Patch File Download Diff File Expand all files Collapse all files

182
lib/mesont/CNTPot.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module CNTPot !*************************************************************************************
!
@ -21,34 +21,34 @@ module CNTPot !*****************************************************************
!
! Carbon nanotubes internal potentials:
! CNTSTRH0, harmonic stretching potential of type 0 with constant Young's modulus
! CNTSTRH1, harmonic stretching potential of type 1 with variable Young's modulus
! CNTSTRH1, harmonic stretching potential of type 1 with variable Young's modulus
! CNTSTRNH0, non-harmonic stretching with fracture potential of type 0
! CNTSTRNH1, non-harmonic stretching with fracture potential of type 1
! CNTBNDH, harmonic bending potential
! CNTBNDHB, harmonic bending-buckling potential
! CNTBNDHBF, harmonic bending-buckling potential with fracture
! CNTBNDHBF, harmonic bending-buckling potential with fracture
! CNTTRS, torsion potential
! CNTBRT, breathing potential
!
! The functional form and force constants of harmonic stretching, bending and
! torsion potentials are taken from:
! torsion potentials are taken from:
! L.V. Zhigilei, Ch. Wei, D. Srivastava, Phys. Rev. B 71, 165417 (2005)
!
! The model of stress-strain curve for the non-harmonic potential with fracture
! The model of stress-strain curve for the non-harmonic potential with fracture
! is developed and parameterized using the following constants
! -- Young's modulus
! -- Young's modulus
! -- maximum linear strain (only for the NH potential of type 1)
! -- tensile strength (or fracture strain)
! -- strain at failure (or fracture strain)
! -- maximum strain.
! All these parameters are assumed to be independent of CNT radius or chriality type.
! In this model, the true strain at failure CNTSTREft and true tensile strength
! In this model, the true strain at failure CNTSTREft and true tensile strength
! CNTSTRSft are slightly different from the imposed values CNTSTREf and CNTSTRSf.
!
! The non-harmonic stretching potentials of types 0 and 1 are different from
! The non-harmonic stretching potentials of types 0 and 1 are different from
! each other by the functional form of the stress-strain curve
!
! Different parameterizations of CNTSTRH0, CNTSTRNH0 and CNTSTRNH1 potentials can be chosen,
! Different parameterizations of CNTSTRH0, CNTSTRNH0 and CNTSTRNH1 potentials can be chosen,
! see subroutine CNTSTRSetParameterization
!
!---------------------------------------------------------------------------------------------------
@ -75,47 +75,47 @@ implicit none
integer(c_int), parameter :: CNTPOT_BFRACTURE = 5
! Harmonic stretching model (constant Young's modulus)
integer(c_int), parameter :: CNTSTRMODEL_H0 = 0
integer(c_int), parameter :: CNTSTRMODEL_H0 = 0
! Harmonic stretching model (Young's modulus depends on radius)
integer(c_int), parameter :: CNTSTRMODEL_H1 = 1
integer(c_int), parameter :: CNTSTRMODEL_H1 = 1
! Non-harmonic stretching with fracture, potential of type 0
integer(c_int), parameter :: CNTSTRMODEL_NH0F = 2
integer(c_int), parameter :: CNTSTRMODEL_NH0F = 2
! Non-harmonic stretching without fracture, potential of type 1
integer(c_int), parameter :: CNTSTRMODEL_NH1 = 3
integer(c_int), parameter :: CNTSTRMODEL_NH1 = 3
! Non-harmonic stretching with fracture, potential of type 1
integer(c_int), parameter :: CNTSTRMODEL_NH1F = 4
! Harmonic stretching model + axial buckling
integer(c_int), parameter :: CNTSTRMODEL_H1B = 5
integer(c_int), parameter :: CNTSTRMODEL_NH1F = 4
! Harmonic stretching model + axial buckling
integer(c_int), parameter :: CNTSTRMODEL_H1B = 5
! Harmonic stretching model + axial buckling + hysteresis
integer(c_int), parameter :: CNTSTRMODEL_H1BH = 6
integer(c_int), parameter :: CNTSTRMODEL_H1BH = 6
integer(c_int), parameter :: CNTBNDMODEL_H = 0 ! Harmonic bending model
integer(c_int), parameter :: CNTBNDMODEL_HB = 1 ! Harmonic bending - buckling model
integer(c_int), parameter :: CNTBNDMODEL_HBF = 2 ! Harmonic bending - buckling - fracture model
integer(c_int), parameter :: CNTBNDMODEL_HBH = 3 ! Harmonic bending - buckling + Hysteresis
integer(c_int), parameter :: CNTPOTNMAX = 4000 ! Maximum number of points in the interpolation tables
!---------------------------------------------------------------------------------------------------
! Parameters of potentials
!---------------------------------------------------------------------------------------------------
! Stretching potential
! Type of the bending model
integer(c_int) :: CNTSTRModel = CNTSTRMODEL_H1
! Type of parameterization
integer(c_int) :: CNTSTRParams = 0
integer(c_int) :: CNTSTRParams = 0
! Type of dependence of the Young's modulus on tube radius
integer(c_int) :: CNTSTRYMT = 0
integer(c_int) :: CNTSTRYMT = 0
! Parameters of non-harmonic potential and fracture model
real(c_double) :: CNTSTRR0 = 6.8d+00 ! Reference radius of nanotubes (A)
! (this parameter is not used for the model
! parametrization, but only for calculation of the
real(c_double) :: CNTSTRR0 = 6.8d+00 ! Reference radius of nanotubes (A)
! (this parameter is not used for the model
! parametrization, but only for calculation of the
! force constant in eV/A)
real(c_double) :: CNTSTRD0 = 3.4d+00 ! CNT wall thickness (A)
real(c_double) :: CNTSTREmin = -0.4d+00 ! Minimum strain in tabulated potential
real(c_double) :: CNTSTRD0 = 3.4d+00 ! CNT wall thickness (A)
real(c_double) :: CNTSTREmin = -0.4d+00 ! Minimum strain in tabulated potential
real(c_double) :: CNTSTREmax = 0.13d+00 ! Maximum strain in tabulated potential.
! Simultaneously, U=0 if E> CNTSTREmax
real(c_double) :: CNTSTREl = 5.0d-02 ! Maximum linear strain
@ -129,17 +129,17 @@ implicit none
real(c_double) :: CNTSTRSi ! Maximum stress (not used in the model) (Pa)
real(c_double) :: CNTSTRDf ! dF/dE at failure
real(c_double) :: CNTSTRAA, CNTSTRBB !
real(c_double) :: CNTSTRAA, CNTSTRBB !
real(c_double) :: CNTSTRAAA, CNTSTRBBB ! Auxiliary constants
real(c_double) :: CNTSTRUl, CNTSTRUf !
real(c_double) :: CNTSTRUl, CNTSTRUf !
! Axial buckling - hysteresis approach
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) :: 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
! Bending potential
integer(c_int) :: CNTBNDModel = CNTBNDMODEL_H ! Type of the bending model
! Buckling model parameters
real(c_double) :: CNTBNDN = 1.0d+00 ! Buckling exponent
@ -149,9 +149,9 @@ implicit none
real(c_double) :: CNTBNDTF = M_PI * 120.0d+00 / 180.0d+00 ! Fracture buckling angle (rad)
real(c_double) :: CNTBNDN1
real(c_double) :: CNTBNDC2
contains !******************************************************************************************
!---------------------------------------------------------------------------------------------------
! Stretching potential
!---------------------------------------------------------------------------------------------------
@ -167,7 +167,7 @@ contains !**********************************************************************
!-------------------------------------------------------------------------------------------
integer(c_int), intent(in) :: PType
!-------------------------------------------------------------------------------------------
select case ( PType )
select case ( PType )
case ( 0 ) ! This parametrization is based on averaged exp. data of Ref. [1]
CNTSTRR0 = 6.8d+00 ! Ref. [1]
CNTSTRD0 = 3.4d+00 ! Ref. [1]
@ -180,7 +180,7 @@ contains !**********************************************************************
case ( 1 ) ! This parameterization is taken from Ref. [2] for (10,10) CNTs.
! These values are obtained in MD simulations with REBO potential.
! Values of Young's modulus, tensile strength and stress here
! are close to those obtained in Ref. [3] for pristine (defectless)
! are close to those obtained in Ref. [3] for pristine (defectless)
! (5,5) CNT in semi-empirical QM calculations based on PM3 model
CNTSTRR0 = 6.785d+00 ! Calculated with the usual formula for (10,10) CNT
CNTSTRD0 = 3.35d+00 ! Ref. [2]
@ -190,7 +190,7 @@ contains !**********************************************************************
CNTSTREf = 27.9d-02 ! Corresponds to maximum strain in Ref. [2]
CNTSTRS0 = 1.031e+12 ! Ref. [2]
CNTSTRSf = 148.5d+09 ! Corresponds to tensile strength in Ref. [2]
case ( 2 ) ! This parametrization is taken from Ref. [3] for (5,5) CNTs
case ( 2 ) ! This parametrization is taken from Ref. [3] for (5,5) CNTs
! with one atom vacancy defect obtained with the semi-empirical QM PM3 model
CNTSTRR0 = 3.43d+00 ! Ref. [3]
CNTSTRD0 = 3.4d+00 ! Ref. [3]
@ -246,14 +246,14 @@ contains !**********************************************************************
!
! Stretching without fracture, harmonic potential, with axial buckling without hysteresis
!
integer(c_int) function CNTSTRH1BCalc ( U, dUdL, L, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function CNTSTRH1BCalc ( U, dUdL, L, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Young's modulus depends on R, see [4].
! Axial buckling without hysteresis.
!-------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: U, dUdL
real(c_double), intent(in) :: L, R0, L0
real(c_double) :: E, K, Kbcl, dUbcl, d, ud
real(c_double) :: E, K, Kbcl, dUbcl, d, ud
!-------------------------------------------------------------------------------------------
E = ( L - L0 ) / L0
K = 86.64d+00 + 100.56d+00 * R0
@ -266,7 +266,7 @@ contains !**********************************************************************
dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc
U = Kbcl * E + dUbcl
dUdL = Kbcl / L0
CNTSTRH1BCalc = CNTPOT_STRETCHING
CNTSTRH1BCalc = CNTPOT_STRETCHING
else ! Return to harmonic potential
d = -0.0142794
dUdL = K * ( d + E - CNTSTREc2 )
@ -274,7 +274,7 @@ contains !**********************************************************************
Ud = 0.5d+00 * L0 * K * d * d
U = 0.5d+00 * L0 * (d+E-CNTSTREc2) * dUdL + dUbcl - Ud
CNTSTRH1BCalc = CNTPOT_STRETCHING
end if
end if
end function CNTSTRH1BCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
@ -290,7 +290,7 @@ contains !**********************************************************************
!-------------------------------------------------------------------------------------------
real(c_double) :: E, K, dUbcl, Ebcl, Kbcl, Edu
real(c_double) :: C, DE, t
!-------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------
E = ( L - L0 ) / L0
K = 86.64d+00 + 100.56d+00 * R0
Kbcl = -10.98d+00 * L0
@ -319,7 +319,7 @@ contains !**********************************************************************
dUdL = Kbcl / L0
CNTSTRH1BHCalc = CNTPOT_SBUCKLING
Ebuc = 0.5d+00 * L0 * K * CNTSTREc1 * CNTSTREc1 - Kbcl * CNTSTREc1 - dUbcl
else ! Already buckled
else ! Already buckled
dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc
U = Kbcl * E + dUbcl
dUdL = Kbcl / L0
@ -331,7 +331,7 @@ contains !**********************************************************************
U = 0.5d+00 * L0 * E * dUdL
CNTSTRH1BHCalc = CNTPOT_STRETCHING
Ebuc = 0.0d+00
end if
end if
end function CNTSTRH1BHCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
@ -345,26 +345,26 @@ contains !**********************************************************************
!-------------------------------------------------------------------------------------------
E = ( L - L0 ) / L0
if ( E < CNTSTREf ) then
dUdL = ( CNTSTRAA - CNTSTRBB * E ) * E
dUdL = ( CNTSTRAA - CNTSTRBB * E ) * E
U = ( CNTSTRAAA - CNTSTRBBB * E ) * E * E
CNTSTRNH0FCalc = CNTPOT_STRETCHING
else
else
dUdL = 0.0d+00
U = 0.0d+00
CNTSTRNH0FCalc = CNTPOT_SFRACTURE
end if
U = L0 * R0 * U
dUdL = R0 * dUdL
dUdL = R0 * dUdL
end function CNTSTRNH0FCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CNTSTRNH0Init () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) :: S
!-------------------------------------------------------------------------------------------
S = M_2PI * CNTSTRD0 * 1.0e-20 / K_MDFU
CNTSTRSl = CNTSTRS0 * CNTSTREl
CNTSTRF0 = CNTSTRS0 * S
CNTSTRFl = CNTSTRSl * S
CNTSTRFf = CNTSTRSf * S
CNTSTRF0 = CNTSTRS0 * S
CNTSTRFl = CNTSTRSl * S
CNTSTRFf = CNTSTRSf * S
CNTSTRAA = CNTSTRF0
CNTSTRBB = ( CNTSTRF0 * CNTSTREf - CNTSTRFf ) / ( CNTSTREf * CNTSTREf )
CNTSTRAAA= CNTSTRAA / 2.0d+00
@ -375,11 +375,11 @@ contains !**********************************************************************
CNTSTRSi = 0.0d+00
CNTSTRDf = 0.0d+00
end subroutine CNTSTRNH0Init !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Stretching without fracture, non-harmonic potential of type 1
!
integer(c_int) function CNTSTRNH1Calc ( U, dUdL, L, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: U, dUdL
real(c_double), intent(in) :: L, R0, L0
@ -390,7 +390,7 @@ contains !**********************************************************************
dUdL = CNTSTRF0 * E
U = 0.5d+00 * E * dUdL
CNTSTRNH1Calc = CNTPOT_STRETCHING
else
else
DE = E - CNTSTREl
C = 1.0 + CNTSTRBB * DE
dUdL = CNTSTRFl + CNTSTRAA * ( 1.0d+00 - 1.0d+00 / C )
@ -398,13 +398,13 @@ contains !**********************************************************************
end if
CNTSTRNH1Calc = CNTPOT_STRETCHING
U = L0 * R0 * U
dUdL = R0 * dUdL
dUdL = R0 * dUdL
end function CNTSTRNH1Calc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Stretching with fracture, non-harmonic potential of type 1
!
integer(c_int) function CNTSTRNH1FCalc ( U, dUdL, L, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: U, dUdL
real(c_double), intent(in) :: L, R0, L0
@ -421,24 +421,24 @@ contains !**********************************************************************
dUdL = CNTSTRFl + CNTSTRAA * ( 1.0d+00 - 1.0d+00 / C )
U = CNTSTRUl + CNTSTRAAA * DE - CNTSTRBBB * dlog ( C )
CNTSTRNH1FCalc = CNTPOT_STRETCHING
else
else
dUdL = 0.0d+00
U = 0.0d+00
CNTSTRNH1FCalc = CNTPOT_SFRACTURE
end if
U = L0 * R0 * U
dUdL = R0 * dUdL
dUdL = R0 * dUdL
end function CNTSTRNH1FCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CNTSTRNH1Init () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) :: S, C, E, t
integer(c_int) :: i, CaseID
!-------------------------------------------------------------------------------------------
S = M_2PI * CNTSTRD0 * 1.0e-20 / K_MDFU
CNTSTRSl = CNTSTRS0 * CNTSTREl
CNTSTRF0 = CNTSTRS0 * S
CNTSTRFl = CNTSTRSl * S
CNTSTRFf = CNTSTRSf * S
CNTSTRF0 = CNTSTRS0 * S
CNTSTRFl = CNTSTRSl * S
CNTSTRFf = CNTSTRSf * S
CNTSTRAA = ( CNTSTRFf - CNTSTRFl ) * ( CNTSTREf * CNTSTRF0 - CNTSTRFl ) / ( CNTSTREf * CNTSTRF0 - CNTSTRFf )
CNTSTRBB = CNTSTRF0 / CNTSTRAA
CNTSTRAAA= CNTSTRFl + CNTSTRAA
@ -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
!
@ -477,7 +477,7 @@ contains !**********************************************************************
CNTSTRCalc = CNTSTRH1BHCalc ( U, dUdL, L, R0, L0, ABF, Ebuc )
end select
end function CNTSTRCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CNTSTRInit ( STRModel, STRParams, YMType, Rref ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int), intent(in) :: STRModel, STRParams, YMType
real(c_double), intent(in) :: Rref
@ -498,13 +498,13 @@ contains !**********************************************************************
else
call CNTSTRNH1Init ()
end if
end if
end if
end subroutine CNTSTRInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Bending potentials
!---------------------------------------------------------------------------------------------------
subroutine BendingGradients ( K, G0, G1, G2, R0, R1, R2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(inout) :: K
real(c_double), dimension(0:2), intent(inout) :: G0, G1, G2
@ -525,7 +525,7 @@ contains !**********************************************************************
G2 = G2 / L2
G1 = - ( G0 + G2 )
end subroutine BendingGradients !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function CNTBNDHCalc ( U, dUdC, C, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Bending model of type 0:Harmonic bending potential.
!-------------------------------------------------------------------------------------------
@ -554,10 +554,10 @@ contains !**********************************************************************
! Check the condition for buckling
if ( C2 .ge. CNTBNDC2 ) then ! Buckling takes place
Theta= M_PI - acos ( C )
Kbnd = 63.8d+00 * R0**2.93d+00
Kbnd = 63.8d+00 * R0**2.93d+00
Kbcl = CNTBNDB * Kbnd / CNTBNDR
DUbcl= Kbnd * ( CNTBNDB * ( M_PI - 2.0d+00 * atan ( 2.0 * CNTBNDR / L0 ) ) - 0.5d+00 * L0 / CNTBNDR ) &
/ CNTBNDR
DUbcl= Kbnd * ( CNTBNDB * ( M_PI - 2.0d+00 * atan ( 2.0 * CNTBNDR / L0 ) ) - 0.5d+00 * L0 / CNTBNDR ) &
/ CNTBNDR
U = Kbcl * abs( Theta )**CNTBNDN - DUbcl
dUdC = Kbcl * CNTBNDN * abs( Theta )**CNTBNDN1 / sqrt ( 1.0d+00 - C * C )
CNTBNDHBCalc = CNTPOT_BBUCKLING
@ -568,7 +568,7 @@ contains !**********************************************************************
CNTBNDHBCalc = CNTPOT_BENDING
end if
end function CNTBNDHBCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function CNTBNDHBFCalc ( U, dUdC, C, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: U, dUdC
real(c_double), intent(in) :: C, R0, L0
@ -586,10 +586,10 @@ contains !**********************************************************************
dUdC = 0.0d+00
CNTBNDHBFCalc = CNTPOT_BFRACTURE
else
Kbnd = 63.8d+00 * R0**2.93d+00
Kbnd = 63.8d+00 * R0**2.93d+00
Kbcl = CNTBNDB * Kbnd / CNTBNDR
DUbcl= Kbnd * ( CNTBNDB * ( M_PI - 2.0d+00 * atan ( 2.0 * CNTBNDR / L0 ) ) - &
0.5d+00 * L0 / CNTBNDR ) / CNTBNDR
0.5d+00 * L0 / CNTBNDR ) / CNTBNDR
U = Kbcl * abs ( Theta )**CNTBNDN - DUbcl
dUdC = Kbcl * CNTBNDN * abs ( Theta )**CNTBNDN1 / sqrt ( 1.0d+00 - C * C )
CNTBNDHBFCalc = CNTPOT_BBUCKLING
@ -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.
!-------------------------------------------------------------------------------------------
@ -611,16 +611,16 @@ contains !**********************************************************************
real(c_double) :: E1, E2, C2, Kbnd, Kbcl,Theta,DUbcl, Ubcl, Cmin,Rmax
!-------------------------------------------------------------------------------------------
Rmax = 340.0d+00
Cmin = 1.0/(Rmax*Rmax)
Cmin = 1.0/(Rmax*Rmax)
E1 = 1.0d+00 - C
E2 = 1.0d+00 + C
E2 = 1.0d+00 + C
! Calculate the square of curvature
C2 = 4.0d+00 * E2 / ( L0 * L0 * E1 )
Theta = M_PI - acos ( C )
if ( C2 .lt. Cmin ) then ! Harmonic bending
if ( C2 .lt. Cmin ) then ! 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 )
dUdC = 2.0d+00 * Kbnd / ( E1 * E1 )
CNTBNDHBHCalc = CNTPOT_BENDING
Ebuc = 0.0
else if ( C2 .ge. Cmin .and. C2 .lt. CNTBNDC2 ) then ! Potential depends on buckling flag of a node
@ -640,7 +640,7 @@ contains !**********************************************************************
dUdC = Kbcl * CNTBNDN * abs( Theta )**CNTBNDN1 / sqrt ( 1.0d+00 - C * C )
Ebuc = 0.0d+00
CNTBNDHBHCalc = CNTPOT_BBUCKLING
end if
end if
else ! Greater than buckling critical point
if ( BBF .eq. 1 ) then ! Already buckled
Theta= M_PI - acos ( C )
@ -653,7 +653,7 @@ contains !**********************************************************************
Ebuc = 0.0d00
CNTBNDHBHCalc = CNTPOT_BBUCKLING
else ! Newly buckled
Theta= M_PI - acos ( C )
Theta= M_PI - acos ( C )
Kbnd = 63.8d+00 * R0**2.93d+00
Kbcl = CNTBNDB * Kbnd / CNTBNDR
DUbcl= 2.0d+00*Kbnd * &
@ -666,11 +666,11 @@ contains !**********************************************************************
end if
end if
end function CNTBNDHBHCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! General
!
integer(c_int) function CNTBNDCalc ( U, dUdC, C, R0, L0, BBF, Ebuc ) !!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: U, dUdC, Ebuc
real(c_double), intent(in) :: C, R0, L0
@ -688,7 +688,7 @@ contains !**********************************************************************
CNTBNDCalc = CNTBNDHBHCalc ( U, dUdC, C, R0, L0, BBF, Ebuc )
end select
end function CNTBNDCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CNTBNDInit ( BNDModel ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int), intent(in) :: BNDModel
real(c_double) :: A, E
@ -710,5 +710,5 @@ contains !**********************************************************************
call CNTSTRInit ( STRModel, STRParams, YMType, Rref )
call CNTBNDInit ( BNDModel )
end subroutine InitCNTPotModule !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module CNTPot !*********************************************************************************

100
lib/mesont/ExportCNT.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Maxim Shugaev (UVA), mvs9t@virginia.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module ExportCNT !**********************************************************************************
@ -28,122 +28,122 @@ contains
bind(c, name = "mesont_lib_InitCNTPotModule")
integer(c_int), intent(in) :: STRModel, STRParams, YMType, BNDModel
real(c_double), intent(in) :: Rref
call InitCNTPotModule(STRModel, STRParams, YMType, BNDModel, Rref)
endsubroutine
subroutine TPBInit_() &
bind(c, name = "mesont_lib_TPBInit")
call TPBInit()
endsubroutine
subroutine TPMInit_(M, N) &
bind(c, name = "mesont_lib_TPMInit")
integer(c_int), intent(in) :: M, N
call TPMInit(M, N)
endsubroutine
subroutine SetTablePath_(TPMFile_, N) &
bind(c, name = "mesont_lib_SetTablePath")
integer(c_int), intent(in) :: N
character(c_char), intent(in), dimension(N) :: TPMFile_
integer :: i
do i = 1, len(TPMFile)
if (i <= N) then
TPMFile(i:i) = TPMFile_(i)
else
else
TPMFile(i:i) = ' '
endif
enddo
endsubroutine
function get_R_ () &
bind(c, name = "mesont_lib_get_R")
real(c_double) :: get_R_
get_R_ = TPMR1
return
endfunction
subroutine TubeStretchingForceField_(U1, U2, F1, F2, S1, S2, X1, X2, R12, L12) &
bind(c, name = "mesont_lib_TubeStretchingForceField")
! Interaction energies associated with nodes X1 and X2
real(c_double), intent(inout) :: U1, U2
! Interaction energies associated with nodes X1 and X2
real(c_double), intent(inout) :: U1, U2
! Forces exerted on nodes X1 and X2
real(c_double), intent(inout), dimension(0:2) :: F1, F2
real(c_double), intent(inout), dimension(0:2) :: F1, F2
! Contributions of nodes X1 and X2 to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2
! Coordinates of the segment nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2
! Coordinates of the segment nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2
! Radius of a nanotube the segment (X1,X2) belongs to
real(c_double), intent(in) :: R12
real(c_double), intent(in) :: R12
! Equilibrium length of segment (X1,X2)
real(c_double), intent(in) :: L12
real(c_double), intent(in) :: L12
call TubeStretchingForceField(U1, U2, F1, F2, S1, S2, X1, X2, R12, L12)
endsubroutine
subroutine TubeBendingForceField_(U1, U2, U3, F1, F2, F3, S1, S2, S3, X1, X2, X3, R123, L123, BBF2) &
bind(c, name = "mesont_lib_TubeBendingForceField")
! Interaction energies associated with nodes X1, X2, and X3
real(c_double), intent(inout) :: U1, U2, U3
real(c_double), intent(inout) :: U1, U2, U3
! Forces exerted on nodes X1, X2, and X3
real(c_double), intent(inout), dimension(0:2) :: F1, F2, F3
real(c_double), intent(inout), dimension(0:2) :: F1, F2, F3
! Contributions of nodes X1, X2, and X3 to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2, S3
! Coordinates of nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2, X3
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2, S3
! Coordinates of nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2, X3
! Radius of nanotube the segment (X1,X2) belongs to
real(c_double), intent(in) :: R123
real(c_double), intent(in) :: R123
! Equilibrium length of segment (X1,X2) and (X2,X3) (It is assumed to be the same for both segments)
real(c_double), intent(in) :: L123
real(c_double), intent(in) :: L123
integer(c_int), intent(inout) :: BBF2
call TubeBendingForceField(U1, U2, U3, F1, F2, F3, S1, S2, S3, X1, X2, X3, R123, L123, BBF2 )
endsubroutine
subroutine SegmentTubeForceField_(U1,U2,U,F1,F2,F,Fe,S1,S2,S,Se,X1,X2,R12,N,X,Xe,BBF,R,E1,E2,Ee,TPMType)&
bind(c, name = "mesont_lib_SegmentTubeForceField")
! Number of nodes in array X
integer(c_int), intent(in) :: N
integer(c_int), intent(in) :: N
! Interaction energies associated with nodes X1 and X2
real(c_double), intent(inout) :: U1, U2
real(c_double), intent(inout) :: U1, U2
! Interaction energies associated with nodes X
real(c_double), intent(inout), dimension(0:N-1) :: U
real(c_double), intent(inout), dimension(0:N-1) :: U
! Forces exerted on nodes X1 and X2
real(c_double), intent(inout), dimension(0:2) :: F1, F2
real(c_double), intent(inout), dimension(0:2) :: F1, F2
! Forces exerted on nodes X
real(c_double), intent(inout), dimension(0:2,0:N-1) :: F
real(c_double), intent(inout), dimension(0:2,0:N-1) :: F
! Force exerted on node Xe (can be updated only if Ee > 0)
real(c_double), intent(inout), dimension(0:2) :: Fe
real(c_double), intent(inout), dimension(0:2) :: Fe
! Contributions of nodes X1 and X2 to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2
! Contributions of nodes X to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2,0:N-1) :: S
real(c_double), intent(inout), dimension(0:2,0:2,0:N-1) :: S
! Contributions of node Xe to the virial stress tensor (can be updated only if Ee > 0)
real(c_double), intent(inout), dimension(0:2,0:2) :: Se
! Coordinates of the segment nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2
real(c_double), intent(inout), dimension(0:2,0:2) :: Se
! Coordinates of the segment nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2
! Radius of nanotube the segment (X1,X2) belongs to
real(c_double), intent(in) :: R12
real(c_double), intent(in) :: R12
! Coordinates of the nanotube nodes
real(c_double), intent(in), dimension(0:2,0:N-1) :: X
real(c_double), intent(in), dimension(0:2,0:N-1) :: X
! Additional node of the extended chain if Ee > 0
real(c_double), intent(in), dimension(0:2) :: Xe
! Bending buckling flags (BBF(i) = 1 in a case of buckling in node i)
integer(c_int), intent(in), dimension(0:N-1) :: BBF
real(c_double), intent(in), dimension(0:2) :: Xe
! Bending buckling flags (BBF(i) = 1 in a case of buckling in node i)
integer(c_int), intent(in), dimension(0:N-1) :: BBF
! Radius of nanotube X
real(c_double), intent(in) :: R
! E1 = 1 if the chain node 0 is a CNT end; E2 = 1 if the chain node N-1 is a CNT end;
integer(c_int), intent(in) :: E1, E2
real(c_double), intent(in) :: R
! E1 = 1 if the chain node 0 is a CNT end; E2 = 1 if the chain node N-1 is a CNT end;
integer(c_int), intent(in) :: E1, E2
! Parameter defining the type of the extended chain (0,1,2)
integer(c_int), intent(in) :: Ee
integer(c_int), intent(in) :: Ee
! Type of the tubular potential (0 or 1)
integer(c_int), intent(in) :: TPMType
integer(c_int), intent(in) :: TPMType
call SegmentTubeForceField(U1, U2, U, F1, F2, F, Fe, S1, S2, S, Se, X1, X2, R12, N, X, Xe, BBF, R, E1, E2, Ee, TPMType)
endsubroutine

20
lib/mesont/LinFun2.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module LinFun2 !************************************************************************************
!
@ -26,7 +26,7 @@ module LinFun2 !****************************************************************
!***************************************************************************************************
use iso_c_binding, only : c_int, c_double, c_char
implicit none
contains !******************************************************************************************
real(c_double) function CalcLinFun1_0 ( i, X, N, P, F ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@ -42,7 +42,7 @@ contains !**********************************************************************
A = 1.0d+00 - A0
CalcLinFun1_0 = A0 * F(i1) + A * F(i)
end function CalcLinFun1_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CalcLinFun1_1 ( S, Sx1, i, X, N, P, F, Fx ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: S, Sx1
integer(c_int), intent(in) :: i, N
@ -78,7 +78,7 @@ contains !**********************************************************************
G0 = B0 * F(i1,j1) + B * F(i1,j)
CalcLinFun2_0 = A0 * G0 + A * G
end function CalcLinFun2_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CalcLinFun2_1 ( S, Sx1, Sy1, i, j, X, Y, N1, N2, P1, P2, F, Fx, Fy ) !!!!!!!!!!!!
real(c_double), intent(out) :: S, Sx1, Sy1
integer(c_int), intent(in) :: i, j, N1, N2
@ -95,19 +95,19 @@ contains !**********************************************************************
A = 1.0d+00 - A0
B0 = ( P2(j) - Y ) / ( P2(j) - P2(j1) )
B = 1.0d+00 - B0
G = B0 * F(i,j1) + B * F(i,j)
G0 = B0 * F(i1,j1) + B * F(i1,j)
S = A0 * G0 + A * G
G = B0 * Fx(i,j1) + B * Fx(i,j)
G0 = B0 * Fx(i1,j1) + B * Fx(i1,j)
Sx1 = A0 * G0 + A * G
G = B0 * Fy(i,j1) + B * Fy(i,j)
G0 = B0 * Fy(i1,j1) + B * Fy(i1,j)
Sy1 = A0 * G0 + A * G
end subroutine CalcLinFun2_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module LinFun2 !********************************************************************************

30
lib/mesont/Spline1.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module Spline1 !************************************************************************************
!
@ -64,7 +64,7 @@ contains !**********************************************************************
!-------------------------------------------------------------------------------------------
X(0) = F(0) / K1(0)
F(0) = - K2(0) / K1(0)
do i = 1, N - 1
do i = 1, N - 1
D = - ( K1(i) + F(i-1) * K0(i) )
X(i) = ( K0(i) * X(i-1) - F(i) ) / D
F(i) = K2(i) / D
@ -85,8 +85,8 @@ contains !**********************************************************************
K0(i) = P(i) - P(i-1)
K1(i) = ( F(i) - F(i-1) ) / K0(i)
end do
select case ( CL )
case (1)
select case ( CL )
case (1)
K1(0) = 2.0d+00 / 3.0d+00
K2(0) = 1.0d+00 / 3.0d+00
D (0) = 2 * ( K1(1) - M(0) ) / K0(1)
@ -98,14 +98,14 @@ contains !**********************************************************************
K1(0) = 1.0d+00
K2(0) = 0.0d+00
D(0) = 0.0d+00
end select
end select
Z = K1(N-1)
do i = 1, N - 2
do i = 1, N - 2
D(i) = 6.0d+00 * ( K1(i+1) - K1(i) )
K2(i) = K0(i+1)
K1(i) = 2.0d+00 * ( K2(i) + K0(i) )
end do
select case ( CR )
select case ( CR )
case (1)
D(N-1) = 2.0d+00 * ( M(N-1) - Z ) / K0(N-1)
K1(N-1) = 2.0d+00 / 3.0d+00
@ -118,14 +118,14 @@ contains !**********************************************************************
K1(N-1) = 1.0d+00
K0(N-1) = 0.0d+00
D(N-1) = 0.0d+00
end select
end select
call sprogonka3 ( N, K0, K1, K2, D, M )
end subroutine CreateSpline1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function CalcSpline1_0 ( i, X, N, P, F, M ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int), intent(in) :: i, N
real(c_double), intent(in) :: X
real(c_double), dimension(0:N-1), intent(in) :: P, F, M
real(c_double), dimension(0:N-1), intent(in) :: P, F, M
integer(c_int) :: j
real(c_double) :: HL, HR, H, H6, H26, HR2, HL2, HRH, HLH
!-------------------------------------------------------------------------------------------
@ -147,7 +147,7 @@ contains !**********************************************************************
real(c_double), intent(out) :: S, S1
integer(c_int), intent(in) :: i, N
real(c_double), intent(in) :: X
real(c_double), dimension(0:N-1), intent(in) :: P, F, M
real(c_double), dimension(0:N-1), intent(in) :: P, F, M
integer(c_int) :: j
real(c_double) :: HL, HR, H, H6, H26, HR2, HL2, HRH, HLH
!-------------------------------------------------------------------------------------------
@ -169,7 +169,7 @@ contains !**********************************************************************
real(c_double), intent(out) :: S, S1, S2
integer(c_int), intent(in) :: i, N
real(c_double), intent(in) :: X
real(c_double), dimension(0:N-1), intent(in) :: P, F, M
real(c_double), dimension(0:N-1), intent(in) :: P, F, M
integer(c_int) :: j
real(c_double) :: HL, HR, H, H6, H26, HR2, HL2, HRH, HLH
!-------------------------------------------------------------------------------------------
@ -185,7 +185,7 @@ contains !**********************************************************************
HRH = HR / H
S = ( M(j) * HR2 * HRH + M(i) * HL2 * HLH ) / 6.0d+00 + ( F(j) - M(j) * H26 ) * HRH + ( F(i) - M(i) * H26 ) * HLH
S1 = ( ( M(i) * HL2 - M(j) * HR2 ) / 2.0d+00 + F(i) - F(j) ) / H - H6 * ( M(i) - M(j) )
S2 = M(j) * HRH + M(i) * HLH
S2 = M(j) * HRH + M(i) * HLH
end subroutine CalcSpline1_2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module Spline1 !********************************************************************************

26
lib/mesont/Spline2.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module Spline2 !************************************************************************************
!
@ -79,7 +79,7 @@ contains !**********************************************************************
Fxx = 0.0d+00
Fyy = 0.0d+00
Fxxyy = 0.0d+00
do II = 0, N2A
FF(0:N1-1) = F(0:N1-1,II)
MM(0) = Fxx(0,II)
@ -87,7 +87,7 @@ contains !**********************************************************************
call CreateSpline1 ( CL, CR, N1, P1, FF, MM, DD, K0, K1, K2 )
Fxx(0:N1-1,II) = MM(0:N1-1)
end do
do II = N2A + 1, N2 - 1
FF(0:N1-N1A-1) = F(N1A:N1-1,II)
MM(0) = Fxx(N1A,II)
@ -103,7 +103,7 @@ contains !**********************************************************************
call CreateSpline1 ( CD, CU, N2A + 1, P2, FF, MM, DD, K0, K1, K2 )
Fyy(II,0:N2A) = MM(0:N2A)
end do
do II = N1A, N1 - 1
MM(0) = Fyy(II,0)
MM(N-1) = Fyy(II,N2-1)
@ -111,19 +111,19 @@ contains !**********************************************************************
call CreateSpline1 ( CD, CU, N2, P2, FF, MM, DD, K0, K1, K2 )
Fyy(II,0:N2-1) = MM(0:N2-1)
end do
FF(0:N1-1) = Fyy(0:N1-1,0)
call CreateSpline1 ( 3, 3, N1, P1, FF, MM, DD, K0, K1, K2 )
Fxxyy(0:N1-1,0) = MM(0:N1-1)
FF(0:N1A) = Fyy(0:N1A,N2A)
call CreateSpline1 ( 3, 3, N1A + 1, P1, FF, MM, DD, K0, K1, K2 )
Fxxyy(0:N1A,N2A) = MM(0:N1A)
FF(0:N1-N1A-1) = Fyy(N1A:N1-1,N2-1 )
call CreateSpline1 ( 3, 3, N1-N1A, P1, FF, MM, DD, K0, K1, K2 )
Fxxyy(N1A:N1-1,N2-1) = MM(0:N1-N1A-1)
do II = 1, N1A
MM(0) = Fxxyy(II,0)
MM(N2A) = Fxxyy(II,N2A)
@ -131,7 +131,7 @@ contains !**********************************************************************
call CreateSpline1 ( 2 , 2, N2A + 1, P2, FF, MM, DD, K0, K1, K2 )
Fxxyy(II,0:N2A) = MM(0:N2A)
end do
do II = N1A + 1, N1 - 2
MM(0) = Fxxyy(II,0)
MM(N-1) = Fxxyy(II,N2-1)
@ -139,9 +139,9 @@ contains !**********************************************************************
call CreateSpline1 ( 2 , 2, N2, P2, FF, MM, DD, K0, K1, K2 )
Fxxyy(II,0:N2-1) = MM(0:N2-1)
end do
end subroutine CreateSpline2Ext !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function CalcSpline2_0 ( i, j, X, Y, N1, N2, P1, P2, F, Fxx, Fyy, Fxxyy ) !!!
integer(c_int), intent(in) :: i, j, N1, N2
real(c_double), intent(in) :: X, Y
@ -182,5 +182,5 @@ contains !**********************************************************************
call ValueSpline1_1 ( S, Sx1, X, P1(i), P1(i1), Gy_0, Gy_1,Gxxy_0, Gxxy_1, P1(i) - P1(i1) )
Sy1 = ValueSpline1_0 ( X, P1(i), P1(i1), Gyy_0, Gyy_1,Gxxyy_0, Gxxyy_1, P1(i) - P1(i1) )
end subroutine CalcSpline2_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module Spline2 !********************************************************************************

140
lib/mesont/TPMForceField.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module TPMForceField !******************************************************************************
!
@ -32,32 +32,32 @@ use iso_c_binding, only : c_int, c_double, c_char
implicit none
contains !******************************************************************************************
subroutine TubeStretchingForceField ( U1, U2, F1, F2, S1, S2, X1, X2, R12, L12 ) !!!!!!!!!!!
! Interaction energies associated with nodes X1 and X2
real(c_double), intent(inout) :: U1, U2
real(c_double), intent(inout) :: U1, U2
! Forces exerted on nodes X1 and X2
real(c_double), intent(inout), dimension(0:2) :: F1, F2
real(c_double), intent(inout), dimension(0:2) :: F1, F2
! Contributions of nodes X1 and X2 to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2
! Coordinates of the segment nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2
! Coordinates of the segment nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2
! Radius of a nanotube the segment (X1,X2) belongs to
real(c_double), intent(in) :: R12
real(c_double), intent(in) :: R12
! Equilibrium length of segment (X1,X2)
real(c_double), intent(in) :: L12
real(c_double), intent(in) :: L12
!-------------------------------------------------------------------------------------------
integer(c_int) :: ii, jj, Event
real(c_double) :: U, F, LL, S, Ubcl
real(c_double), dimension(0:2) :: DX, FF
!-------------------------------------------------------------------------------------------
DX = X2 - X1
DX = X2 - X1
LL = S_V3norm3 ( DX )
Event = CNTSTRCalc ( U, F, LL, R12, L12, 0, Ubcl )
U = U / 2.0d+00
FF = DX * F / LL
F1 = F1 + FF
U1 = U1 + U
@ -68,25 +68,25 @@ contains !**********************************************************************
do ii = 0, 2
do jj = 0, 2
S = - 0.5d+00 * DX(ii) * FF(jj)
S1(ii,jj) = S1(ii,jj) + S
S2(ii,jj) = S2(ii,jj) + S
S1(ii,jj) = S1(ii,jj) + S
S2(ii,jj) = S2(ii,jj) + S
end do
end do
end subroutine TubeStretchingForceField !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TubeBendingForceField ( U1, U2, U3, F1, F2, F3, S1, S2, S3, X1, X2, X3, R123, L123, BBF2 )
subroutine TubeBendingForceField ( U1, U2, U3, F1, F2, F3, S1, S2, S3, X1, X2, X3, R123, L123, BBF2 )
! Interaction energies associated with nodes X1, X2, and X3
real(c_double), intent(inout) :: U1, U2, U3
real(c_double), intent(inout) :: U1, U2, U3
! Forces exerted on nodes X1, X2, and X3
real(c_double), intent(inout), dimension(0:2) :: F1, F2, F3
real(c_double), intent(inout), dimension(0:2) :: F1, F2, F3
! Contributions of nodes X1, X2, and X3 to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2, S3
! Coordinates of nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2, X3
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2, S3
! Coordinates of nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2, X3
! Radius of nanotube the segment (X1,X2) belongs to
real(c_double), intent(in) :: R123
real(c_double), intent(in) :: R123
! Equilibrium length of segment (X1,X2) and (X2,X3) (It is assumed to be the same for both segments)
real(c_double), intent(in) :: L123
real(c_double), intent(in) :: L123
integer(c_int), intent(inout) :: BBF2
!-------------------------------------------------------------------------------------------
integer(c_int) :: ii, jj, Event
@ -105,9 +105,9 @@ contains !**********************************************************************
U = U / 4.0d+00
F = - F
F1 = F1 + G0 * F
F2 = F2 + G1 * F
F3 = F3 + G2 * F
F1 = F1 + G0 * F
F2 = F2 + G1 * F
F3 = F3 + G2 * F
U1 = U1 + U
U2 = U2 + 2.0d+00 * U
@ -117,16 +117,16 @@ contains !**********************************************************************
do ii = 0, 2
do jj = 0, 2
S = 0.5d+00 * ( X1(ii) - X2(ii) ) * G0(jj)
S1(ii,jj) = S1(ii,jj) + S
S2(ii,jj) = S2(ii,jj) + S
S1(ii,jj) = S1(ii,jj) + S
S2(ii,jj) = S2(ii,jj) + S
S = 0.5d+00 * ( X3(ii) - X2(ii) ) * G2(jj)
S3(ii,jj) = S3(ii,jj) + S
S2(ii,jj) = S2(ii,jj) + S
S3(ii,jj) = S3(ii,jj) + S
S2(ii,jj) = S2(ii,jj) + S
end do
end do
end subroutine TubeBendingForceField !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The purpose of subroutine SegmentTubeForceField is to calculate interaction forces
! The purpose of subroutine SegmentTubeForceField is to calculate interaction forces
! (as well potential energies and components of the virial stress tensor) between a segment
! (X1,X2) and a sequence of segments which belongs to a single CNT.
@ -134,7 +134,7 @@ contains !**********************************************************************
! neighbor list of segment (X1,X2).
! The nodes in X are assumed to be ordered according to their physical appearance in the nanotube.
! It means that (X(i),X(i+1)) are either correspond to a real segment or divided by segments
! It means that (X(i),X(i+1)) are either correspond to a real segment or divided by segments
! that do not belong to a nanotube.
! Concept of the extended chain:
@ -143,45 +143,45 @@ contains !**********************************************************************
! If node Xe follows XN and Xe is the nanotube end, then the extended chain is (X1,...,XN,Xe) and Ee = 2.
! In all other cases, the extended chain coincides with (X1,...,XN) and Ee = 0.
! If the extended chain contains additional node, then non-zero force is exerted on this node.
subroutine SegmentTubeForceField ( U1, U2, U, F1, F2, F, Fe, S1, S2, S, Se, X1, X2, R12, N, X, Xe,&
BBF, R, E1, E2, Ee, TPMType )
BBF, R, E1, E2, Ee, TPMType )
! Number of nodes in array X
integer(c_int), intent(in) :: N
! Interaction energies associated with nodes X1 and X2
real(c_double), intent(inout) :: U1, U2
integer(c_int), intent(in) :: N
! Interaction energies associated with nodes X1 and X2
real(c_double), intent(inout) :: U1, U2
! Interaction energies associated with nodes X
real(c_double), intent(inout), dimension(0:N-1) :: U
real(c_double), intent(inout), dimension(0:N-1) :: U
! Forces exerted on nodes X1 and X2
real(c_double), intent(inout), dimension(0:2) :: F1, F2
real(c_double), intent(inout), dimension(0:2) :: F1, F2
! Forces exerted on nodes X
real(c_double), intent(inout), dimension(0:2,0:N-1) :: F
real(c_double), intent(inout), dimension(0:2,0:N-1) :: F
! Force exerted on node Xe (can be updated only if Ee > 0)
real(c_double), intent(inout), dimension(0:2) :: Fe
real(c_double), intent(inout), dimension(0:2) :: Fe
! Contributions of nodes X1 and X2 to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2
! Contributions of nodes X to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2,0:N-1) :: S
real(c_double), intent(inout), dimension(0:2,0:2,0:N-1) :: S
! Contributions of node Xe to the virial stress tensor (can be updated only if Ee > 0)
real(c_double), intent(inout), dimension(0:2,0:2) :: Se
! Coordinates of the segment nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2
real(c_double), intent(inout), dimension(0:2,0:2) :: Se
! Coordinates of the segment nodes
real(c_double), intent(in), dimension(0:2) :: X1, X2
! Radius of a nanotube the segment (X1,X2) belongs to
real(c_double), intent(in) :: R12
real(c_double), intent(in) :: R12
! Coordinates of the nanotube nodes
real(c_double), intent(in), dimension(0:2,0:N-1) :: X
real(c_double), intent(in), dimension(0:2,0:N-1) :: X
! Additional node of the extended chain if Ee > 0
real(c_double), intent(in), dimension(0:2) :: Xe
real(c_double), intent(in), dimension(0:2) :: Xe
! Bending buckling flags (BBF(i) = 1 in a case of buckling in node i)
integer(c_int), intent(in), dimension(0:N-1) :: BBF
integer(c_int), intent(in), dimension(0:N-1) :: BBF
! Radius of nanotube X
real(c_double), intent(in) :: R
real(c_double), intent(in) :: R
! E1 = 1 if the chain node 0 is a CNT end; E1 = 2 if the chain node N-1 is a CNT end
integer(c_int), intent(in) :: E1, E2
integer(c_int), intent(in) :: E1, E2
! Parameter defining the type of the extended chain (0,1,2)
integer(c_int), intent(in) :: Ee
integer(c_int), intent(in) :: Ee
! Type of the tubular potential (0 or 1)
integer(c_int), intent(in) :: TPMType
integer(c_int), intent(in) :: TPMType
!-------------------------------------------------------------------------------------------
integer(c_int) :: k, ii, jj, IntSign
integer(c_int) :: BType, EType, LocalTPMType
@ -191,18 +191,18 @@ contains !**********************************************************************
real(c_double), dimension(0:2) :: G, DG, DQ, XX
real(c_double) :: UT, DR, DS, DS1
! Interaction energies associated with nodes X1 and X2
real(c_double) :: xU1, xU2
real(c_double) :: xU1, xU2
! Interaction energies associated with nodes X
real(c_double), dimension(0:N-1) :: xU
real(c_double), dimension(0:N-1) :: xU
! Forces exerted on nodes X1 and X2
real(c_double), dimension(0:2) :: xF1, xF2
real(c_double), dimension(0:2) :: xF1, xF2
! Forces exerted on nodes X
real(c_double), dimension(0:2,0:N-1) :: xF
real(c_double), dimension(0:2,0:N-1) :: xF
! Force exerted on node Xe (can be updated only if Ee > 0)
real(c_double), dimension(0:2) :: xFe
real(c_double), dimension(0:2) :: xFe
!-------------------------------------------------------------------------------------------
! Looking for a buckling point
! Looking for a buckling point
BType = 0
do k = 0, N - 1
if ( BBF(k) == 1 ) then
@ -212,7 +212,7 @@ contains !**********************************************************************
end do
! Choosing the LocalTPMType and Etype.
! LocalTPMType is set to 0 if both ends of the chain are nanotube ends or the chain contains a buckling point.
! LocalTPMType is set to 0 if both ends of the chain are nanotube ends or the chain contains a buckling point.
! Overwise, LocalTPMType = TPMType.
if ( BType == 1 ) then
LocalTPMType = 0
@ -231,7 +231,7 @@ contains !**********************************************************************
if ( EType1 .and. EType2 ) then
LocalTPMType = 0
else
LocalTPMType = TPMType
LocalTPMType = TPMType
if ( EType1 ) then
EType = 1
else if ( EType2 ) then
@ -241,7 +241,7 @@ contains !**********************************************************************
end if
end if
end if
if ( LocalTPMType == 0 ) then
IntSign = TPMInteractionFW0 ( QQ, UT, xU1, xU2, xU, xF1, xF2, xF, G1, G2, X1, X2, N, N, X )
else
@ -263,10 +263,10 @@ contains !**********************************************************************
! Contributions to the virial stresses tensor
do ii = 0, 2
DR = 0.125d+00 * ( X2(ii) - X1(ii) )
DR = 0.125d+00 * ( X2(ii) - X1(ii) )
do jj = 0, 2
DS = DR * ( xF2(jj) - xF1(jj) )
S1(ii,jj) = S1(ii,jj) + DS
S1(ii,jj) = S1(ii,jj) + DS
S2(ii,jj) = S2(ii,jj) + DS
end do
end do
@ -278,24 +278,24 @@ contains !**********************************************************************
do ii = 0, 2
do jj = 0, 2
DS = DQ(ii) * xFe(jj)
S1(ii,jj) = S1(ii,jj) + DS
S1(ii,jj) = S1(ii,jj) + DS
S2(ii,jj) = S1(ii,jj) + DS
Se(ii,jj) = Se(ii,jj) + DS
end do
end do
end if
do k = 0, N - 2
DQ = 0.5d+00 * ( X(0:2,k+1) + X(0:2,k) ) - XX
DQ = 0.5d+00 * ( X(0:2,k+1) + X(0:2,k) ) - XX
call ApplyPeriodicBC ( DQ )
DQ = 0.125d+00 * DQ
G = G1(0:2,k+1) + G2(0:2,k)
DQ = 0.125d+00 * DQ
G = G1(0:2,k+1) + G2(0:2,k)
DG = G1(0:2,k+1) - G2(0:2,k)
do ii = 0, 2
DR = 0.125d+00 * ( X(ii,k+1) - X(ii,k) )
DR = 0.125d+00 * ( X(ii,k+1) - X(ii,k) )
do jj = 0, 2
DS = DQ(ii) * G(jj)
DS1 = DS + DR * DG(jj)
S1(ii,jj) = S1(ii,jj) + DS
S1(ii,jj) = S1(ii,jj) + DS
S2(ii,jj) = S2(ii,jj) + DS
S(ii,jj,k) = S(ii,jj,k) + DS1
S(ii,jj,k+1) = S(ii,jj,k+1) + DS1

18
lib/mesont/TPMGeom.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module TPMGeom !************************************************************************************
!
@ -39,20 +39,20 @@ implicit none
!---------------------------------------------------------------------------------------------------
! Global variables
!---------------------------------------------------------------------------------------------------
! Coordinates of the whole domain
real(c_double) :: DomXmin, DomXmax, DomYmin, DomYmax, DomZmin, DomZmax
real(c_double) :: DomLX, DomLY, DomLZ
real(c_double) :: DomLXhalf, DomLYhalf, DomLZhalf
! Boundary conditions
! Boundary conditions
integer(c_int) :: BC_X = 0
integer(c_int) :: BC_Y = 0
integer(c_int) :: BC_Z = 0
! Skin parameter in NBL and related algorithms
real(c_double) :: Rskin = 1.0d+00
contains !******************************************************************************************
subroutine ApplyPeriodicBC ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@ -111,12 +111,12 @@ contains !**********************************************************************
! cosA, cosine of the angle between lines.
! D1, D2, displacements.
! L12, unit vector directed along the closest distance.
!-------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------
real(c_double), intent(inout) :: H, cosA, D1, D2
real(c_double), dimension(0:2), intent(out) :: L12
real(c_double), dimension(0:2), intent(in) :: R1, L1, R2, L2
!-------------------------------------------------------------------------------------------
real(c_double), intent(in) :: Prec
real(c_double), intent(in) :: Prec
real(c_double), dimension(0:2) :: Q1, Q2, R
real(c_double) :: C, DD1, DD2, C1, C2
!-------------------------------------------------------------------------------------------
@ -151,5 +151,5 @@ contains !**********************************************************************
L12 = L12 / H
end if
end function LineLine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TPMGeom !********************************************************************************

18
lib/mesont/TPMLib.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module TPMLib !*************************************************************************************
!
@ -61,7 +61,7 @@ contains !**********************************************************************
!---------------------------------------------------------------------------------------------------
! Simple mathematical functions
!---------------------------------------------------------------------------------------------------
real(c_double) function rad ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: X
!-------------------------------------------------------------------------------------------
@ -76,7 +76,7 @@ contains !**********************************************************************
integer(c_int) function signum ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: X
!-------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------
if ( X > 0 ) then
signum = 1
else if ( X < 0 ) then
@ -111,7 +111,7 @@ contains !**********************************************************************
subroutine V3_ort ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), dimension(0:2), intent(inout) :: V
!-------------------------------------------------------------------------------------------
real(c_double) :: Vabs
real(c_double) :: Vabs
!-------------------------------------------------------------------------------------------
Vabs = S_V3norm3 ( V )
V(0) = V(0) / Vabs
@ -131,7 +131,7 @@ contains !**********************************************************************
!---------------------------------------------------------------------------------------------------
! Handling of spherical and Euler angles
!---------------------------------------------------------------------------------------------------
subroutine RotationMatrix3 ( M, Psi, Tet, Phi ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Ksi, Tet and Phi are Euler angles
!-------------------------------------------------------------------------------------------
@ -156,7 +156,7 @@ contains !**********************************************************************
M(2,1) = - sT * cK
M(2,2) = cT
end subroutine RotationMatrix3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine EulerAngles ( Psi, Tet, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: Tet, Psi
real(c_double), dimension(0:2), intent(in) :: L
@ -165,7 +165,7 @@ contains !**********************************************************************
Psi = atan2 ( L(1), L(0) )
if ( Psi > M_3PI_2 ) then
Psi = Psi - M_3PI_2
else
else
Psi = Psi + M_PI_2
end if
end subroutine EulerAngles !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@ -211,5 +211,5 @@ contains !**********************************************************************
close ( unit = Fuid )
Fuid = -1
end subroutine CloseFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TPMLib !*********************************************************************************

34
lib/mesont/TPMM0.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module TPMM0 !**************************************************************************************
!
@ -32,7 +32,7 @@ use iso_c_binding, only : c_int, c_double, c_char
implicit none
contains !******************************************************************************************
integer(c_int) function TPMInteractionFSS ( Q, U, F1_1, F1_2, F2_1, F2_2, R1_1, R1_2, R2_1, R2_2, EType )
real(c_double), intent(inout) :: Q, U
real(c_double), dimension(0:2), intent(inout) :: F1_1, F1_2, F2_1, F2_2
@ -48,9 +48,9 @@ contains !**********************************************************************
L2 = S_V3norm3 ( Laxis2 )
Laxis2 = Laxis2 / L2
if ( EType < 2 ) then
TPMInteractionFSS = TPMInteractionF ( Q, U, F1_1, F1_2, F2_1, F2_2, Fd, R1_1, R1_2, R2_1, R2_2, 1 )
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 )
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 )
@ -59,7 +59,7 @@ contains !**********************************************************************
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 )
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
@ -74,13 +74,13 @@ contains !**********************************************************************
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
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 )
real(c_double), intent(inout) :: U, U1, U2
integer(c_int), intent(in) :: N, NMAX
@ -142,7 +142,7 @@ contains !**********************************************************************
if ( D < Dmina ) Dmina = D
end if
end if
DR = R1 - R(0:2,i+1)
call ApplyPeriodicBC ( DR )
Dminb = sqr ( DR(0) ) + sqr ( DR(1) ) + sqr ( DR(2) )
@ -165,20 +165,20 @@ contains !**********************************************************************
if ( D < Dminb ) Dminb = D
end if
end if
if ( Dmina < Dminb ) then
EType = 1
else
EType = 2
end if
if ( TPMInteractionFSS ( QQ(i), Ua, F1_1a, F1_2a, F2_1a, F2_2a, R1, R2, R(0:2,i), R(0:2,i+1), &
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
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
@ -190,5 +190,5 @@ contains !**********************************************************************
end if
end do
end function TPMInteractionFW0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TPMM0 !**********************************************************************************

62
lib/mesont/TPMM1.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module TPMM1 !**************************************************************************************
!
@ -56,9 +56,9 @@ implicit none
real(c_double), dimension(0:TPM_MAX_CHAIN-1) :: W, C
real(c_double), dimension(0:2) :: RR, E10
real(c_double) :: L10, D10
contains !******************************************************************************************
subroutine PairWeight1 ( W, E1_1, E1_2, E2_1, E2_2, R2_1, R2_2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: W
real(c_double), dimension(0:2), intent(out) :: E1_1, E1_2, E2_1, E2_2
@ -78,7 +78,7 @@ contains !**********************************************************************
E2_2 = 0.0d+00
return
end if
E20 = 0.5d+00 * ( R2_2 - R2_1 )
E20 = 0.5d+00 * ( R2_2 - R2_1 )
L20 = sqrt ( S_V3xx ( E20 ) + sqr ( TPMR2 ) )
D20 = L10 + L20 + TPBRcutoff + RSkin
if ( D > D20 * D20 ) then
@ -96,12 +96,12 @@ contains !**********************************************************************
t = ( D - D10 ) / D20
W = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
dWdD = 3.0d+00 * t * ( t - 1.0d+00 ) / D20
E1_1 = dWdD * ( t * E10 - E )
E1_2 = dWdD * ( - t * E10 - E )
E2_1 = dWdD * ( E + t * E20 )
E2_2 = dWdD * ( E - t * E20 )
E1_1 = dWdD * ( t * E10 - E )
E1_2 = dWdD * ( - t * E10 - E )
E2_1 = dWdD * ( E + t * E20 )
E2_2 = dWdD * ( E - t * E20 )
end subroutine PairWeight1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function EndWeight1 ( W, E1_1, E1_2, E2_1, E2_2, R1_1, R1_2, R2_1, R2_2 ) !!!
real(c_double), intent(out) :: W
real(c_double), dimension(0:2), intent(out) :: E1_1, E1_2, E2_1, E2_2
@ -114,7 +114,7 @@ contains !**********************************************************************
E = 0.5d+00 * ( R2_1 + R2_2 - ( R1_1 + R1_2 ) )
call ApplyPeriodicBC ( E )
D = S_V3norm3 ( E )
E20 = 0.5d+00 * ( R2_2 - R2_1 )
E20 = 0.5d+00 * ( R2_2 - R2_1 )
L20 = sqrt ( S_V3xx ( E20 ) + sqr ( TPMR2 ) )
D1 = L10 + L20 + TPBRcutoff + RSkin
if ( D < D1 ) then
@ -126,7 +126,7 @@ contains !**********************************************************************
E2_2 = 0.0d+00
return
end if
D2 = D1 + TPMC3
D2 = D1 + TPMC3
if ( D > D2 ) then
EndWeight1 = 2
W = 0.0d+00
@ -142,13 +142,13 @@ contains !**********************************************************************
t = ( D - D1 ) / TPMC3
W = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
dWdD = 3.0d+00 * t * ( t - 1.0d+00 ) / TPMC3
E1_1 = dWdD * ( E10 - E )
E1_2 = dWdD * ( - E10 - E )
E2_1 = dWdD * ( E + E20 )
E2_2 = dWdD * ( E - E20 )
E1_1 = dWdD * ( E10 - E )
E1_2 = dWdD * ( - E10 - E )
E2_1 = dWdD * ( E + E20 )
E2_2 = dWdD * ( E - E20 )
end function EndWeight1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMInteractionFC1 ( Q, U, F1, F2, P1, P2, Pe, Pe1, R1, R2, Q1, Q2, Qe, Qe1, EType )
integer(c_int) function TPMInteractionFC1 ( Q, U, F1, F2, P1, P2, Pe, Pe1, R1, R2, Q1, Q2, Qe, Qe1, EType )
real(c_double), intent(out) :: Q, U
real(c_double), dimension(0:2), intent(out) :: F1, F2, P1, P2, Pe, Pe1
real(c_double), dimension(0:2), intent(in) :: R1, R2, Q1, Q2, Qe, Qe1
@ -196,11 +196,11 @@ contains !**********************************************************************
end if
call TPMSegmentForces ( P1a, P2a, F1a, F2a, R1, R2, QX, M, L )
end if
if ( CaseID > 0 ) then
IntSignb = TPMInteractionF ( Qb, Ub, F1b, F2b, P1b, P2b, Peeb, R1, R2, Q1, Q2, 0 )
end if
if ( CaseID == 0 ) then
TPMInteractionFC1 = IntSigna
Q = Qa
@ -226,7 +226,7 @@ contains !**********************************************************************
TPMInteractionFC1 = 0
if ( IntSigna > 0 .or. IntSignb > 0 ) TPMInteractionFC1 = 1
W1 = 1.0d+00 - W
DU = Ub - Ua
DU = Ub - Ua
Q = W * Qa + W1 * Qb
U = W * Ua + W1 * Ub
Pe = ( W * Peea / D ) * Me
@ -240,7 +240,7 @@ contains !**********************************************************************
end if
end if
end function TPMInteractionFC1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMInteractionFW1 ( QQ, U, U1, U2, UU, F1, F2, F, Fe, G1, G2, R1, R2, N, NMAX, R, Re, EType )
real(c_double), intent(out) :: U, U1, U2
integer(c_int), intent(in) :: N, NMAX, EType
@ -251,7 +251,7 @@ contains !**********************************************************************
real(c_double), dimension(0:2,0:NMAX-1), intent(in) :: R
!-------------------------------------------------------------------------------------------
integer(c_int) :: i, j
real(c_double) :: Q, WW, DD
real(c_double) :: Q, WW, DD
!-------------------------------------------------------------------------------------------
Q1 = 0.0d+00
Q2 = 0.0d+00
@ -259,7 +259,7 @@ contains !**********************************************************************
Z1 = 0.0d+00
Z2 = 0.0d+00
TPMInteractionFW1 = 0
E10 = 0.5d+00 * ( R2 - R1 )
E10 = 0.5d+00 * ( R2 - R1 )
L10 = sqrt ( S_V3xx ( E10 ) + sqr ( TPMR1 ) )
D10 = TPMR1 + TPMR2 + TPMC123 * TPBRcutoff + RSkin
E10 = E10 / L10
@ -293,10 +293,10 @@ contains !**********************************************************************
Qe = 0.0d+00
Qe1 = 0.0d+00
end if
TPMInteractionFW1 = TPMInteractionFC1 ( Q, U, F1, F2, S1, S2, Pe, Pe1, R1, R2, Q1, Q2, Qe, Qe1, EType )
TPMInteractionFW1 = TPMInteractionFC1 ( Q, U, F1, F2, S1, S2, Pe, Pe1, R1, R2, Q1, Q2, Qe, Qe1, EType )
if ( TPMInteractionFW1 == 0 ) return
W(0:N-2) = W(0:N-2) / WW
E1(0:2,0:N-2) = E1(0:2,0:N-2) / WW
E2(0:2,0:N-2) = E2(0:2,0:N-2) / WW
@ -322,7 +322,7 @@ contains !**********************************************************************
do j = 0, N - 2
if ( j == 0 ) then
DR = EE1(0:2,0) * ( 1.0d+00 - W(0) )
else
else
DR = - W(j) * EE1(0:2,0)
end if
F(0:2,0) = F(0:2,0) + C(j) * DR
@ -337,8 +337,8 @@ contains !**********************************************************************
else if ( j == i - 1 ) then
G1(0:2,i) = G1(0:2,i) + C(j) * ( EE2(0:2,j) - W(j) * EE2(0:2,i-1) )
G2(0:2,i) = G2(0:2,i) - C(j) * W(j) * EE1(0:2,i)
else
G1(0:2,i) = G1(0:2,i) - C(j) * W(j) * EE2(0:2,i-1)
else
G1(0:2,i) = G1(0:2,i) - C(j) * W(j) * EE2(0:2,i-1)
G2(0:2,i) = G2(0:2,i) - C(j) * W(j) * EE1(0:2,i)
end if
end do
@ -348,7 +348,7 @@ contains !**********************************************************************
do j = 0, N - 2
if ( j == N - 2 ) then
DR = EE2(0:2,N-2) * ( 1.0d+00 - W(N-2) )
else
else
DR = - W(j) * EE2(0:2,N-2)
end if
F(0:2,N-1) = F(0:2,N-1) + C(j) * DR
@ -368,5 +368,5 @@ contains !**********************************************************************
G1(0:2,N-1) = F(0:2,N-1)
G2(0:2,0) = F(0:2,0)
end function TPMInteractionFW1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TPMM1 !**********************************************************************************

88
lib/mesont/TubePotBase.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module TubePotBase !********************************************************************************
!
@ -25,10 +25,10 @@ module TubePotBase !************************************************************
!
!---------------------------------------------------------------------------------------------------
!
! This module contains basic parameters for all modules involved into calculations of tubular
! This module contains basic parameters for all modules involved into calculations of tubular
! potentials.
!
! It includes definitions of
!
! It includes definitions of
! -- TPBU, Lennard-Jones (12-6) potential
! -- TPBQ, Transfer function
!
@ -43,7 +43,7 @@ implicit none
!---------------------------------------------------------------------------------------------------
! Constants
!---------------------------------------------------------------------------------------------------
! Types of the potential with respect to the breathing mode
integer(c_int), parameter :: TP_POT_MONO_R = 0
integer(c_int), parameter :: TP_POT_POLY_R = 1
@ -56,7 +56,7 @@ implicit none
! Mass of C atom
real(c_double), parameter :: TPBMc = 12.0107d+00 ! (Da)
! Parameters of the Van der Waals interaction between carbon atoms in graphene sheets, see
! Stuart S.J., Tutein A.B., Harrison J.A., J. Chem. Phys. 112(14), 2000
real(c_double), parameter :: TPBEcc = 0.00284d+00 ! (eV)
@ -66,73 +66,73 @@ implicit none
! Dresselhaus et al, Carbon 33(7), 1995
real(c_double), parameter :: TPBAcc = 1.421d+00 ! (A)
real(c_double), parameter :: TPBDcc = 4.0d+00 / ( TPBConstD * TPBAcc * TPBAcc ) ! (1/A^2)
! Specific heat of carbon nanotubes
real(c_double), parameter :: TPBSHcc = 600.0d+00 / K_MDCU ! (eV/(Da*K))
! Cutoff distances for the interactomic potential and transfer function.
! Changes in these parameters can result in necessity to change some numerical parameters too.
real(c_double), parameter :: TPBRmincc = 0.001d+00 * TPBScc ! (A)
real(c_double), parameter :: TPBRcutoffcc = 3.0d+00 * TPBScc ! (A)
real(c_double), parameter :: TPBRcutoff1cc = 2.16d+00 * TPBScc ! (A)
! Parameters of the transfer function for non-bonded interaction between carbon atoms
real(c_double), parameter :: TPBQScc = 7.0d+00 ! (A)
real(c_double), parameter :: TPBQRcutoff1cc = 8.0d+00 ! (A)
!---------------------------------------------------------------------------------------------------
! Global variables
!---------------------------------------------------------------------------------------------------
! Set to .true. to generate diagnostic and warning messages
logical :: TPErrCheck = .true.
character*512 :: TPErrMsg = ''
real(c_double) :: TPGeomPrec = 1.0d-06 ! Geometric precision, see TPInt
logical :: TPErrCheck = .true.
character*512 :: TPErrMsg = ''
real(c_double) :: TPGeomPrec = 1.0d-06 ! Geometric precision, see TPInt
integer(c_int) :: TPPotType = TP_POT_MONO_R ! Type of the potential with respect to the breathing mode
! Parameters of the interatomic potential and atoms distribution at the surface
! of the tube
real(c_double) :: TPBM = TPBMc ! Mass of an atom (Da)
real(c_double) :: TPBE = TPBEcc ! Depth of the energy well in (12-6) LJ interatomic potential (eV)
real(c_double) :: TPBS = TPBScc ! Sigma parameter of (12-6) LJ interatomic potential (A)
real(c_double) :: TPBD = TPBDcc ! Numerical density of atoms at the tube surface (1/A^2)
real(c_double) :: TPBSH = TPBSHcc ! Specific heat (eV/(Da*K))
real(c_double) :: TPBRmin = TPBRmincc ! (A)
real(c_double) :: TPBRcutoff = TPBRcutoffcc ! (A)
real(c_double) :: TPBRcutoff1 = TPBRcutoff1cc ! (A)
! Parameters of the transfer function
real(c_double) :: TPBQS = TPBQScc ! Sigma parameter of the transfer function (A)
real(c_double) :: TPBQRcutoff1 = TPBQRcutoff1cc! (A)
! Auxiliary variables
real(c_double) :: TPBE4, TPBE24, TPBDRcutoff, TPBQDRcutoff
real(c_double) :: TPBQR0 ! Constant-value distance for the transfer function (A)
! Table of inter-particle potential, force, and transfer function
integer(c_int) :: TPBN = TPBNMAX
real(c_double) :: TPBDR
real(c_double), dimension(0:TPBNMAX-1) :: TPBQ
real(c_double), dimension(0:TPBNMAX-1) :: TPBU, TPBdUdR
contains !******************************************************************************************
integer(c_int) function TPBsizeof () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
TPBsizeof = 8 * ( size ( TPBQ ) + size ( TPBU ) + size ( TPBdUdR ) )
end function TPBsizeof !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Interpolation
!---------------------------------------------------------------------------------------------------
real(c_double) function TPBQInt0 ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: R
real(c_double), intent(in) :: R
!-------------------------------------------------------------------------------------------
real(c_double) :: Z, RR
integer(c_int) :: i
@ -140,7 +140,7 @@ contains !**********************************************************************
if ( R < TPBRmin ) then
!call PrintStdLogMsg ( TPErrMsg )
!write ( TPErrMsg, '(a,e20.10,a,e20.10)' ) ': R < Rmin: R=', R, ', Rmin=', TPBRmin
!call Error ( 'TPBQInt0', TPErrMsg )
!call Error ( 'TPBQInt0', TPErrMsg )
elseif ( R > TPBRcutoff ) then
TPBQInt0 = 0.0d+00
return
@ -149,11 +149,11 @@ contains !**********************************************************************
i = int ( RR )
RR = RR - i
Z = 1.0d+00 - RR
TPBQInt0 = TPBQ(i) * Z + TPBQ(i+1) * RR
TPBQInt0 = TPBQ(i) * Z + TPBQ(i+1) * RR
end function TPBQInt0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function TPBUInt0 ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: R
real(c_double), intent(in) :: R
!-------------------------------------------------------------------------------------------
real(c_double) :: Z, RR
integer(c_int) :: i
@ -161,7 +161,7 @@ contains !**********************************************************************
if ( R < TPBRmin ) then
!call PrintStdLogMsg ( TPErrMsg )
!write ( TPErrMsg, '(a,e20.10,a,e20.10)' ) ': R < Rmin: R=', R, ', Rmin=', TPBRmin
!call Error ( 'TPBUInt0', TPErrMsg )
!call Error ( 'TPBUInt0', TPErrMsg )
elseif ( R > TPBRcutoff ) then
TPBUInt0 = 0.0d+00
return
@ -170,12 +170,12 @@ contains !**********************************************************************
i = int ( RR )
RR = RR - i
Z = 1.0d+00 - RR
TPBUInt0 = TPBU(i) * Z + TPBU(i+1) * RR
TPBUInt0 = TPBU(i) * Z + TPBU(i+1) * RR
end function TPBUInt0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPBUInt1 ( U, dUdR, R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: U, dUdR
real(c_double), intent(in) :: R
real(c_double), intent(in) :: R
!-------------------------------------------------------------------------------------------
real(c_double) :: Z, RR
integer(c_int) :: i
@ -183,7 +183,7 @@ contains !**********************************************************************
if ( R < TPBRmin ) then
!call PrintStdLogMsg ( TPErrMsg )
!write ( TPErrMsg, '(a,e20.10,a,e20.10)' ) ': R < Rmin: R=', R, ', Rmin=', TPBRmin
!call Error ( 'TPBUInt1', TPErrMsg )
!call Error ( 'TPBUInt1', TPErrMsg )
elseif ( R > TPBRcutoff ) then
TPBU = 0.0d+00
TPBdUdR = 0.0d+00
@ -193,14 +193,14 @@ contains !**********************************************************************
i = int ( RR )
RR = RR - i
Z = 1.0d+00 - RR
U = TPBU(i) * Z + TPBU(i+1) * RR
U = TPBU(i) * Z + TPBU(i+1) * RR
dUdR = TPBdUdR(i) * Z + TPBdUdR(i+1) * RR
end subroutine TPBUInt1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Calculation
!---------------------------------------------------------------------------------------------------
real(c_double) function TPBQCalc0 ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: R
!-------------------------------------------------------------------------------------------
@ -210,7 +210,7 @@ contains !**********************************************************************
TPBQCalc0 = 0.0d+00
else if ( R < TPBQR0 ) then
TPBQCalc0 = 1.0d+00
else
else
Z = TPBQS / R
Z = Z * Z * Z
Z = Z * Z
@ -222,7 +222,7 @@ contains !**********************************************************************
endif
endif
end function TPBQCalc0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function TPBUCalc0 ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: R
!-------------------------------------------------------------------------------------------
@ -230,7 +230,7 @@ contains !**********************************************************************
!-------------------------------------------------------------------------------------------
if ( R > TPBRcutoff ) then
TPBUCalc0 = 0.0d+00
else
else
Z = TPBS / R
Z = Z * Z * Z
Z = Z * Z
@ -251,7 +251,7 @@ contains !**********************************************************************
if ( R > TPBRcutoff ) then
U = 0.0d+00
dUdR = 0.0d+00
else
else
Z = TPBS / R
Z = Z * Z * Z
Z = Z * Z
@ -266,7 +266,7 @@ contains !**********************************************************************
endif
endif
end subroutine TPBUCalc1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPBSegmentForces ( F1, F2, F, M, Laxis, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), dimension(0:2), intent(out) :: F1, F2
real(c_double), dimension(0:2), intent(in) :: F, M, Laxis
@ -284,7 +284,7 @@ contains !**********************************************************************
!---------------------------------------------------------------------------------------------------
! Initialization
!---------------------------------------------------------------------------------------------------
subroutine TPBInit () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) :: R
integer(c_int) :: i

348
lib/mesont/TubePotMono.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module TubePotMono !********************************************************************************
!
@ -27,22 +27,22 @@ module TubePotMono !************************************************************
!
! Four potentials and transfer functions are calculated in this module:
!
! 1. SSTP (segment - semi-infinite tube parallel): Linear density of the potential along
! the segment axis which is produced by a parallel semi-infinite tube. 2D tables for this potential
! 1. SSTP (segment - semi-infinite tube parallel): Linear density of the potential along
! the segment axis which is produced by a parallel semi-infinite tube. 2D tables for this potential
! are generated at initialization or can be loaded from a file.
!
! 2. STP (segment - tube parallel): Linear density of the potential along the segment axis
! which is produced by a parallel infinite tubes. This is only a particular case of the SSTP potential,
! but it is considered separately for computational efficiency. 1D tables of this potential are taken
! 2. STP (segment - tube parallel): Linear density of the potential along the segment axis
! which is produced by a parallel infinite tubes. This is only a particular case of the SSTP potential,
! but it is considered separately for computational efficiency. 1D tables of this potential are taken
! from 2D tables of SSTP potential.
!
! 3. SST (segment - semi-infinite tube): Potential for a segment produced by an arbitrary-
! oriented semi-infinite tube. This potential can not be kept in 2D tables, therefore, all
! data are calculated 'on fly' with the help of SSTP potential and numerical integration along the
! oriented semi-infinite tube. This potential can not be kept in 2D tables, therefore, all
! data are calculated 'on fly' with the help of SSTP potential and numerical integration along the
! segment axis
!
! 4. ST (segment - tube): Potential for a segment produced by an arbitrary-oriented
! infinitely long tube. 2D tables for this potential are generated at initialization or can be
! 4. ST (segment - tube): Potential for a segment produced by an arbitrary-oriented
! infinitely long tube. 2D tables for this potential are generated at initialization or can be
! loaded from a file.
!
!***************************************************************************************************
@ -62,77 +62,77 @@ implicit none
integer(c_int), parameter :: TPMNZMAX = 129
integer(c_int), parameter :: TPMNEMAX = 128
integer(c_int), parameter :: TPMNHMAX = 1001
integer(c_int), parameter :: TPMNXMAX = 1001
integer(c_int), parameter :: TPMNMAX = 1001
integer(c_int), parameter :: TPMNHMAX = 1001
integer(c_int), parameter :: TPMNXMAX = 1001
integer(c_int), parameter :: TPMNMAX = 1001
!---------------------------------------------------------------------------------------------------
! Global variables
!---------------------------------------------------------------------------------------------------
integer(c_int) :: TPMStartMode = 1
character(len=512) :: TPMFile = 'MESONT-TABTP.xrs'
integer(c_int) :: TPMUnitID ! Unit for the tabulated potential file
integer(c_int) :: TPMNZ = TPMNZMAX
integer(c_int) :: TPMNZ1 = TPMNZMAX - 1
integer(c_int) :: TPMNE = TPMNEMAX
integer(c_int) :: TPMNE1 = TPMNEMAX - 1
integer(c_int) :: TPMNH = TPMNHMAX
integer(c_int) :: TPMNH1 = TPMNHMAX - 1
integer(c_int) :: TPMNX = TPMNXMAX
integer(c_int) :: TPMNX1 = TPMNXMAX - 1
integer :: TPMChiIndM ! Chirality index M
integer :: TPMChiIndM ! Chirality index M
integer :: TPMChiIndN ! Chirality index N
real(c_double) :: TPMR1
real(c_double) :: TPMR2
real(c_double) :: TPMHmax
real(c_double) :: TPMDH
! Parameters of empirical correction functions
integer(c_int) :: TPMAN = 20
real(c_double) :: TPMAHmin
real(c_double) :: TPMAHmax
real(c_double) :: TPMADH
real(c_double), dimension(0:TPMNHMAX-1) :: TPMAH, TPMAF, TPMAFxx
! Fitting parameters that depend on the SWCNT chirality
real(c_double) :: TPMCaA = 0.22d+00 ! 0.22 for (10,10) CNTs
real(c_double) :: TPMCeA = 0.35d+00 ! 0.35 for (10,10) CNTs
real(c_double) :: TPMAHmin0 = 10.0d+00 ! 10.0 A for (10,10) CNTs
! Parameters of SSTP integrator
real(c_double) :: TPMDE
real(c_double), dimension(0:TPMNEMAX-1) :: TPMCE, TPMSE
! Additional parameters for SSTP potential
real(c_double) :: TPMSSTPDelta = 0.25d+00
integer(c_int) :: TPMSSTPNH
integer(c_int) :: TPMSSTPNX
real(c_double) :: TPMSSTPX1
real(c_double) :: TPMSSTPXmax
real(c_double) :: TPMSSTPDX
real(c_double), dimension(0:TPMNHMAX-1,0:TPMNXMAX-1) :: TPMSSTPG
real(c_double), dimension(0:TPMNHMAX-1,0:TPMNXMAX-1) :: TPMSSTPF, TPMSSTPFxx, TPMSSTPFyy, TPMSSTPFxxyy
real(c_double), dimension(0:TPMNHMAX-1) :: TPMSSTPH
real(c_double), dimension(0:TPMNXMAX-1) :: TPMSSTPX
! Additional parameters for STP potential
! In calculations of this potential, some parameters of SSTP potential are also used.
! In particular, STP potential has no its own integrator. All data comes from SSTP integrator.
! It does not result in any computational inefficiency unless the STP potential is used without SSTP one.
integer(c_int) :: TPMNN = 10
real(c_double), dimension(0:TPMNHMAX-1) :: TPMSTPG
real(c_double), dimension(0:TPMNHMAX-1) :: TPMSTPF, TPMSTPFxx
@ -140,25 +140,25 @@ implicit none
! Parameters for ST potential
! Minimum gap dh for ST-potential
real(c_double) :: TPMSTDelta = 1.0d+00
real(c_double) :: TPMSTDelta = 1.0d+00
! Number of subdivisions for every grid step in ST-integrator
integer(c_int) :: TPMSTNXS = 10
integer(c_int) :: TPMSTNXS = 10
real(c_double) :: TPMSTXmax
real(c_double) :: TPMSTH1
real(c_double) :: TPMSTH2
real(c_double) :: TPMSTDH12
real(c_double), dimension(0:TPMNHMAX-1,0:TPMNXMAX-1) :: TPMSTG
real(c_double), dimension(0:TPMNHMAX-1,0:TPMNXMAX-1) :: TPMSTF, TPMSTFxx, TPMSTFyy, TPMSTFxxyy
real(c_double), dimension(0:TPMNHMAX-1) :: TPMSTH
real(c_double), dimension(0:TPMNXMAX-1) :: TPMSTX
! Switching parameters
! Height switch (at H=0 in SST-potential)
integer(c_int) :: TPMHSwitch = 0 ! 1, use h-switch; 0, do not use the switch
real(c_double) :: TPMHS = 3.0d+00 ! Switch height, Angstrom
! Angle switch
integer(c_int) :: TPMASwitch = 0 ! 1, use a-switch; 0, do not use the switch
real(c_double) :: TPMAS = 3.0d+00 ! Switch angle, degree
@ -220,7 +220,7 @@ contains !**********************************************************************
R2Y = TPMR1 * TPMSE(j)
R2Z = Zmin - D
do k = 0, TPMNZ1 ! Integration over the length of the second tube
R = sqr ( R2X - R1X ) + sqr ( R2Y - R1Y ) + sqr ( R2Z )
R = sqr ( R2X - R1X ) + sqr ( R2Y - R1Y ) + sqr ( R2Z )
if ( R < Rcutoff2 ) then
R = dsqrt ( R )
if ( k == 0 .or. k == TPMNZ1 ) then
@ -239,7 +239,7 @@ contains !**********************************************************************
Q = Q * C
U = U * sqr ( TPBD ) * C
end subroutine TPMSSTPIntegrator !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMSSTPInt0 ( Q, U, H, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function returns the transfer function Q and potential U for the SSTP potential
! calculated by interpolation in the table without switch.
@ -253,7 +253,7 @@ contains !**********************************************************************
i = 1 + int ( H / TPMDH )
j = 1 + int ( ( X + TPMSSTPXMax ) / TPMSSTPDX )
if ( ( i < TPMSSTPNH ) .and. ( j > TPMSSTPNX ) ) then
!call PrintTPErrMsg ()
!call PrintTPErrMsg ()
!!call PrintStdLogMsg (TPErrMsg )
!write ( TPErrMsg, '(a,e20.10,a,e20.10)' ) 'Tubes intersect each other: H=', H, ', X=', X
!call Error ( 'TPMSSTPInt0', TPErrMsg )
@ -273,7 +273,7 @@ contains !**********************************************************************
else
XX = X
end if
Q = CalcLinFun2_0 ( i, j, H, XX, TPMNH, TPMNX, TPMSSTPH, TPMSSTPX, TPMSSTPG )
Q = CalcLinFun2_0 ( i, j, H, XX, TPMNH, TPMNX, TPMSSTPH, TPMSSTPX, TPMSSTPG )
U = CalcSpline2_0 ( i, j, H, XX, TPMNH, TPMNX, TPMSSTPH, TPMSSTPX, TPMSSTPF, TPMSSTPFxx, TPMSSTPFyy, TPMSSTPFxxyy )
TPMSSTPInt0 = 1
end function TPMSSTPInt0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@ -303,7 +303,7 @@ contains !**********************************************************************
end if
end if
end function TPMSSTPInt0S !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMSSTPInt1 ( Q, U, Uh, Ux, H, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function returns the transfer function Q, potential U, and derivatives Uh=dU/dH and
! Ux=dU/dX for the SSTP potential calculated by interpolation in the table without switch
@ -317,7 +317,7 @@ contains !**********************************************************************
i = 1 + int ( H / TPMDH )
j = 1 + int ( ( X + TPMSSTPXMax ) / TPMSSTPDX )
if ( ( i < TPMSSTPNH ) .and. ( j > TPMSSTPNX ) ) then
!call PrintTPErrMsg ()
!call PrintTPErrMsg ()
!!call PrintStdLogMsg ( TPErrMsg )
!write ( TPErrMsg, '(a,e20.10,a,e20.10)' ) 'Tubes intersect each other: H=', H, ', X=', X
!call Error ( 'TPMSSTPInt1', TPErrMsg )
@ -339,15 +339,15 @@ contains !**********************************************************************
else
XX = X
end if
Q = CalcLinFun2_0 ( i, j, H, XX, TPMNH, TPMNX, TPMSSTPH, TPMSSTPX, TPMSSTPG )
Q = CalcLinFun2_0 ( i, j, H, XX, TPMNH, TPMNX, TPMSSTPH, TPMSSTPX, TPMSSTPG )
call CalcSpline2_1 ( U, Uh, Ux, i, j, H, XX, TPMNH, TPMNX, TPMSSTPH, TPMSSTPX, TPMSSTPF, &
TPMSSTPFxx, TPMSSTPFyy, TPMSSTPFxxyy )
TPMSSTPFxx, TPMSSTPFyy, TPMSSTPFxxyy )
TPMSSTPInt1 = 1
end function TPMSSTPInt1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMSSTPInt1S ( Q, U, Uh, Ux, H, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function returns the transfer function Q, potential U, and derivatives Uh=dU/dH and
! Ux=dU/dX for the SSTP potential calculated by interpolation in the table and switch to
! Ux=dU/dX for the SSTP potential calculated by interpolation in the table and switch to
! the case of zero H.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Q, U, Uh, Ux
@ -375,7 +375,7 @@ contains !**********************************************************************
end if
end if
end function TPMSSTPInt1S !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPMSSTPWrite () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function writes the table of the SSTP potential to a disk file.
!-------------------------------------------------------------------------------------------
@ -412,7 +412,7 @@ contains !**********************************************************************
end do
end do
end subroutine TPMSSTPRead !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPMSSTPInit () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function calculates the table of the SSTP potential.
!-------------------------------------------------------------------------------------------
@ -424,7 +424,7 @@ contains !**********************************************************************
TPMDE = M_2PI / TPMNE
E = 0.0d+00
do i = 0, TPMNE1
TPMCE(i) = cos ( E )
TPMCE(i) = cos ( E )
TPMSE(i) = sin ( E )
E = E + TPMDE
end do
@ -444,7 +444,7 @@ contains !**********************************************************************
do j = 0, TPMNX1
if ( ( i .ge. TPMSSTPNH ) .or. ( j .le. TPMSSTPNX ) ) then
call TPMSSTPIntegrator ( TPMSSTPG(i,j), TPMSSTPF(i,j), TPMSSTPH(i), TPMSSTPX(j) )
print '(2i5,a,e20.10,a,e20.10,a,e20.10,a,e20.10)', i, j, ' H=', TPMSSTPH(i), &
print '(2i5,a,e20.10,a,e20.10,a,e20.10,a,e20.10)', i, j, ' H=', TPMSSTPH(i), &
', X=', TPMSSTPX(j), ', Q=', TPMSSTPG(i,j), ', U=', TPMSSTPF(i,j)
end if
end do
@ -454,11 +454,11 @@ contains !**********************************************************************
call TPMSSTPRead ()
end if
call CreateSpline2Ext ( 3, 3, 3, 3, TPMNH, TPMSSTPNH, TPMNX, TPMSSTPNX, TPMNMAX, TPMSSTPH, TPMSSTPX, &
TPMSSTPF, TPMSSTPFxx, TPMSSTPFyy, TPMSSTPFxxyy, FF, MM, DD, K0, K1, K2 )
TPMSSTPF, TPMSSTPFxx, TPMSSTPFyy, TPMSSTPFxxyy, FF, MM, DD, K0, K1, K2 )
end subroutine TPMSSTPInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! STP potential for an infinite tube interacting with a parallel segment. No actual initialization
! STP potential for an infinite tube interacting with a parallel segment. No actual initialization
! is necessary for this potential, since the data are taken from the table for the SSTP potential.
!---------------------------------------------------------------------------------------------------
@ -473,7 +473,7 @@ contains !**********************************************************************
!-------------------------------------------------------------------------------------------
i = 1 + int ( H / TPMDH )
if ( i < TPMSSTPNH ) then
!call PrintTPErrMsg ()
!call PrintTPErrMsg ()
!!call PrintStdLogMsg ( TPErrMsg )
!write ( TPErrMsg, '(a,e20.10)' ) 'Tubes intersect each other: H=', H
!call Error ( 'TPMSTPInt0', TPErrMsg )
@ -485,13 +485,13 @@ contains !**********************************************************************
return
end if
if ( i == TPMNH ) i = TPMNH - 1
Q = CalcLinFun1_0 ( i, H, TPMNH, TPMSSTPH, TPMSTPG )
U = CalcSpline1_0 ( i, H, TPMNH, TPMSSTPH, TPMSTPF, TPMSTPFxx )
Q = CalcLinFun1_0 ( i, H, TPMNH, TPMSSTPH, TPMSTPG )
U = CalcSpline1_0 ( i, H, TPMNH, TPMSSTPH, TPMSTPF, TPMSTPFxx )
TPMSTPInt0 = 1
end function TPMSTPInt0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMSTPInt1 ( Q, U, dUdH, H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function returns the transfer function Q, potential U, and derivative dUdH for
! This function returns the transfer function Q, potential U, and derivative dUdH for
! the STP potential calculated by interpolation in the table.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Q, U, dUdH
@ -500,7 +500,7 @@ contains !**********************************************************************
!-------------------------------------------------------------------------------------------
i = 1 + int ( H / TPMDH )
if ( i < TPMSSTPNH ) then
!call PrintTPErrMsg ()
!call PrintTPErrMsg ()
!!call PrintStdLogMsg ( TPErrMsg )
!write ( TPErrMsg, '(a,e20.10)' ) 'Tubes intersect each other: H=', H
!call Error ( 'TPMSTPInt0', TPErrMsg )
@ -512,11 +512,11 @@ contains !**********************************************************************
TPMSTPInt1 = 0
return
end if
Q = CalcLinFun1_0 ( i, H, TPMNH, TPMSSTPH, TPMSTPG )
call CalcSpline1_1 ( U, dUdH, i, H, TPMNH, TPMSSTPH, TPMSTPF, TPMSTPFxx )
Q = CalcLinFun1_0 ( i, H, TPMNH, TPMSSTPH, TPMSTPG )
call CalcSpline1_1 ( U, dUdH, i, H, TPMNH, TPMSSTPH, TPMSTPF, TPMSTPFxx )
TPMSTPInt1 = 1
end function TPMSTPInt1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPMSTPInit () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function initializes the table of the STP potential
!-------------------------------------------------------------------------------------------
@ -524,7 +524,7 @@ contains !**********************************************************************
TPMSTPF(0:TPMNH1) = TPMSSTPF(0:TPMNH1,TPMNX1)
TPMSTPFxx(0:TPMNH1) = TPMSSTPFxx(0:TPMNH1,TPMNX1)
end subroutine TPMSTPInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Fitting functions for the SST and ST potentials.
! This correction functions are chosen empirically to improve accuracy of the SST and ST potentials.
@ -567,14 +567,14 @@ contains !**********************************************************************
X = - ( X1_2 - X1_1 ) / 2.0d+00
do j = 0, ( TPTNXMAX - 1 ) / 2
HH = sqrt ( TPMAH(i) * TPMAH(i) + sqr ( X * sin ( M_PI_2 ) ) )
IntSign = TPMSTPInt0 ( Qb, Ub, HH )
IntSign = TPMSTPInt0 ( Qb, Ub, HH )
call TPTSetSegPosition2 ( TPTSeg1, R1_1, R1_2 )
call TPTSetSegPosition2 ( TPTSeg2, R2_1, R2_2 )
IntSign = TPTSectionPotential ( Qa, Ua, Fa, Ma, TPTSeg1, j, TPTSeg2 )
IntSign = TPTSectionPotential ( Qa, Ua, Fa, Ma, TPTSeg1, j, TPTSeg2 )
if ( j == 0 ) then
Uamin = Ua
Ubmin = Ub
else
else
if ( Ua < Uamin ) then
Uamin = Ua
end if
@ -614,7 +614,7 @@ contains !**********************************************************************
return
end if
i = 1 + int ( ( H - TPMAHmin ) / TPMADH )
TPMA0 = CalcSpline1_0 ( i, H, TPMAN, TPMAH, TPMAF, TPMAFxx )
TPMA0 = CalcSpline1_0 ( i, H, TPMAN, TPMAH, TPMAF, TPMAFxx )
end function TPMA0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPMA1 ( A, Ah, H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@ -646,9 +646,9 @@ contains !**********************************************************************
i = 1 + int ( ( H - TPMAHmin ) / TPMADH )
call CalcSpline1_1 ( A, Ah, i, H, TPMAN, TPMAH, TPMAF, TPMAFxx )
end subroutine TPMA1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function TPMCu0 ( H, cosA, sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function returns the correction function for the magnitude of the potential.
! This function returns the correction function for the magnitude of the potential.
!-------------------------------------------------------------------------------------------
real(c_double), intent(in) :: H, cosA, sinA
!-------------------------------------------------------------------------------------------
@ -656,7 +656,7 @@ contains !**********************************************************************
end function TPMCu0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPMCu1 ( Cu, CuH, CuA, H, cosA, sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The subroutine calculates the correction function Cu for the magnitude of the potential and
! The subroutine calculates the correction function Cu for the magnitude of the potential and
! its derivatives CuH and CuA.
!-------------------------------------------------------------------------------------------
real(c_double), intent(ouT) :: Cu, CuH, CuA
@ -664,7 +664,7 @@ contains !**********************************************************************
real(c_double) :: AA, AAh, D
!-------------------------------------------------------------------------------------------
call TPMA1 ( AA, AAh, H )
D = sqr ( sinA )
D = sqr ( sinA )
AA = AA - 1.0d+00
Cu = 1.0d+00 + AA * D
CuH = AAh * D
@ -672,7 +672,7 @@ contains !**********************************************************************
end subroutine TPMCu1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function TPMCa0 ( cosA, sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function returns the correction function for the argument of the potential.
! This function returns the correction function for the argument of the potential.
! If correction is not necessary, it should return sinA.
!-------------------------------------------------------------------------------------------
real(c_double), intent(in) :: cosA, sinA
@ -681,21 +681,21 @@ contains !**********************************************************************
end function TPMCa0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPMCa1 ( Ca, CaA, Ka, KaA, cosA, sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This subroutine calculates the correction function Cu for the depth of the potential well
! and its derivatives CuH and CuA. If correction is not necessary, it returns Ca = sinA
! This subroutine calculates the correction function Cu for the depth of the potential well
! and its derivatives CuH and CuA. If correction is not necessary, it returns Ca = sinA
! and CaA = cosA.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Ca, CaA, Ka, KaA
real(c_double), intent(in) :: cosA, sinA
!-------------------------------------------------------------------------------------------
Ka = 1.0d+00 / ( 1.0d+00 - TPMCaA * sqr ( sinA ) )
Ca = sinA * Ka
Ca = sinA * Ka
KaA = 2.0d+00 * TPMCaA * sinA * cosA * sqr ( Ka )
CaA = cosA * Ka + sinA * KaA
end subroutine TPMCa1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function TPMCe0 ( sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function returns the correction function for the argument of the potential.
! This function returns the correction function for the argument of the potential.
! If correction is not necessary, it returns sinA.
!-------------------------------------------------------------------------------------------
real(c_double), intent(in) :: sinA
@ -713,19 +713,19 @@ contains !**********************************************************************
CeA = - 2.0d+00 * TPMCeA * sinA * cosA
Ke = TPMCeA
end subroutine TPMCe1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! SST potential for the semi-infinite tube interacting with segment.
! This potential does not need any initialization. All necessary data is taken from tables of the
! SST potential for the semi-infinite tube interacting with segment.
! This potential does not need any initialization. All necessary data is taken from tables of the
! SSTP potential.
!---------------------------------------------------------------------------------------------------
integer(c_int) function TPMSSTPotential ( Q, U, X1, X2, H, cosA, D, N ) !!!!!!!!!!!!!!!!!!!!
! This function calculates the transfer function Q and potential U applied to a segment
! from a semi-infinite tube based on the numerical integration (trapezoidal rule) along the segment
! axis for non-parallel objects.
! This function calculates the transfer function Q and potential U applied to a segment
! from a semi-infinite tube based on the numerical integration (trapezoidal rule) along the segment
! axis for non-parallel objects.
! Relative position of the nanotube and segment is given by axial positions of the segment
! ends X1 and X2, height H, cosA= cos(A), where A is the cross-axis angle, and the displacement
! ends X1 and X2, height H, cosA= cos(A), where A is the cross-axis angle, and the displacement
! D of a nanotube end.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Q, U
@ -761,7 +761,7 @@ contains !**********************************************************************
end function TPMSSTPotential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMSSTPotentialPar ( Q, U, R1_1, Laxis1, R2_1, Laxis2, L1, N ) !!!!!
! Potential for a segment and a semi-infinite tube is calculated by the numerical
! Potential for a segment and a semi-infinite tube is calculated by the numerical
! integration (trapezoidal rule) along the segment axis for parallel objects.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Q, U
@ -799,9 +799,9 @@ contains !**********************************************************************
Q = Q * DX
U = U * DX
end function TPMSSTPotentialPar !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMSSTForces ( Q, U, F1, F2, Fd, X1, X2, H, cosA, D, N ) !!!!!!!!!!!
! Potential and forces applied to the segment from a semi-infinite tube are calculated
! Potential and forces applied to the segment from a semi-infinite tube are calculated
! by the numerical integration (trapezoidal rule) along the segment axis.
! Non-parallel case.
!-------------------------------------------------------------------------------------------
@ -859,7 +859,7 @@ contains !**********************************************************************
X = X + DX
end do
if ( TPMSSTForces == 0 ) return
C = DX * Cu
I0 = I0 * C
Ih = Ih * C
@ -868,7 +868,7 @@ contains !**********************************************************************
Ix = Ix * C
Ix1 = Ix1 * C
DX = X2 - X1
DX = X2 - X1
Q = Q * C
U = I0
@ -877,23 +877,23 @@ contains !**********************************************************************
C1 = Ka * Ka * Ih1
C = h * ( C1 + cosA * Ke * Ix ) / DX
F1(0) = X2 * Uh - C
F1(0) = X2 * Uh - C
F2(0) = C - X1 * Uh
C = ( cosA * sinA * C1 + ( Ke * sinA - sinA ) * Ix ) / DX
F1(1) = Ua - X2 * C
F1(1) = Ua - X2 * C
F2(1) = X1 * C - Ua
C1 = Ca * Ca * Ih1 + cosA * Ix
C2 = Ca * Ca * Ih2 + cosA * Ix1
F1(2) = ( U - X2 * C1 + C2 ) / DX
F2(2) = ( X1 * C1 - C2 - U ) / DX
Fd = Ce * Ix
end function TPMSSTForces !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMSSTForcesPar ( Q, U, F1, F2, Fd, R1_1, Laxis1, R2_1, Laxis2, L1, N )
! Potential and forces applied to the segment from a semi-infinite tube are calculated by
integer(c_int) function TPMSSTForcesPar ( Q, U, F1, F2, Fd, R1_1, Laxis1, R2_1, Laxis2, L1, N )
! Potential and forces applied to the segment from a semi-infinite tube are calculated by
! the numerical integration (trapezoidal rule) along the segment axis.
! Parallel case
!-------------------------------------------------------------------------------------------
@ -942,8 +942,8 @@ contains !**********************************************************************
U = U + Us
Fd = Fd + Usx
end if
F1 = F1 + Gamma * Fs
F2 = F2 + Beta * Fs
F1 = F1 + Gamma * Fs
F2 = F2 + Beta * Fs
end if
X = X + DX
end do
@ -954,15 +954,15 @@ contains !**********************************************************************
F1 = DX * F1 + Fs
F2 = DX * F2 - Fs
end function TPMSSTForcesPar !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! ST: Potential for a infinite tube interacting with a segment
!--------------------------------------------------------------------------------------------------
!
! These functions are used to smooth boundaries in (H,X) domain for ST potential
!
real(c_double) function TPMSTXMin0 ( H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: H
!-------------------------------------------------------------------------------------------
@ -979,7 +979,7 @@ contains !**********************************************************************
TPMSTXMin0 = sqrt ( TPMSTH2 * TPMSTH2 - H * H ) &
* ( 1.0d+00 - X * X * X * ( 3.0d+00 * X * ( 2.0d+00 * X - 5.0d+00 ) + 10.0d+00 ) )
end function TPMSTXMin0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function TPMSTXMax0 ( H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: H
!-------------------------------------------------------------------------------------------
@ -1005,11 +1005,11 @@ contains !**********************************************************************
F = 1.0d+00 - X * X * X * ( 3.0d+00 * X * ( 2.0d+00 * X - 5.0d+00 ) + 10.0d+00 )
X = X * ( X - 1.0d+00 )
dFdX = - 30.0d+00 * X * X
XMin = sqrt ( TPMSTH2 * TPMSTH2 - H * H )
XMin = sqrt ( TPMSTH2 * TPMSTH2 - H * H )
dXMindH = - H * F / XMin + XMin * dFdX / TPMSTDH12
XMin = XMin * F
end subroutine TPMSTXMin1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPMSTXMax1 ( XMax, dXMaxdH, H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: XMax, dXMaxdH
real(c_double), intent(in) :: H
@ -1017,11 +1017,11 @@ contains !**********************************************************************
XMax = sqrt ( TPMSTXMax * TPMSTXMax - H * H )
dXMaxdH = - H / XMax
end subroutine TPMSTXMax1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! ST (segment-tube) table
!
subroutine TPMSTIntegrator ( G, F, Q, U, H, X, DX ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(inout) :: G, F, Q, U
real(c_double), intent(in) :: H, X, DX
@ -1036,7 +1036,7 @@ contains !**********************************************************************
U = 0.0d+00
HH = dsqrt ( sqr ( H ) + sqr ( X ) )
if ( HH > TPMHmax ) return
IntSign = TPMSTPInt0 ( Q, U, HH )
IntSign = TPMSTPInt0 ( Q, U, HH )
if ( IntSign == 1 ) then
G = G + Q * DDX
F = F + U * DDX
@ -1059,10 +1059,10 @@ contains !**********************************************************************
G = 0.0d+00
F = 0.0d+00
TPMSTInt0 = 2
!call PrintTPErrMsg ()
!call PrintTPErrMsg ()
!!call PrintStdLogMsg ( TPErrMsg )
!all Error ( 'TPMSTInt0', 'H < 0' )
!!return
!!return
end if
S = signum ( X )
XA = dabs ( X )
@ -1074,7 +1074,7 @@ contains !**********************************************************************
if ( XXX < 0.0d+00 ) then
j = 1
XXXX = 0.0d+00
!call PrintTPErrMsg ()
!call PrintTPErrMsg ()
!write ( TPErrMsg, '(a,e20.10,a,e20.10,a,e20.10)' ) 'X < XMin, X=', XA, ', XMin=', XMin, ', H=', H
!call Error ( 'TPMSTInt0', TPErrMsg )
else if ( XXX .ge. 1.0d+00 ) then
@ -1084,8 +1084,8 @@ contains !**********************************************************************
XXXX = XXX
j = 1 + int ( XXXX * TPMNX1 )
end if
G = S * CalcLinFun2_0 ( i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, TPMSTG )
F = S * CalcSpline2_0 ( i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, TPMSTF, TPMSTFxx, TPMSTFyy, TPMSTFxxyy )
G = S * CalcLinFun2_0 ( i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, TPMSTG )
F = S * CalcSpline2_0 ( i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, TPMSTF, TPMSTFxx, TPMSTFyy, TPMSTFxxyy )
TPMSTInt0 = 1
end function TPMSTInt0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@ -1109,10 +1109,10 @@ contains !**********************************************************************
Fh = 0.0d+00
Fx = 0.0d+00
TPMSTInt1 = 2
!call PrintTPErrMsg ()
!call PrintTPErrMsg ()
!!call PrintStdLogMsg ( trim ( TPErrMsg ) )
!call Error ( 'TPMSTInt1', 'H < 0' )
!!return
!!return
end if
S = signum ( X )
XA = dabs ( X )
@ -1126,7 +1126,7 @@ contains !**********************************************************************
j = 1
XXX = 0.0d+00
XXXX = 0.0d+00
!call PrintTPErrMsg ()
!call PrintTPErrMsg ()
!write ( TPErrMsg, '(a,e20.10,a,e20.10,a,e20.10)' ) 'X < XMin, X=', XA, ', XMin=', XMin, ', H=', H
!call Error ( 'TPMSTInt', TPErrMsg )
else if ( XXX .ge. 1.0d+00 ) then
@ -1134,19 +1134,19 @@ contains !**********************************************************************
XXX = 1.0d+00
XXXX = 1.0d+00
else
XXXX = XXX
XXXX = XXX
j = 1 + int ( XXXX * TPMNX1 )
end if
G = S * CalcLinFun2_0 ( i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, TPMSTG )
G = S * CalcLinFun2_0 ( i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, TPMSTG )
call CalcSpline2_1 ( F, Fh, Fx, i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, &
TPMSTF, TPMSTFxx, TPMSTFyy, TPMSTFxxyy )
TPMSTF, TPMSTFxx, TPMSTFyy, TPMSTFxxyy )
Fx = Fx / DX
Fh = Fh - Fx * ( dXMaxdH * XXX + dXMindH * ( 1.0d+00 - XXX ) )
F = F * S
Fh = Fh * S
TPMSTInt1 = 1
end function TPMSTInt1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMSTPotential ( Q, U, X1, X2, H, cosA, CaseID ) !!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: Q, U
real(c_double), intent(in) :: X1, X2, H, cosA
@ -1157,7 +1157,7 @@ contains !**********************************************************************
if ( CaseID == MD_LINES_PAR ) then
TPMSTPotential = TPMSTPInt0 ( Q, U, H )
U = U * ( X2 - X1 )
return
return
end if
TPMSTPotential = 0
sinA = dsqrt ( 1.0d+00 - cosA * cosA )
@ -1192,12 +1192,12 @@ contains !**********************************************************************
F2(2) = - U
Q = Q * DX
U = U * DX
return
return
end if
sinA = dsqrt ( 1.0d+00 - cosA * cosA )
call TPMCa1 ( Ca, CaA, Ka, KaA, cosA, sinA )
IntSign1 = TPMSTInt1 ( GG1, FF1, Fh1, Fx1, H, X1 * Ca )
IntSign1 = TPMSTInt1 ( GG1, FF1, Fh1, Fx1, H, X1 * Ca )
IntSign2 = TPMSTInt1 ( GG2, FF2, Fh2, Fx2, H, X2 * Ca )
if ( ( IntSign1 .ne. 1 ) .and. ( IntSign2 .ne. 1 ) ) then
Q = 0.0d+00
@ -1205,11 +1205,11 @@ contains !**********************************************************************
F1 = 0.0d+00
F2 = 0.0d+00
TPMSTForces = 0
return
return
end if
call TPMCu1 ( Cu, CuH, CuA, H, cosA, sinA )
Q = Cu * ( GG2 - GG1 ) / Ca
U = Cu * ( FF2 - FF1 ) / Ca
@ -1218,19 +1218,19 @@ contains !**********************************************************************
D = CuH * U / Cu + Cu * ( Fh2 - Fh1 ) / Ca
F1(0) = ( X2 * D - C ) / DX
F2(0) = ( C - X1 * D ) / DX
C = cosA * B
D = ( CuA / Cu - CaA / Ca ) * U + CaA * Cu * ( X2 * Fx2 - X1 * Fx1 ) / Ca
F1(1) = ( D - X2 * C ) / DX
F2(1) = ( X1 * C - D ) / DX
F1(2) = Cu * Fx1
F2(2) = - Cu * Fx2
TPMSTForces = 1
end function TPMSTForces !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMSTForceTorque( Qi, Ui, Fi, Ti, Q, U, F, T, Psi, PsiA, Cap, L, H, cosA, CaseID )
integer(c_int) function TPMSTForceTorque( Qi, Ui, Fi, Ti, Q, U, F, T, Psi, PsiA, Cap, L, H, cosA, CaseID )
real(c_double), intent(out) :: Qi, Ui, Fi, Ti, Q, U, F, T, Psi, PsiA, Cap
real(c_double), intent(in) :: L, H, cosA
integer(c_int), intent(in) :: CaseID
@ -1254,14 +1254,14 @@ contains !**********************************************************************
Psi = 0.0d+00
PsiA = 0.0d+00
Cap = 0.0d+00
return
return
end if
L2 = 0.5d+00 * L
sinA = dsqrt ( 1.0d+00 - cosA * cosA )
call TPMCa1 ( Ca, CaA, Ka, KaA, cosA, sinA )
IntSign = TPMSTInt1 ( GG, FF, Fh, Fx, H, L2 * Ca )
IntSign = TPMSTInt1 ( GGi, FFi, Fhi, Fxi, H, TPMSTXmax )
IntSign = TPMSTInt1 ( GG, FF, Fh, Fx, H, L2 * Ca )
IntSign = TPMSTInt1 ( GGi, FFi, Fhi, Fxi, H, TPMSTXmax )
if ( IntSign .ne. 1 ) then
Qi = 0.0d+00
Ui = 0.0d+00
@ -1275,11 +1275,11 @@ contains !**********************************************************************
PsiA = 0.0d+00
Cap = 0.0d+00
TPMSTForceTorque = 0
return
return
end if
call TPMCu1 ( Cu, CuH, CuA, H, cosA, sinA )
Psi = Cu / Ca
PsiA = ( CuA * Ca - Cu * CaA ) / Ca / Ca
Cap = CuA / Cu - KaA / Ka - cosA / sinA
@ -1287,16 +1287,16 @@ contains !**********************************************************************
Ui = 2.0d+00 * Psi * FFi
Fi = - 2.0d+00 * ( CuH * FFi / Ca + Psi * Fhi )
Ti = - Cap * Ui
Q = 2.0d+00 * Cu * GG / Ca
U = 2.0d+00 * Cu * FF / Ca
F = - 2.0d+00 * ( CuH * FF / Ca + Psi * Fh )
T = - 2.0d+00 * ( ( CuA * Ka - Cu * KaA ) / ( Ka * Ka * sinA ) - Cu * cosA / ( Ka * sinA * sinA ) ) * FF &
- 2.0d+00 * Cu / ( Ka * sinA ) * Fx * L2 * ( KaA * sinA + Ka * cosA )
TPMSTForceTorque = 1
end function TPMSTForceTorque !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPMSTInit () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) :: X, Q, U, DX, DDX, XMin, XMax
integer(c_int) :: i, j, k
@ -1307,7 +1307,7 @@ contains !**********************************************************************
TPMSTDH12 = TPMSTH2 - TPMSTH1
TPMSTXmax = TPMHMax + TPMSTDelta
DX = 1.0 / TPMNX1
do j = 0, TPMNX1
do j = 0, TPMNX1
TPMSTX(j) = DX * j
end do
do i = 0, TPMNH1
@ -1323,7 +1323,7 @@ contains !**********************************************************************
TPMSTF(i,0) = 0.0d+00
TPMSTFyy(i,0) = U
TPMSTFyy(i,TPMNX1) = 0.0d+00
do j = 1, TPMNX1
do j = 1, TPMNX1
TPMSTG(i,j) = TPMSTG(i,j-1)
TPMSTF(i,j) = TPMSTF(i,j-1)
do k = 0, TPMSTNXS - 1
@ -1334,18 +1334,18 @@ contains !**********************************************************************
end do
end do
call CreateSpline2 ( 3, 3, 3, 3, TPMNH, TPMNX, TPMNMAX, TPMSTH, TPMSTX, TPMSTF, TPMSTFxx, &
TPMSTFyy, TPMSTFxxyy, FF, MM, DD, K0, K1, K2 )
TPMSTFyy, TPMSTFxxyy, FF, MM, DD, K0, K1, K2 )
end subroutine TPMSTInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Interaction functions: They can be used for calculation of the potential and forces between a
! Interaction functions: They can be used for calculation of the potential and forces between a
! segment and an infinite or semi-infinite nanotube.
!---------------------------------------------------------------------------------------------------
subroutine TPMSegmentForces ( F2_1, F2_2, F1_1, F1_2, R1_1, R1_2, R2, Laxis2, L2 ) !!!!!!!!!
real(c_double), dimension(0:2), intent(out) :: F2_1, F2_2
real(c_double), dimension(0:2), intent(in) :: F1_1, F1_2, R1_1, R1_2, R2, Laxis2
real(c_double), intent(in) :: L2
real(c_double), intent(in) :: L2
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2) :: F, M, RR
!-------------------------------------------------------------------------------------------
@ -1389,11 +1389,11 @@ contains !**********************************************************************
L2 = S_V3norm3 ( Laxis2 )
Laxis1 = Laxis1 / L1
Laxis2 = Laxis2 / L2
L1 = 0.5d+00 * L1
L1 = 0.5d+00 * L1
L2 = 0.5d+00 * L2
if ( SType2 == 2 ) Laxis2 = - Laxis2
GeomID = LineLine ( H, cosA, D1, D2, L12, R1, Laxis1, R2, Laxis2, TPGeomPrec )
! Angle switch
if ( TPMASwitch == 0 ) then
if ( GeomID == MD_LINES_PAR ) then
@ -1414,7 +1414,7 @@ contains !**********************************************************************
SwitchID = 1
end if
end if
if ( SwitchID < 2 ) then
D2 = D2 - L2
if ( SType2 == 0 ) then
@ -1429,14 +1429,14 @@ contains !**********************************************************************
Ly = Ly * S
if ( IntSigna > 0 ) then
F1_1a = F1(0) * L12 + F1(1) * Ly + F1(2) * Laxis1
F1_2a = F2(0) * L12 + F2(1) * Ly + F2(2) * Laxis1
F1_2a = F2(0) * L12 + F2(1) * Ly + F2(2) * Laxis1
else
F1_1a = 0.0d+00
F1_2a = 0.0d+00
end if
end if
if ( SwitchID > 0 ) then
if ( SwitchID > 0 ) then
if ( SType2 == 0 ) then
call LinePoint ( H, L12, R2, Laxis2, R1 )
L12 = L12 - R1
@ -1451,7 +1451,7 @@ contains !**********************************************************************
else
L12 = L12 / H
F1_1b = F1(0) * L12 + F1(2) * Laxis1
F1_2b = F2(0) * L12 + F2(2) * Laxis1
F1_2b = F2(0) * L12 + F2(2) * Laxis1
end if
else
F1_1b = 0.0d+00
@ -1480,7 +1480,7 @@ contains !**********************************************************************
F1_2 = F1_2b
Fd = Fdb
TPMInteractionF = IntSignb
else
else
W1 = 1.0d+00 - W
Q = W * Qa + W1 * Qb
U = W * Ua + W1 * Ub
@ -1490,7 +1490,7 @@ contains !**********************************************************************
Fd = W * Fda + W1 * Fdb
TPMInteractionF = 0
if ( IntSigna > 0 .or. IntSignb > 0 ) TPMInteractionF = 1
end if
end if
! Calculation of forces for the complimentary tube
if ( SType2 == 2 ) Laxis2 = - Laxis2
@ -1510,7 +1510,7 @@ contains !**********************************************************************
F2_1 = F2_1 - DR
end if
end function TPMInteractionF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPMInteractionU ( Q, U, R1_1, R1_2, R2_1, R2_2, SType2 ) !!!!!!!!!!!
real(c_double), intent(inout) :: Q, U
real(c_double), dimension(0:2), intent(in) :: R1_1, R1_2, R2_1, R2_2
@ -1560,7 +1560,7 @@ contains !**********************************************************************
IntSigna = TPMSSTPotential ( Qa, Ua, D1 - L1, D1 + L1, H, cosA, D2 - L2, TPMNN )
end if
end if
if ( SwitchID > 0 ) then
if ( Stype2 == 0 ) then
call LinePoint ( H, L12, R2, Laxis2, R1 )
@ -1573,7 +1573,7 @@ contains !**********************************************************************
IntSignb = TPMSSTPotentialPar ( Qb, Ub, R1_1, Laxis1, R2_2, Laxis2, 2.0d+00 * L1, TPMNN )
end if
end if
if ( SwitchID == 0 ) then
Q = Qa
U = Ua
@ -1612,14 +1612,14 @@ contains !**********************************************************************
D2 = 2.0d+00 * Delta
do i = 0, 2
DD = - Delta
do j = 0 , 1
RR = R1_1
do j = 0 , 1
RR = R1_1
RR(i) = RR(i) + DD
IntSign = TPMInteractionU ( QQ, U1_1(j,i), RR, R1_2, R2_1, R2_2, SType2 )
RR = R1_2
RR = R1_2
RR(i) = RR(i) + DD
IntSign = TPMInteractionU ( QQ, U1_2(j,i), R1_1, RR, R2_1, R2_2, SType2 )
RR = R2_1
RR = R2_1
RR(i) = RR(i) + DD;
IntSign = TPMInteractionU ( QQ, U2_1(j,i), R1_1, R1_2, RR, R2_2, SType2 )
RR = R2_2
@ -1635,7 +1635,7 @@ contains !**********************************************************************
F2_2(i) = ( U2_2(0,i) - U2_2(1,i) ) / D2
end do
end function TPMInteractionFNum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Initialization
!---------------------------------------------------------------------------------------------------
@ -1646,7 +1646,7 @@ contains !**********************************************************************
character(len=512) :: PDate
!-------------------------------------------------------------------------------------------
TPPotType = TP_POT_MONO_R
! Here we calculate the radius of nanotubes
RT = TPBAcc * sqrt ( 3.0d+00 * ( ChiIndM * ChiIndM + ChiIndN * ChiIndN + ChiIndM * ChiIndN ) ) / M_2PI
@ -1661,7 +1661,7 @@ contains !**********************************************************************
TPMHmax = TPMR1 + TPMR2 + TPBRcutoff
TPMDH = TPMHmax / TPMNH1
! Parameters of the angle switch
TPMASMin = sqr ( cos ( rad ( TPMAS ) ) )
TPMASMax = 1.0d+00 - TPGeomPrec
@ -1669,30 +1669,30 @@ contains !**********************************************************************
if ( TPMStartMode == 1 ) then
TPMUnitID = OpenFile ( TPMFile, 'rt', '' )
read ( unit = TPMUnitID, fmt = '()' )
read ( unit = TPMUnitID, fmt = '()' )
read ( unit = TPMUnitID, fmt = '()' )
read ( unit = TPMUnitID, fmt = '()' )
read ( unit = TPMUnitID, fmt = '()' )
read ( unit = TPMUnitID, fmt = '()' )
else
TPMUnitID = OpenFile ( TPMFile, 'wt', '' )
call fdate( PDate )
write ( unit = TPMUnitID, fmt = '(a,a)' ) 'DATE ', 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'
write ( unit = TPMUnitID, fmt = '(a)' ) &
'A. N. Volkov, L. V. Zhigilei, J. Phys. Chem. C 114, 5513-5531, 2010. doi: 10.1021/jp906142h'
end if
call TPMSSTPInit ()
call TPMSTPInit ()
DX = TPMR1 + TPMR2 + TPBRcutoff
call TPMAInit ( - DX, DX, - DX, DX )
call TPMAInit ( - DX, DX, - DX, DX )
call TPMSTInit ()
call CloseFile ( TPMUnitID )
end subroutine TPMInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TubePotMono !****************************************************************************

60
lib/mesont/TubePotTrue.f90
View File

@ -9,9 +9,9 @@
! the GNU General Public License.
!
! See the README file in the top-level LAMMPS directory.
!
!
! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
module TubePotTrue !********************************************************************************
!
@ -25,8 +25,8 @@ module TubePotTrue !************************************************************
!
!---------------------------------------------------------------------------------------------------
!
! This module implements calculation of the true potential and transfer functions for interaction
! between two cylinder segments of nanotubes by direct integration over the surfaces of both
! This module implements calculation of the true potential and transfer functions for interaction
! between two cylinder segments of nanotubes by direct integration over the surfaces of both
! segments.
!
!***************************************************************************************************
@ -55,15 +55,15 @@ implicit none
integer(c_int) :: NX, NE ! Number of nodes for numerical integration
real(c_double) :: DX, DE ! Spacings
real(c_double), dimension(0:2,0:2) :: M ! Transformation matrix
real(c_double), dimension(0:TPTNXMAX-1,0:TPTNXMAX-1,0:2) :: Rtab! Node coordinates
real(c_double), dimension(0:TPTNXMAX-1,0:TPTNXMAX-1,0:2) :: Rtab! Node coordinates
end type TPTSEG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Global variables
!---------------------------------------------------------------------------------------------------
type(TPTSEG) :: TPTSeg1, TPTSeg2 ! Two segments
contains !******************************************************************************************
subroutine TPTSegAxisVector ( S, Laxis ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@ -96,11 +96,11 @@ contains !**********************************************************************
!-------------------------------------------------------------------------------------------
call TPTSegAxisVector ( S, Laxis )
call TPTSegRadVector ( S, Lrad, Eps )
R(0) = S%X + X * Laxis(0) + S%R * Lrad(0)
R(0) = S%X + X * Laxis(0) + S%R * Lrad(0)
R(1) = S%Y + X * Laxis(1) + S%R * Lrad(1)
R(2) = S%Z + X * Laxis(2) + S%R * Lrad(2)
end subroutine TPTRadiusVector !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPTCalcSegNodeTable ( S ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
type(TPTSEG), intent(inout) :: S
!-------------------------------------------------------------------------------------------
@ -111,7 +111,7 @@ contains !**********************************************************************
call RotationMatrix3 ( S%M, S%Psi, S%Theta, S%Phi )
do i = 0, S%NX - 1
Eps = 0.0d+00
do j = 0, S%NE - 1
do j = 0, S%NE - 1
call TPTRadiusVector ( S, S%Rtab(i,j,0:2), X, Eps )
Eps = Eps + S%DE
end do
@ -133,12 +133,12 @@ contains !**********************************************************************
S%Z = Rcenter(2)
call TPTCalcSegNodeTable ( S )
end subroutine TPTSetSegPosition1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPTSetSegPosition2 ( S, R1, R2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
type(TPTSEG), intent(inout) :: S
real(c_double), dimension(0:2), intent(in) :: R1, R2
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2) :: R, Laxis
real(c_double), dimension(0:2) :: R, Laxis
real(c_double) :: L
!-------------------------------------------------------------------------------------------
R = 0.5 * ( R1 + R2 )
@ -147,7 +147,7 @@ contains !**********************************************************************
Laxis = Laxis / L
call TPTSetSegPosition1 ( S, R, Laxis, L )
end subroutine TPTSetSegPosition2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPTCheckIntersection ( S1, S2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
type(TPTSEG), intent(in) :: S1, S2
!-------------------------------------------------------------------------------------------
@ -193,7 +193,7 @@ contains !**********************************************************************
Xmin = 0.0d+00
Xmax = 0.0d+00
TPTCalcPointRange = 0
return
return
end if
Distance = sqrt ( TPBRcutoff * TPBRcutoff - Distance * Distance )
Xmin = Displacement - Distance
@ -218,13 +218,13 @@ contains !**********************************************************************
R2_2(1) = - X2_2 * sin ( A )
R2_2(2) = X2_2 * cos ( A )
end subroutine TPTGetEnds !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Tubular potential
!---------------------------------------------------------------------------------------------------
integer(c_int) function TPTPointPotential ( Q, U, F, R, S ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function returns the potential U and force F applied to an atom in position R and
! This function returns the potential U and force F applied to an atom in position R and
! produced by the segment S.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Q, U
@ -243,7 +243,7 @@ contains !**********************************************************************
F = 0.0d+00
if ( TPTCalcPointRange ( S, Xmin, Xmax, R ) == 0 ) return
X = - S%L / 2.0
do i = 0, S%NX - 1
do i = 0, S%NX - 1
if ( X > Xmin .and. X < Xmax ) then
QQ = 0.0d+00
UU = 0.0d+00
@ -264,7 +264,7 @@ contains !**********************************************************************
Q = Q + 0.5d+00 * QQ
U = U + 0.5d+00 * UU
F = F + 0.5d+00 * FF
else
else
Q = Q + QQ
U = U + UU
F = F + FF
@ -279,7 +279,7 @@ contains !**********************************************************************
end function TPTPointPotential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPTSectionPotential ( Q, U, F, M, S, i, Ssource ) !!!!!!!!!!!!!!!!!!
! This function returns the potential U, force F and torque M produced by the segment Ssource
! This function returns the potential U, force F and torque M produced by the segment Ssource
! and applied to the i-th circular cross-section of the segment S.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Q, U
@ -290,7 +290,7 @@ contains !**********************************************************************
integer(c_int) :: j
real(c_double), dimension(0:2) :: R, Fp, Mp, Lrad
real(c_double) :: Qp, Up, Eps
real(c_double) :: Coeff
real(c_double) :: Coeff
!-------------------------------------------------------------------------------------------
TPTSectionPotential = 0
Q = 0.0d+00
@ -321,7 +321,7 @@ contains !**********************************************************************
end function TPTSectionPotential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPTSegmentPotential ( Q, U, F, M, S, Ssource ) !!!!!!!!!!!!!!!!!!!!!!
! This function returns the potential U, force F and torque M produced by the segment
! This function returns the potential U, force F and torque M produced by the segment
! Ssource and applied to the segment S.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Q, U
@ -340,7 +340,7 @@ contains !**********************************************************************
TPTSegmentPotential = 2
return
end if
do i = 0, S%NX - 1
do i = 0, S%NX - 1
if ( TPTSectionPotential ( Qc, Uc, Fc, Mc, S, i, Ssource ) == 1 ) then
if ( i == 0 .or. i == S%NX - 1 ) then
Q = Q + 0.5d+00 * Qc
@ -355,17 +355,17 @@ contains !**********************************************************************
end if
TPTSegmentPotential = 1
end if
end do
end do
Q = Q * S%DX
U = U * S%DX
F = F * S%DX
M = M * S%DX
end function TPTSegmentPotential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Forces
!---------------------------------------------------------------------------------------------------
subroutine TPTSegmentForces ( F1, F2, F, M, Laxis, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), dimension(0:2), intent(out) :: F1, F2
real(c_double), dimension(0:2), intent(in) :: F, M, Laxis
@ -376,10 +376,10 @@ contains !**********************************************************************
FF = 0.5d+00 * F
MM = M / L
call V3_V3xxV3 ( FFF, MM, Laxis )
F1 = FF - FFF
F2 = FF + FFF
F1 = FF - FFF
F2 = FF + FFF
end subroutine TPTSegmentForces !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPTInteractionF ( Q, U, F1_1, F1_2, F2_1, F2_2, R1_1, R1_2, R2_1, R2_2 )
! This function returns the potential and forces applied to the ends of segments.
!-------------------------------------------------------------------------------------------
@ -413,7 +413,7 @@ contains !**********************************************************************
!---------------------------------------------------------------------------------------------------
! Initialization
!---------------------------------------------------------------------------------------------------
subroutine TPTInit ( R1, R2, NX, NE ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: R1, R2
integer(c_int), intent(in) :: NX, NE

1
tools/coding_standard/whitespace.py
View File

@ -51,6 +51,7 @@ patterns:
- "*.py"
- "*.rst"
- "*.sh"
- "*.f90"
- ".gitignore"
- "README"
- "requirements.txt"