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

fix code related issues

Browse Source
This commit is contained in:
iafoss
2020-03-20 14:21:24 -04:00
parent bdbbe00ec7
commit 53ccc4c607
98 changed files with 4210 additions and 3657 deletions
Download Patch File Download Diff File Expand all files Collapse all files

734
lib/mesont/CNTPot.f90 Normal file
View File

@ -0,0 +1,734 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !*************************************************************************************
!
! TMD Library: Mesoscopic potential for internal modes in CNTs
!
!---------------------------------------------------------------------------------------------------
!
! Implementation of 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 Youngs 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
! CNTTRS, torsion potential
! CNTBRT, breathing potential
!
! The functional form and force constants of harmonic streatching, bending and
! 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 non-harmonic potential with fracture
! is developed and parameterized with the help of constant
! -- Young's modulus (Pa),
! -- maximal linear strain (only for the NH potential of type 1)
! -- tensile strength (or fracture strain, Pa),
! -- strain at failure (or fracture strain)
! -- maximal strain.
! All these parameters are assumed to be independent of SWCNT radius or type.
! In this model true strain at failure CNTSTREft and true tensile strength
! CNTSTRSft are slightly different from imposed values CNTSTREf and CNTSTRSf.
! This difference is really small and is not taken into account.
!
! 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, see subroutine CNTSTRSetParameterization
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 08.02.m.m.2.m, 2017
!
!***************************************************************************************************
use TPMLib
use iso_c_binding, only : c_int, c_double, c_char
implicit none
!---------------------------------------------------------------------------------------------------
! Constants
!---------------------------------------------------------------------------------------------------
integer(c_int), parameter :: CNTPOT_STRETCHING = 0
integer(c_int), parameter :: CNTPOT_SBUCKLING = 1
integer(c_int), parameter :: CNTPOT_SFRACTURE = 2
integer(c_int), parameter :: CNTPOT_BENDING = 3
integer(c_int), parameter :: CNTPOT_BBUCKLING = 4
integer(c_int), parameter :: CNTPOT_BFRACTURE = 5
integer(c_int), parameter :: CNTSTRMODEL_H0 = 0 ! Harmonic stetching model (constant Young's modulus)
integer(c_int), parameter :: CNTSTRMODEL_H1 = 1 ! Harmonic stretching model (Young's modulus depends on radius)
integer(c_int), parameter :: CNTSTRMODEL_NH0F = 2 ! Non-harmonic stretching with fracture, potential of type 0
integer(c_int), parameter :: CNTSTRMODEL_NH1 = 3 ! Non-harmonic stretching without fracture, potential of type 1
integer(c_int), parameter :: CNTSTRMODEL_NH1F = 4 ! Non-harmonic stretching with fracture, potential of type 1
integer(c_int), parameter :: CNTSTRMODEL_H1B = 5 ! Harmonic stetching model + axial buckling
integer(c_int), parameter :: CNTSTRMODEL_H1BH = 6 ! Harmonic stetching model + axial buckling + hysteresis
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 ! Maximal number of points in interpolation tables
!---------------------------------------------------------------------------------------------------
! Parameters of potentials
!---------------------------------------------------------------------------------------------------
! Stretching potential
integer(c_int) :: CNTSTRModel = CNTSTRMODEL_H1! Type of the bending model
integer(c_int) :: CNTSTRParams = 0 ! Type of parameterization
integer(c_int) :: CNTSTRYMT = 0 ! Type of dependence of the Young's modulus on tube radius
! 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
! paramerization, but only for calcuation of the
! force constant in eV/A)
real(c_double) :: CNTSTRD0 = 3.4d+00 ! CNT wall thickness (diameter of carbon atom), A
real(c_double) :: CNTSTREmin = -0.4d+00 ! Minimal strain in tabulated potential
real(c_double) :: CNTSTREmax = 0.13d+00 ! Maximal strain in tabulated potential. Simultaneously, U=0 if E> CNTSTREmax
real(c_double) :: CNTSTREl = 5.0d-02 ! Maximal linear strain
real(c_double) :: CNTSTREf = 12.0d-02 ! Strain at failure
real(c_double) :: CNTSTRS0 = 0.850e+12 ! Young's modulus, Pa
real(c_double) :: CNTSTRSl ! Maximal linear strees, Pa
real(c_double) :: CNTSTRSf = 75.0d+09 ! Tensile strength, Pa
real(c_double) :: CNTSTRF0 ! Elastic force constant, eV/A**2
real(c_double) :: CNTSTRFl ! Maximal linear force, eV/A**2
real(c_double) :: CNTSTRFf ! Tensile force at failure, eV/A**2
real(c_double) :: CNTSTRSi ! Maximal available stress (reference parameter, not used in the model), Pa
real(c_double) :: CNTSTRDf ! dF/dE at failure
real(c_double) :: CNTSTRAA, CNTSTRBB !
real(c_double) :: CNTSTRAAA, CNTSTRBBB ! | Auxilary constants
real(c_double) :: CNTSTRUl, CNTSTRUf ! /
! Axial buckling - hysteresis approch
real(c_double) :: CNTSTREc = -0.0142d+00 ! The minimal buckling strain
real(c_double) :: CNTSTREc1 = -0.04d+00 ! Critical axial buckling strain
real(c_double) :: CNTSTREc2 = -0.45d+00 ! Maximal buckling strain (the pot is harmonic for larger strains(in abs val))
!real(c_double) :: CNTSTRAmin
!real(c_double) :: CNTSTRAmax
!real(c_double) :: CNTSTRDA
! Bending potential
integer(c_int) :: CNTBNDModel = CNTBNDMODEL_H ! Type of the bending model
!real(c_double) :: CNTBNDAmin
!real(c_double) :: CNTBNDAmax
!real(c_double) :: CNTBNDDA
! Buckling model parameters
real(c_double) :: CNTBNDN = 1.0d+00 ! Buckling exponent
real(c_double) :: CNTBNDB = 0.68d+00 ! Buckling number
real(c_double) :: CNTBNDR = 275.0d+00 ! Critical radius of curvarure, A
! This is mean value for (10,10) SWCNT
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
!---------------------------------------------------------------------------------------------------
subroutine CNTSTRSetParameterization ( PType ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup parameters for further parameterization of streatching models
! References:
! [1] Yu M.-F. et al., Phys. Rev. Lett. 84(24), 5552 (2000)
! [2] Liew K.M. et al., Acta Materialia 52, 2521 (2004)
! [3] Mielke S.L. et al., Chem. Phys. Lett. 390, 413 (2004)
! [4] Zhigilei L.V. et al., Phys. Rev. B 71, 165417 (2005)
! [5] Kelly B.T., Physics of graphite, 1981
!-------------------------------------------------------------------------------------------
integer(c_int), intent(in) :: 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]
CNTSTREmin = -0.4d+00 ! Chosen arbitrary
CNTSTREmax = 3.64d-02 ! = CNTSTREf + 0.005
CNTSTREl = 2.0d-02 ! Chosen arbitrary
CNTSTREf = 3.14d-02 ! Ref. [1]
CNTSTRS0 = 1.002e+12 ! Ref. [1]
CNTSTRSf = 30.0d+09 ! Ref. [1]
case ( 1 ) ! This parameterization is taken from Ref. [2] for (10,10) SWCNT
! These values are obtained in MD simulatuions with REBO potential
! Values of Young's modulus, Tensile strenght and stress here
! are close to those obtained in Ref. [3] for pristine (defectless)
! (5,5) SWCNT in semiempirical QM calcuilations based on PM3 model
CNTSTRR0 = 6.785d+00 ! Calculated with usual formula for (10,10) CNT
CNTSTRD0 = 3.35d+00 ! Ref. [2]
CNTSTREmin = -0.4d+00 ! Chosen arbitrary
CNTSTREmax = 28.4d-02 ! = CNTSTREf + 0.005
CNTSTREl = 5.94d-02 ! Ref. [2]
CNTSTREf = 27.9d-02 ! Corresponds to Maximal 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) SWCNT
! with one atom vacancy defect obtained by semiempirical QM PM3 model
CNTSTRR0 = 3.43d+00 ! Ref. [3]
CNTSTRD0 = 3.4d+00 ! Ref. [3]
CNTSTREmin = -0.4d+00 ! Chosen arbitrary
CNTSTREmax = 15.8d-02 ! = CNTSTREf + 0.005
CNTSTREl = 6.00d-02 ! Chosed similar to Ref. [2]
CNTSTREf = 15.3d-02 ! Ref. [3]
CNTSTRS0 = 1.100e+12 ! Ref. [3]
CNTSTRSf = 100.0d+09 ! Ref. [3]
case ( 3 ) ! This special parameterization changes the only value of Young's modulus
! with accordance with the stretching constant in Ref. [4]
CNTSTRS0 = ( 86.64d+00 + 100.56d+00 * CNTSTRR0 ) * K_MDFU / ( M_2PI * CNTSTRR0 * CNTSTRD0 * 1.0e-20 ) ! Ref. [4]
case ( 4 ) ! This special parameterization changes the only value of Young's modulus
! making it equal to the in-plane Young's modulus of graphite
CNTSTRR0 = 6.785d+00 ! Calculated with usual formula for (10,10) CNT
CNTSTRD0 = 3.4d+00 ! Ref. [1]
CNTSTRS0 = 1.06e+12 ! Ref. [5]
end select
end subroutine CNTSTRSetParameterization !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Stretching without fracture, harmonic potential
!
integer(c_int) function CNTSTRH0Calc ( U, dUdL, L, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Young's modulus is independent of R
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: U, dUdL
real(c_double), intent(in) :: L, R0, L0
real(c_double) :: E
!-------------------------------------------------------------------------------------------
E = ( L - L0 ) / L0
dUdL = R0 * CNTSTRF0 * E
U = 0.5d+00 * L0 * E * dUdL
CNTSTRH0Calc = CNTPOT_STRETCHING
end function CNTSTRH0Calc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function CNTSTRH1Calc ( U, dUdL, L, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Young's modulus depends on R, see [4]
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: U, dUdL
real(c_double), intent(in) :: L, R0, L0
real(c_double) :: E, K
!-------------------------------------------------------------------------------------------
E = ( L - L0 ) / L0
K = 86.64d+00 + 100.56d+00 * R0
dUdL = K * E
U = 0.5d+00 * L0 * E * dUdL
CNTSTRH1Calc = CNTPOT_STRETCHING
end function CNTSTRH1Calc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Stretching without fracture, harmonic potential, with axial buckling without hysteresis
!
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
!-------------------------------------------------------------------------------------------
E = ( L - L0 ) / L0
K = 86.64d+00 + 100.56d+00 * R0
Kbcl = -10.98d+00 * L0
if ( E .gt. CNTSTREc ) then !Harmonic stretching
dUdL = K * E
U = 0.5d+00 * L0 * E * dUdL
CNTSTRH1BCalc = CNTPOT_STRETCHING
else if ( E .gt. CNTSTREc2 ) then !Axial buckling
dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc
U = Kbcl * E + dUbcl
dUdL = Kbcl / L0
CNTSTRH1BCalc = CNTPOT_STRETCHING !should be buckling, but doesn't work for some reason...
else !Return to harmonic potential
d = -0.0142794
dUdL = K * ( d + E - CNTSTREc2 )
dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc + Kbcl * CNTSTREc2
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 function CNTSTRH1BCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Stretching without fracture, harmonic potential, with axial buckling with hysteresis
!
integer(c_int) function CNTSTRH1BHCalc ( U, dUdL, L, R0, L0, ABF, Ebuc ) !!!!!!!!!!!!!!!!!!!!!!!!
! Young's modulus depends on R, see [4]
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: U, dUdL, Ebuc
real(c_double), intent(in) :: L, R0, L0
integer(c_int), intent(in) :: ABF
!-------------------------------------------------------------------------------------------
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
if ( E .gt. CNTSTREc ) then ! harmonic potential - no buckling
dUdL = K * E
U = 0.5d+00 * L0 * E * dUdL
CNTSTRH1BHCalc = CNTPOT_STRETCHING
Ebuc = 0.0d+00
else if ( E .gt. CNTSTREc1 ) then !above minimal buckling strain, but not at critical strain
if ( ABF .eq. 0 ) then ! not buckled. Continue harmonic potential
dUdL = K * E
U = 0.5d+00 * L0 * E * dUdL
CNTSTRH1BHCalc = CNTPOT_STRETCHING
Ebuc = 0.0d+00
else ! relaxing from buckled state. Use buckling potential
dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc
U = Kbcl * E + dUbcl
dUdL = Kbcl / L0
CNTSTRH1BHCalc = CNTPOT_SBUCKLING
Ebuc = 0.0d+00
end if
else if( E .gt. CNTSTREc2 ) then ! Axial buckling strain region
if ( ABF .eq. 0 ) then !newly buckled
dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc
U = Kbcl * E + dUbcl
dUdL = Kbcl / L0
CNTSTRH1BHCalc = CNTPOT_SBUCKLING
Ebuc = 0.5d+00 * L0 * K * CNTSTREc1 * CNTSTREc1 - Kbcl * CNTSTREc1 - dUbcl
else ! already buckled
dUbcl = 0.5d+00 * L0 * K * CNTSTREc * CNTSTREc - Kbcl * CNTSTREc
U = Kbcl * E + dUbcl
dUdL = Kbcl / L0
CNTSTRH1BHCalc = CNTPOT_SBUCKLING
Ebuc = 0.0d+00
end if
else ! Maximum strain and return to harmonic potential
dUdL = K * E
U = 0.5d+00 * L0 * E * dUdL
CNTSTRH1BHCalc = CNTPOT_STRETCHING
Ebuc = 0.0d+00
end if
end function CNTSTRH1BHCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Stretching with fracture, non-harmonic potential of type 0
!
integer(c_int) function CNTSTRNH0FCalc ( U, dUdL, L, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: U, dUdL
real(c_double), intent(in) :: L, R0, L0
real(c_double) :: E, DE, t
!-------------------------------------------------------------------------------------------
E = ( L - L0 ) / L0
if ( E < CNTSTREf ) then
dUdL = ( CNTSTRAA - CNTSTRBB * E ) * E
U = ( CNTSTRAAA - CNTSTRBBB * E ) * E * E
CNTSTRNH0FCalc = CNTPOT_STRETCHING
else
dUdL = 0.0d+00
U = 0.0d+00
CNTSTRNH0FCalc = CNTPOT_SFRACTURE
end if
U = L0 * R0 * U
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
CNTSTRAA = CNTSTRF0
CNTSTRBB = ( CNTSTRF0 * CNTSTREf - CNTSTRFf ) / ( CNTSTREf * CNTSTREf )
CNTSTRAAA= CNTSTRAA / 2.0d+00
CNTSTRBBB= CNTSTRAA / 3.0d+00
CNTSTRUl = 0.0d+00
CNTSTRUf = ( CNTSTRAAA - CNTSTRBBB * CNTSTREf ) * CNTSTREf * CNTSTREf
! These two values are not defined yet
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
real(c_double) :: E, C, DE, t
!-------------------------------------------------------------------------------------------
E = ( L - L0 ) / L0
if ( E < CNTSTREl ) then
dUdL = CNTSTRF0 * E
U = 0.5d+00 * E * dUdL
CNTSTRNH1Calc = CNTPOT_STRETCHING
else
DE = E - CNTSTREl
C = 1.0 + CNTSTRBB * DE
dUdL = CNTSTRFl + CNTSTRAA * ( 1.0d+00 - 1.0d+00 / C )
U = CNTSTRUl + CNTSTRAAA * DE - CNTSTRBBB * dlog ( C )
end if
CNTSTRNH1Calc = CNTPOT_STRETCHING
U = L0 * R0 * U
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
real(c_double) :: E, C, DE, t
!character(c_char)*512 :: Msg
!-------------------------------------------------------------------------------------------
E = ( L - L0 ) / L0
if ( E < CNTSTREl ) then
dUdL = CNTSTRF0 * E
U = 0.5d+00 * E * dUdL
CNTSTRNH1FCalc = CNTPOT_STRETCHING
else if ( E < CNTSTREf ) then
DE = E - CNTSTREl
C = 1.0 + CNTSTRBB * DE
dUdL = CNTSTRFl + CNTSTRAA * ( 1.0d+00 - 1.0d+00 / C )
U = CNTSTRUl + CNTSTRAAA * DE - CNTSTRBBB * dlog ( C )
CNTSTRNH1FCalc = CNTPOT_STRETCHING
else
!write ( Msg, * ) 'F Strains', E, CNTSTREf
!call PrintStdLogMsg ( Msg )
dUdL = 0.0d+00
U = 0.0d+00
CNTSTRNH1FCalc = CNTPOT_SFRACTURE
end if
U = L0 * R0 * U
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
CNTSTRAA = ( CNTSTRFf - CNTSTRFl ) * ( CNTSTREf * CNTSTRF0 - CNTSTRFl ) / ( CNTSTREf * CNTSTRF0 - CNTSTRFf )
CNTSTRBB = CNTSTRF0 / CNTSTRAA
CNTSTRAAA= CNTSTRFl + CNTSTRAA
CNTSTRBBB= CNTSTRAA / CNTSTRBB
CNTSTRSi = CNTSTRSl + CNTSTRAA / S
C = 1.0 + CNTSTRBB * ( CNTSTREf - CNTSTREl )
CNTSTRDf = CNTSTRF0 / C / C
CNTSTRUl = 0.5d+00 * CNTSTRFl * CNTSTREl
CNTSTRUf = CNTSTRUl + ( CNTSTRFl + CNTSTRAA ) * ( CNTSTREf - CNTSTREl ) - CNTSTRAA * dlog ( C ) / CNTSTRBB
end subroutine CNTSTRNH1Init !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! General
!
!integer(c_int) function CNTSTRCalc ( U, dUdL, L, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function CNTSTRCalc ( U, dUdL, L, R0, L0 , ABF, Ebuc ) !!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: U, dUdL, Ebuc
real(c_double), intent(in) :: L, R0, L0
integer(c_int), intent(in) :: ABF
!-------------------------------------------------------------------------------------------
Ebuc = 0.0d+00
select case ( CNTSTRModel )
case ( CNTSTRMODEL_H0 )
CNTSTRCalc = CNTSTRH0Calc ( U, dUdL, L, R0, L0 )
case ( CNTSTRMODEL_H1 )
CNTSTRCalc = CNTSTRH1Calc ( U, dUdL, L, R0, L0 )
case ( CNTSTRMODEL_NH0F )
CNTSTRCalc = CNTSTRNH0FCalc ( U, dUdL, L, R0, L0 )
case ( CNTSTRMODEL_NH1 )
CNTSTRCalc = CNTSTRNH1Calc ( U, dUdL, L, R0, L0 )
case ( CNTSTRMODEL_NH1F )
CNTSTRCalc = CNTSTRNH1FCalc ( U, dUdL, L, R0, L0 )
case ( CNTSTRMODEL_H1B )
CNTSTRCalc = CNTSTRH1BCalc ( U, dUdL, L, R0, L0 )
case ( CNTSTRMODEL_H1BH )
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
!real(c_double) :: A
!integer(c_int) :: i
!-------------------------------------------------------------------------------------------
CNTSTRModel = STRModel
CNTSTRParams = STRParams
CNTSTRYMT = YMType
if ( STRModel .ne. CNTSTRMODEL_H1 ) then
call CNTSTRSetParameterization ( STRParams )
if ( YMType == 2 ) then
call CNTSTRSetParameterization ( 4 )
else if ( YMType == 1 ) then
CNTSTRR0 = Rref
call CNTSTRSetParameterization ( 3 )
end if
if ( STRModel == CNTSTRMODEL_NH0F ) then
call CNTSTRNH0Init ()
else
call CNTSTRNH1Init ()
end if
end if
!CNTSTRAmin = -0.4d+00
!CNTSTRAmax = 0.4d+00
!CNTSTRDA = ( CNTSTRAmax - CNTSTRAmin ) / ( CNTPOTN - 1 )
!A = CNTSTRAmin
!do i = 0, CNTPOTN - 1
! CNTSTRU(i) = 0.5d+00 * A * A
! CNTSTRdUdA(i) = A
! A = A + CNTSTRDA
!end do
end subroutine CNTSTRInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Bending potentials
!---------------------------------------------------------------------------------------------------
subroutine BendingGradients ( K, G0, G1, G2, R0, R1, R2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This functions calculates degreeiest for bending forces
!-------------------------------------------------------------------------------------------
real(c_double), intent(inout) :: K
real(c_double), dimension(0:2), intent(inout) :: G0, G1, G2
real(c_double), dimension(0:2), intent(in) :: R0, R1, R2
real(c_double), dimension(0:2) :: DR0, DR2
real(c_double) :: L0, L2
!-------------------------------------------------------------------------------------------
DR0 = R0 - R1
DR2 = R2 - R1
L0 = S_V3norm3 ( DR0 )
L2 = S_V3norm3 ( DR2 )
DR0 = DR0 / L0
DR2 = DR2 / L2
K = S_V3xV3 ( DR0, DR2 )
G0 = DR2 - K * DR0
G2 = DR0 - K * DR2
G0 = G0 / L0
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
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: U, dUdC
real(c_double), intent(in) :: C, R0, L0
real(c_double) :: E, K
!-------------------------------------------------------------------------------------------
E = 1.0d+00 - C
K = 2.0d+00 * ( 63.8d+00 * R0**2.93d+00 ) / L0
U = K * ( 1.0d+00 + C ) / E
dUdC = 2.0d+00 * K / ( E * E )
CNTBNDHCalc = CNTPOT_BENDING
end function CNTBNDHCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function CNTBNDHBCalc ( U, dUdC, C, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Bending model of type 1:
! Harmonic bending potential with buckling
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: U, dUdC
real(c_double), intent(in) :: C, R0, L0
real(c_double) :: E1, E2, C2, Kbnd, Kbcl, Theta, DUbcl
!-------------------------------------------------------------------------------------------
E1 = 1.0d+00 - C
E2 = 1.0d+00 + C
! Calculate the square of curvature
C2 = 4.0d+00 * E2 / ( L0 * L0 * E1 )
! 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
Kbcl = CNTBNDB * Kbnd / 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
else ! 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 )
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
real(c_double) :: E1, E2, C2, Kbnd, Kbcl, Theta, DUbcl
!-------------------------------------------------------------------------------------------
E1 = 1.0d+00 - C
E2 = 1.0d+00 + C
! Calculate the square of curvature
C2 = 4.0d+00 * E2 / ( L0 * L0 * E1 )
! Check the condition for buckling
if ( C2 .ge. CNTBNDC2 ) then ! Buckling takes place
Theta= M_PI - acos ( C )
if ( Theta > CNTBNDTF ) then ! Fracture takes place
U = 0.0d+00
dUdC = 0.0d+00
CNTBNDHBFCalc = CNTPOT_BFRACTURE
else
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
U = Kbcl * abs ( Theta )**CNTBNDN - DUbcl
dUdC = Kbcl * CNTBNDN * abs ( Theta )**CNTBNDN1 / sqrt ( 1.0d+00 - C * C )
CNTBNDHBFCalc = CNTPOT_BBUCKLING
end if
else ! 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 )
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 approch.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: U, dUdC, Ebuc
real(c_double), intent(in) :: C , R0, L0
integer(c_int), intent(in) :: BBF
real(c_double) :: E1, E2, C2, Kbnd, Kbcl,Theta,DUbcl, Ubcl, Cmin,Rmax
!-------------------------------------------------------------------------------------------
Rmax = 340.0d+00
Cmin = 1.0/(Rmax*Rmax)
E1 = 1.0d+00 - C
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
Kbnd = 2.0d+00 * ( 63.8d+00 * R0**2.93d+00 ) / L0
U = Kbnd * E2 / E1
dUdC = 2.0d+00 * Kbnd / ( E1 * E1 )
CNTBNDHBHCalc = CNTPOT_BENDING
Ebuc = 0.0
else if ( C2 .ge. Cmin .and. C2 .lt. CNTBNDC2 ) then !Potential here depends on buckling flag of node
if ( BBF .eq. 0 ) then ! Not buckled yet. Continue harmonic bending
Kbnd = 2.0d+00 * ( 63.8d+00 * R0**2.93d+00 ) / L0
U = Kbnd * E2 / E1
dUdC = 2.0d+00 * Kbnd / ( E1 * E1 )
CNTBNDHBHCalc = CNTPOT_BENDING
Ebuc = 0.0d+00
else ! Already has been buckled or is buckled. Use buckling potential until Cmin.
Theta= M_PI - acos ( C )
Kbnd = 63.8d+00 * R0**2.93d+00
Kbcl = CNTBNDB * Kbnd / CNTBNDR
DUbcl= 2.0d+00*Kbnd * (1.0d+00+cos(l0/Rmax+M_PI))/(1.0d+00-cos(l0/Rmax+M_PI))/L0-Kbcl*abs(l0/Rmax)**CNTBNDN
U = Kbcl * abs( Theta )**CNTBNDN + DUbcl
dUdC = Kbcl * CNTBNDN * abs( Theta )**CNTBNDN1 / sqrt ( 1.0d+00 - C * C )
Ebuc = 0.0d+00
CNTBNDHBHCalc = CNTPOT_BBUCKLING
end if
else ! Greater than buckling critical point
if ( BBF .eq. 1 ) then ! Already buckled
Theta= M_PI - acos ( C )
Kbnd = 63.8d+00 * R0**2.93d+00
Kbcl = CNTBNDB * Kbnd / CNTBNDR
DUbcl= 2.0d+00*Kbnd * (1.0d+00+cos(l0/Rmax+M_PI))/(1.0d+00-cos(l0/Rmax+M_PI))/L0-Kbcl*abs(l0/Rmax)**CNTBNDN
U = Kbcl * abs( Theta )**CNTBNDN + DUbcl
dUdC = Kbcl * CNTBNDN * abs( Theta )**CNTBNDN1 / sqrt ( 1.0d+00 - C * C )
Ebuc = 0.0d00
CNTBNDHBHCalc = CNTPOT_BBUCKLING
else ! Newly buckled
Theta= M_PI - acos ( C )
Kbnd = 63.8d+00 * R0**2.93d+00
Kbcl = CNTBNDB * Kbnd / CNTBNDR
DUbcl= 2.0d+00*Kbnd * (1.0d+00+cos(l0/Rmax+M_PI))/(1.0d+00-cos(l0/Rmax+M_PI))/L0-Kbcl*abs(l0/Rmax)**CNTBNDN
U = Kbcl * abs( Theta )**CNTBNDN + DUbcl
dUdC = Kbcl * CNTBNDN * abs( Theta )**CNTBNDN1 / sqrt ( 1.0d+00 - C * C )
Ebuc = 2.0d+00*Kbnd * (1.0d+00+cos(l0/CNTBNDR+M_PI)) / (1.0d+00-cos(l0/CNTBNDR+M_PI))/L0- Kbcl*abs(l0/CNTBNDR)**CNTBNDN-dUbcl
CNTBNDHBHCalc = CNTPOT_BBUCKLING
end if
end if
end function CNTBNDHBHCalc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! General
!
! integer(c_int) function CNTBNDCalc ( U, dUdC, C, R0, L0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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
integer(c_int), intent(in) :: BBF
!-------------------------------------------------------------------------------------------
Ebuc = 0.0d+00
select case ( CNTBNDModel )
case ( CNTBNDMODEL_H )
CNTBNDCalc = CNTBNDHCalc ( U, dUdC, C, R0, L0 )
case ( CNTBNDMODEL_HB )
CNTBNDCalc = CNTBNDHBCalc ( U, dUdC, C, R0, L0 )
case ( CNTBNDMODEL_HBF )
CNTBNDCalc = CNTBNDHBFCalc ( U, dUdC, C, R0, L0 )
case ( CNTBNDMODEL_HBH )
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
integer(c_int) :: i
!-------------------------------------------------------------------------------------------
CNTBNDModel= BNDModel
CNTBNDN1 = CNTBNDN - 1.0d+00
CNTBNDC2 = 1.0d+00 / ( CNTBNDR * CNTBNDR )
!CNTBNDAmin = -1.0d+00
!CNTBNDAmax = 0.99d+00
!CNTBNDDA = ( CNTBNDAmax - CNTBNDAmin ) / ( CNTPOTN - 1 )
!A = CNTBNDAmin
!do i = 0, CNTPOTN - 1
! E = 1.0d+00 - A
! CNTBNDU(i) = 2.0d+00 * ( 1.0d+00 + A ) / E
! CNTBNDdUdA(i) = 4.0d+00 / E / E
! A = A + CNTBNDDA
!end do
end subroutine CNTBNDInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Module initialization
!---------------------------------------------------------------------------------------------------
subroutine InitCNTPotModule ( STRModel, STRParams, YMType, BNDModel, Rref ) !!!!!!!!!!!!!!!!
integer(c_int), intent(in) :: STRModel, STRParams, YMType, BNDModel
real(c_double), intent(in) :: Rref
!-------------------------------------------------------------------------------------------
call CNTSTRInit ( STRModel, STRParams, YMType, Rref )
call CNTBNDInit ( BNDModel )
end subroutine InitCNTPotModule !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module CNTPot !*********************************************************************************

125
lib/mesont/ExportCNT.f90 Normal file
View File

@ -0,0 +1,125 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !*******************************************************************************
use iso_c_binding
use CNTPot
use TPMLib
use TubePotMono
use TPMForceField
use iso_c_binding, only : c_int, c_double, c_char
implicit none
contains
subroutine InitCNTPotModule_(STRModel, STRParams, YMType, BNDModel, Rref) &
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_(TPMSSTPFile_, N1, TPMAFile_, N2) &
bind(c, name = "mesont_lib_SetTablePath")
integer(c_int), intent(in) :: N1, N2
character(c_char), intent(in), dimension(N1) :: TPMSSTPFile_
character(c_char), intent(in), dimension(N2) :: TPMAFile_
integer :: i
do i = 1, len(TPMSSTPFile)
if (i <= N1) then
TPMSSTPFile(i:i) = TPMSSTPFile_(i)
else
TPMSSTPFile(i:i) = ' '
endif
enddo
do i = 1, len(TPMAFile)
if (i <= N2) then
TPMAFile(i:i) = TPMAFile_(i)
else
TPMAFile(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")
real(c_double), intent(inout) :: U1, U2 ! Interaction energies associated with nodes X1 and X2
real(c_double), intent(inout), dimension(0:2) :: F1, F2 ! Forces exerted on nodes X1 and X2
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2 ! Contributions of nodes X1 and X2 to the virial stress tensor
real(c_double), intent(in), dimension(0:2) :: X1, X2 ! Coordinates of the segmnet nodes
real(c_double), intent(in) :: R12 ! Radius of nanotube the segment (X1,X2) belongs to
real(c_double), intent(in) :: L12 ! Equilubrium length of segment (X1,X2)
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")
real(c_double), intent(inout) :: U1, U2, U3 ! Interaction energies associated with nodes X1, X2, and X3
real(c_double), intent(inout), dimension(0:2) :: F1, F2, F3 ! Forces exerted on nodes X1, X2, and X3
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2, S3 ! Contributions of nodes X1, X2, and X3 to the virial stress tensor
real(c_double), intent(in), dimension(0:2) :: X1, X2, X3 ! Coordinates of nodes
real(c_double), intent(in) :: R123 ! Radius of nanotube the segment (X1,X2) belongs to
real(c_double), intent(in) :: L123 ! Equilubrium length of segment (X1,X2) and (X2,X3) (It is assumed to be the same for both segments)
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")
integer(c_int), intent(in) :: N ! Number of nodes in array X
real(c_double), intent(inout) :: U1, U2 ! Interaction energies associated with nodes X1 and X2
real(c_double), intent(inout), dimension(0:N-1) :: U ! Interaction energies associated with nodes X
real(c_double), intent(inout), dimension(0:2) :: F1, F2 ! Forces exerted on nodes X1 and X2
real(c_double), intent(inout), dimension(0:2,0:N-1) :: F ! Forces exerted on nodes X
real(c_double), intent(inout), dimension(0:2) :: Fe ! Force exerted on node Xe (can be updated only if Ee > 0)
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2 ! Contributions of nodes X1 and X2 to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2,0:N-1) :: S ! Contributions of nodes X to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2) :: Se ! Contributions of node Xe to the virial stress tensor (can be updated only if Ee > 0)
real(c_double), intent(in), dimension(0:2) :: X1, X2 ! Coordinates of the segmnet nodes
real(c_double), intent(in) :: R12 ! Radius of nanotube the segment (X1,X2) belongs to
real(c_double), intent(in), dimension(0:2,0:N-1) :: X ! Coordinates of the nanotube nodes
real(c_double), intent(in), dimension(0:2) :: Xe ! Additional node of the extended chain if Ee > 0
integer(c_int), intent(in), dimension(0:N-1) :: BBF ! Bending buckling flags (BBF(i) = 1 in a case of buckling in node i)
real(c_double), intent(in) :: R ! Radius of nanotube X
integer(c_int), intent(in) :: E1, E2 ! E1 = 1 if the chnane node 0 is a CNT end; E2 = 1 if the chnane node N-1 is a CNT end;
integer(c_int), intent(in) :: Ee ! Parameter defining the type of the extended chain (0,1,2)
integer(c_int), intent(in) :: TPMType ! Type of the tubular potential (0 or 1)
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
endmodule ExportCNT !**************************************************************************

113
lib/mesont/LinFun2.f90 Normal file
View File

@ -0,0 +1,113 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !************************************************************************************
!
! TMD Library: Bi-linear functions and their derivatives
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
!
!***************************************************************************************************
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 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int), intent(in) :: i, N
real(c_double), intent(in) :: X
real(c_double), dimension(0:N-1), intent(in) :: P
real(c_double), dimension(0:N-1), intent(inout) :: F
integer(c_int) :: i1
real(c_double) :: A, A0
!-------------------------------------------------------------------------------------------
i1 = i - 1
A0 = ( P(i) - X ) / ( P(i) - P(i1) )
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
real(c_double), intent(in) :: X
real(c_double), dimension(0:N-1), intent(in) :: P
real(c_double), dimension(0:N-1), intent(inout) :: F, Fx
integer(c_int) :: i1
real(c_double) :: A, A0
!-------------------------------------------------------------------------------------------
i1 = i - 1
A0 = ( P(i) - X ) / ( P(i) - P(i1) )
A = 1.0d+00 - A0
S = A0 * F(i1) + A * F(i)
Sx1 = A0 * Fx(i1) + A * Fx(i)
end subroutine CalcLinFun1_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function CalcLinFun2_0 ( i, j, X, Y, N1, N2, P1, P2, F ) !!
integer(c_int), intent(in) :: i, j, N1, N2
real(c_double), intent(in) :: X, Y
real(c_double), dimension(0:N1-1), intent(in) :: P1
real(c_double), dimension(0:N2-1), intent(in) :: P2
real(c_double), dimension(0:N1-1,0:N2-1), intent(inout) :: F
integer(c_int) :: i1, j1
real(c_double) :: A, A0, B, B0, G, G0
!-------------------------------------------------------------------------------------------
i1 = i - 1
j1 = j - 1
A0 = ( P1(i) - X ) / ( P1(i) - P1(i1) )
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)
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
real(c_double), intent(in) :: X, Y
real(c_double), dimension(0:N1-1), intent(in) :: P1
real(c_double), dimension(0:N2-1), intent(in) :: P2
real(c_double), dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fx, Fy
integer(c_int) :: i1, j1
real(c_double) :: A, A0, B, B0, G, G0
!-------------------------------------------------------------------------------------------
i1 = i - 1
j1 = j - 1
A0 = ( P1(i) - X ) / ( P1(i) - P1(i1) )
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 !********************************************************************************

56
lib/mesont/Makefile.gfortran Normal file
View File

@ -0,0 +1,56 @@
SHELL = /bin/sh
# which file will be copied to Makefile.lammps
EXTRAMAKE = Makefile.lammps.gfortran
# ------ FILES ------
SRC = LinFun2.f90 Spline1.f90 Spline2.f90 TPMLib.f90 TPMGeom.f90 TubePotBase.f90 TubePotTrue.f90 \
TubePotMono.f90 TPMM0.f90 TPMM1.f90 CNTPot.f90 TPMForceField.f90 ExportCNT.f90
FILES = $(SRC) Makefile
# ------ DEFINITIONS ------
LIB = libmesont.a
OBJ = $(SRC:.f90=.o)
# ------ SETTINGS ------
F90 = gfortran
CC = gcc
F90FLAGS = -O3 -fPIC -ffast-math -ftree-vectorize -fexpensive-optimizations -fno-second-underscore -g -ffree-line-length-none -cpp
#F90FLAGS = -O
ARCHIVE = ar
ARCHFLAG = -rc
LINK = g++
LINKFLAGS = -O
USRLIB =
SYSLIB =
# ------ MAKE PROCEDURE ------
lib: $(OBJ)
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
@cp $(EXTRAMAKE) Makefile.lammps
# ------ COMPILE RULES ------
%.o:%.F
$(F90) $(F90FLAGS) -c $<
%.o:%.f90
$(F90) $(F90FLAGS) -c $<
%.o:%.c
$(CC) $(F90FLAGS) -c $<
#include .depend
# ------ CLEAN ------
clean:
-rm *.o *.mod $(LIB)
tar:
-tar -cvf ../MESONT.tar $(FILES)

52
lib/mesont/Makefile.ifort Normal file
View File

@ -0,0 +1,52 @@
SHELL = /bin/sh
# which file will be copied to Makefile.lammps
EXTRAMAKE = Makefile.lammps.ifort
# ------ FILES ------
SRC = LinFun2.f90 Spline1.f90 Spline2.f90 TPMLib.f90 TPMGeom.f90 TubePotBase.f90 TubePotTrue.f90 \
TubePotMono.f90 TPMM0.f90 TPMM1.f90 CNTPot.f90 TPMForceField.f90 ExportCNT.f90
FILES = $(SRC) Makefile
# ------ DEFINITIONS ------
LIB = libmesont.a
OBJ = $(SRC:.f90=.o)
# ------ SETTINGS ------
F90 = ifort
F90FLAGS = -Ofast -fPIC -ipo -fpp
ARCHIVE = ar
ARCHFLAG = -rc
USRLIB =
SYSLIB =
# ------ MAKE PROCEDURE ------
lib: $(OBJ)
$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
@cp $(EXTRAMAKE) Makefile.lammps
# ------ COMPILE RULES ------
%.o:%.F
$(F90) $(F90FLAGS) -c $<
%.o:%.f90
$(F90) $(F90FLAGS) -c $<
%.o:%.c
$(CC) $(F90FLAGS) -c $<
#include .depend
# ------ CLEAN ------
clean:
-rm *.o *.mod $(LIB)
tar:
-tar -cvf ../MESONT.tar $(FILES)

5
lib/mesont/Makefile.lammps.gfortran Normal file
View File

@ -0,0 +1,5 @@
# Settings that the LAMMPS build will import when this package library is used
mesont_SYSINC =
mesont_SYSLIB = -lgfortran
mesont_SYSPATH =

5
lib/mesont/Makefile.lammps.ifort Normal file
View File

@ -0,0 +1,5 @@
# Settings that the LAMMPS build will import when this package library is used
mesont_SYSINC =
mesont_SYSLIB = -lifcore -lsvml -limf -ldl -lstdc++ -lgfortran
mesont_SYSPATH = -L/opt/intel/fce/10.0.023/lib

1
lib/mesont/Makefile.serial Normal file
View File

@ -0,0 +1 @@
Makefile.gfortran

67
lib/mesont/README Normal file
View File

@ -0,0 +1,67 @@
USER-MESONT is a LAMMPS package for simulation of nanomechanics of carbon
nanotubes (CNTs). The model is based on a coarse-grained representation
of CNTs as "flexible cylinders" consisting of a variable number of
segments. Internal interactions within a CNT and the van der Waals
interaction between the tubes are described by a mesoscopic force
field designed and parameterized based on the results of atomic-level
molecular dynamics simulations. The description of the force field
is provided in the papers listed below.
This folder contains a Fortran library implementing basic level functions
describing stretching, bending, and intertube components of the CNT tubular
potential model (TPM) mesoscopic force field.
This library was created by Alexey N. Volkov, University of Alabama,
avolkov1@ua.edu.
--
References:
L. V. Zhigilei, C. Wei, and D. Srivastava, Mesoscopic model for dynamic
simulations of carbon nanotubes, Phys. Rev. B 71, 165417, 2005.
A. N. Volkov and L. V. Zhigilei, Structural stability of carbon nanotube
films: The role of bending buckling, ACS Nano 4, 6187-6195, 2010.
A. N. Volkov, K. R. Simov, and L. V. Zhigilei, Mesoscopic model for simulation
of CNT-based materials, Proceedings of the ASME International Mechanical
Engineering Congress and Exposition (IMECE2008), ASME paper IMECE2008-68021,
2008.
A. N. Volkov and L. V. Zhigilei, Mesoscopic interaction potential for carbon
nanotubes of arbitrary length and orientation, J. Phys. Chem. C 114, 5513-5531,
2010.
B. K. Wittmaack, A. H. Banna, A. N. Volkov, L. V. Zhigilei, Mesoscopic
modeling of structural self-organization of carbon nanotubes into vertically
aligned networks of nanotube bundles, Carbon 130, 69-86, 2018.
B. K. Wittmaack, A. N. Volkov, L. V. Zhigilei, Mesoscopic modeling of the
uniaxial compression and recovery of vertically aligned carbon nanotube
forests, Compos. Sci. Technol. 166, 66-85, 2018.
B. K. Wittmaack, A. N. Volkov, L. V. Zhigilei, Phase transformation as the
mechanism of mechanical deformation of vertically aligned carbon nanotube
arrays: Insights from mesoscopic modeling, Carbon 143, 587-597, 2019.
A. N. Volkov and L. V. Zhigilei, Scaling laws and mesoscopic modeling of
thermal conductivity in carbon nanotube materials, Phys. Rev. Lett. 104,
215902, 2010.
A. N. Volkov, T. Shiga, D. Nicholson, J. Shiomi, and L. V. Zhigilei, Effect
of bending buckling of carbon nanotubes on thermal conductivity of carbon
nanotube materials, J. Appl. Phys. 111, 053501, 2012.
A. N. Volkov and L. V. Zhigilei, Heat conduction in carbon nanotube materials:
Strong effect of intrinsic thermal conductivity of carbon nanotubes, Appl.
Phys. Lett. 101, 043113, 2012.
W. M. Jacobs, D. A. Nicholson, H. Zemer, A. N. Volkov, and L. V. Zhigilei,
Acoustic energy dissipation and thermalization in carbon nanotubes: Atomistic
modeling and mesoscopic description, Phys. Rev. B 86, 165414, 2012.
A. N. Volkov and A. H. Banna, Mesoscopic computational model of covalent
cross-links and mechanisms of load transfer in cross-linked carbon nanotube
films with continuous networks of bundles, Comp. Mater. Sci. 176, 109410, 2020.

192
lib/mesont/Spline1.f90 Normal file
View File

@ -0,0 +1,192 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !************************************************************************************
!
! TMD Library: One-dimensional cubic spline function
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
!
!***************************************************************************************************
use iso_c_binding, only : c_int, c_double, c_char
implicit none
contains !******************************************************************************************
real(c_double) function ValueSpline1_0 ( X, Xi, Xi_1, Yi, Yi_1, Mi, Mi_1, Hi_1 ) !!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: X, Xi, Xi_1, Yi, Yi_1, Mi, Mi_1, Hi_1
real(c_double) :: H26, HL, HR
!-------------------------------------------------------------------------------------------
H26 = Hi_1 * Hi_1 / 6.0
Hl = X - Xi_1
Hr = Xi - X
ValueSpline1_0 = ( ( Mi_1 * Hr * Hr * Hr + Mi * Hl * Hl * Hl ) / 6.0 + ( Yi_1 - Mi_1 * H26 ) * Hr + ( Yi - Mi * H26 ) * Hl ) / Hi_1
end function ValueSpline1_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine ValueSpline1_1 ( S, S1, X, Xi, Xi_1, Yi, Yi_1, Mi, Mi_1, Hi_1 ) !!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: S, S1
real(c_double), intent(in) :: X, Xi, Xi_1, Yi, Yi_1, Mi, Mi_1, Hi_1
real(c_double) :: H6, H26, HL, HR, HL2, HR2
!-------------------------------------------------------------------------------------------
H6 = Hi_1 / 6.0d+00
H26 = Hi_1 * H6
HL = X - Xi_1
HR = Xi - X
HL2 = HL * HL
HR2 = HR * HR
S = ( ( Mi_1 * HR2 * Hr + Mi * HL2 * Hl ) / 6.0 + ( Yi_1 - Mi_1 * H26 ) * HR + ( Yi - Mi * H26 ) * HL ) / Hi_1
S1 = ( ( Mi * HL2 - Mi_1 * HR2 ) / 2.0d+00 + Yi - Yi_1 ) / Hi_1 - H6 * ( Mi - Mi_1 )
end subroutine ValueSpline1_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine sprogonka3 ( N, K0, K1, K2, F, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! K0[i] * X[i-1] + K1[i] * X[I] + K2[i] * X[i+1] = F[i]
! i = 0..(N-1)
!-------------------------------------------------------------------------------------------
integer(c_int), intent(in) :: N
real(c_double), dimension(0:N-1), intent(in) :: K0, K1, K2
real(c_double), dimension(0:N-1), intent(inout) :: F, X
real(c_double) :: D
integer(c_int) :: i
!-------------------------------------------------------------------------------------------
X(0) = F(0) / K1(0)
F(0) = - K2(0) / K1(0)
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
end do
do i = N - 2, 0, -1
X(i) = X(i) + F(i) * X(i+1)
end do
end subroutine sprogonka3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CreateSpline1 ( CL, CR, N, P, F, M, D, K0, K1, K2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int), intent(in) :: CL, CR, N
real(c_double), dimension (0:N-1), intent(in) :: P, F
real(c_double), dimension (0:N-1), intent(inout):: M, D, K0, K1, K2
integer(c_int) :: i
real(c_double) :: Z
!-------------------------------------------------------------------------------------------
do i = 1, N - 1
K0(i) = P(i) - P(i-1)
K1(i) = ( F(i) - F(i-1) ) / K0(i)
end do
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)
case (2)
K1(0) = 1.0d+00
K2(0) = 0.0d+00
D(0) = M(0)
case (3)
K1(0) = 1.0d+00
K2(0) = 0.0d+00
D(0) = 0.0d+00
end select
Z = K1(N-1)
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 )
case (1)
D(N-1) = 2.0d+00 * ( M(N-1) - Z ) / K0(N-1)
K1(N-1) = 2.0d+00 / 3.0d+00
K0(N-1) = 1.0d+00 / 3.0d+00
case (2)
K1(N-1) = 1.0d+00
K0(N-1) = 0.0d+00
D(N-1) = M(N-1)
case (3)
K1(N-1) = 1.0d+00
K0(N-1) = 0.0d+00
D(N-1) = 0.0d+00
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
integer(c_int) :: j
real(c_double) :: HL, HR, H, H6, H26, HR2, HL2, HRH, HLH
!-------------------------------------------------------------------------------------------
j = i - 1
HL = X - P(j)
HR = P(i) - X
H = P(i) - P(j)
H6 = H / 6.0d+00
H26 = H * H6
HL2 = HL * HL
HR2 = HR * HR
HLH = HL / H
HRH = HR / H
CalcSpline1_0 = ( M(j) * HR2 * HRH + M(i) * HL2 * HLH ) / 6.0d+00 + ( F(j) - M(j) * H26 ) * HRH + ( F(i) - M(i) * H26 ) * HLH
end function CalcSpline1_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CalcSpline1_1 ( S, S1, i, X, N, P, F, M ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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
integer(c_int) :: j
real(c_double) :: HL, HR, H, H6, H26, HR2, HL2, HRH, HLH
!-------------------------------------------------------------------------------------------
j = i - 1
HL = X - P(j)
HR = P(i) - X
H = P(i) - P(j)
H6 = H / 6.0d+00
H26 = H * H6
HL2 = HL * HL
HR2 = HR * HR
HLH = HL / H
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) )
end subroutine CalcSpline1_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CalcSpline1_2 ( S, S1, S2, i, X, N, P, F, M ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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
integer(c_int) :: j
real(c_double) :: HL, HR, H, H6, H26, HR2, HL2, HRH, HLH
!-------------------------------------------------------------------------------------------
j = i - 1
HL = X - P(j)
HR = P(i) - X
H = P(i) - P(j)
H6 = H / 6.0d+00
H26 = H * H6
HL2 = HL * HL
HR2 = HR * HR
HLH = HL / H
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
end subroutine CalcSpline1_2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module Spline1 !********************************************************************************

186
lib/mesont/Spline2.f90 Normal file
View File

@ -0,0 +1,186 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !************************************************************************************
!
! TMD Library: Two-dimensional cubic spline function
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
!
!***************************************************************************************************
use Spline1
use iso_c_binding, only : c_int, c_double, c_char
implicit none
contains !******************************************************************************************
subroutine CreateSpline2 ( CL, CD, CR, CU, N1, N2, N, P1, P2, F, Fxx, Fyy, Fxxyy, FF, MM, DD, K0, K1, K2 )
integer(c_int), intent(in) :: CL, CD, CR, CU, N1, N2, N
real(c_double), dimension(0:N1-1), intent(in) :: P1
real(c_double), dimension(0:N2-1), intent(in) :: P2
real(c_double), dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy
real(c_double), dimension(0:N-1), intent(inout) :: FF, MM, DD, K0, K1, K2
integer(c_int) :: II
!-------------------------------------------------------------------------------------------
do II = 0, N2 - 1
FF(0:N1-1) = F(0:N1-1,II)
MM(0) = Fxx(0,II)
MM(N1-1) = Fxx(N1-1,II)
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 = 0, N1 - 1
MM(0) = Fyy(II,0)
MM(N-1) = Fyy(II,N2-1)
FF(0:N2-1) = F(II,0:N2-1)
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:N1-1) = Fyy(0:N1-1,N2-1 )
call CreateSpline1 ( 3, 3, N1, P1, FF, MM, DD, K0, K1, k2 )
Fxxyy(0:N1-1,N2-1) = MM(0:N1-1)
do II = 1, N1 - 2
MM(0) = Fxxyy(II,0)
MM(N-1) = Fxxyy(II,N2-1)
FF(0:N2-1) = Fxx(II,0:N2-1)
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 CreateSpline2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CreateSpline2Ext ( CL, CD, CR, CU, N1, N1A, N2, N2A, N, P1, P2, F, Fxx, Fyy, Fxxyy, FF, MM, DD, K0, K1, K2 )
integer(c_int), intent(in) :: CL, CD, CR, CU, N1, N1A, N2, N2A, N
real(c_double), dimension(0:N1-1), intent(in) :: P1
real(c_double), dimension(0:N2-1), intent(in) :: P2
real(c_double), dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy
real(c_double), dimension(0:N-1), intent(inout) :: FF, MM, DD, K0, K1, K2
integer(c_int) :: II
!-------------------------------------------------------------------------------------------
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)
MM(N1-1) = Fxx(N1-1,II)
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)
MM(N1-N1A-1) = Fxx(N1-1,II)
call CreateSpline1 ( CL, CR, N1 - N1A, P1, FF, MM, DD, K0, K1, K2 )
Fxx(N1A:N1-1,II) = MM(0:N1-N1A-1)
end do
do II = 0, N1A - 1
MM(0) = Fyy(II,0)
MM(N2A) = Fyy(II,N2A)
FF(0:N2A) = F(II,0:N2A)
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)
FF(0:N2-1) = F(II,0:N2-1)
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)
FF(0:N2A) = Fxx(II,0:N2A)
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)
FF(0:N2-1) = Fxx(II,0:N2-1)
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
real(c_double), dimension(0:N1-1), intent(in) :: P1
real(c_double), dimension(0:N2-1), intent(in) :: P2
real(c_double), dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy
integer(c_int) :: i1, j1
real(c_double) :: T, Gy_0, Gy_1, Gxxy_0, Gxxy_1
!-------------------------------------------------------------------------------------------
i1 = i - 1
j1 = j - 1
T = P2(j) - P2(j1)
Gy_0 = ValueSpline1_0 ( Y, P2(j), P2(j1), F(i,j), F(i,j1), Fyy(i,j), Fyy(i,j1), T )
Gy_1 = ValueSpline1_0 ( Y, P2(j), P2(j1), F(i1,j), F(i1,j1), Fyy(i1,j), Fyy(i1,j1), T )
Gxxy_0 = ValueSpline1_0 ( Y, P2(j), P2(j1), Fxx(i,j), Fxx(i,j1), Fxxyy(i,j), Fxxyy(i,j1), T )
Gxxy_1 = ValueSpline1_0 ( Y, P2(j), P2(j1), Fxx(i1,j), Fxx(i1,j1), Fxxyy(i1,j), Fxxyy(i1,j1), T )
CalcSpline2_0 = ValueSpline1_0 ( X, P1(i), P1(i1), Gy_0, Gy_1,Gxxy_0, Gxxy_1, P1(i) - P1(i1) )
end function CalcSpline2_0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CalcSpline2_1 ( S, Sx1, Sy1, i, j, X, Y, N1, N2, P1, P2, F, Fxx, Fyy, Fxxyy ) !!!
real(c_double), intent(out) :: S, Sx1, Sy1
integer(c_int), intent(in) :: i, j, N1, N2
real(c_double), intent(in) :: X, Y
real(c_double), dimension(0:N1-1), intent(in) :: P1
real(c_double), dimension(0:N2-1), intent(in) :: P2
real(c_double), dimension(0:N1-1,0:N2-1), intent(inout) :: F, Fxx, Fyy, Fxxyy
integer(c_int) :: i1, j1
real(c_double) :: T, Gy_0, Gy_1, Gxxy_0, Gxxy_1
real(c_double) :: Gyy_0, Gyy_1, Gxxyy_0, Gxxyy_1
!-------------------------------------------------------------------------------------------
i1 = i - 1
j1 = j - 1
T = P2(j) - P2(j1)
call ValueSpline1_1 ( Gy_0, Gyy_0, Y, P2(j), P2(j1), F(i,j), F(i,j1), Fyy(i,j), Fyy(i,j1), T )
call ValueSpline1_1 ( Gy_1, Gyy_1, Y, P2(j), P2(j1), F(i1,j), F(i1,j1), Fyy(i1,j), Fyy(i1,j1), T )
call ValueSpline1_1 ( Gxxy_0, Gxxyy_0, Y, P2(j), P2(j1), Fxx(i,j), Fxx(i,j1), Fxxyy(i,j), Fxxyy(i,j1), T )
call ValueSpline1_1 ( Gxxy_1, Gxxyy_1, Y, P2(j), P2(j1), Fxx(i1,j), Fxx(i1,j1), Fxxyy(i1,j), Fxxyy(i1,j1), T )
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 !********************************************************************************

