317 lines
16 KiB
Fortran
317 lines
16 KiB
Fortran
! ------------ ----------------------------------------------------------
|
|
! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
|
! https://www.lammps.org/ 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 !******************************************************************************
|
|
!
|
|
! Calculation of the TMD force field
|
|
!
|
|
!---------------------------------------------------------------------------------------------------
|
|
!
|
|
! PGI Fortran, Intel Fortran
|
|
!
|
|
! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, version 09.01, 2020
|
|
!
|
|
!***************************************************************************************************
|
|
|
|
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 ) !!!!!!!!!!!
|
|
! Interaction energies associated with nodes X1 and X2
|
|
real(c_double), intent(inout) :: U1, U2
|
|
! Forces exerted on nodes X1 and X2
|
|
real(c_double), intent(inout), dimension(0:2) :: F1, F2
|
|
! Contributions of nodes X1 and X2 to the virial stress tensor
|
|
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2
|
|
! Coordinates of the segment nodes
|
|
real(c_double), intent(in), dimension(0:2) :: X1, X2
|
|
! Radius of a nanotube the segment (X1,X2) belongs to
|
|
real(c_double), intent(in) :: R12
|
|
! Equilibrium length of segment (X1,X2)
|
|
real(c_double), intent(in) :: L12
|
|
!-------------------------------------------------------------------------------------------
|
|
integer(c_int) :: ii, jj, Event
|
|
real(c_double) :: U, F, LL, S, Ubcl
|
|
real(c_double), dimension(0:2) :: DX, FF
|
|
!-------------------------------------------------------------------------------------------
|
|
DX = X2 - X1
|
|
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 )
|
|
! Interaction energies associated with nodes X1, X2, and X3
|
|
real(c_double), intent(inout) :: U1, U2, U3
|
|
! Forces exerted on nodes X1, X2, and X3
|
|
real(c_double), intent(inout), dimension(0:2) :: F1, F2, F3
|
|
! Contributions of nodes X1, X2, and X3 to the virial stress tensor
|
|
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2, S3
|
|
! Coordinates of nodes
|
|
real(c_double), intent(in), dimension(0:2) :: X1, X2, X3
|
|
! Radius of nanotube the segment (X1,X2) belongs to
|
|
real(c_double), intent(in) :: R123
|
|
! Equilibrium length of segment (X1,X2) and (X2,X3) (It is assumed to be the same for both segments)
|
|
real(c_double), intent(in) :: L123
|
|
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 energies and components of the virial stress tensor) between a segment
|
|
! (X1,X2) and a sequence of segments which 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 appearance in the nanotube.
|
|
! It means that (X(i),X(i+1)) are either correspond to a real segment or divided by segments
|
|
! that do not belong to a nanotube.
|
|
|
|
! Concept of the extended chain:
|
|
! Let's consider a sequence of nodes (X1,X2,...,XN) forming continuous part of a nanotube.
|
|
! If node Xe precedes 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, the extended chain coincides with (X1,...,XN) and Ee = 0.
|
|
! If the extended chain contains additional node, then non-zero force is exerted on this node.
|
|
|
|
subroutine SegmentTubeForceField ( U1, U2, U, F1, F2, F, Fe, S1, S2, S, Se, X1, X2, R12, N, X, Xe,&
|
|
BBF, R, E1, E2, Ee, TPMType )
|
|
! Number of nodes in array X
|
|
integer(c_int), intent(in) :: N
|
|
! Interaction energies associated with nodes X1 and X2
|
|
real(c_double), intent(inout) :: U1, U2
|
|
! Interaction energies associated with nodes X
|
|
real(c_double), intent(inout), dimension(0:N-1) :: U
|
|
! Forces exerted on nodes X1 and X2
|
|
real(c_double), intent(inout), dimension(0:2) :: F1, F2
|
|
! Forces exerted on nodes X
|
|
real(c_double), intent(inout), dimension(0:2,0:N-1) :: F
|
|
! Force exerted on node Xe (can be updated only if Ee > 0)
|
|
real(c_double), intent(inout), dimension(0:2) :: Fe
|
|
! Contributions of nodes X1 and X2 to the virial stress tensor
|
|
real(c_double), intent(inout), dimension(0:2,0:2) :: S1, S2
|
|
! Contributions of nodes X to the virial stress tensor
|
|
real(c_double), intent(inout), dimension(0:2,0:2,0:N-1) :: S
|
|
! Contributions of node Xe to the virial stress tensor (can be updated only if Ee > 0)
|
|
real(c_double), intent(inout), dimension(0:2,0:2) :: Se
|
|
! Coordinates of the segment nodes
|
|
real(c_double), intent(in), dimension(0:2) :: X1, X2
|
|
! Radius of a nanotube the segment (X1,X2) belongs to
|
|
real(c_double), intent(in) :: R12
|
|
! Coordinates of the nanotube nodes
|
|
real(c_double), intent(in), dimension(0:2,0:N-1) :: X
|
|
! Additional node of the extended chain if Ee > 0
|
|
real(c_double), intent(in), dimension(0:2) :: Xe
|
|
! Bending buckling flags (BBF(i) = 1 in a case of buckling in node i)
|
|
integer(c_int), intent(in), dimension(0:N-1) :: BBF
|
|
! Radius of nanotube X
|
|
real(c_double), intent(in) :: R
|
|
! E1 = 1 if the chain node 0 is a CNT end; E1 = 2 if the chain node N-1 is a CNT end
|
|
integer(c_int), intent(in) :: E1, E2
|
|
! Parameter defining the type of the extended chain (0,1,2)
|
|
integer(c_int), intent(in) :: Ee
|
|
! Type of the tubular potential (0 or 1)
|
|
integer(c_int), intent(in) :: TPMType
|
|
!-------------------------------------------------------------------------------------------
|
|
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
|
|
! Interaction energies associated with nodes X1 and X2
|
|
real(c_double) :: xU1, xU2
|
|
! Interaction energies associated with nodes X
|
|
real(c_double), dimension(0:N-1) :: xU
|
|
! Forces exerted on nodes X1 and X2
|
|
real(c_double), dimension(0:2) :: xF1, xF2
|
|
! Forces exerted on nodes X
|
|
real(c_double), dimension(0:2,0:N-1) :: xF
|
|
! Force exerted on node Xe (can be updated only if Ee > 0)
|
|
real(c_double), dimension(0:2) :: xFe
|
|
!-------------------------------------------------------------------------------------------
|
|
|
|
! 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 ! The first node in the chain is the tube end
|
|
EType1 = .true.
|
|
else
|
|
EType1 = .false.
|
|
end if
|
|
if ( E2 == 1 ) then ! The 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 ! The first node in the extended chain is the tube end
|
|
EType = 3
|
|
else if ( Ee == 2 ) then ! The 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 !**************************************************************************
|