forked from lijiext/lammps
1699 lines
93 KiB
Fortran
1699 lines
93 KiB
Fortran
! ------------ ----------------------------------------------------------
|
|
! 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 TubePotMono !********************************************************************************
|
|
!
|
|
! Approximate tubular potentials and transfer functions for mono-radius tubes.
|
|
!
|
|
!---------------------------------------------------------------------------------------------------
|
|
!
|
|
! Intel Fortran
|
|
!
|
|
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 13.00, 2020
|
|
!
|
|
!---------------------------------------------------------------------------------------------------
|
|
!
|
|
! Four potentials and transfer functions are calculated in this module:
|
|
!
|
|
! 1. SSTP (segment - semi-infinite tube parallel): Linear density of the potential along
|
|
! the segment axis which is produced by a parallel semi-infinite tube. 2D tables for this potential
|
|
! are generated at initialization or can be loaded from a file.
|
|
!
|
|
! 2. STP (segment - tube parallel): Linear density of the potential along the segment axis
|
|
! which is produced by a parallel infinite tubes. This is only a particular case of the SSTP potential,
|
|
! but it is considered separately for computational efficiency. 1D tables of this potential are taken
|
|
! from 2D tables of SSTP potential.
|
|
!
|
|
! 3. SST (segment - semi-infinite tube): Potential for a segment produced by an arbitrary-
|
|
! oriented semi-infinite tube. This potential can not be kept in 2D tables, therefore, all
|
|
! data are calculated 'on fly' with the help of SSTP potential and numerical integration along the
|
|
! segment axis
|
|
!
|
|
! 4. ST (segment - tube): Potential for a segment produced by an arbitrary-oriented
|
|
! infinitely long tube. 2D tables for this potential are generated at initialization or can be
|
|
! loaded from a file.
|
|
!
|
|
!***************************************************************************************************
|
|
|
|
use TPMLib
|
|
use TPMGeom
|
|
use TubePotBase
|
|
use TubePotTrue
|
|
use LinFun2
|
|
use Spline2
|
|
use iso_c_binding, only : c_int, c_double, c_char
|
|
implicit none
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
! Constants
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
integer(c_int), parameter :: TPMNZMAX = 129
|
|
integer(c_int), parameter :: TPMNEMAX = 128
|
|
|
|
integer(c_int), parameter :: TPMNHMAX = 1001
|
|
integer(c_int), parameter :: TPMNXMAX = 1001
|
|
integer(c_int), parameter :: TPMNMAX = 1001
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
! Global variables
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
integer(c_int) :: TPMStartMode = 1
|
|
character*512 :: TPMFile = 'MESONT-TABTP.xrs'
|
|
integer(c_int) :: TPMUnitID ! Unit for the tabulated potential file
|
|
|
|
integer(c_int) :: TPMNZ = TPMNZMAX
|
|
integer(c_int) :: TPMNZ1 = TPMNZMAX - 1
|
|
integer(c_int) :: TPMNE = TPMNEMAX
|
|
integer(c_int) :: TPMNE1 = TPMNEMAX - 1
|
|
|
|
integer(c_int) :: TPMNH = TPMNHMAX
|
|
integer(c_int) :: TPMNH1 = TPMNHMAX - 1
|
|
integer(c_int) :: TPMNX = TPMNXMAX
|
|
integer(c_int) :: TPMNX1 = TPMNXMAX - 1
|
|
|
|
integer :: TPMChiIndM ! Chirality index M
|
|
integer :: TPMChiIndN ! Chirality index N
|
|
real(c_double) :: TPMR1
|
|
real(c_double) :: TPMR2
|
|
|
|
real(c_double) :: TPMHmax
|
|
real(c_double) :: TPMDH
|
|
|
|
! Parameters of empirical correction functions
|
|
|
|
integer(c_int) :: TPMAN = 20
|
|
real(c_double) :: TPMAHmin
|
|
real(c_double) :: TPMAHmax
|
|
real(c_double) :: TPMADH
|
|
real(c_double), dimension(0:TPMNHMAX-1) :: TPMAH, TPMAF, TPMAFxx
|
|
|
|
! Fitting parameters that depend on the SWCNT chirality
|
|
|
|
real(c_double) :: TPMCaA = 0.22d+00 ! 0.22 for (10,10) CNTs
|
|
real(c_double) :: TPMCeA = 0.35d+00 ! 0.35 for (10,10) CNTs
|
|
real(c_double) :: TPMAHmin0 = 10.0d+00 ! 10.0 A for (10,10) CNTs
|
|
|
|
! Parameters of SSTP integrator
|
|
|
|
real(c_double) :: TPMDE
|
|
real(c_double), dimension(0:TPMNEMAX-1) :: TPMCE, TPMSE
|
|
|
|
! Additional parameters for SSTP potential
|
|
|
|
real(c_double) :: TPMSSTPDelta = 0.25d+00
|
|
integer(c_int) :: TPMSSTPNH
|
|
integer(c_int) :: TPMSSTPNX
|
|
|
|
real(c_double) :: TPMSSTPX1
|
|
real(c_double) :: TPMSSTPXmax
|
|
real(c_double) :: TPMSSTPDX
|
|
|
|
real(c_double), dimension(0:TPMNHMAX-1,0:TPMNXMAX-1) :: TPMSSTPG
|
|
real(c_double), dimension(0:TPMNHMAX-1,0:TPMNXMAX-1) :: TPMSSTPF, TPMSSTPFxx, TPMSSTPFyy, TPMSSTPFxxyy
|
|
real(c_double), dimension(0:TPMNHMAX-1) :: TPMSSTPH
|
|
real(c_double), dimension(0:TPMNXMAX-1) :: TPMSSTPX
|
|
|
|
! Additional parameters for STP potential
|
|
|
|
! In calculations of this potential, some parameters of SSTP potential are also used.
|
|
! In particular, STP potential has no its own integrator. All data comes from SSTP integrator.
|
|
! It does not result in any computational inefficiency unless the STP potential is used without SSTP one.
|
|
|
|
integer(c_int) :: TPMNN = 10
|
|
real(c_double), dimension(0:TPMNHMAX-1) :: TPMSTPG
|
|
real(c_double), dimension(0:TPMNHMAX-1) :: TPMSTPF, TPMSTPFxx
|
|
|
|
! Parameters for ST potential
|
|
|
|
! Minimum gap dh for ST-potential
|
|
real(c_double) :: TPMSTDelta = 1.0d+00
|
|
! Number of subdivisions for every grid step in ST-integrator
|
|
integer(c_int) :: TPMSTNXS = 10
|
|
real(c_double) :: TPMSTXmax
|
|
real(c_double) :: TPMSTH1
|
|
real(c_double) :: TPMSTH2
|
|
real(c_double) :: TPMSTDH12
|
|
|
|
real(c_double), dimension(0:TPMNHMAX-1,0:TPMNXMAX-1) :: TPMSTG
|
|
real(c_double), dimension(0:TPMNHMAX-1,0:TPMNXMAX-1) :: TPMSTF, TPMSTFxx, TPMSTFyy, TPMSTFxxyy
|
|
real(c_double), dimension(0:TPMNHMAX-1) :: TPMSTH
|
|
real(c_double), dimension(0:TPMNXMAX-1) :: TPMSTX
|
|
|
|
! Switching parameters
|
|
|
|
! Height switch (at H=0 in SST-potential)
|
|
integer(c_int) :: TPMHSwitch = 0 ! 1, use h-switch; 0, do not use the switch
|
|
real(c_double) :: TPMHS = 3.0d+00 ! Switch height, Angstrom
|
|
|
|
! Angle switch
|
|
integer(c_int) :: TPMASwitch = 0 ! 1, use a-switch; 0, do not use the switch
|
|
real(c_double) :: TPMAS = 3.0d+00 ! Switch angle, degree
|
|
real(c_double) :: TPMASMin
|
|
real(c_double) :: TPMASMax
|
|
real(c_double) :: TPMASDelta
|
|
|
|
! These variables are used to print error message if intertube force filed fails
|
|
integer(c_int) :: Err_CNT1 = 0, Err_CNT1_Node = 0, Err_CNT2 = 0, &
|
|
Err_CNT2_Node1 = 0, Err_CNT2_Node2 = 0, Err_EType = 0
|
|
|
|
contains !******************************************************************************************
|
|
|
|
integer(c_int) function TPMsizeof () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
TPMsizeof = 8 * ( size ( TPMAH ) + size ( TPMAF ) + size ( TPMAFxx ) &
|
|
+ size ( TPMCE ) + size ( TPMSE ) + size ( TPMSSTPG ) + size ( TPMSSTPF ) &
|
|
+ size ( TPMSSTPFxx ) + size ( TPMSSTPFyy ) + size ( TPMSSTPFxxyy ) &
|
|
+ size ( TPMSSTPH ) + size ( TPMSSTPX ) + size ( TPMSTPG ) + size ( TPMSTPF ) &
|
|
+ size ( TPMSTPFxx ) + size ( TPMSTG ) + size ( TPMSTF ) + size ( TPMSTFxx ) &
|
|
+ size ( TPMSTFyy ) + size ( TPMSTFxxyy ) + size ( TPMSTH ) + size ( TPMSTX ) )
|
|
end function TPMsizeof !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
! Printing error message if intertube force field fails
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
subroutine PrintTPErrMsg () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
!write ( TPErrMsg, fmt = '(a,i8,a,i8,a,i8,a,i8,a,i8,a,i1)' ) 'CNT ', Err_CNT1, ' [', &
|
|
! Err_CNT1_Node,'] with CNT ', Err_CNT2, ' [', Err_CNT2_Node1, ', ', Err_CNT2_Node2, '] E=', Err_EType
|
|
!call PrintStdLogMsg ( TPErrMsg )
|
|
end subroutine PrintTPErrMsg !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
! SSTP: Linear potential density for the tube interacting with parallel semi-infinite tube
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
subroutine TPMSSTPIntegrator ( Q, U, H, D ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function calculates the transfer function Q and potential U between an infinitely long
|
|
! tube and a cross-section of another parallel tube for given height H and displacement D.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U
|
|
real(c_double), intent(in) :: H, D
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i, j, k
|
|
real(c_double) :: C, Zmin, Zmax, DZ, R1X, R1Y, R2X, R2Y, R2Z, R, Rcutoff2
|
|
!-------------------------------------------------------------------------------------------
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
Zmin = D - TPBRcutoff
|
|
Zmax = D + TPBRcutoff
|
|
Rcutoff2 = TPBRcutoff * TPBRcutoff
|
|
if ( Zmin < 0.0d+00 ) Zmin = 0.0d+00
|
|
DZ = ( Zmax - Zmin ) / TPMNZ1
|
|
do i = 0, TPMNE1 ! Integration over the section of the first tube
|
|
R1X = TPMR1 * TPMCE(i)
|
|
R1Y = TPMR1 * TPMSE(i)
|
|
do j = 0, TPMNE1 ! !Integration over the section of the second tube
|
|
R2X = TPMR1 * TPMCE(j) + H
|
|
R2Y = TPMR1 * TPMSE(j)
|
|
R2Z = Zmin - D
|
|
do k = 0, TPMNZ1 ! Integration over the length of the second tube
|
|
R = sqr ( R2X - R1X ) + sqr ( R2Y - R1Y ) + sqr ( R2Z )
|
|
if ( R < Rcutoff2 ) then
|
|
R = dsqrt ( R )
|
|
if ( k == 0 .or. k == TPMNZ1 ) then
|
|
Q = Q + 0.5d+00 * TPBQCalc0 ( R )
|
|
U = U + 0.5d+00 * TPBUCalc0 ( R )
|
|
else
|
|
Q = Q + TPBQCalc0 ( R )
|
|
U = U + TPBUCalc0 ( R )
|
|
end if
|
|
end if
|
|
R2Z = R2Z + DZ
|
|
end do
|
|
end do
|
|
end do
|
|
C = sqr ( TPMDE ) * TPMR1 * TPMR2 * DZ
|
|
Q = Q * C
|
|
U = U * sqr ( TPBD ) * C
|
|
end subroutine TPMSSTPIntegrator !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSSTPInt0 ( Q, U, H, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function returns the transfer function Q and potential U for the SSTP potential
|
|
! calculated by interpolation in the table without switch.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U
|
|
real(c_double), intent(in) :: H, X
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i, j
|
|
real(c_double) :: XX
|
|
!-------------------------------------------------------------------------------------------
|
|
i = 1 + int ( H / TPMDH )
|
|
j = 1 + int ( ( X + TPMSSTPXMax ) / TPMSSTPDX )
|
|
if ( ( i < TPMSSTPNH ) .and. ( j > TPMSSTPNX ) ) then
|
|
!call PrintTPErrMsg ()
|
|
!!call PrintStdLogMsg (TPErrMsg )
|
|
!write ( TPErrMsg, '(a,e20.10,a,e20.10)' ) 'Tubes intersect each other: H=', H, ', X=', X
|
|
!call Error ( 'TPMSSTPInt0', TPErrMsg )
|
|
end if
|
|
if ( i > TPMNH1 ) then
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
TPMSSTPInt0 = 0
|
|
return
|
|
end if
|
|
if ( X .le. - TPMSSTPXmax ) then
|
|
j = 1
|
|
XX = - TPMSSTPXmax
|
|
else if ( X .ge. TPMSSTPXmax ) then
|
|
j = TPMNX1
|
|
XX = TPMSSTPXmax
|
|
else
|
|
XX = X
|
|
end if
|
|
Q = CalcLinFun2_0 ( i, j, H, XX, TPMNH, TPMNX, TPMSSTPH, TPMSSTPX, TPMSSTPG )
|
|
U = CalcSpline2_0 ( i, j, H, XX, TPMNH, TPMNX, TPMSSTPH, TPMSSTPX, TPMSSTPF, TPMSSTPFxx, TPMSSTPFyy, TPMSSTPFxxyy )
|
|
TPMSSTPInt0 = 1
|
|
end function TPMSSTPInt0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSSTPInt0S ( Q, U, H, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function returns the transfer function Q and potential U for the SSTP potential
|
|
! calculated by interpolation in the table and switch to the case of zero H.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U
|
|
real(c_double), intent(in) :: H, X
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: IntSign
|
|
real(c_double) :: t, W, Qa, Ua
|
|
!-------------------------------------------------------------------------------------------
|
|
if ( TPMHSwitch == 0 ) then
|
|
TPMSSTPInt0S = TPMSSTPInt0 ( Q, U, H, X )
|
|
else
|
|
if ( H > TPMHS ) then
|
|
TPMSSTPInt0S = TPMSSTPInt0 ( Q, U, H, X )
|
|
else
|
|
t = H / TPMHS
|
|
W = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
|
|
TPMSSTPInt0S = TPMSSTPInt0 ( Q, U, H, X )
|
|
IntSign = TPMSSTPInt0 ( Qa, Ua, 0.0d+00, X )
|
|
Q = W * Qa + ( 1.0d+00 - W ) * Q
|
|
U = W * Ua + ( 1.0d+00 - W ) * U
|
|
end if
|
|
end if
|
|
end function TPMSSTPInt0S !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSSTPInt1 ( Q, U, Uh, Ux, H, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function returns the transfer function Q, potential U, and derivatives Uh=dU/dH and
|
|
! Ux=dU/dX for the SSTP potential calculated by interpolation in the table without switch
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U, Uh, Ux
|
|
real(c_double), intent(in) :: H, X
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i, j
|
|
real(c_double) :: XX
|
|
!-------------------------------------------------------------------------------------------
|
|
i = 1 + int ( H / TPMDH )
|
|
j = 1 + int ( ( X + TPMSSTPXMax ) / TPMSSTPDX )
|
|
if ( ( i < TPMSSTPNH ) .and. ( j > TPMSSTPNX ) ) then
|
|
!call PrintTPErrMsg ()
|
|
!!call PrintStdLogMsg ( TPErrMsg )
|
|
!write ( TPErrMsg, '(a,e20.10,a,e20.10)' ) 'Tubes intersect each other: H=', H, ', X=', X
|
|
!call Error ( 'TPMSSTPInt1', TPErrMsg )
|
|
end if
|
|
if ( i > TPMNH1 ) then
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
Uh = 0.0d+00
|
|
Ux = 0.0d+00
|
|
TPMSSTPInt1 = 0
|
|
return
|
|
end if
|
|
if ( X .le. - TPMSSTPXmax ) then
|
|
j = 1
|
|
XX = - TPMSSTPXmax
|
|
else if ( X .ge. TPMSSTPXmax ) then
|
|
j = TPMNX1
|
|
XX = TPMSSTPXmax
|
|
else
|
|
XX = X
|
|
end if
|
|
Q = CalcLinFun2_0 ( i, j, H, XX, TPMNH, TPMNX, TPMSSTPH, TPMSSTPX, TPMSSTPG )
|
|
call CalcSpline2_1 ( U, Uh, Ux, i, j, H, XX, TPMNH, TPMNX, TPMSSTPH, TPMSSTPX, TPMSSTPF, &
|
|
TPMSSTPFxx, TPMSSTPFyy, TPMSSTPFxxyy )
|
|
TPMSSTPInt1 = 1
|
|
end function TPMSSTPInt1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSSTPInt1S ( Q, U, Uh, Ux, H, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function returns the transfer function Q, potential U, and derivatives Uh=dU/dH and
|
|
! Ux=dU/dX for the SSTP potential calculated by interpolation in the table and switch to
|
|
! the case of zero H.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U, Uh, Ux
|
|
real(c_double), intent(in) :: H, X
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: IntSign
|
|
real(c_double) :: t, W, W1, dWdH, Qa, Ua, Uha, Uxa
|
|
!-------------------------------------------------------------------------------------------
|
|
if ( TPMHSwitch == 0 ) then
|
|
TPMSSTPInt1S = TPMSSTPInt1 ( Q, U, Uh, Ux, H, X )
|
|
else
|
|
if ( H > TPMHS ) then
|
|
TPMSSTPInt1S = TPMSSTPInt1 ( Q, U, Uh, Ux, H, X )
|
|
else
|
|
t = H / TPMHS
|
|
W = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
|
|
dWdH = 6.0d+00 * t * ( t - 1.0d+00 ) / TPMHS
|
|
TPMSSTPInt1S = TPMSSTPInt1 ( Q, U, Uh, Ux, H, X )
|
|
IntSign = TPMSSTPInt1 ( Qa, Ua, Uha, Uxa, 0.0d+00, X )
|
|
W1 = 1.0d+00 - W
|
|
Q = W * Qa + W1 * Q
|
|
U = W * Ua + W1 * U
|
|
Uh = W1 * Uh + ( Ua - U ) * dWdH
|
|
Ux = W * Uxa + W1 * Ux
|
|
end if
|
|
end if
|
|
end function TPMSSTPInt1S !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMSSTPWrite () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function writes the table of the SSTP potential to a disk file.
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i, j
|
|
!-------------------------------------------------------------------------------------------
|
|
write ( unit = TPMUnitID, fmt = '(4i8)' ) TPMChiIndM, TPMChiIndN, TPMNH1, TPMNX1
|
|
do i = 0, TPMNH1
|
|
do j = 0, TPMNX1
|
|
if ( ( i .ge. TPMSSTPNH ) .or. ( j .le. TPMSSTPNX ) ) &
|
|
write ( unit = TPMUnitID, fmt = '(2e26.17)' ) TPMSSTPG(i,j), TPMSSTPF(i,j)
|
|
end do
|
|
end do
|
|
end subroutine TPMSSTPWrite !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMSSTPRead () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function reads the table of the SSTP potential from a disk file.
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i, j
|
|
integer(c_int) :: iTPMChiIndM, iTPMChiIndN, iTPMNH1, iTPMNX1
|
|
!-------------------------------------------------------------------------------------------
|
|
read ( unit = TPMUnitID, fmt = '(4i8)' ) iTPMChiIndM, iTPMChiIndN, iTPMNH1, iTPMNX1
|
|
if ( iTPMChiIndM .NE. TPMChiIndM .OR. iTPMChiIndN .NE. TPMChiIndN ) then
|
|
print *, 'ERROR in [TPMSSTPRead]: iTPMChiIndM .NE. TPMChiIndM .OR. iTPMChiIndN .NE. TPMChiIndN'
|
|
stop
|
|
end if
|
|
if ( iTPMNH1 .NE. TPMNH1 .OR. iTPMNX1 .NE. TPMNX1 ) then
|
|
print *, 'ERROR in [TPMSSTPRead]: iTPMNH1 .NE. TPMNH1 .OR. iTPMNX1 .NE. TPMNX1'
|
|
stop
|
|
end if
|
|
do i = 0, TPMNH1
|
|
do j = 0, TPMNX1
|
|
if ( ( i .ge. TPMSSTPNH ) .or. ( j .le. TPMSSTPNX ) ) &
|
|
read ( unit = TPMUnitID, fmt = '(2e26.17)' ) TPMSSTPG(i,j), TPMSSTPF(i,j)
|
|
end do
|
|
end do
|
|
end subroutine TPMSSTPRead !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMSSTPInit () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function calculates the table of the SSTP potential.
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i, j
|
|
real(c_double) :: E
|
|
character(c_char) :: Msg
|
|
real(c_double), dimension(0:TPMNMAX-1) :: FF, DD, MM, K0, K1, K2
|
|
!-------------------------------------------------------------------------------------------
|
|
TPMDE = M_2PI / TPMNE
|
|
E = 0.0d+00
|
|
do i = 0, TPMNE1
|
|
TPMCE(i) = cos ( E )
|
|
TPMSE(i) = sin ( E )
|
|
E = E + TPMDE
|
|
end do
|
|
do i = 0, TPMNH1
|
|
TPMSSTPH(i) = TPMDH * i
|
|
end do
|
|
TPMSSTPX1 = - 2.0d+00 * TPMSSTPDelta
|
|
TPMSSTPXmax = TPBRcutoff + TPMSSTPDelta
|
|
TPMSSTPDX = 2.0 * TPMSSTPXmax / TPMNX1
|
|
do j = 0, TPMNX1
|
|
TPMSSTPX(j) = - TPMSSTPXmax + TPMSSTPDX * j
|
|
end do
|
|
TPMSSTPNH = 1 + int ( ( TPMR1 + TPMR2 + TPMSSTPDelta ) / TPMDH )
|
|
TPMSSTPNX = int ( ( TPMSSTPXMax + TPMSSTPX1 ) / TPMSSTPDX ) - 1
|
|
if ( TPMStartMode == 0 ) then
|
|
do i = 0, TPMNH1
|
|
do j = 0, TPMNX1
|
|
if ( ( i .ge. TPMSSTPNH ) .or. ( j .le. TPMSSTPNX ) ) then
|
|
call TPMSSTPIntegrator ( TPMSSTPG(i,j), TPMSSTPF(i,j), TPMSSTPH(i), TPMSSTPX(j) )
|
|
print '(2i5,a,e20.10,a,e20.10,a,e20.10,a,e20.10)', i, j, ' H=', TPMSSTPH(i), &
|
|
', X=', TPMSSTPX(j), ', Q=', TPMSSTPG(i,j), ', U=', TPMSSTPF(i,j)
|
|
end if
|
|
end do
|
|
end do
|
|
call TPMSSTPWrite ()
|
|
else
|
|
call TPMSSTPRead ()
|
|
end if
|
|
call CreateSpline2Ext ( 3, 3, 3, 3, TPMNH, TPMSSTPNH, TPMNX, TPMSSTPNX, TPMNMAX, TPMSSTPH, TPMSSTPX, &
|
|
TPMSSTPF, TPMSSTPFxx, TPMSSTPFyy, TPMSSTPFxxyy, FF, MM, DD, K0, K1, K2 )
|
|
end subroutine TPMSSTPInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
! STP potential for an infinite tube interacting with a parallel segment. No actual initialization
|
|
! is necessary for this potential, since the data are taken from the table for the SSTP potential.
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
integer(c_int) function TPMSTPInt0 ( Q, U, H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function returns the transfer function Q and potential U for the STP potential
|
|
! calculated by interpolation in the table.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U
|
|
real(c_double), intent(in) :: H
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i
|
|
!-------------------------------------------------------------------------------------------
|
|
i = 1 + int ( H / TPMDH )
|
|
if ( i < TPMSSTPNH ) then
|
|
!call PrintTPErrMsg ()
|
|
!!call PrintStdLogMsg ( TPErrMsg )
|
|
!write ( TPErrMsg, '(a,e20.10)' ) 'Tubes intersect each other: H=', H
|
|
!call Error ( 'TPMSTPInt0', TPErrMsg )
|
|
end if
|
|
if ( H > TPMHmax ) then
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
TPMSTPInt0 = 0
|
|
return
|
|
end if
|
|
if ( i == TPMNH ) i = TPMNH - 1
|
|
Q = CalcLinFun1_0 ( i, H, TPMNH, TPMSSTPH, TPMSTPG )
|
|
U = CalcSpline1_0 ( i, H, TPMNH, TPMSSTPH, TPMSTPF, TPMSTPFxx )
|
|
TPMSTPInt0 = 1
|
|
end function TPMSTPInt0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSTPInt1 ( Q, U, dUdH, H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function returns the transfer function Q, potential U, and derivative dUdH for
|
|
! the STP potential calculated by interpolation in the table.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U, dUdH
|
|
real(c_double), intent(in) :: H
|
|
integer(c_int) :: i
|
|
!-------------------------------------------------------------------------------------------
|
|
i = 1 + int ( H / TPMDH )
|
|
if ( i < TPMSSTPNH ) then
|
|
!call PrintTPErrMsg ()
|
|
!!call PrintStdLogMsg ( TPErrMsg )
|
|
!write ( TPErrMsg, '(a,e20.10)' ) 'Tubes intersect each other: H=', H
|
|
!call Error ( 'TPMSTPInt0', TPErrMsg )
|
|
end if
|
|
if ( H > TPMHmax ) then
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
dUdH = 0.0d+00
|
|
TPMSTPInt1 = 0
|
|
return
|
|
end if
|
|
Q = CalcLinFun1_0 ( i, H, TPMNH, TPMSSTPH, TPMSTPG )
|
|
call CalcSpline1_1 ( U, dUdH, i, H, TPMNH, TPMSSTPH, TPMSTPF, TPMSTPFxx )
|
|
TPMSTPInt1 = 1
|
|
end function TPMSTPInt1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMSTPInit () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function initializes the table of the STP potential
|
|
!-------------------------------------------------------------------------------------------
|
|
TPMSTPG(0:TPMNH1) = TPMSSTPG(0:TPMNH1,TPMNX1)
|
|
TPMSTPF(0:TPMNH1) = TPMSSTPF(0:TPMNH1,TPMNX1)
|
|
TPMSTPFxx(0:TPMNH1) = TPMSSTPFxx(0:TPMNH1,TPMNX1)
|
|
end subroutine TPMSTPInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
! Fitting functions for the SST and ST potentials.
|
|
! This correction functions are chosen empirically to improve accuracy of the SST and ST potentials.
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
subroutine TPMAInit ( X1_1, X1_2, X2_1, X2_2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(in) :: X1_1, X1_2, X2_1, X2_2
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), dimension(0:2) :: R1_1, R1_2, R2_1, R2_2
|
|
real(c_double), dimension(0:2) :: Fa, Ma
|
|
real(c_double) :: Qa, Ua, Qb, Ub, X, H, HH, Ucoeff, Uamin, Ubmin
|
|
integer(c_int) :: i, j, IntSign
|
|
real(c_double), dimension(0:TPMNHMAX-1) :: D, K0, K1, K2
|
|
integer(c_int) :: iTPMChiIndM, iTPMChiIndN, iTPMAN
|
|
!-------------------------------------------------------------------------------------------
|
|
TPMAHmin = TPMR1 + TPMR2 + TPMSTDelta
|
|
TPMAHmax = TPMR1 + TPMR2 + 0.95d+00 * TPBRcutoff
|
|
TPMADH = ( TPMAHmax - TPMAHmin ) / ( TPMAN - 1 )
|
|
if ( TPMStartMode == 1 ) then
|
|
read ( unit = TPMUnitID, fmt = '(4i8)' ) iTPMChiIndM, iTPMChiIndN, iTPMAN
|
|
if ( iTPMChiIndM .NE. TPMChiIndM .OR. iTPMChiIndN .NE. TPMChiIndN ) then
|
|
print *, 'ERROR in [TPMAInit]: iTPMChiIndM .NE. TPMChiIndM .OR. iTPMChiIndN .NE. TPMChiIndN'
|
|
stop
|
|
end if
|
|
if ( iTPMAN .NE. TPMAN ) then
|
|
print *, 'ERROR in [TPMAInit]: iTPMAN .NE. TPMAN'
|
|
stop
|
|
end if
|
|
do i = 0, TPMAN - 1
|
|
TPMAH(i) = TPMAHmin + i * TPMADH
|
|
read ( unit = TPMUnitID, fmt = * ) TPMAF(i)
|
|
end do
|
|
call CreateSpline1 ( 3, 3, TPMAN, TPMAH, TPMAF, TPMAFxx, D, K0, K1, K2 )
|
|
return
|
|
end if
|
|
call TPTInit ( TPMR1, TPMR2, TPTNXMAX, TPTNEMAX )
|
|
do i = 0, TPMAN - 1
|
|
TPMAH(i) = TPMAHmin + i * TPMADH
|
|
call TPTGetEnds ( R1_1, R1_2, R2_1, R2_2, X1_1, X1_2, X2_1, X2_2, TPMAH(i), M_PI_2 )
|
|
X = - ( X1_2 - X1_1 ) / 2.0d+00
|
|
do j = 0, ( TPTNXMAX - 1 ) / 2
|
|
HH = sqrt ( TPMAH(i) * TPMAH(i) + sqr ( X * sin ( M_PI_2 ) ) )
|
|
IntSign = TPMSTPInt0 ( Qb, Ub, HH )
|
|
call TPTSetSegPosition2 ( TPTSeg1, R1_1, R1_2 )
|
|
call TPTSetSegPosition2 ( TPTSeg2, R2_1, R2_2 )
|
|
IntSign = TPTSectionPotential ( Qa, Ua, Fa, Ma, TPTSeg1, j, TPTSeg2 )
|
|
if ( j == 0 ) then
|
|
Uamin = Ua
|
|
Ubmin = Ub
|
|
else
|
|
if ( Ua < Uamin ) then
|
|
Uamin = Ua
|
|
end if
|
|
if ( Ub < Ubmin ) then
|
|
Ubmin = Ub
|
|
end if
|
|
end if
|
|
X = X + TPTSeg1%DX
|
|
end do
|
|
TPMAF(i) = Uamin / Ubmin
|
|
end do
|
|
write ( unit = TPMUnitID, fmt = '(4i8)' ) TPMChiIndM, TPMChiIndN, TPMAN
|
|
do i = 0, TPMAN - 1
|
|
write ( unit = TPMUnitID, fmt = * ) TPMAF(i)
|
|
end do
|
|
call CreateSpline1 ( 3, 3, TPMAN, TPMAH, TPMAF, TPMAFxx, D, K0, K1, K2 )
|
|
end subroutine TPMAInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
real(c_double) function TPMA0 ( H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(in) :: H
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i
|
|
real(c_double) :: A0, t, S
|
|
!-------------------------------------------------------------------------------------------
|
|
if ( H > TPMAHmax ) then
|
|
TPMA0 = 1.0d+00
|
|
return
|
|
else if ( H < TPMAHmin ) then
|
|
if ( H < TPMAHmin0 ) then
|
|
TPMA0 = 1.5d+00
|
|
return
|
|
end if
|
|
A0 = 1.5d+00
|
|
t = ( H - TPMAHmin0 ) / TPMAHmin
|
|
S = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
|
|
TPMA0 = ( 1.0d+00 - S ) * CalcSpline1_0 ( 1, H, TPMAN, TPMAH, TPMAF, TPMAFxx ) + A0 * S
|
|
return
|
|
end if
|
|
i = 1 + int ( ( H - TPMAHmin ) / TPMADH )
|
|
TPMA0 = CalcSpline1_0 ( i, H, TPMAN, TPMAH, TPMAF, TPMAFxx )
|
|
end function TPMA0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMA1 ( A, Ah, H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(out) :: A, Ah
|
|
real(c_double), intent(in) :: H
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i
|
|
real(c_double) :: A0, t, S, dSdH
|
|
!-------------------------------------------------------------------------------------------
|
|
if ( H > TPMAHmax ) then
|
|
A = 1.0d+00
|
|
Ah = 0.0d+00
|
|
return
|
|
else if ( H < TPMAHmin ) then
|
|
if ( H < TPMAHmin0 ) then
|
|
A = 1.5d+00
|
|
Ah = 0.0d+00
|
|
return
|
|
end if
|
|
A0 = 1.5d+00
|
|
t = ( H - TPMAHmin0 ) / TPMAHmin
|
|
S = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
|
|
dSdH = 6.0d+00 * t * ( t - 1.0d+00 ) / TPMAHmin
|
|
call CalcSpline1_1 ( A, Ah, 1, H, TPMAN, TPMAH, TPMAF, TPMAFxx )
|
|
Ah = Ah * ( 1.0d+00 - S ) + dSdH * ( A0 - A )
|
|
A = A * ( 1.0d+00 - S ) + A0 * S
|
|
return
|
|
end if
|
|
i = 1 + int ( ( H - TPMAHmin ) / TPMADH )
|
|
call CalcSpline1_1 ( A, Ah, i, H, TPMAN, TPMAH, TPMAF, TPMAFxx )
|
|
end subroutine TPMA1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
real(c_double) function TPMCu0 ( H, cosA, sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function returns the correction function for the magnitude of the potential.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(in) :: H, cosA, sinA
|
|
!-------------------------------------------------------------------------------------------
|
|
TPMCu0 = 1.0d+00 + ( TPMA0 ( H ) - 1.0d+00 ) * sqr ( sinA )
|
|
end function TPMCu0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMCu1 ( Cu, CuH, CuA, H, cosA, sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! The subroutine calculates the correction function Cu for the magnitude of the potential and
|
|
! its derivatives CuH and CuA.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(ouT) :: Cu, CuH, CuA
|
|
real(c_double), intent(in) :: H, cosA, sinA
|
|
real(c_double) :: AA, AAh, D
|
|
!-------------------------------------------------------------------------------------------
|
|
call TPMA1 ( AA, AAh, H )
|
|
D = sqr ( sinA )
|
|
AA = AA - 1.0d+00
|
|
Cu = 1.0d+00 + AA * D
|
|
CuH = AAh * D
|
|
CuA = AA * 2.0d+0 * cosA * sinA
|
|
end subroutine TPMCu1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
real(c_double) function TPMCa0 ( cosA, sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function returns the correction function for the argument of the potential.
|
|
! If correction is not necessary, it should return sinA.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(in) :: cosA, sinA
|
|
!-------------------------------------------------------------------------------------------
|
|
TPMCa0 = sinA / ( 1.0d+00 - TPMCaA * sqr ( sinA ) )
|
|
end function TPMCa0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMCa1 ( Ca, CaA, Ka, KaA, cosA, sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This subroutine calculates the correction function Cu for the depth of the potential well
|
|
! and its derivatives CuH and CuA. If correction is not necessary, it returns Ca = sinA
|
|
! and CaA = cosA.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Ca, CaA, Ka, KaA
|
|
real(c_double), intent(in) :: cosA, sinA
|
|
!-------------------------------------------------------------------------------------------
|
|
Ka = 1.0d+00 / ( 1.0d+00 - TPMCaA * sqr ( sinA ) )
|
|
Ca = sinA * Ka
|
|
KaA = 2.0d+00 * TPMCaA * sinA * cosA * sqr ( Ka )
|
|
CaA = cosA * Ka + sinA * KaA
|
|
end subroutine TPMCa1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
real(c_double) function TPMCe0 ( sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! This function returns the correction function for the argument of the potential.
|
|
! If correction is not necessary, it returns sinA.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(in) :: sinA
|
|
!-------------------------------------------------------------------------------------------
|
|
TPMCe0 = 1.0d+00 - TPMCeA * sinA * sinA
|
|
end function TPMCe0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMCe1 ( Ce, CeA, Ke, cosA, sinA ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! If correction is not necessary, it returns Ce = 1 and CeA = 0.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Ce, CeA, Ke
|
|
real(c_double), intent(in) :: cosA, sinA
|
|
!-------------------------------------------------------------------------------------------
|
|
Ce = 1.0d+00 - TPMCeA * sinA * sinA
|
|
CeA = - 2.0d+00 * TPMCeA * sinA * cosA
|
|
Ke = TPMCeA
|
|
end subroutine TPMCe1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
! SST potential for the semi-infinite tube interacting with segment.
|
|
! This potential does not need any initialization. All necessary data is taken from tables of the
|
|
! SSTP potential.
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
integer(c_int) function TPMSSTPotential ( Q, U, X1, X2, H, cosA, D, N ) !!!!!!!!!!!!!!!!!!!!
|
|
! This function calculates the transfer function Q and potential U applied to a segment
|
|
! from a semi-infinite tube based on the numerical integration (trapezoidal rule) along the segment
|
|
! axis for non-parallel objects.
|
|
! Relative position of the nanotube and segment is given by axial positions of the segment
|
|
! ends X1 and X2, height H, cosA= cos(A), where A is the cross-axis angle, and the displacement
|
|
! D of a nanotube end.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U
|
|
real(c_double), intent(in) :: X1, X2, H, cosA, D
|
|
integer(c_int), intent(in) :: N ! Number of nodes for numerical integration
|
|
real(c_double) :: sinA, Qs, Us, DX, X, XX, HH, Cu, Ca, Ce
|
|
integer(c_int) :: i
|
|
!-------------------------------------------------------------------------------------------
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
DX = ( X2 - X1 ) / ( N - 1 )
|
|
X = X1
|
|
sinA = dsqrt ( 1.0d+00 - cosA * cosA )
|
|
Cu = TPMCu0 ( H, cosA, sinA )
|
|
Ca = TPMCa0 ( cosA, sinA )
|
|
Ce = TPMCe0 ( sinA )
|
|
TPMSSTPotential = 0
|
|
do i = 0, N - 1
|
|
XX = X * cosA - Ce * D
|
|
HH = sqrt ( H * H + sqr ( X * Ca ) )
|
|
if ( TPMSSTPInt0S ( Qs, Us, HH, XX ) > 0 ) TPMSSTPotential = 1
|
|
if ( i == 0 .or. i == N - 1 ) then
|
|
Q = Q + 0.5d+00 * Qs
|
|
U = U + 0.5d+00 * Us
|
|
else
|
|
Q = Q + Qs
|
|
U = U + Us
|
|
end if
|
|
X = X + DX
|
|
end do
|
|
Q = Cu * Q * DX
|
|
U = Cu * U * DX
|
|
end function TPMSSTPotential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSSTPotentialPar ( Q, U, R1_1, Laxis1, R2_1, Laxis2, L1, N ) !!!!!
|
|
! Potential for a segment and a semi-infinite tube is calculated by the numerical
|
|
! integration (trapezoidal rule) along the segment axis for parallel objects.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U
|
|
real(c_double), dimension(0:2), intent(in) :: R1_1, Laxis1, R2_1, Laxis2
|
|
real(c_double), intent(in) :: L1
|
|
integer(c_int), intent(in) :: N ! Number of nodes for numerical integration
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double) :: Qs, Us, DX, X, S, H
|
|
real(c_double), dimension(0:2) :: R1, L12
|
|
integer(c_int) :: i
|
|
!-------------------------------------------------------------------------------------------
|
|
DX = L1 / ( N - 1 )
|
|
X = 0.0d+00
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
TPMSSTPotentialPar = 0
|
|
do i = 0, N - 1
|
|
R1 = R1_1 + X * Laxis1
|
|
call LinePoint ( S, L12, R2_1, Laxis2, R1 )
|
|
L12 = L12 - R1
|
|
call ApplyPeriodicBC ( L12 )
|
|
H = S_V3norm3 ( L12 )
|
|
if ( TPMSSTPInt0S ( Qs, Us, H, S ) > 0 ) then
|
|
TPMSSTPotentialPar = 1
|
|
if ( i == 0 .or. i == TPMNN - 1 ) then
|
|
Q = Q + 0.5d+00 * Qs
|
|
U = U + 0.5d+00 * Us
|
|
else
|
|
Q = Q + Qs
|
|
U = U + Us
|
|
end if
|
|
X = X + DX
|
|
end if
|
|
end do
|
|
Q = Q * DX
|
|
U = U * DX
|
|
end function TPMSSTPotentialPar !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSSTForces ( Q, U, F1, F2, Fd, X1, X2, H, cosA, D, N ) !!!!!!!!!!!
|
|
! Potential and forces applied to the segment from a semi-infinite tube are calculated
|
|
! by the numerical integration (trapezoidal rule) along the segment axis.
|
|
! Non-parallel case.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U, Fd
|
|
real(c_double), dimension(0:2), intent(out) :: F1, F2
|
|
real(c_double), intent(in) :: X1, X2, H, cosA, D
|
|
integer(c_int), intent(in) :: N ! Number of nodes for numerical integration
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double) :: DX, sinA
|
|
real(c_double) :: Qs, Us, Ush, Usx, Fx, Fy, Fz
|
|
real(c_double) :: C, C1, C2, I0, Ih, Ih1, Ih2, Ix, Ix1, X, XX, HH
|
|
real(c_double) :: Ca, CaA, Ka, KaA, Cu, CuH, CuA, Ce, CeA, Ke, Uh, Ua
|
|
integer(c_int) :: IntSign, i
|
|
!-------------------------------------------------------------------------------------------
|
|
I0 = 0.0d+00
|
|
Ih = 0.0d+00
|
|
Ih1 = 0.0d+00
|
|
Ih2 = 0.0d+00
|
|
Ix = 0.0d+00
|
|
Ix1 = 0.0d+00
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
F1 = 0.0d+00
|
|
F2 = 0.0d+00
|
|
Fd = 0.0d+00
|
|
sinA = dsqrt ( 1.0d+00 - cosA * cosA )
|
|
X = X1
|
|
DX = ( X2 - X1 ) / ( N - 1 )
|
|
TPMSSTForces = 0
|
|
call TPMCa1 ( Ca, CaA, Ka, KaA, cosA, sinA )
|
|
call TPMCu1 ( Cu, CuH, CuA, H, cosA, sinA )
|
|
call TPMCe1 ( Ce, CeA, Ke, cosA, sinA )
|
|
do i = 0, N - 1
|
|
XX = X * cosA - Ce * D
|
|
HH = sqrt ( H * H + sqr ( X * Ca ) )
|
|
if ( TPMSSTPInt1S ( Qs, Us, Ush, Usx, HH, XX ) > 0 ) TPMSSTForces = 1
|
|
if ( i == 0 .or. i == N - 1 ) then
|
|
Qs = 0.5d+0 * Qs
|
|
Us = 0.5d+0 * Us
|
|
Ush = 0.5d+0 * Ush
|
|
Usx = 0.5d+0 * Usx
|
|
end if
|
|
if ( HH .le. TPGeomPrec ) then
|
|
Ush = 0.0d+00
|
|
else
|
|
Ush = Ush / HH
|
|
end if
|
|
Q = Q + Qs
|
|
I0 = I0 + Us
|
|
Ih = Ih + Ush
|
|
Ih1 = Ih1 + X * Ush
|
|
Ih2 = Ih2 + X * X * Ush
|
|
Ix = Ix + Usx
|
|
Ix1 = Ix1 + X * Usx
|
|
X = X + DX
|
|
end do
|
|
if ( TPMSSTForces == 0 ) return
|
|
|
|
C = DX * Cu
|
|
I0 = I0 * C
|
|
Ih = Ih * C
|
|
Ih1 = Ih1 * C
|
|
Ih2 = Ih2 * C
|
|
Ix = Ix * C
|
|
Ix1 = Ix1 * C
|
|
|
|
DX = X2 - X1
|
|
|
|
Q = Q * C
|
|
U = I0
|
|
Uh = ( CuH * U / Cu + h * Ih ) / DX
|
|
Ua = ( CuA * I0 / Cu + Ca * CaA * Ih2 - sinA * Ix1 - CeA * D * Ix ) / DX
|
|
|
|
C1 = Ka * Ka * Ih1
|
|
C = h * ( C1 + cosA * Ke * Ix ) / DX
|
|
F1(0) = X2 * Uh - C
|
|
F2(0) = C - X1 * Uh
|
|
|
|
C = ( cosA * sinA * C1 + ( Ke * sinA - sinA ) * Ix ) / DX
|
|
F1(1) = Ua - X2 * C
|
|
F2(1) = X1 * C - Ua
|
|
|
|
C1 = Ca * Ca * Ih1 + cosA * Ix
|
|
C2 = Ca * Ca * Ih2 + cosA * Ix1
|
|
F1(2) = ( U - X2 * C1 + C2 ) / DX
|
|
F2(2) = ( X1 * C1 - C2 - U ) / DX
|
|
|
|
Fd = Ce * Ix
|
|
end function TPMSSTForces !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSSTForcesPar ( Q, U, F1, F2, Fd, R1_1, Laxis1, R2_1, Laxis2, L1, N )
|
|
! Potential and forces applied to the segment from a semi-infinite tube are calculated by
|
|
! the numerical integration (trapezoidal rule) along the segment axis.
|
|
! Parallel case
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(out) :: Q, U, Fd
|
|
real(c_double), dimension(0:2), intent(out) :: F1, F2
|
|
real(c_double), dimension(0:2), intent(in) :: R1_1, Laxis1, R2_1, Laxis2
|
|
real(c_double), intent(in) :: L1
|
|
integer(c_int), intent(in) :: N ! Number of nodes for numerical integration
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double) :: Qs, Us, Ush, Usx, DX, X, S, H, Beta, Gamma
|
|
real(c_double), dimension(0:2) :: R1, L12, Fs
|
|
integer(c_int) :: i, N1
|
|
!-------------------------------------------------------------------------------------------
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
F1 = 0.0d+00
|
|
F2 = 0.0d+00
|
|
Fd = 0.0d+00
|
|
X = 0.0d+00
|
|
N1 = N - 1
|
|
DX = L1 / N1
|
|
TPMSSTForcesPar = 0
|
|
do i = 0, N1
|
|
R1 = R1_1 + X * Laxis1
|
|
call LinePoint ( S, L12, R2_1, Laxis2, R1 )
|
|
L12 = L12 - R1
|
|
call ApplyPeriodicBC ( L12 )
|
|
H = S_V3norm3 ( L12 )
|
|
if ( TPMSSTPInt1S ( Qs, Us, Ush, Usx, H, S ) > 0 ) then
|
|
TPMSSTForcesPar = 1
|
|
if ( H .ge. TPGeomPrec ) then
|
|
Fs = Ush * L12 / H - Usx * Laxis2
|
|
else
|
|
Fs = - Usx * Laxis2
|
|
end if
|
|
Beta = real ( i ) / N1
|
|
Gamma = 1.0d+00 - Beta
|
|
if ( i == 0 .or. i == N1 ) then
|
|
Q = Q + 0.5d+00 * Qs
|
|
U = U + 0.5d+00 * Us
|
|
Fd = Fd + 0.5d+00 * Usx
|
|
Gamma = 0.5d+00 * Gamma
|
|
Beta = 0.5d+000 * Beta
|
|
else
|
|
Q = Q + Qs
|
|
U = U + Us
|
|
Fd = Fd + Usx
|
|
end if
|
|
F1 = F1 + Gamma * Fs
|
|
F2 = F2 + Beta * Fs
|
|
end if
|
|
X = X + DX
|
|
end do
|
|
Q = Q * DX
|
|
U = U * DX
|
|
Fd = Fd * DX
|
|
Fs = U * Laxis1 / L1
|
|
F1 = DX * F1 + Fs
|
|
F2 = DX * F2 - Fs
|
|
end function TPMSSTForcesPar !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
! ST: Potential for a infinite tube interacting with a segment
|
|
!--------------------------------------------------------------------------------------------------
|
|
|
|
!
|
|
! These functions are used to smooth boundaries in (H,X) domain for ST potential
|
|
!
|
|
|
|
real(c_double) function TPMSTXMin0 ( H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(in) :: H
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double) :: X
|
|
!-------------------------------------------------------------------------------------------
|
|
if ( H < TPMSTH1 ) then
|
|
TPMSTXMin0 = sqrt ( TPMSTH2 * TPMSTH2 - H * H )
|
|
return
|
|
else if ( H > TPMSTH2 ) then
|
|
TPMSTXMin0 = 0.0d+00
|
|
return
|
|
end if
|
|
X = ( H - TPMSTH1 ) / TPMSTDH12
|
|
TPMSTXMin0 = sqrt ( TPMSTH2 * TPMSTH2 - H * H ) &
|
|
* ( 1.0d+00 - X * X * X * ( 3.0d+00 * X * ( 2.0d+00 * X - 5.0d+00 ) + 10.0d+00 ) )
|
|
end function TPMSTXMin0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
real(c_double) function TPMSTXMax0 ( H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(in) :: H
|
|
!-------------------------------------------------------------------------------------------
|
|
TPMSTXMax0 = sqrt ( TPMSTXMax * TPMSTXMax - H * H )
|
|
end function TPMSTXMax0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMSTXMin1 ( XMin, dXMindH, H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(out) :: XMin, dXMindH
|
|
real(c_double), intent(in) :: H
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double) :: X, F, dFdX
|
|
!-------------------------------------------------------------------------------------------
|
|
if ( H < TPMSTH1 ) then
|
|
XMin = sqrt ( TPMSTH2 * TPMSTH2 - H * H )
|
|
dXMindH = - H / XMin
|
|
return
|
|
else if ( H > TPMSTH2 ) then
|
|
XMin = 0.0d+00
|
|
dXMindH = 0.0d+00
|
|
return
|
|
end if
|
|
X = ( H - TPMSTH1 ) / TPMSTDH12
|
|
F = 1.0d+00 - X * X * X * ( 3.0d+00 * X * ( 2.0d+00 * X - 5.0d+00 ) + 10.0d+00 )
|
|
X = X * ( X - 1.0d+00 )
|
|
dFdX = - 30.0d+00 * X * X
|
|
XMin = sqrt ( TPMSTH2 * TPMSTH2 - H * H )
|
|
dXMindH = - H * F / XMin + XMin * dFdX / TPMSTDH12
|
|
XMin = XMin * F
|
|
end subroutine TPMSTXMin1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMSTXMax1 ( XMax, dXMaxdH, H ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(out) :: XMax, dXMaxdH
|
|
real(c_double), intent(in) :: H
|
|
!-------------------------------------------------------------------------------------------
|
|
XMax = sqrt ( TPMSTXMax * TPMSTXMax - H * H )
|
|
dXMaxdH = - H / XMax
|
|
end subroutine TPMSTXMax1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
!
|
|
! ST (segment-tube) table
|
|
!
|
|
|
|
subroutine TPMSTIntegrator ( G, F, Q, U, H, X, DX ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(inout) :: G, F, Q, U
|
|
real(c_double), intent(in) :: H, X, DX
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double) :: FFx, HH, DDX
|
|
integer(c_int) :: IntSign
|
|
!-------------------------------------------------------------------------------------------
|
|
DDX = 0.5 * DX
|
|
G = G + Q * DDX
|
|
F = F + U * DDX
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
HH = dsqrt ( sqr ( H ) + sqr ( X ) )
|
|
if ( HH > TPMHmax ) return
|
|
IntSign = TPMSTPInt0 ( Q, U, HH )
|
|
if ( IntSign == 1 ) then
|
|
G = G + Q * DDX
|
|
F = F + U * DDX
|
|
end if
|
|
end subroutine TPMSTIntegrator !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSTInt0 ( G, F, H, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(out) :: G, F
|
|
real(c_double), intent(in) :: H, X
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i, j
|
|
real(c_double) :: S, XA, XXX, XXXX, XMin, XMax
|
|
!-------------------------------------------------------------------------------------------
|
|
if ( H > TPMHmax ) then
|
|
G = 0.0d+00
|
|
F = 0.0d+00
|
|
TPMSTInt0 = 0
|
|
return
|
|
else if ( H < 0.0d+00 ) then
|
|
G = 0.0d+00
|
|
F = 0.0d+00
|
|
TPMSTInt0 = 2
|
|
!call PrintTPErrMsg ()
|
|
!!call PrintStdLogMsg ( TPErrMsg )
|
|
!all Error ( 'TPMSTInt0', 'H < 0' )
|
|
!!return
|
|
end if
|
|
S = signum ( X )
|
|
XA = dabs ( X )
|
|
i = 1 + int ( H / TPMDH )
|
|
if ( i > TPMNH1 ) i = TPMNH1
|
|
XMin = TPMSTXMin0 ( H )
|
|
XMax = TPMSTXMax0 ( H )
|
|
XXX = ( XA - XMin ) / ( XMax - XMin )
|
|
if ( XXX < 0.0d+00 ) then
|
|
j = 1
|
|
XXXX = 0.0d+00
|
|
!call PrintTPErrMsg ()
|
|
!write ( TPErrMsg, '(a,e20.10,a,e20.10,a,e20.10)' ) 'X < XMin, X=', XA, ', XMin=', XMin, ', H=', H
|
|
!call Error ( 'TPMSTInt0', TPErrMsg )
|
|
else if ( XXX .ge. 1.0d+00 ) then
|
|
j = TPMNX1
|
|
XXXX = 1.0d+00
|
|
else
|
|
XXXX = XXX
|
|
j = 1 + int ( XXXX * TPMNX1 )
|
|
end if
|
|
G = S * CalcLinFun2_0 ( i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, TPMSTG )
|
|
F = S * CalcSpline2_0 ( i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, TPMSTF, TPMSTFxx, TPMSTFyy, TPMSTFxxyy )
|
|
TPMSTInt0 = 1
|
|
end function TPMSTInt0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSTInt1 ( G, F, Fh, Fx, H, X ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(inout) :: G, F, Fh, Fx
|
|
real(c_double), intent(in) :: H, X
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i, j
|
|
real(c_double) :: S, XA, DX, XXX, XXXX, XMin, XMax, dXMindH, dXMaxdH
|
|
!-------------------------------------------------------------------------------------------
|
|
if ( H > TPMHmax ) then
|
|
G = 0.0d+00
|
|
F = 0.0d+00
|
|
Fh = 0.0d+00
|
|
Fx = 0.0d+00
|
|
TPMSTInt1 = 0
|
|
return
|
|
else if ( H < 0.0d+00 ) then
|
|
G = 0.0d+00
|
|
F = 0.0d+00
|
|
Fh = 0.0d+00
|
|
Fx = 0.0d+00
|
|
TPMSTInt1 = 2
|
|
!call PrintTPErrMsg ()
|
|
!!call PrintStdLogMsg ( trim ( TPErrMsg ) )
|
|
!call Error ( 'TPMSTInt1', 'H < 0' )
|
|
!!return
|
|
end if
|
|
S = signum ( X )
|
|
XA = dabs ( X )
|
|
i = 1 + int ( H / TPMDH )
|
|
if ( i > TPMNH1 ) i = TPMNH1
|
|
call TPMSTXMin1 ( XMin, dXMindH, H )
|
|
call TPMSTXMax1 ( XMax, dXMaxdH, H )
|
|
DX = XMax - XMin
|
|
XXX = ( XA - XMin ) / DX
|
|
if ( XXX < 0.0d+00 ) then
|
|
j = 1
|
|
XXX = 0.0d+00
|
|
XXXX = 0.0d+00
|
|
!call PrintTPErrMsg ()
|
|
!write ( TPErrMsg, '(a,e20.10,a,e20.10,a,e20.10)' ) 'X < XMin, X=', XA, ', XMin=', XMin, ', H=', H
|
|
!call Error ( 'TPMSTInt', TPErrMsg )
|
|
else if ( XXX .ge. 1.0d+00 ) then
|
|
j = TPMNX1
|
|
XXX = 1.0d+00
|
|
XXXX = 1.0d+00
|
|
else
|
|
XXXX = XXX
|
|
j = 1 + int ( XXXX * TPMNX1 )
|
|
end if
|
|
G = S * CalcLinFun2_0 ( i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, TPMSTG )
|
|
call CalcSpline2_1 ( F, Fh, Fx, i, j, H, XXXX, TPMNH, TPMNX, TPMSTH, TPMSTX, &
|
|
TPMSTF, TPMSTFxx, TPMSTFyy, TPMSTFxxyy )
|
|
Fx = Fx / DX
|
|
Fh = Fh - Fx * ( dXMaxdH * XXX + dXMindH * ( 1.0d+00 - XXX ) )
|
|
F = F * S
|
|
Fh = Fh * S
|
|
TPMSTInt1 = 1
|
|
end function TPMSTInt1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSTPotential ( Q, U, X1, X2, H, cosA, CaseID ) !!!!!!!!!!!!!!!!!!!
|
|
real(c_double), intent(out) :: Q, U
|
|
real(c_double), intent(in) :: X1, X2, H, cosA
|
|
integer(c_int), intent(in) :: CaseID
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double) :: sinA, GG1, GG2, FF1, FF2, Ca, Cu
|
|
!-------------------------------------------------------------------------------------------
|
|
if ( CaseID == MD_LINES_PAR ) then
|
|
TPMSTPotential = TPMSTPInt0 ( Q, U, H )
|
|
U = U * ( X2 - X1 )
|
|
return
|
|
end if
|
|
TPMSTPotential = 0
|
|
sinA = dsqrt ( 1.0d+00 - cosA * cosA )
|
|
Cu = TPMCu0 ( H, cosA, sinA )
|
|
Ca = TPMCa0 ( cosA, sinA )
|
|
if ( TPMSTInt0 ( GG1, FF1, H, X1 * Ca ) > 0 ) TPMSTPotential = 1
|
|
if ( TPMSTInt0 ( GG2, FF2, H, X2 * Ca ) > 0 ) TPMSTPotential = 1
|
|
Q = Cu * ( GG2 - GG1 ) / Ca
|
|
U = Cu * ( FF2 - FF1 ) / Ca
|
|
end function TPMSTPotential !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSTForces ( Q, U, F1, F2, X1, X2, H, cosA, CaseID ) !!!!!!!!!!!!!!
|
|
real(c_double), intent(out) :: Q, U
|
|
real(c_double), dimension(0:2), intent(out) :: F1, F2
|
|
real(c_double), intent(in) :: X1, X2, H, cosA
|
|
integer(c_int), intent(in) :: CaseID
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double) :: DX, sinA
|
|
real(c_double) :: GG1, GG2, FF1, FF2, Fh1, Fh2, Fx1, Fx2
|
|
real(c_double) :: B, C, D
|
|
real(c_double) :: Ca, CaA, Ka, KaA, Cu, CuH, CuA
|
|
integer(c_int) :: IntSign1, IntSign2
|
|
!-------------------------------------------------------------------------------------------
|
|
DX = X2 - X1
|
|
if ( CaseID == MD_LINES_PAR ) then
|
|
F1 = 0.0d+00
|
|
F2 = 0.0d+00
|
|
TPMSTForces = TPMSTPInt1 ( Q, U, F1(0), H )
|
|
F1(0) = F1(0) * 0.5 * DX
|
|
F2(0) = F1(0)
|
|
F1(2) = U
|
|
F2(2) = - U
|
|
Q = Q * DX
|
|
U = U * DX
|
|
return
|
|
end if
|
|
|
|
sinA = dsqrt ( 1.0d+00 - cosA * cosA )
|
|
call TPMCa1 ( Ca, CaA, Ka, KaA, cosA, sinA )
|
|
IntSign1 = TPMSTInt1 ( GG1, FF1, Fh1, Fx1, H, X1 * Ca )
|
|
IntSign2 = TPMSTInt1 ( GG2, FF2, Fh2, Fx2, H, X2 * Ca )
|
|
if ( ( IntSign1 .ne. 1 ) .and. ( IntSign2 .ne. 1 ) ) then
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
F1 = 0.0d+00
|
|
F2 = 0.0d+00
|
|
TPMSTForces = 0
|
|
return
|
|
end if
|
|
|
|
call TPMCu1 ( Cu, CuH, CuA, H, cosA, sinA )
|
|
|
|
Q = Cu * ( GG2 - GG1 ) / Ca
|
|
U = Cu * ( FF2 - FF1 ) / Ca
|
|
|
|
B = Cu * ( Fx2 - Fx1 ) / sinA
|
|
C = H * B / sinA
|
|
D = CuH * U / Cu + Cu * ( Fh2 - Fh1 ) / Ca
|
|
F1(0) = ( X2 * D - C ) / DX
|
|
F2(0) = ( C - X1 * D ) / DX
|
|
|
|
C = cosA * B
|
|
D = ( CuA / Cu - CaA / Ca ) * U + CaA * Cu * ( X2 * Fx2 - X1 * Fx1 ) / Ca
|
|
F1(1) = ( D - X2 * C ) / DX
|
|
F2(1) = ( X1 * C - D ) / DX
|
|
|
|
F1(2) = Cu * Fx1
|
|
F2(2) = - Cu * Fx2
|
|
|
|
TPMSTForces = 1
|
|
end function TPMSTForces !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMSTForceTorque( Qi, Ui, Fi, Ti, Q, U, F, T, Psi, PsiA, Cap, L, H, cosA, CaseID )
|
|
real(c_double), intent(out) :: Qi, Ui, Fi, Ti, Q, U, F, T, Psi, PsiA, Cap
|
|
real(c_double), intent(in) :: L, H, cosA
|
|
integer(c_int), intent(in) :: CaseID
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double) :: L2, sinA
|
|
real(c_double) :: GG, FF, Fh, Fx, GGi, FFi, Fhi, Fxi
|
|
real(c_double) :: B, C, D
|
|
real(c_double) :: Ca, CaA, Ka, KaA, Cu, CuH, CuA
|
|
integer(c_int) :: IntSign
|
|
!-------------------------------------------------------------------------------------------
|
|
if ( CaseID == MD_LINES_PAR ) then
|
|
TPMSTForceTorque = TPMSTPInt1 ( Q, U, F, H )
|
|
Q = Q * L
|
|
U = U * L
|
|
F = - F * L
|
|
T = 0.0d+00
|
|
Qi = 0.0d+00
|
|
Ui = 0.0d+00
|
|
Fi = 0.0d+00
|
|
Ti = 0.0d+00
|
|
Psi = 0.0d+00
|
|
PsiA = 0.0d+00
|
|
Cap = 0.0d+00
|
|
return
|
|
end if
|
|
|
|
L2 = 0.5d+00 * L
|
|
sinA = dsqrt ( 1.0d+00 - cosA * cosA )
|
|
call TPMCa1 ( Ca, CaA, Ka, KaA, cosA, sinA )
|
|
IntSign = TPMSTInt1 ( GG, FF, Fh, Fx, H, L2 * Ca )
|
|
IntSign = TPMSTInt1 ( GGi, FFi, Fhi, Fxi, H, TPMSTXmax )
|
|
if ( IntSign .ne. 1 ) then
|
|
Qi = 0.0d+00
|
|
Ui = 0.0d+00
|
|
Fi = 0.0d+00
|
|
Ti = 0.0d+00
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
F = 0.0d+00
|
|
T = 0.0d+00
|
|
Psi = 0.0d+00
|
|
PsiA = 0.0d+00
|
|
Cap = 0.0d+00
|
|
TPMSTForceTorque = 0
|
|
return
|
|
end if
|
|
|
|
call TPMCu1 ( Cu, CuH, CuA, H, cosA, sinA )
|
|
|
|
Psi = Cu / Ca
|
|
PsiA = ( CuA * Ca - Cu * CaA ) / Ca / Ca
|
|
Cap = CuA / Cu - KaA / Ka - cosA / sinA
|
|
Qi = 2.0d+00 * Psi * GGi
|
|
Ui = 2.0d+00 * Psi * FFi
|
|
Fi = - 2.0d+00 * ( CuH * FFi / Ca + Psi * Fhi )
|
|
Ti = - Cap * Ui
|
|
|
|
Q = 2.0d+00 * Cu * GG / Ca
|
|
U = 2.0d+00 * Cu * FF / Ca
|
|
F = - 2.0d+00 * ( CuH * FF / Ca + Psi * Fh )
|
|
T = - 2.0d+00 * ( ( CuA * Ka - Cu * KaA ) / ( Ka * Ka * sinA ) - Cu * cosA / ( Ka * sinA * sinA ) ) * FF &
|
|
- 2.0d+00 * Cu / ( Ka * sinA ) * Fx * L2 * ( KaA * sinA + Ka * cosA )
|
|
|
|
TPMSTForceTorque = 1
|
|
end function TPMSTForceTorque !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
subroutine TPMSTInit () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
real(c_double) :: X, Q, U, DX, DDX, XMin, XMax
|
|
integer(c_int) :: i, j, k
|
|
real(c_double), dimension(0:TPMNMAX-1) :: FF, DD, MM, K0, K1, K2
|
|
!-------------------------------------------------------------------------------------------
|
|
TPMSTH1 = TPMR1 + TPMR2
|
|
TPMSTH2 = TPMSTH1 + TPMSTDelta
|
|
TPMSTDH12 = TPMSTH2 - TPMSTH1
|
|
TPMSTXmax = TPMHMax + TPMSTDelta
|
|
DX = 1.0 / TPMNX1
|
|
do j = 0, TPMNX1
|
|
TPMSTX(j) = DX * j
|
|
end do
|
|
do i = 0, TPMNH1
|
|
TPMSTH(i) = TPMDH * i
|
|
XMin = TPMSTXmin0 ( TPMSTH(i) )
|
|
XMax = TPMSTXMax0 ( TPMSTH(i) )
|
|
Q = 0.0d+00
|
|
U = 0.0d+00
|
|
DX = ( XMax - XMin ) * TPMSTX(1) / TPMSTNXS
|
|
X = XMin
|
|
call TPMSTIntegrator ( TPMSTG(i,0), TPMSTF(i,0), Q, U, TPMSTH(i), X, DX )
|
|
TPMSTG(i,0) = 0.0d+00
|
|
TPMSTF(i,0) = 0.0d+00
|
|
TPMSTFyy(i,0) = U
|
|
TPMSTFyy(i,TPMNX1) = 0.0d+00
|
|
do j = 1, TPMNX1
|
|
TPMSTG(i,j) = TPMSTG(i,j-1)
|
|
TPMSTF(i,j) = TPMSTF(i,j-1)
|
|
do k = 0, TPMSTNXS - 1
|
|
X = X + DX
|
|
call TPMSTIntegrator ( TPMSTG(i,j), TPMSTF(i,j), Q, U, TPMSTH(i), X, DX )
|
|
end do
|
|
if ( j < TPMNX1 ) DX = ( XMax - XMin ) * ( TPMSTX(j+1) - TPMSTX(j) ) / TPMSTNXS
|
|
end do
|
|
end do
|
|
call CreateSpline2 ( 3, 3, 3, 3, TPMNH, TPMNX, TPMNMAX, TPMSTH, TPMSTX, TPMSTF, TPMSTFxx, &
|
|
TPMSTFyy, TPMSTFxxyy, FF, MM, DD, K0, K1, K2 )
|
|
end subroutine TPMSTInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
! Interaction functions: They can be used for calculation of the potential and forces between a
|
|
! segment and an infinite or semi-infinite nanotube.
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
subroutine TPMSegmentForces ( F2_1, F2_2, F1_1, F1_2, R1_1, R1_2, R2, Laxis2, L2 ) !!!!!!!!!
|
|
real(c_double), dimension(0:2), intent(out) :: F2_1, F2_2
|
|
real(c_double), dimension(0:2), intent(in) :: F1_1, F1_2, R1_1, R1_2, R2, Laxis2
|
|
real(c_double), intent(in) :: L2
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), dimension(0:2) :: F, M, RR
|
|
!-------------------------------------------------------------------------------------------
|
|
RR = R1_1 - R2
|
|
! Taking into account periodic boundary conditions
|
|
call ApplyPeriodicBC ( RR )
|
|
call V3_V3xxV3 ( M, RR, F1_1 )
|
|
RR = R1_2 - R2
|
|
! Taking into account periodic boundary conditions
|
|
call ApplyPeriodicBC ( RR )
|
|
call V3_V3xxV3 ( F, RR, F1_2 )
|
|
M = - ( M + F )
|
|
F = - ( F1_1 + F1_2 )
|
|
call TPBSegmentForces ( F2_1, F2_2, F, M, Laxis2, L2 )
|
|
end subroutine TPMSegmentForces !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
!
|
|
! Interaction of a segment with a semi-infinite or infinite tube
|
|
!
|
|
|
|
integer(c_int) function TPMInteractionF ( Q, U, F1_1, F1_2, F2_1, F2_2, Fd, R1_1, R1_2, R2_1, R2_2, SType2 )
|
|
! SType2 in the type of the second segment:
|
|
! SType2 == 0, internal segment;
|
|
! Stype2 == 1, point R2_1 is the end of the tube;
|
|
! Stype2 == 2, point R2_2 in the end of the tube.
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), intent(inout) :: Q, U, Fd
|
|
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) :: SType2
|
|
real(c_double), dimension(0:2) :: R1, R2, Laxis1, Laxis2, F1, F2, L12, Ly, DR, F1_1a, F1_2a, F1_1b, F1_2b
|
|
real(c_double) :: H, cosA, D1, D2, L1, L2, cosA2, t, W, W1, dWdt, Qa, Ua, Qb, Ub, Fda, Fdb, FF
|
|
integer(c_int) :: GeomID, SwitchID, S, IntSigna, IntSignb
|
|
!-------------------------------------------------------------------------------------------
|
|
R1 = 0.5d+00 * ( R1_1 + R1_2 )
|
|
R2 = 0.5d+00 * ( 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
|
|
L1 = 0.5d+00 * L1
|
|
L2 = 0.5d+00 * L2
|
|
if ( SType2 == 2 ) Laxis2 = - Laxis2
|
|
GeomID = LineLine ( H, cosA, D1, D2, L12, R1, Laxis1, R2, Laxis2, TPGeomPrec )
|
|
|
|
! Angle switch
|
|
if ( TPMASwitch == 0 ) then
|
|
if ( GeomID == MD_LINES_PAR ) then
|
|
SwitchID = 2
|
|
else
|
|
SwitchID = 0
|
|
end if
|
|
else
|
|
cosA2 = cosA * cosA
|
|
if ( cosA2 .ge. TPMASMax .or. GeomID == MD_LINES_PAR ) then
|
|
SwitchID = 2
|
|
else if ( cosA2 .le. TPMASMin ) then
|
|
SwitchID = 0
|
|
else
|
|
t = ( cosA2 - TPMASMin ) / TPMASDelta
|
|
W = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
|
|
dWdt = 6.0d+00 * t * ( t - 1.0d+00 ) / TPMASDelta
|
|
SwitchID = 1
|
|
end if
|
|
end if
|
|
|
|
if ( SwitchID < 2 ) then
|
|
D2 = D2 - L2
|
|
if ( SType2 == 0 ) then
|
|
IntSigna = TPMSTForces ( Qa, Ua, F1, F2, D1 - L1, D1 + L1, H, cosA, MD_LINES_NONPAR )
|
|
Fda = 0.0d+00
|
|
else
|
|
IntSigna = TPMSSTForces ( Qa, Ua, F1, F2, Fda, D1 - L1, D1 + L1, H, cosA, D2, TPMNN )
|
|
end if
|
|
call V3_V3xxV3 ( Ly, Laxis1, Laxis2 )
|
|
S = signum ( S_V3xV3 ( Ly, L12 ) )
|
|
call V3_V3xxV3 ( Ly, Laxis1, L12 )
|
|
Ly = Ly * S
|
|
if ( IntSigna > 0 ) then
|
|
F1_1a = F1(0) * L12 + F1(1) * Ly + F1(2) * Laxis1
|
|
F1_2a = F2(0) * L12 + F2(1) * Ly + F2(2) * Laxis1
|
|
else
|
|
F1_1a = 0.0d+00
|
|
F1_2a = 0.0d+00
|
|
end if
|
|
end if
|
|
|
|
if ( SwitchID > 0 ) then
|
|
if ( SType2 == 0 ) then
|
|
call LinePoint ( H, L12, R2, Laxis2, R1 )
|
|
L12 = L12 - R1
|
|
call ApplyPeriodicBC ( L12 )
|
|
H = S_V3norm3 ( L12 )
|
|
IntSignb = TPMSTForces ( Qb, Ub, F1, F2, - L1, L1, H, cosA, MD_LINES_PAR )
|
|
Fdb = 0.0d+00
|
|
if ( IntSignb > 0 ) then
|
|
if ( H .le. TPGeomPrec ) then
|
|
F1_1b = F1(2) * Laxis1
|
|
F1_2b = F2(2) * Laxis1
|
|
else
|
|
L12 = L12 / H
|
|
F1_1b = F1(0) * L12 + F1(2) * Laxis1
|
|
F1_2b = F2(0) * L12 + F2(2) * Laxis1
|
|
end if
|
|
else
|
|
F1_1b = 0.0d+00
|
|
F1_2b = 0.0d+00
|
|
end if
|
|
else if ( Stype2 == 1 ) then
|
|
IntSignb = TPMSSTForcesPar ( Qb, Ub, F1_1b, F1_2b, Fdb, R1_1, Laxis1, R2_1, Laxis2, &
|
|
2.0d+00 * L1, TPMNN )
|
|
else
|
|
IntSignb = TPMSSTForcesPar ( Qb, Ub, F1_1b, F1_2b, Fdb, R1_1, Laxis1, R2_2, Laxis2, &
|
|
2.0d+00 * L1, TPMNN )
|
|
end if
|
|
end if
|
|
|
|
if ( SwitchID == 0 ) then
|
|
Q = Qa
|
|
U = Ua
|
|
F1_1 = F1_1a
|
|
F1_2 = F1_2a
|
|
Fd = Fda
|
|
TPMInteractionF = IntSigna
|
|
else if ( SwitchID == 2 ) then
|
|
Q = Qb
|
|
U = Ub
|
|
F1_1 = F1_1b
|
|
F1_2 = F1_2b
|
|
Fd = Fdb
|
|
TPMInteractionF = IntSignb
|
|
else
|
|
W1 = 1.0d+00 - W
|
|
Q = W * Qa + W1 * Qb
|
|
U = W * Ua + W1 * Ub
|
|
Ly = Ly * ( Ua - Ub ) * dWdt * cosA * sqrt ( 1.0d+00 - sqr ( cosA ) ) / L1
|
|
F1_1 = W * F1_1a + W1 * F1_1b - Ly
|
|
F1_2 = W * F1_2a + W1 * F1_2b + Ly
|
|
Fd = W * Fda + W1 * Fdb
|
|
TPMInteractionF = 0
|
|
if ( IntSigna > 0 .or. IntSignb > 0 ) TPMInteractionF = 1
|
|
end if
|
|
|
|
! Calculation of forces for the complimentary tube
|
|
if ( SType2 == 2 ) Laxis2 = - Laxis2
|
|
call TPMSegmentForces ( F2_1, F2_2, F1_1, F1_2, R1_1, R1_2, R2, Laxis2, 2.0d+00 * L2 )
|
|
! After the previous subroutine call, F2_1*Laxis2 = F2_2*Laxis2, but this is not true for a semi-infinite tube.
|
|
! The force along the tube should be applied to the end of the tube, while for the
|
|
! another point corresponding force is equal to zero.
|
|
if ( SType2 == 1 ) then
|
|
FF = S_V3xV3 ( F2_1, Laxis2 )
|
|
DR = ( Fd - FF ) * Laxis2
|
|
F2_1 = F2_1 + DR
|
|
F2_2 = F2_2 - DR
|
|
else if ( SType2 == 2 ) then
|
|
FF = S_V3xV3 ( F2_2, Laxis2 )
|
|
DR = ( - Fd - FF ) * Laxis2
|
|
F2_2 = F2_2 + DR
|
|
F2_1 = F2_1 - DR
|
|
end if
|
|
end function TPMInteractionF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMInteractionU ( Q, U, R1_1, R1_2, R2_1, R2_2, SType2 ) !!!!!!!!!!!
|
|
real(c_double), intent(inout) :: Q, U
|
|
real(c_double), dimension(0:2), intent(in) :: R1_1, R1_2, R2_1, R2_2
|
|
integer(c_int), intent(in) :: SType2
|
|
!-------------------------------------------------------------------------------------------
|
|
real(c_double), dimension(0:2) :: R1, R2, Laxis1, Laxis2, F1, F2, L12, DR
|
|
real(c_double) :: H, cosA, D1, D2, L1, L2, cosA2, t, W, Qa, Ua, Qb, Ub
|
|
integer(c_int) :: GeomID, SwitchID, IntSigna, IntSignb
|
|
!-------------------------------------------------------------------------------------------
|
|
R1 = 0.5d+00 * ( R1_1 + R1_2 )
|
|
R2 = 0.5d+00 * ( 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
|
|
if ( SType2 == 2 ) Laxis2 = - Laxis2
|
|
GeomID = LineLine ( H, cosA, D1, D2, L12, R1, Laxis1, R2, Laxis2, TPGeomPrec )
|
|
L1 = 0.5d+00 * L1
|
|
L2 = 0.5d+00 * L2
|
|
|
|
! Angle switch
|
|
if ( TPMASwitch == 0 ) then
|
|
if ( GeomID == MD_LINES_PAR ) then
|
|
SwitchID = 2
|
|
else
|
|
SwitchID = 0
|
|
end if
|
|
else
|
|
cosA2 = cosA * cosA
|
|
if ( cosA2 .ge. TPMASMax .or. GeomID == MD_LINES_PAR ) then
|
|
SwitchID = 2
|
|
else if ( cosA2 .le. TPMASMin ) then
|
|
SwitchID = 0
|
|
else
|
|
t = ( cosA2 - TPMASMin ) / TPMASDelta
|
|
W = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t )
|
|
SwitchID = 1
|
|
end if
|
|
end if
|
|
|
|
if ( SwitchID < 2 ) then
|
|
if ( Stype2 == 0 ) then
|
|
IntSigna = TPMSTPotential ( Qa, Ua, D1 - L1, D1 + L1, H, cosA, MD_LINES_NONPAR )
|
|
else
|
|
IntSigna = TPMSSTPotential ( Qa, Ua, D1 - L1, D1 + L1, H, cosA, D2 - L2, TPMNN )
|
|
end if
|
|
end if
|
|
|
|
if ( SwitchID > 0 ) then
|
|
if ( Stype2 == 0 ) then
|
|
call LinePoint ( H, L12, R2, Laxis2, R1 )
|
|
L12 = L12 - R1
|
|
call ApplyPeriodicBC ( L12 )
|
|
IntSignb = TPMSTPotential ( Qb, Ub, - L1, L1, S_V3norm3 ( L12 ), cosA, MD_LINES_PAR )
|
|
else if ( Stype2 == 1 ) then
|
|
IntSignb = TPMSSTPotentialPar ( Qb, Ub, R1_1, Laxis1, R2_1, Laxis2, 2.0d+00 * L1, TPMNN )
|
|
else
|
|
IntSignb = TPMSSTPotentialPar ( Qb, Ub, R1_1, Laxis1, R2_2, Laxis2, 2.0d+00 * L1, TPMNN )
|
|
end if
|
|
end if
|
|
|
|
if ( SwitchID == 0 ) then
|
|
Q = Qa
|
|
U = Ua
|
|
TPMInteractionU = IntSigna
|
|
else if ( SwitchID == 2 ) then
|
|
Q = Qb
|
|
U = Ub
|
|
TPMInteractionU = IntSignb
|
|
else
|
|
Q = W * Qa + ( 1.0d+00 - W ) * Qb
|
|
U = W * Ua + ( 1.0d+00 - W ) * Ub
|
|
TPMInteractionU = 0
|
|
if ( IntSigna > 0 .or. IntSignb > 0 ) TPMInteractionU = 1
|
|
end if
|
|
end function TPMInteractionU !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
integer(c_int) function TPMInteractionFNum ( Q, U, F1_1, F1_2, F2_1, F2_2, R1_1, R1_2, R2_1, R2_2, Stype2, Delta )
|
|
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) :: SType2
|
|
real(c_double), intent(in) :: Delta
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: i, j, IntSign
|
|
real(c_double) :: QQ, DD, D2
|
|
real(c_double), dimension(0:1,0:2) :: U1_1, U1_2, U2_1, U2_2
|
|
real(c_double), dimension(0:2) :: RR
|
|
!-------------------------------------------------------------------------------------------
|
|
U = 0.0d+00
|
|
F1_1 = 0.0d+00
|
|
F1_2 = 0.0d+00
|
|
F2_1 = 0.0d+00
|
|
F2_2 = 0.0d+00
|
|
TPMInteractionFNum = TPMInteractionU ( Q, U, R1_1, R1_2, R2_1, R2_2, SType2 )
|
|
if ( TPMInteractionFNum == 0 ) return
|
|
D2 = 2.0d+00 * Delta
|
|
do i = 0, 2
|
|
DD = - Delta
|
|
do j = 0 , 1
|
|
RR = R1_1
|
|
RR(i) = RR(i) + DD
|
|
IntSign = TPMInteractionU ( QQ, U1_1(j,i), RR, R1_2, R2_1, R2_2, SType2 )
|
|
RR = R1_2
|
|
RR(i) = RR(i) + DD
|
|
IntSign = TPMInteractionU ( QQ, U1_2(j,i), R1_1, RR, R2_1, R2_2, SType2 )
|
|
RR = R2_1
|
|
RR(i) = RR(i) + DD;
|
|
IntSign = TPMInteractionU ( QQ, U2_1(j,i), R1_1, R1_2, RR, R2_2, SType2 )
|
|
RR = R2_2
|
|
RR(i) = RR(i) + DD
|
|
IntSign = TPMInteractionU ( QQ, U2_2(j,i), R1_1, R1_2, R2_1, RR, SType2 )
|
|
DD = DD + D2
|
|
end do
|
|
end do
|
|
do i = 0, 2
|
|
F1_1(i) = ( U1_1(0,i) - U1_1(1,i) ) / D2
|
|
F1_2(i) = ( U1_2(0,i) - U1_2(1,i) ) / D2
|
|
F2_1(i) = ( U2_1(0,i) - U2_1(1,i) ) / D2
|
|
F2_2(i) = ( U2_2(0,i) - U2_2(1,i) ) / D2
|
|
end do
|
|
end function TPMInteractionFNum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
!---------------------------------------------------------------------------------------------------
|
|
! Initialization
|
|
!---------------------------------------------------------------------------------------------------
|
|
|
|
subroutine TPMInit ( ChiIndM, ChiIndN ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
integer(c_int), intent(in) :: ChiIndM, ChiIndN
|
|
real(c_double) :: RT, DX
|
|
character*512 :: PDate
|
|
!-------------------------------------------------------------------------------------------
|
|
TPPotType = TP_POT_MONO_R
|
|
|
|
! Here we calculate the radius of nanotubes
|
|
RT = TPBAcc * sqrt ( 3.0d+00 * ( ChiIndM * ChiIndM + ChiIndN * ChiIndN + ChiIndM * ChiIndN ) ) / M_2PI
|
|
|
|
TPMChiIndM = ChiIndM
|
|
TPMChiIndN = ChiIndN
|
|
TPMR1 = RT
|
|
TPMR2 = RT
|
|
|
|
TPMCaA = 0.275d+00 * ( 1.0d+00 - 1.0d+00 / ( 1.0d+00 + 0.59d+00 * RT ) )
|
|
TPMCeA = 0.35d+00 + 0.0226d+00 * ( RT - 6.785d+00 )
|
|
TPMAHmin0 = 10.0d+00 * ( RT / 6.785d+00 )**1.5
|
|
|
|
TPMHmax = TPMR1 + TPMR2 + TPBRcutoff
|
|
TPMDH = TPMHmax / TPMNH1
|
|
|
|
! Parameters of the angle switch
|
|
TPMASMin = sqr ( cos ( rad ( TPMAS ) ) )
|
|
TPMASMax = 1.0d+00 - TPGeomPrec
|
|
TPMASDelta = TPMASMax - TPMASMin
|
|
|
|
if ( TPMStartMode == 1 ) then
|
|
TPMUnitID = OpenFile ( TPMFile, 'rt', '' )
|
|
read ( unit = TPMUnitID, fmt = '()' )
|
|
read ( unit = TPMUnitID, fmt = '()' )
|
|
read ( unit = TPMUnitID, fmt = '()' )
|
|
else
|
|
TPMUnitID = OpenFile ( TPMFile, 'wt', '' )
|
|
call fdate( PDate )
|
|
write ( unit = TPMUnitID, fmt = '(a,a)' ) 'DATE ', PDate
|
|
write ( unit = TPMUnitID, fmt = '(a,i3,a,i3,a)' ) &
|
|
'Tabulated data of the tubular potential for (', ChiIndM, ',', ChiIndN, ') CNTs'
|
|
write ( unit = TPMUnitID, fmt = '(a)' ) &
|
|
'A. N. Volkov, L. V. Zhigilei, J. Phys. Chem. C 114, 5513-5531, 2010. doi: 10.1021/jp906142h'
|
|
end if
|
|
|
|
call TPMSSTPInit ()
|
|
|
|
call TPMSTPInit ()
|
|
|
|
DX = TPMR1 + TPMR2 + TPBRcutoff
|
|
call TPMAInit ( - DX, DX, - DX, DX )
|
|
|
|
call TPMSTInit ()
|
|
|
|
call CloseFile ( TPMUnitID )
|
|
|
|
end subroutine TPMInit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
end module TubePotMono !****************************************************************************
|