290
lib/mesont/TPMForceField.f90 Normal file
View File

@ -0,0 +1,290 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !************************************************************************************
!
! TMD Library: Calculation of the TMD force field
!
!---------------------------------------------------------------------------------------------------
!
! PGI Fortran, Intel Fortran
!
! Alexey N. Volkov, University of Alabama (avolkov1@ua.edu), Version 09.01.33, 2018
!
!***************************************************************************************************
use CNTPot
use TPMM0
use TPMM1
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 ) !!!!!!!!!!!
real(c_double), intent(inout) :: U1, U2 ! Interaction energies associated with nodes X1 and X2
real(c_double), intent(inout), dimension(0:2) :: F1, F2 ! Forces exerted on nodes X1 and X2
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2 ! Contributions of nodes X1 and X2 to the virial stress tensor
real(c_double), intent(in), dimension(0:2) :: X1, X2 ! Coordinates of the segmnet nodes
real(c_double), intent(in) :: R12 ! Radius of nanotube the segment (X1,X2) belongs to
real(c_double), intent(in) :: L12 ! Equilubrium length of segment (X1,X2)
!-------------------------------------------------------------------------------------------
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
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
F2 = F2 - FF
U2 = U2 + U
! Stress
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
end do
end do
end subroutine TubeStretchingForceField !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TubeBendingForceField ( U1, U2, U3, F1, F2, F3, S1, S2, S3, X1, X2, X3, R123, L123, BBF2 )
real(c_double), intent(inout) :: U1, U2, U3 ! Interaction energies associated with nodes X1, X2, and X3
real(c_double), intent(inout), dimension(0:2) :: F1, F2, F3 ! Forces exerted on nodes X1, X2, and X3
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2, S3 ! Contributions of nodes X1, X2, and X3 to the virial stress tensor
real(c_double), intent(in), dimension(0:2) :: X1, X2, X3 ! Coordinates of nodes
real(c_double), intent(in) :: R123 ! Radius of nanotube the segment (X1,X2) belongs to
real(c_double), intent(in) :: L123 ! Equilubrium length of segment (X1,X2) and (X2,X3) (It is assumed to be the same for both segments)
integer(c_int), intent(inout) :: BBF2
!-------------------------------------------------------------------------------------------
integer(c_int) :: ii, jj, Event
real(c_double) :: U, F, K, S, Ubcl
real(c_double), dimension(0:2) :: G0, G1, G2
!-------------------------------------------------------------------------------------------
call BendingGradients ( K, G0, G1, G2, X1, X2, X3 )
Event = CNTBNDCalc ( U, F, K, R123, L123, BBF2, Ubcl )
if ( Event == CNTPOT_BBUCKLING ) then
BBF2 = 1
else
BBF2 = 0
end if
U = U / 4.0d+00
F = - F
F1 = F1 + G0 * F
F2 = F2 + G1 * F
F3 = F3 + G2 * F
U1 = U1 + U
U2 = U2 + 2.0d+00 * U
U3 = U3 + U
! Stress
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
S = 0.5d+00 * ( X3(ii) - X2(ii) ) * G2(jj)
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
! (as well potential nergies and componets of the virial stress tensor) between a segment
! (X1,X2) and a sequence of segments with node coordinates that belongs to a single CNT
! It is assumed that X contains ALL nodes of a single CNT that are included into the
! neighbor list of segment (X1,X2)
! The nodes in X are assumed to be ordered according to their physical appearence in the nanotube
! It means that (X(i),X(i+1)) are either correspond to a real segment or divided by a segments
! that do not belong to a nanotube.
! Concept of the extendend chain:
! Let's consider a sequant of nodes (X1,X2,...,XN) forming continuous part of a nanotube.
! If node Xe preceeds X1 and Xe is the nanotube end, then the extended chain is (Xe,X1,...,XN) and Ee = 1.
! 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, extended chain coincides with (X1,...,XN) and Ee = 0
! If the extended chain contains additional node, then non-zero force is exterted 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 )
integer(c_int), intent(in) :: N ! Number of nodes in array X
real(c_double), intent(inout) :: U1, U2 ! Interaction energies associated with nodes X1 and X2
real(c_double), intent(inout), dimension(0:N-1) :: U ! Interaction energies associated with nodes X
real(c_double), intent(inout), dimension(0:2) :: F1, F2 ! Forces exerted on nodes X1 and X2
real(c_double), intent(inout), dimension(0:2,0:N-1) :: F ! Forces exerted on nodes X
real(c_double), intent(inout), dimension(0:2) :: Fe ! Force exerted on node Xe (can be updated only if Ee > 0)
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2 ! Contributions of nodes X1 and X2 to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2,0:N-1) :: S ! Contributions of nodes X to the virial stress tensor
real(c_double), intent(inout), dimension(0:2,0:2) :: Se ! Contributions of node Xe to the virial stress tensor (can be updated only if Ee > 0)
real(c_double), intent(in), dimension(0:2) :: X1, X2 ! Coordinates of the segmnet nodes
real(c_double), intent(in) :: R12 ! Radius of nanotube the segment (X1,X2) belongs to
real(c_double), intent(in), dimension(0:2,0:N-1) :: X ! Coordinates of the nanotube nodes
real(c_double), intent(in), dimension(0:2) :: Xe ! Additiona node of the extended chain if Ee > 0
integer(c_int), intent(in), dimension(0:N-1) :: BBF ! Bending buckling flags (BBF(i) = 1 in a case of buckling in node i)
real(c_double), intent(in) :: R ! Radius of nanotube X
integer(c_int), intent(in) :: E1, E2 ! E1 = 1 if the chnane node 0 is a CNT end; E1 = 2 if the chnane node N-1 is a CNT end;
integer(c_int), intent(in) :: Ee ! Parameter defining the type of the extended chain (0,1,2)
integer(c_int), intent(in) :: TPMType ! Type of the tubular potential (0 or 1)
!-------------------------------------------------------------------------------------------
integer(c_int) :: k, ii, jj, IntSign
integer(c_int) :: BType, EType, LocalTPMType
real(c_double), dimension(0:2,0:N-1) :: G1, G2
real(c_double), dimension(0:N-1) :: QQ
logical :: EType1, EType2
real(c_double), dimension(0:2) :: G, DG, DQ, XX
real(c_double) :: UT, DR, DS, DS1
real(c_double) :: xU1, xU2 ! Interaction energies associated with nodes X1 and X2
real(c_double), dimension(0:N-1) :: xU ! Interaction energies associated with nodes X
real(c_double), dimension(0:2) :: xF1, xF2 ! Forces exerted on nodes X1 and X2
real(c_double), dimension(0:2,0:N-1) :: xF ! Forces exerted on nodes X
real(c_double), dimension(0:2) :: xFe ! Force exerted on node Xe (can be updated only if Ee > 0)
!-------------------------------------------------------------------------------------------
!U1 = 0.0d+00
!U2 = 0.0d+00
!U = 0.0d+00
!F1 = 0.0d+00
!F2 = 0.0d+00
!F = 0.0d+00
!S1 = 0.0d+00
!S2 = 0.0d+00
!S = 0.0d+00
! Looking for a buckling point
BType = 0
do k = 0, N - 1
if ( BBF(k) == 1 ) then
BType = 1
exit
end if
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.
! Overwise, LocalTPMType = TPMType.
if ( BType == 1 ) then
LocalTPMType = 0
EType = 0
else
if ( E1 == 1 ) then ! First node in the chain is the tube end
EType1 = .true.
else
EType1 = .false.
end if
if ( E2 == 1 ) then ! Last node in the chain is the tube end
EType2 = .true.
else
EType2 = .false.
end if
if ( EType1 .and. EType2 ) then
LocalTPMType = 0
else
LocalTPMType = TPMType
if ( EType1 ) then
EType = 1
else if ( EType2 ) then
EType = 2
else ! No tube ends in the chain
EType = 0
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
if ( EType == 0 ) then
if ( Ee == 1 ) then ! First node in the extended chain is the tube end
EType = 3
else if ( Ee == 2 ) then ! Last node in the extended chain is the tube end
EType = 4
end if
end if
IntSign = TPMInteractionFW1 ( QQ, UT, xU1, xU2, xU, xF1, xF2, xF, xFe, G1, G2, X1, X2, N, N, X, Xe, EType )
end if
if ( IntSign == 0 ) return ! No interaction
! Final potential energies
U1 = U1 + 0.5d+00 * xU1
U2 = U2 + 0.5d+00 * xU2
U(0:N-1) = U(0:N-1) + 0.5d+00 * xU(0:N-1)
! Contributions to the virial stresses tensor
do ii = 0, 2
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
S2(ii,jj) = S2(ii,jj) + DS
end do
end do
XX = 0.5d+00 * ( X2 + X1 )
if ( EType > 2 ) then
DQ = Xe - XX
call ApplyPeriodicBC ( DQ )
DQ = DQ / 6.0d+00
do ii = 0, 2
do jj = 0, 2
DS = DQ(ii) * xFe(jj)
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
call ApplyPeriodicBC ( DQ )
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) )
do jj = 0, 2
DS = DQ(ii) * G(jj)
DS1 = DS + DR * DG(jj)
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
end do
end do
end do
! Final forces
F1 = F1 + 0.5d+00 * xF1
F2 = F2 + 0.5d+00 * xF2
F(0:2,0:N-1) = F(0:2,0:N-1) + 0.5d+00 * xF(0:2,0:N-1)
if ( EType > 2 ) then
Fe = Fe + 0.5d+00 * xFe
end if
end subroutine SegmentTubeForceField !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TPMForceField !**************************************************************************

