From 61855c5058a0e88fb2a445fdf1ae549c8d57e894 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 25 Aug 2021 23:46:48 -0400 Subject: [PATCH] apply whitespace checking/fixing also to free-format Fortran files --- lib/mesont/CNTPot.f90 | 182 +++++++-------- lib/mesont/ExportCNT.f90 | 100 ++++---- lib/mesont/LinFun2.f90 | 20 +- lib/mesont/Spline1.f90 | 30 +-- lib/mesont/Spline2.f90 | 26 +-- lib/mesont/TPMForceField.f90 | 140 +++++------ lib/mesont/TPMGeom.f90 | 18 +- lib/mesont/TPMLib.f90 | 18 +- lib/mesont/TPMM0.f90 | 34 +-- lib/mesont/TPMM1.f90 | 62 ++--- lib/mesont/TubePotBase.f90 | 88 +++---- lib/mesont/TubePotMono.f90 | 348 ++++++++++++++-------------- lib/mesont/TubePotTrue.f90 | 60 ++--- tools/coding_standard/whitespace.py | 1 + 14 files changed, 564 insertions(+), 563 deletions(-) diff --git a/lib/mesont/CNTPot.f90 b/lib/mesont/CNTPot.f90 index 4a130074fd..5c83ce63e5 100644 --- a/lib/mesont/CNTPot.f90 +++ b/lib/mesont/CNTPot.f90 @@ -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 !********************************************************************************* diff --git a/lib/mesont/ExportCNT.f90 b/lib/mesont/ExportCNT.f90 index 20166a956e..039503d51e 100644 --- a/lib/mesont/ExportCNT.f90 +++ b/lib/mesont/ExportCNT.f90 @@ -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 diff --git a/lib/mesont/LinFun2.f90 b/lib/mesont/LinFun2.f90 index 8d01a25eda..d0c10088e9 100644 --- a/lib/mesont/LinFun2.f90 +++ b/lib/mesont/LinFun2.f90 @@ -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 !******************************************************************************** diff --git a/lib/mesont/Spline1.f90 b/lib/mesont/Spline1.f90 index bd796e35cc..95d7317049 100644 --- a/lib/mesont/Spline1.f90 +++ b/lib/mesont/Spline1.f90 @@ -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 !******************************************************************************** diff --git a/lib/mesont/Spline2.f90 b/lib/mesont/Spline2.f90 index b449f25c07..2dd7fb131b 100644 --- a/lib/mesont/Spline2.f90 +++ b/lib/mesont/Spline2.f90 @@ -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 !******************************************************************************** diff --git a/lib/mesont/TPMForceField.f90 b/lib/mesont/TPMForceField.f90 index a7312e09f8..5fa0514eea 100644 --- a/lib/mesont/TPMForceField.f90 +++ b/lib/mesont/TPMForceField.f90 @@ -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 diff --git a/lib/mesont/TPMGeom.f90 b/lib/mesont/TPMGeom.f90 index c5b8fcfaa7..638f2cd27e 100644 --- a/lib/mesont/TPMGeom.f90 +++ b/lib/mesont/TPMGeom.f90 @@ -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 !******************************************************************************** diff --git a/lib/mesont/TPMLib.f90 b/lib/mesont/TPMLib.f90 index 3a186917c0..c3c9f5eb4d 100644 --- a/lib/mesont/TPMLib.f90 +++ b/lib/mesont/TPMLib.f90 @@ -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 !********************************************************************************* diff --git a/lib/mesont/TPMM0.f90 b/lib/mesont/TPMM0.f90 index bd2e598b69..d2b38fe558 100644 --- a/lib/mesont/TPMM0.f90 +++ b/lib/mesont/TPMM0.f90 @@ -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 !********************************************************************************** diff --git a/lib/mesont/TPMM1.f90 b/lib/mesont/TPMM1.f90 index d6c8a193db..d8186683ca 100644 --- a/lib/mesont/TPMM1.f90 +++ b/lib/mesont/TPMM1.f90 @@ -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 !********************************************************************************** diff --git a/lib/mesont/TubePotBase.f90 b/lib/mesont/TubePotBase.f90 index dbb15c1457..d278a2390f 100644 --- a/lib/mesont/TubePotBase.f90 +++ b/lib/mesont/TubePotBase.f90 @@ -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 diff --git a/lib/mesont/TubePotMono.f90 b/lib/mesont/TubePotMono.f90 index 8c19043a67..fc66d15c65 100644 --- a/lib/mesont/TubePotMono.f90 +++ b/lib/mesont/TubePotMono.f90 @@ -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 !**************************************************************************** diff --git a/lib/mesont/TubePotTrue.f90 b/lib/mesont/TubePotTrue.f90 index 2cb300cc32..98e7450a99 100644 --- a/lib/mesont/TubePotTrue.f90 +++ b/lib/mesont/TubePotTrue.f90 @@ -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 diff --git a/tools/coding_standard/whitespace.py b/tools/coding_standard/whitespace.py index 2f84e1aa6d..18fb3308df 100644 --- a/tools/coding_standard/whitespace.py +++ b/tools/coding_standard/whitespace.py @@ -51,6 +51,7 @@ patterns: - "*.py" - "*.rst" - "*.sh" + - "*.f90" - ".gitignore" - "README" - "requirements.txt"