159
lib/mesont/TPMGeom.f90 Normal file
View File

@ -0,0 +1,159 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !************************************************************************************
!
! TMD Library: Geometry functions
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
!
!***************************************************************************************************
use TPMLib
use iso_c_binding, only : c_int, c_double, c_char
implicit none
!---------------------------------------------------------------------------------------------------
! Constants
!---------------------------------------------------------------------------------------------------
integer(c_int), parameter :: MD_LINES_NONPAR = 0
integer(c_int), parameter :: MD_LINES_PAR = 1
!---------------------------------------------------------------------------------------------------
! 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
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 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This subroutine changes coortinates of the point accorning to periodic boundary conditions
! it order to makesure that the point is inside the computational cell
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2), intent(inout) :: R
!-------------------------------------------------------------------------------------------
! These commented lines implemment the more general, but less efficient algorithm
!if ( BC_X == 1 ) R(0) = R(0) - DomLX * roundint ( R(0) / DomLX )
!if ( BC_Y == 1 ) R(1) = R(1) - DomLY * roundint ( R(1) / DomLY )
!if ( BC_Z == 1 ) R(2) = R(2) - DomLZ * roundint ( R(2) / DomLZ )
if ( BC_X == 1 ) then
if ( R(0) .GT. DomLXHalf ) then
R(0) = R(0) - DomLX
else if ( R(0) .LT. - DomLXHalf ) then
R(0) = R(0) + DomLX
end if
end if
if ( BC_Y == 1 ) then
if ( R(1) .GT. DomLYHalf ) then
R(1) = R(1) - DomLY
else if ( R(1) .LT. - DomLYHalf ) then
R(1) = R(1) + DomLY
end if
end if
if ( BC_Z == 1 ) then
if ( R(2) .GT. DomLZHalf ) then
R(2) = R(2) - DomLZ
else if ( R(2) .LT. - DomLZHalf ) then
R(2) = R(2) + DomLZ
end if
end if
end subroutine ApplyPeriodicBC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine LinePoint ( Displacement, Q, R1, L1, R0 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This function calculates the point Q of projection of point R0 on line (R1,L1)
! Q = R1 + Disaplacement * L1
!-------------------------------------------------------------------------------------------
real(c_double), intent(inout) :: Displacement
real(c_double), dimension(0:2), intent(inout) :: Q
real(c_double), dimension(0:2), intent(in) :: R1, L1, R0
!--------------------------------------------------------------------------------------------
Q = R0 - R1
! Here we take into account periodic boundaries
call ApplyPeriodicBC ( Q )
Displacement = S_V3xV3 ( Q, L1 )
Q = R1 + Displacement * L1
end subroutine LinePoint !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function LineLine ( H, cosA, D1, D2, L12, R1, L1, R2, L2, Prec ) !!!!!!!!!!!!!!!!!
! This function determines the neares distance H between two lines (R1,L1) and (R2,L2)
!-------------------------------------------------------------------------------------------
! Input values:
! R1, L1, point and direction of line 1
! R2, L2, point and direction of line 2
! Prec, precision for the case L1 * L2 = 0 (parallel lines)
! Return values:
! H, minimal distance between lines
! cosA, cosine of angle between lines
! D1, D2, displacemets
! L12, unit vector directed along the closes 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), dimension(0:2) :: Q1, Q2, R
real(c_double) :: C, DD1, DD2, C1, C2
!-------------------------------------------------------------------------------------------
cosA = S_V3xV3 ( L1, L2 )
C = 1.0 - sqr ( cosA )
if ( C < Prec ) then ! Lines are parallel to each other
LineLine = MD_LINES_PAR
return
end if
LineLine = MD_LINES_NONPAR
R = R2 - R1
! Here we take into account periodic boundaries
call ApplyPeriodicBC ( R )
DD1 = S_V3xV3 ( R, L1 )
DD2 = S_V3xV3 ( R, L2 )
D1 = ( cosA * DD2 - DD1 ) / C
D2 = ( DD2 - cosA * DD1 ) / C
Q1 = R1 - D1 * L1
Q2 = R2 - D2 * L2
L12 = Q2 - Q1
! Here we take into account periodic boundaries
call ApplyPeriodicBC ( L12 )
H = S_V3norm3 ( L12 )
if ( H < Prec ) then ! Lines intersect each other
C1 = signum ( D1 )
C2 = signum ( D1 ) * signum ( cosA )
Q1 = C1 * L1
Q2 = C2 * L2
call V3_V3xxV3 ( L12, Q1, Q2 )
call V3_ort ( L12 )
else ! No intersection
L12 = L12 / H
end if
end function LineLine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TPMGeom !********************************************************************************

219
lib/mesont/TPMLib.f90 Normal file
View File

@ -0,0 +1,219 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !*************************************************************************************
!
! TMD Library: Basic constants, types, and mathematical functions
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
!
!***************************************************************************************************
use iso_c_binding, only : c_int, c_double, c_char
implicit none
!---------------------------------------------------------------------------------------------------
! Mathematical constants
!---------------------------------------------------------------------------------------------------
real(c_double), parameter :: M_PI_2 = 1.57079632679489661923
real(c_double), parameter :: M_PI = 3.14159265358979323846
real(c_double), parameter :: M_3PI_2 = 4.71238898038468985769
real(c_double), parameter :: M_2PI = 6.28318530717958647692
real(c_double), parameter :: M_PI_180 = 0.017453292519943295769
!---------------------------------------------------------------------------------------------------
! Physical unit constants
!---------------------------------------------------------------------------------------------------
real(c_double), parameter :: K_AMU = 1.66056E-27 ! a.m.u. (atomic mass unit, Dalton)
real(c_double), parameter :: K_EV = 1.60217646e-19 ! eV (electron-volt)
real(c_double), parameter :: K_MDLU = 1.0E-10 ! MD length unit (m)
real(c_double), parameter :: K_MDEU = K_EV ! MD energy unit (J)
real(c_double), parameter :: K_MDMU = K_AMU ! MD mass unit (kg)
real(c_double), parameter :: K_MDFU = K_MDEU / K_MDLU ! MD force unit (N)
real(c_double), parameter :: K_MDCU = K_MDEU / K_MDMU ! MD specific heat unit (J/(kg*K))
!---------------------------------------------------------------------------------------------------
! Global variables
!---------------------------------------------------------------------------------------------------
integer(c_int) :: StdUID = 31
contains !******************************************************************************************
!---------------------------------------------------------------------------------------------------
! Simple mathematical functions
!---------------------------------------------------------------------------------------------------
real(c_double) function rad ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: X
!-------------------------------------------------------------------------------------------
rad = X * M_PI_180
end function rad !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function sqr ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: X
!-------------------------------------------------------------------------------------------
sqr = X * X
end function sqr !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function signum ( X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: X
!-------------------------------------------------------------------------------------------
if ( X > 0 ) then
signum = 1
else if ( X < 0 ) then
signum = -1
else
signum = 0
end if
end function signum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Vector & matrix functions
!---------------------------------------------------------------------------------------------------
real(c_double) function S_V3xx ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), dimension(0:2), intent(in) :: V
!-------------------------------------------------------------------------------------------
S_V3xx = V(0) * V(0) + V(1) * V(1) + V(2) * V(2)
end function S_V3xx !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function S_V3xV3 ( V1, V2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), dimension(0:2), intent(in) :: V1, V2
!-------------------------------------------------------------------------------------------
S_V3xV3 = V1(0) * V2(0) + V1(1) * V2(1) + V1(2) * V2(2)
end function S_V3xV3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function S_V3norm3 ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), dimension(0:2), intent(in) :: V
!-------------------------------------------------------------------------------------------
S_V3norm3 = dsqrt ( V(0) * V(0) + V(1) * V(1) + V(2) * V(2) )
end function S_V3norm3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine V3_ort ( V ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Vector production
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2), intent(inout) :: V
!-------------------------------------------------------------------------------------------
real(c_double) :: Vabs
!-------------------------------------------------------------------------------------------
Vabs = S_V3norm3 ( V )
V(0) = V(0) / Vabs
V(1) = V(1) / Vabs
V(2) = V(2) / Vabs
end subroutine V3_ort !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine V3_V3xxV3 ( V, V1, V2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Vector production
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2), intent(out) :: V
real(c_double), dimension(0:2), intent(in) :: V1, V2
!-------------------------------------------------------------------------------------------
V(0) = V1(1) * V2(2) - V1(2) * V2(1)
V(1) = V1(2) * V2(0) - V1(0) * V2(2)
V(2) = V1(0) * V2(1) - V1(1) * V2(0)
end subroutine V3_V3xxV3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Handling of spherical and Euler angles
!---------------------------------------------------------------------------------------------------
subroutine RotationMatrix3 ( M, Psi, Tet, Phi ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Ksi, Tet and Phi are Euler angles
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2,0:2), intent(out) :: M
real(c_double), intent(in) :: Psi, Tet, Phi
!-------------------------------------------------------------------------------------------
real(c_double) :: cK, sK, cT, sT, cP, sP
!-------------------------------------------------------------------------------------------
cK = dcos ( Psi )
sK = dsin ( Psi )
cT = dcos ( Tet )
sT = dsin ( Tet )
cP = dcos ( Phi )
sP = dsin ( Phi )
M(0,0) = cP * cK - sK * sP * cT
M(0,1) = cP * sK + sP * cT * cK
M(0,2) = sP * sT
M(1,0) = - sP * cK - cP * cT * sK
M(1,1) = - sP * sK + cP * cT * cK
M(1,2) = cP * sT
M(2,0) = sT * sK
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
!-------------------------------------------------------------------------------------------
Tet = acos ( L(2) )
Psi = atan2 ( L(1), L(0) )
if ( Psi > M_3PI_2 ) then
Psi = Psi - M_3PI_2
else
Psi = Psi + M_PI_2
end if
end subroutine EulerAngles !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! File inout and output
!---------------------------------------------------------------------------------------------------
integer(c_int) function OpenFile ( Name, Params, Path ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
character*(*), intent(in) :: Name, Params, Path
!-------------------------------------------------------------------------------------------
integer(c_int) :: Fuid
character*512 :: FullName, Msg, Name1, Action1, Status1, Form1, Position1
!-------------------------------------------------------------------------------------------
OpenFile = StdUID + 1
if ( Params(1:1) == 'r' ) then
Action1 = 'read'
Status1 = 'old'
Position1 = 'rewind'
else if ( Params(1:1) == 'w' ) then
Action1 = 'write'
Status1 = 'replace'
Position1 = 'rewind'
else if ( Params(1:1) == 'a' ) then
Action1 = 'write'
Status1 = 'old'
Position1 = 'append'
endif
if ( Params(2:2) == 'b' ) then
Form1 = 'binary'
else
Form1 = 'formatted'
endif
open ( unit = OpenFile, file = Name, form = Form1, action = Action1, status = Status1, position = Position1 )
StdUID = StdUID + 1
return
end function OpenFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine CloseFile ( Fuid ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int), intent(inout) :: Fuid
!-------------------------------------------------------------------------------------------
if ( Fuid < 0 ) return
close ( unit = Fuid )
Fuid = -1
end subroutine CloseFile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TPMLib !*********************************************************************************

195
lib/mesont/TPMM0.f90 Normal file
View File

@ -0,0 +1,195 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !**************************************************************************************
!
! TMD Library: Combined/Weighted potential of type 0
!
! Direct application of SST potential to calculation of segment-segment interaction
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
!
!***************************************************************************************************
!use TMDCounters
use TubePotMono
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
real(c_double), dimension(0:2), intent(in) :: R1_1, R1_2, R2_1, R2_2
integer(c_int), intent(in) :: EType
!-------------------------------------------------------------------------------------------
real(c_double) :: Qa, Ua, Fd, L2
real(c_double), dimension(0:2) :: F1_1a, F1_2a, F2_1a, F2_2a, R2_3, R2, Laxis2, F
integer(c_int) :: IntSign
!-------------------------------------------------------------------------------------------
! C_TPM_4 = C_TPM_4 + 1
R2 = 0.5d+00 * ( R2_1 + R2_2 )
Laxis2 = R2_2 - R2_1
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 )
R2_3 = R2_2 + R2_2 - R2_1
IntSign = TPMInteractionF ( Qa, Ua, F1_1a, F1_2a, F2_1a, F2_2a, Fd, R1_1, R1_2, R2_2, R2_3, 1 )
if ( IntSign > 0 ) then
TPMInteractionFSS = 1
call TPMSegmentForces ( F2_1a, F2_2a, F1_1a, F1_2a, R1_1, R1_2, R2, Laxis2, L2 )
F = ( Fd - S_V3xV3 ( F2_2a, Laxis2 ) ) * Laxis2
F2_2a = F2_2a + F
F2_1a = F2_1a - F
end if
else
TPMInteractionFSS = TPMInteractionF ( Q, U, F1_1, F1_2, F2_1, F2_2, Fd, R1_1, R1_2, R2_1, R2_2, 2 )
R2_3 = R2_1 + R2_1 - R2_2
IntSign = TPMInteractionF ( Qa, Ua, F1_1a, F1_2a, F2_1a, F2_2a, Fd, R1_1, R1_2, R2_1, R2_3, 1 )
if ( IntSign > 0 ) then
TPMInteractionFSS = 1
call TPMSegmentForces ( F2_1a, F2_2a, F1_1a, F1_2a, R1_1, R1_2, R2, Laxis2, L2 )
F = ( - Fd - S_V3xV3 ( F2_1a, Laxis2 ) ) * Laxis2
F2_1a = F2_1a + F
F2_2a = F2_2a - F
end if
end if
if ( IntSign > 0 ) then
Q = Q - Qa
if ( Q < 0.0d+00 ) Q = 0.0d+00
U = U - Ua
F2_1 = F2_1 - F2_1a
F2_2 = F2_2 - F2_2a
F1_1 = F1_1 - F1_1a
F1_2 = F1_2 - F1_2a
end if
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
real(c_double), dimension(0:NMAX-1), intent(out) :: QQ, UU
real(c_double), dimension(0:2), intent(out) :: F1, F2
real(c_double), dimension(0:2,0:NMAX-1), intent(out) :: F, G1, G2
real(c_double), dimension(0:2), intent(in) :: R1, R2
real(c_double), dimension(0:2,0:NMAX-1), intent(in) :: R
!-------------------------------------------------------------------------------------------
integer(c_int) :: i, SType2, GeomID, EType
real(c_double) :: Ua
real(c_double), dimension(0:2) :: F1_1a, F1_2a, F2_1a, F2_2a
real(c_double), dimension(0:2) :: R1a, R2a, Laxis1, Laxis2, L12, DR
real(c_double) :: L1, L2, D1, D2, H, cosA, D, Dmina, Dminb
!-------------------------------------------------------------------------------------------
QQ = 0.0d+00
U = 0.0d+00
U1 = 0.0d+00
U2 = 0.0d+00
UU = 0.0d+00
F1 = 0.0d+00
F2 = 0.0d+00
F = 0.0d+00
G1 = 0.0d+00
G2 = 0.0d+00
TPMInteractionFW0 = 0
do i = 0, N - 2
R1a = 0.5d+00 * ( R1 + R2 )
R2a = 0.5d+00 * ( R(0:2,i+1) + R(0:2,i) )
Laxis1 = R2 - R1
Laxis2 = R(0:2,i+1) - R(0:2,i)
L1 = S_V3norm3 ( Laxis1 )
L2 = S_V3norm3 ( Laxis2 )
Laxis1 = Laxis1 / L1
Laxis2 = Laxis2 / L2
L2 = 0.5d+00 * L2
L1 = 0.5d+00 * L1
GeomID = LineLine ( H, cosA, D1, D2, L12, R1a, Laxis1, R2a, Laxis2, TPGeomPrec )
DR = R1 - R(0:2,i)
call ApplyPeriodicBC ( DR )
Dmina = sqr ( DR(0) ) + sqr ( DR(1) ) + sqr ( DR(2) )
DR = R2 - R(0:2,i)
call ApplyPeriodicBC ( DR )
D = sqr ( DR(0) ) + sqr ( DR(1) ) + sqr ( DR(2) )
if ( D < Dmina ) Dmina = D
if ( GeomID == MD_LINES_NONPAR ) then
D = ( D2 - L2 ) * cosA
if ( D > D1 - L1 .and. D < D1 + L1 ) then
D = sqr ( D2 - L2 ) * ( 1.0d+00 - sqr ( cosA ) ) + sqr ( H )
if ( D < Dmina ) Dmina = D
end if
else
call LinePoint ( D, DR, R1, Laxis1, R(0:2,i) )
if ( D > 0.0d+00 .and. D < 2.0d+00 * L1 ) then
DR = DR - R(0:2,i)
call ApplyPeriodicBC ( DR )
D = sqr ( DR(0) ) + sqr ( DR(1) ) + sqr ( DR(2) )
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) )
DR = R2 - R(0:2,i+1)
call ApplyPeriodicBC ( DR )
D = sqr ( DR(0) ) + sqr ( DR(1) ) + sqr ( DR(2) )
if ( D < Dminb ) Dminb = D
if ( GeomID == MD_LINES_NONPAR ) then
D = ( D2 + L2 ) * cosA
if ( D > D1 - L1 .and. D < D1 + L1 ) then
D = sqr ( D2 + L2 ) * ( 1.0d+00 - sqr ( cosA ) ) + sqr ( H )
if ( D < Dminb ) Dminb = D
end if
else
call LinePoint ( D, DR, R1, Laxis1, R(0:2,i+1) )
if ( D > 0.0d+00 .and. D < 2.0d+00 * L1 ) then
DR = DR - R(0:2,i+1)
call ApplyPeriodicBC ( DR )
D = sqr ( DR(0) ) + sqr ( DR(1) ) + sqr ( DR(2) )
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), EType ) > 0 ) then
TPMInteractionFW0 = 1
U = U + Ua
Ua = 0.25d+00 * Ua
U1 = U1 + Ua
U2 = U2 + Ua
UU(i) = UU(i) + Ua
UU(i+1) = UU(i+1) + Ua
F1 = F1 + F1_1a
F2 = F2 + F1_2a
F(0:2,i) = F(0:2,i) + F2_1a
F(0:2,i+1) = F(0:2,i+1) + F2_2a
G2(0:2,i) = F2_1a
G1(0:2,i+1) = F2_2a
end if
end do
end function TPMInteractionFW0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TPMM0 !**********************************************************************************

379
lib/mesont/TPMM1.f90 Normal file
View File

@ -0,0 +1,379 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !**************************************************************************************
!
! TMD Library: Combined/Weighted potential of type 3
!
! Weighting functions are the same as in potential of type 2.
! Calculation of the combined potential is based on the 'extended' chain.
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran.
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
!
!***************************************************************************************************
!use TMDCounters
use TubePotMono
use iso_c_binding, only : c_int, c_double, c_char
implicit none
!---------------------------------------------------------------------------------------------------
! Constants
!---------------------------------------------------------------------------------------------------
! Maximal length of a segment chain
integer(c_int), parameter :: TPM_MAX_CHAIN = 100
!---------------------------------------------------------------------------------------------------
! Numerical parameters
!---------------------------------------------------------------------------------------------------
! Switching parameters
real(c_double) :: TPMC123 = 1.0d+00 ! Non-dimensional
real(c_double) :: TPMC3 = 10.0d+00 ! (A)
!---------------------------------------------------------------------------------------------------
! Global variables
!---------------------------------------------------------------------------------------------------
! These global variables are used to speedup calculations
real(c_double), dimension(0:2,0:TPM_MAX_CHAIN-1) :: E1, E2, EE1, EE2
real(c_double), dimension(0:2) :: Q1, Q2, Qe, Qe1, DR, Z1, Z2, S1, S2, Pe, Pe1
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
real(c_double), dimension(0:2), intent(in) :: R2_1, R2_2
!-------------------------------------------------------------------------------------------
real(c_double) :: D, L20, D20, t, dWdD
real(c_double), dimension(0:2) :: E, E20
!-------------------------------------------------------------------------------------------
E = 0.5d+00 * ( R2_1 + R2_2 ) - RR
call ApplyPeriodicBC ( E )
D = E(0) * E(0) + E(1) * E(1) + E(2) * E(2)
if ( D < D10 * D10 ) then
W = 1.0d+00
E1_1 = 0.0d+00
E1_2 = 0.0d+00
E2_1 = 0.0d+00
E2_2 = 0.0d+00
return
end if
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
W = 0.0d+00
E1_1 = 0.0d+00
E1_2 = 0.0d+00
E2_1 = 0.0d+00
E2_2 = 0.0d+00
return
end if
D = sqrt ( D )
E = E / D
E20 = E20 / L20
D20 = D20 - D10
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 )
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
real(c_double), dimension(0:2), intent(in) :: R1_1, R1_2, R2_1, R2_2
!-------------------------------------------------------------------------------------------
real(c_double) :: D, L20
real(c_double) :: D1, D2, t, dWdD
real(c_double), dimension(0:2) :: RR, E, E20
!-------------------------------------------------------------------------------------------
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 )
L20 = sqrt ( S_V3xx ( E20 ) + sqr ( TPMR2 ) )
D1 = L10 + L20 + TPBRcutoff + RSkin
if ( D < D1 ) then
EndWeight1 = 0
W = 1.0d+00
E1_1 = 0.0d+00
E1_2 = 0.0d+00
E2_1 = 0.0d+00
E2_2 = 0.0d+00
return
end if
D2 = D1 + TPMC3
if ( D > D2 ) then
EndWeight1 = 2
W = 0.0d+00
E1_1 = 0.0d+00
E1_2 = 0.0d+00
E2_1 = 0.0d+00
E2_2 = 0.0d+00
return
end if
EndWeight1 = 1
E = E / D
E20 = E20 / L20
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 )
end function EndWeight1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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
integer(c_int), intent(in) :: EType
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2) :: M, QX, Me, F1a, F2a, P1a, P2a, F1b, F2b, P1b, P2b, ER1, ER2, EQe, EQe1
real(c_double) :: W, W1, D, Qa, Qb, Ua, Ub, L, Pee, Peea, Peeb, DU
integer(c_int) :: IntSigna, IntSignb, CaseID
!-------------------------------------------------------------------------------------------
if ( EType == 0 ) then
! C_TPM_0 = C_TPM_0 + 1
TPMInteractionFC1 = TPMInteractionF ( Q, U, F1, F2, P1, P2, Pee, R1, R2, Q1, Q2, 0 )
Pe = 0.0d+00
Pe1 = 0.0d+00
else if ( EType < 3 ) then
! C_TPM_1 = C_TPM_1 + 1
QX = 0.5d+00 * ( Q1 + Q2 )
M = Q2 - Q1
L = S_V3norm3 ( M )
M = M / L
Me = Qe - QX
D = S_V3norm3 ( Me )
if ( EType == 1 ) then
TPMInteractionFC1 = TPMInteractionF ( Q, U, F1, F2, P1, P2, Pee, R1, R2, QX - D * M, QX, 1 )
else
TPMInteractionFC1 = TPMInteractionF ( Q, U, F1, F2, P1, P2, Pee, R1, R2, QX, QX + D * M, 2 )
end if
call TPMSegmentForces ( P1, P2, F1, F2, R1, R2, QX, M, L )
Pe = ( Pee / D ) * Me
Pe1 = 0.0d+00
QX = 0.5d+00 * Pe
P1 = P1 + QX
P2 = P2 + QX
else
CaseID = EndWeight1 ( W, ER1, ER2, EQe, Eqe1, R1, R2, Qe, Qe1 )
if ( CaseID < 2 ) then
QX = 0.5d+00 * ( Q1 + Q2 )
M = Q2 - Q1
L = S_V3norm3 ( M )
M = M / L
Me = Qe - QX
D = S_V3norm3 ( Me )
if ( EType == 3 ) then
IntSigna = TPMInteractionF ( Qa, Ua, F1a, F2a, P1a, P2a, Peea, R1, R2, QX - D * M, QX, 1 )
else
IntSigna = TPMInteractionF ( Qa, Ua, F1a, F2a, P1a, P2a, Peea, R1, R2, QX, QX + D * M, 2 )
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
! C_TPM_1 = C_TPM_1 + 1
TPMInteractionFC1 = IntSigna
Q = Qa
U = Ua
F1 = F1a
F2 = F2a
Pe = ( Peea / D ) * Me
Pe1 = 0.0d+00
QX = 0.5d+00 * Pe
P1 = P1a + QX
P2 = P2a + QX
else if ( CaseID == 2 ) then
! C_TPM_0 = C_TPM_0 + 1
TPMInteractionFC1 = IntSignb
Q = Qb
U = Ub
F1 = F1b
F2 = F2b
P1 = P1b
P2 = P2b
Pe = 0.0d+00
Pe1 = 0.0d+00
else
! C_TPM_2 = C_TPM_2 + 1
TPMInteractionFC1 = 0
if ( IntSigna > 0 .or. IntSignb > 0 ) TPMInteractionFC1 = 1
W1 = 1.0d+00 - W
DU = Ub - Ua
Q = W * Qa + W1 * Qb
U = W * Ua + W1 * Ub
Pe = ( W * Peea / D ) * Me
QX = 0.5d+00 * Pe
F1 = W * F1a + W1 * F1b + DU * ER1
F2 = W * F2a + W1 * F2b + DU * ER2
P1 = W * P1a + W1 * P1b + QX
P2 = W * P2a + W1 * P2b + QX
Pe = Pe - DU * EQe
Pe1 = - DU * EQe1
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
real(c_double), dimension(0:NMAX-1), intent(out) :: QQ, UU
real(c_double), dimension(0:2), intent(out) :: F1, F2, Fe
real(c_double), dimension(0:2,0:NMAX-1), intent(out) :: F, G1, G2
real(c_double), dimension(0:2), intent(in) :: R1, R2, Re
real(c_double), dimension(0:2,0:NMAX-1), intent(in) :: R
!-------------------------------------------------------------------------------------------
integer(c_int) :: i, j
real(c_double) :: Q, WW, DD
!-------------------------------------------------------------------------------------------
Q1 = 0.0d+00
Q2 = 0.0d+00
WW = 0.0d+00
Z1 = 0.0d+00
Z2 = 0.0d+00
TPMInteractionFW1 = 0
E10 = 0.5d+00 * ( R2 - R1 )
L10 = sqrt ( S_V3xx ( E10 ) + sqr ( TPMR1 ) )
D10 = TPMR1 + TPMR2 + TPMC123 * TPBRcutoff + RSkin
E10 = E10 / L10
RR = 0.5d+00 * ( R1 + R2 )
do i = 0, N - 2
call PairWeight1 ( W(i), E1(0:2,i), E2(0:2,i), EE1(0:2,i), EE2(0:2,i), R(0:2,i), R(0:2,i+1) )
Q1 = Q1 + W(i) * R(0:2,i)
Q2 = Q2 + W(i) * R(0:2,i+1)
WW = WW + W(i)
Z1 = Z1 + E1(0:2,i)
Z2 = Z2 + E2(0:2,i)
end do
if ( WW .le. TPGeomPrec ) return
Q1 = Q1 / WW
Q2 = Q2 / WW
Z1 = Z1 / WW
Z2 = Z2 / WW
if ( EType == 1 ) then
Qe = R(0:2,0)
Qe1 = R(0:2,1)
else if ( EType == 2 ) then
Qe = R(0:2,N-1)
Qe1 = R(0:2,N-2)
else if ( EType == 3 ) then
Qe = Re
Qe1 = R(0:2,0)
else if ( EType == 4 ) then
Qe = Re
Qe1 = R(0:2,N-1)
else
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 )
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
EE1(0:2,0:N-2) = EE1(0:2,0:N-2) / WW
EE2(0:2,0:N-2) = EE2(0:2,0:N-2) / WW
G1(0:2,0:N-1) = 0.0d+00
G2(0:2,0:N-1) = 0.0d+00
U1 = 0.25d+00 * U
U2 = U1
UU = 0.0d+00
do i = 0, N - 2
QQ(i) = W(i) * Q
DD = W(i) * U1
UU(i) = UU(i) + DD
UU(i+1) = UU(i+1) + DD
end do
do i = 0, N - 2
C(i) = S_V3xV3 ( S1, R(0:2,i) ) + S_V3xV3 ( S2, R(0:2,i+1) )
F1 = F1 + C(i) * ( E1(0:2,i) - W(i) * Z1 )
F2 = F2 + C(i) * ( E2(0:2,i) - W(i) * Z2 )
end do
F(0:2,0) = W(0) * S1
do j = 0, N - 2
if ( j == 0 ) then
DR = EE1(0:2,0) * ( 1.0d+00 - W(0) )
else
DR = - W(j) * EE1(0:2,0)
end if
F(0:2,0) = F(0:2,0) + C(j) * DR
end do
do i = 1, N - 2
G1(0:2,i) = W(i-1) * S2
G2(0:2,i) = W(i) * S1
do j = 0, N - 2
if ( j == i ) then
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) * ( EE1(0:2,j) - W(j) * EE1(0:2,i) )
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)
G2(0:2,i) = G2(0:2,i) - C(j) * W(j) * EE1(0:2,i)
end if
end do
F(0:2,i) = G1(0:2,i) + G2(0:2,i)
end do
F(0:2,N-1) = W(N-2) * S2
do j = 0, N - 2
if ( j == N - 2 ) then
DR = EE2(0:2,N-2) * ( 1.0d+00 - W(N-2) )
else
DR = - W(j) * EE2(0:2,N-2)
end if
F(0:2,N-1) = F(0:2,N-1) + C(j) * DR
end do
Fe = 0.0d+00
if ( EType == 1 ) then
F(0:2,0) = F(0:2,0) - Pe
else if ( EType == 2 ) then
F(0:2,N-1) = F(0:2,N-1) - Pe
else if ( EType == 3 ) then
F(0:2,0) = F(0:2,0) - Pe1
Fe = - Pe
else if ( EType == 4 ) then
F(0:2,N-1) = F(0:2,N-1) - Pe1
Fe = - Pe
end if
G1(0:2,N-1) = F(0:2,N-1)
G2(0:2,0) = F(0:2,0)
end function TPMInteractionFW1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TPMM1 !**********************************************************************************

332
lib/mesont/TubePotBase.f90 Normal file
View File

@ -0,0 +1,332 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !********************************************************************************
!
! TMD Library: Non-Bonded pair interaction potential and transfer functions for atoms composing
! nanotubes.
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
!
!---------------------------------------------------------------------------------------------------
!
! This module contains basic parameters for all modules involved into calculations of tubular
! potentials.
!
! It includes definitions of
! -- TPBU, Lennard-Jones (12-6) potential
! -- TPBQ, Transfer function
!
! All default values are adjusted for non-bonded carbob-carbon interaction in carbon nanotubes.
!
!***************************************************************************************************
use TPMLib
use iso_c_binding, only : c_int, c_double, c_char
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
! Maximal number of elements in corresponding tables
integer(c_int), parameter :: TPBNMAX = 2001
! Numerical constants
real(c_double), parameter :: TPbConstD = 5.196152422706632d+00 ! = 3.0**1.5
! Mass of C atom
real(c_double), parameter :: TPBMc = 12.0107d+00 ! (Da)
! Parameters of the Van der Waals inteaction 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)
real(c_double), parameter :: TPBScc = 3.4d+00 ! (A)
! Lattice parameter and numerical density of atoms for a graphene sheet, see
! 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 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 = TPBScc ! (A)
!real(c_double), parameter :: TPBQRcutoff1cc = 2.16d+00 * TPBScc ! (A)
real(c_double), parameter :: TPBQScc = 7.0d+00 ! (A)
real(c_double), parameter :: TPBQRcutoff1cc = 8.0d+00 ! (A)
!---------------------------------------------------------------------------------------------------
! Global variables
!---------------------------------------------------------------------------------------------------
logical :: TPErrCheck = .true. ! Set to .true. to generate diagnostic and warning messages
character*512 :: TPErrMsg = '' ! Typically, this variable is set up in F_tt ()
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
! Physical 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 LJ (12-6) interatomic potential (eV)
real(c_double) :: TPBS = TPBScc ! Sigma parameter of LJ (12-6) 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)
! Physical parameters of the transfer function
real(c_double) :: TPBQS = TPBQScc ! Sigma parameter of the transfer function (A)
real(c_double) :: TPBQRcutoff1 = TPBQRcutoff1cc ! (A)
! Auxilary 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 = sizeof ( TPBU ) + sizeof ( TPBdUdR )
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) :: Z, RR
integer(c_int) :: i
!-------------------------------------------------------------------------------------------
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 )
elseif ( R > TPBRcutoff ) then
TPBQInt0 = 0.0d+00
return
endif
RR = ( R - TPBRmin ) / TPBDR
i = int ( RR )
RR = RR - i
Z = 1.0d+00 - 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) :: Z, RR
integer(c_int) :: i
!-------------------------------------------------------------------------------------------
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 )
elseif ( R > TPBRcutoff ) then
TPBUInt0 = 0.0d+00
return
endif
RR = ( R - TPBRmin ) / TPBDR
i = int ( RR )
RR = RR - i
Z = 1.0d+00 - 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) :: Z, RR
integer(c_int) :: i
!-------------------------------------------------------------------------------------------
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 )
elseif ( R > TPBRcutoff ) then
TPBU = 0.0d+00
TPBdUdR = 0.0d+00
return
endif
RR = ( R - TPBRmin ) / TPBDR
i = int ( RR )
RR = RR - i
Z = 1.0d+00 - 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
!-------------------------------------------------------------------------------------------
real(c_double) :: Z, t, S
!-------------------------------------------------------------------------------------------
if ( R > TPBRcutoff ) then
TPBQCalc0 = 0.0d+00
else if ( R < TPBQR0 ) then
TPBQCalc0 = 1.0d+00
else
Z = TPBQS / R
Z = Z * Z * Z
Z = Z * Z
TPBQCalc0 = 4.0d+00 * ( 1.0d+00 - Z ) * Z
if ( R > TPBQRcutoff1 ) then
t = ( R - TPBQRcutoff1 ) / TPBQDRcutoff
S = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
TPBQCalc0 = TPBQCalc0 * S
endif
endif
end function TPBQCalc0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) function TPBUCalc0 ( R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: R
!-------------------------------------------------------------------------------------------
real(c_double) :: Z, t, S
!-------------------------------------------------------------------------------------------
if ( R > TPBRcutoff ) then
TPBUCalc0 = 0.0d+00
else
Z = TPBS / R
Z = Z * Z * Z
Z = Z * Z
TPBUCalc0 = TPBE4 * ( Z - 1.0d+00 ) * Z
if ( R > TPBRcutoff1 ) then
t = ( R - TPBRcutoff1 ) / TPBDRcutoff
S = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
TPBUCalc0 = TPBUCalc0 * S
endif
endif
end function TPBUCalc0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPBUCalc1 ( U, dUdR, R ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(out) :: U, dUdR
real(c_double), intent(in) :: R
real(c_double) :: Z, t, S, dSdR
!-------------------------------------------------------------------------------------------
if ( R > TPBRcutoff ) then
U = 0.0d+00
dUdR = 0.0d+00
else
Z = TPBS / R
Z = Z * Z * Z
Z = Z * Z
U = TPBE4 * ( Z - 1.0d+00 ) * Z
dUdR = TPBE24 * ( 2.0d+00 * Z - 1.0d+00 ) * Z / R
if ( R > TPBRcutoff1 ) then
t = ( R - TPBRcutoff1 ) / TPBDRcutoff
S = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
dSdR = 6.0d+00 * t * ( t - 1.0d+00 ) / TPBDRcutoff
dUdR = dUdR * S + U * dSdR
U = U * S
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
real(c_double), intent(in) :: L
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2) :: FF, MM, FFF
!-------------------------------------------------------------------------------------------
FF = 0.5d+00 * F
MM = M / L
call V3_V3xxV3 ( FFF, MM, Laxis )
F1 = FF - FFF
F2 = FF + FFF
end subroutine TPBSegmentForces !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Printing
!---------------------------------------------------------------------------------------------------
! subroutine TPBPrint ( FileName ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! character(c_char)*(*), intent(in) :: FileName
! !-------------------------------------------------------------------------------------------
! integer(c_int) :: Fuid
! integer(c_int) :: i
! real(c_double) :: R
! !-------------------------------------------------------------------------------------------
! Fuid = OpenFile ( FileName, "wt", outputpath )
! write ( Fuid, '(a)' ) 'TITLE="TPB Potentials"'
! write ( Fuid, '(a)' ) 'VARIABLES="R" "Q" "U" "dUdR"'
! write ( Fuid, '(a)' ) 'ZONE'
! R = TPBRmin
! do i = 0, TPBN - 1
! write ( Fuid, '(4e22.12)' ) R, TPBQ(i), TPBU(i), TPBDUDR(i)
! R = R + TPBDR
! end do
! call CloseFile ( Fuid )
! end subroutine TPBPrint !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Initialization
!---------------------------------------------------------------------------------------------------
subroutine TPBInit () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) :: R
integer(c_int) :: i
!-------------------------------------------------------------------------------------------
TPBE4 = 4.0d+00 * TPBE
TPBE24 = - 24.0d+00 * TPBE
TPBDRcutoff = TPBRcutoff - TPBRcutoff1
TPBQDRcutoff = TPBRcutoff - TPBQRcutoff1
TPBQR0 = TPBQS * 2.0d+00 ** ( 1.0d+00 / 6.0d+00 )
TPBDR = ( TPBRcutoff - TPBRmin ) / ( TPBN - 1 )
R = TPBRmin
do i = 0, TPBN - 1
TPBQ(i) = TPBQCalc0 ( R )
call TPBUCalc1 ( TPBU(i), TPBdUdR(i), R )
R = R + TPBDR
enddo
end subroutine TPBInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TubePotBase !****************************************************************************

1675
lib/mesont/TubePotMono.f90 Normal file
View File

File diff suppressed because it is too large Load Diff

442
lib/mesont/TubePotTrue.f90 Normal file
View File

@ -0,0 +1,442 @@
! ------------ ----------------------------------------------------------
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
! http://lammps.sandia.gov, Sandia National Laboratories
! Steve Plimpton, sjplimp@sandia.gov
!
! Copyright (2003) Sandia Corporation. Under the terms of Contract
! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
! certain rights in this software. This software is distributed under
! 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 !********************************************************************************
!
! TMD Library: True tubular potential and transfer function
!
!---------------------------------------------------------------------------------------------------
!
! Intel Fortran
!
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017
!
!---------------------------------------------------------------------------------------------------
!
! This module implements calculation of true potential and transfer functions for interaction
! between two cylinder segments of nanotubes by direct integration over the surfaces of both
! segments.
!
!***************************************************************************************************
use TPMGeom
use TubePotBase
use iso_c_binding, only : c_int, c_double, c_char
implicit none
!---------------------------------------------------------------------------------------------------
! Constants
!---------------------------------------------------------------------------------------------------
integer(c_int), parameter :: TPTNXMAX = 257
integer(c_int), parameter :: TPTNEMAX = 128
!---------------------------------------------------------------------------------------------------
! Types
!---------------------------------------------------------------------------------------------------
type TPTSEG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double) :: X, Y, Z
real(c_double) :: Psi, Theta, Phi ! Euler's angles
real(c_double) :: R ! Segment radius
real(c_double) :: L ! Segment length
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
end type TPTSEG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Global variables
!---------------------------------------------------------------------------------------------------
type(TPTSEG) :: TPTSeg1, TPTSeg2 ! Two segments
contains !******************************************************************************************
subroutine TPTSegAxisVector ( S, Laxis ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
type(TPTSEG), intent(in) :: S
real(c_double), dimension(0:2), intent(out) :: Laxis
!-------------------------------------------------------------------------------------------
Laxis(0:2) = S%M(2,0:2)
end subroutine TPTSegAxisVector !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPTSegRadVector ( S, Lrad, Eps ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
type(TPTSEG), intent(in) :: S
real(c_double), dimension(0:2), intent(out) :: Lrad
real(c_double), intent(in) :: Eps
!-------------------------------------------------------------------------------------------
real(c_double) :: Ce, Se
!-------------------------------------------------------------------------------------------
Ce = cos ( Eps )
Se = sin ( Eps )
Lrad(0) = Ce * S%M(0,0) + Se * S%M(1,0)
Lrad(1) = Ce * S%M(0,1) + Se * S%M(1,1)
Lrad(2) = Ce * S%M(0,2) + Se * S%M(1,2)
end subroutine TPTSegRadVector !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPTRadiusVector ( S, R, X, Eps ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
type(TPTSEG), intent(in) :: S
real(c_double), dimension(0:2), intent(out) :: R
real(c_double), intent(in) :: X, Eps
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2) :: Laxis, Lrad
!-------------------------------------------------------------------------------------------
call TPTSegAxisVector ( S, Laxis )
call TPTSegRadVector ( S, Lrad, Eps )
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
!-------------------------------------------------------------------------------------------
real(c_double) :: X, Eps
integer(c_int) :: i, j
!-------------------------------------------------------------------------------------------
X = - S%L / 2.0
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
call TPTRadiusVector ( S, S%Rtab(i,j,0:2), X, Eps )
Eps = Eps + S%DE
end do
X = X + S%DX
end do
end subroutine TPTCalcSegNodeTable !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPTSetSegPosition1 ( S, Rcenter, Laxis, L ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
type(TPTSEG), intent(inout) :: S
real(c_double), dimension(0:2), intent(in) :: Rcenter, Laxis
real(c_double), intent(in) :: L
!-------------------------------------------------------------------------------------------
S%L = L
S%DX = L / ( S%NX - 1 )
call EulerAngles ( S%Psi, S%Theta, Laxis )
S%Phi= 0.0d+00
S%X = Rcenter(0)
S%Y = Rcenter(1)
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) :: L
!-------------------------------------------------------------------------------------------
R = 0.5 * ( R1 + R2 )
Laxis = R2 - R1
L = S_V3norm3 ( Laxis )
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
!-------------------------------------------------------------------------------------------
integer(c_int) :: i, j
real(c_double) :: L1, L2, Displacement, D
real(c_double), dimension(0:2) :: Laxis, Q, R
!-------------------------------------------------------------------------------------------
L2 = S1%L / 2.0
L1 = - L2
call TPTSegAxisVector ( S1, Laxis )
R(0) = S1%X
R(1) = S1%Y
R(2) = S1%Z
do i = 0, S2%NX - 1
do j = 0, S2%NE - 1
call LinePoint ( Displacement, Q, R, Laxis, S2%Rtab(i,j,0:2) )
D = sqrt ( sqr ( Q(0) - S2%Rtab(i,j,0) ) + sqr ( Q(1) - S2%Rtab(i,j,1) ) + sqr ( Q(2) - S2%Rtab(i,j,2) ) )
if ( Displacement > L1 .and. Displacement < L2 .and. D < S1%R ) then
TPTCheckIntersection = 1
return
end if
end do
end do
TPTCheckIntersection = 0
end function TPTCheckIntersection !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPTCalcPointRange ( S, Xmin, Xmax, Re ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
type(TPTSEG), intent(in) :: S
real(c_double), intent(out) :: Xmin, Xmax
real(c_double), dimension(0:2), intent(in) :: Re
!-------------------------------------------------------------------------------------------
real(c_double) :: Displacement, Distance
real(c_double), dimension(0:2) :: Laxis, Q, R
!-------------------------------------------------------------------------------------------
call TPTSegAxisVector ( S, Laxis )
R(0) = S%X
R(1) = S%Y
R(2) = S%Z
call LinePoint ( Displacement, Q, R, Laxis, Re )
Distance = sqrt ( sqr ( Q(0) - Re(0) ) + sqr ( Q(1) - Re(1) ) + sqr ( Q(2) - Re(2) ) ) - S%R
if ( TPBRcutoff < Distance ) then
Xmin = 0.0d+00
Xmax = 0.0d+00
TPTCalcPointRange = 0
return
end if
Distance = sqrt ( TPBRcutoff * TPBRcutoff - Distance * Distance )
Xmin = Displacement - Distance
Xmax = Displacement + Distance
TPTCalcPointRange = 1
end function TPTCalcPointRange !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine TPTGetEnds ( R1_1, R1_2, R2_1, R2_2, X1_1, X1_2, X2_1, X2_2, H, A ) !!!!!!!!!!!!!
real(c_double), dimension(0:2), intent(out) :: R1_1, R1_2, R2_1, R2_2
real(c_double), intent(in) :: X1_1, X1_2, X2_1, X2_2, H, A
!-------------------------------------------------------------------------------------------
R1_1(0) = 0.0d+00
R1_1(1) = 0.0d+00
R1_1(2) = X1_1
R1_2(0) = 0.0d+00
R1_2(1) = 0.0d+00
R1_2(2) = X1_2
R2_1(0) = H
R2_1(1) = - X2_1 * sin ( A )
R2_1(2) = X2_1 * cos ( A )
R2_2(0) = H
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 the atom in position R and
! produced by the segment S.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Q, U
real(c_double), dimension(0:2), intent(out) :: F
real(c_double), dimension(0:2), intent(in) :: R
type(TPTSEG), intent(in) :: S
!-------------------------------------------------------------------------------------------
integer(c_int) :: i, j
real(c_double), dimension(0:2) :: RR, FF
real(c_double) :: QQ, UU, UUU, FFF, Rabs
real(c_double) :: Coeff, Xmin, Xmax, X
!-------------------------------------------------------------------------------------------
TPTPointPotential = 0
Q = 0.0d+00
U = 0.0d+00
F = 0.0d+00
if ( TPTCalcPointRange ( S, Xmin, Xmax, R ) == 0 ) return
X = - S%L / 2.0
do i = 0, S%NX - 1
if ( X > Xmin .and. X < Xmax ) then
QQ = 0.0d+00
UU = 0.0d+00
FF = 0.0d+00
do j = 0, S%NE - 1
RR(0:2) = S%Rtab(i,j,0:2) - R(0:2)
Rabs = S_V3norm3 ( RR )
if ( Rabs < TPBRcutoff ) then
QQ = QQ + TPBQCalc0 ( Rabs )
call TPBUCalc1 ( UUU, FFF, Rabs )
UU = UU + UUU
FFF = FFF / Rabs
FF = FF + FFF * RR
TPTPointPotential = 1
end if
end do
if ( i == 0 .or. i == S%NX - 1 ) then
Q = Q + 0.5d+00 * QQ
U = U + 0.5d+00 * UU
F = F + 0.5d+00 * FF
else
Q = Q + QQ
U = U + UU
F = F + FF
end if
end if
X = X + S%DX
end do
Coeff = TPBD * S%DX * S%R * S%DE
Q = Q * S%DX * S%R * S%DE
U = U * Coeff
F = F * Coeff
end function TPTPointPotential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer(c_int) function TPTSectionPotential ( Q, U, F, M, S, i, Ssource ) !!!!!!!!!!!!!!!!!!!!!!!
! This funcion 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
real(c_double), dimension(0:2), intent(out) :: F, M
type(TPTSEG), intent(in) :: S, Ssource
integer(c_int), intent(in) :: i
!-------------------------------------------------------------------------------------------
integer(c_int) :: j
real(c_double), dimension(0:2) :: R, Fp, Mp, Lrad
real(c_double) :: Qp, Up, Eps
real(c_double) :: Coeff
!-------------------------------------------------------------------------------------------
TPTSectionPotential = 0
Q = 0.0d+00
U = 0.0d+00
F = 0.0d+00
M = 0.0d+00
Eps = 0.0d+00
do j = 0, S%NE - 1
call TPTSegRadVector ( S, Lrad, Eps )
if ( TPTPointPotential ( Qp, Up, Fp, S%Rtab(i,j,0:2), Ssource ) == 1 ) then
Q = Q + Qp
U = U + Up
F = F + Fp
R(0) = S%Rtab(i,j,0) - S%X
R(1) = S%Rtab(i,j,1) - S%Y
R(2) = S%Rtab(i,j,2) - S%Z
call V3_V3xxV3 ( Mp, R, Fp )
M = M + Mp
TPTSectionPotential = 1
end if
Eps = Eps + S%DE
end do
Coeff = TPBD * S%R * S%DE
Q = Q * S%R * S%DE
U = U * Coeff
F = F * Coeff
M = M * Coeff
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
! Ssource and applied to the segment S.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Q, U
real(c_double), dimension(0:2), intent(out) :: F, M
type(TPTSEG), intent(in) :: S, Ssource
integer(c_int) :: i
real(c_double), dimension(0:2) :: Fc, Mc
real(c_double) :: Qc, Uc
!-------------------------------------------------------------------------------------------
TPTSegmentPotential = 0
Q = 0.0d+00
U = 0.0d+00
F = 0.0d+00
M = 0.0d+00
if ( TPTCheckIntersection ( S, Ssource ) == 1 ) then
TPTSegmentPotential = 2
return
end if
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
U = U + 0.5d+00 * Uc
F = F + 0.5d+00 * Fc
M = M + 0.5d+00 * Mc
else
Q = Q + Qc
U = U + Uc
F = F + Fc
M = M + Mc
end if
TPTSegmentPotential = 1
end if
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
real(c_double), intent(in) :: L
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2) :: MM, FF, FFF
!-------------------------------------------------------------------------------------------
FF = 0.5d+00 * F
MM = M / L
call V3_V3xxV3 ( FFF, MM, Laxis )
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 appliend to the ends of segments.
!-------------------------------------------------------------------------------------------
real(c_double), intent(out) :: Q, U
real(c_double), dimension(0:2), intent(out) :: F1_1, F1_2, F2_1, F2_2
real(c_double), dimension(0:2), intent(in) :: R1_1, R1_2, R2_1, R2_2
!-------------------------------------------------------------------------------------------
real(c_double), dimension(0:2) :: R1, R2, Laxis1, Laxis2, DR, F1, M1, F2, M2
real(c_double) :: L1, L2
!-------------------------------------------------------------------------------------------
R1 = 0.5 * ( R1_1 + R1_2 )
R2 = 0.5 * ( R2_1 + R2_2 )
Laxis1 = R1_2 - R1_1
Laxis2 = R2_2 - R2_1
L1 = S_V3norm3 ( Laxis1 )
L2 = S_V3norm3 ( Laxis2 )
Laxis1 = Laxis1 / L1
Laxis2 = Laxis2 / L2
DR = R2 - R1
call TPTSetSegPosition1 ( TPTSeg1, R1, Laxis1, L1 )
call TPTSetSegPosition1 ( TPTSeg2, R2, Laxis2, L2 )
TPTInteractionF = TPTSegmentPotential ( Q, U, F1, M1, TPTSeg1, TPTSeg2 )
if ( TPTInteractionF .ne. 1 ) return
call V3_V3xxV3 ( M2, DR, F1 )
F2 = - F1
M2 = - M1 - M2
call TPTSegmentForces ( F1_1, F1_2, F1, M1, Laxis1, L1 )
call TPTSegmentForces ( F2_1, F2_2, F2, M2, Laxis2, L2 )
end function TPTInteractionF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------------------------------
! Initialization
!---------------------------------------------------------------------------------------------------
subroutine TPTInit ( R1, R2, NX, NE ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(c_double), intent(in) :: R1, R2
integer(c_int), intent(in) :: NX, NE
!-------------------------------------------------------------------------------------------
TPTSeg1%X = 0.0d+00
TPTSeg1%Y = 0.0d+00
TPTSeg1%Z = 0.0d+00
TPTSeg1%Psi = 0.0d+00
TPTSeg1%Theta = 0.0d+00
TPTSeg1%Phi = 0.0d+00
TPTSeg1%R = R1
TPTSeg1%NX = NX
TPTSeg1%NE = NE
TPTSeg1%DE = M_2PI / NE
TPTSeg2%X = 0.0d+00
TPTSeg2%Y = 0.0d+00
TPTSeg2%Z = 0.0d+00
TPTSeg2%Psi = 0.0d+00
TPTSeg2%Theta = 0.0d+00
TPTSeg2%Phi = 0.0d+00
TPTSeg2%R = R2
TPTSeg2%NX = NX
TPTSeg2%NE = NE
TPTSeg2%DE = M_2PI / NE
end subroutine TPTInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module TubePotTrue !****************************************************************************