📄 conrot.f90
字号:
Subroutine conrot4( Nbs, Xb, Yb, Zb, &
Nvib, Nbend, Ntor, &
maxvib, maxbend, maxtor, &
IVIB, IBEND, ITOR, &
RVIB, RBEND, RTOR, &
mapf, BoxSize, dUintra, &
beta, success, Seed, th_h, phi_h )
implicit none
! Nbs is the number of beads in the molecule.
integer, intent(in) :: Nbs
! Xb, Yb, and Zb are the coordinates of the beads.
real, dimension(Nbs), intent(inout) :: Xb, Yb, Zb
integer, intent(in) :: Nvib, Nbend, Ntor
integer, intent(in) :: maxvib, maxbend, maxtor
integer, dimension(2,maxvib) :: IVIB
integer, dimension(3,maxbend) :: IBEND
integer, dimension(4,maxtor) :: ITOR
real, dimension(2,maxvib) :: RVIB
real, dimension(2,maxbend) :: RBEND
real, dimension(4,maxtor) :: RTOR
integer, dimension(-2:4) :: mapf
! BoxSize is the length of the simulation box.
real, intent(in) :: BoxSize
real :: dUintra
real, intent(in) :: beta
logical :: success
! Seed is the current random number generator seed value.
integer, intent(inout) :: Seed
integer, dimension(0:180) :: th_h
integer, dimension(-180:180) :: phi_h
! Local Variables
integer :: h, i, j, k, m
integer :: gate, gate_tmp
integer :: nsol, sol, nb, ibd
integer :: sol_n, sol_m
integer :: nblocks = 10
integer :: ntry
integer, dimension( 1 ) :: isol
integer, dimension( 3 ) :: indx
logical :: F3a1_sol, F3a2_sol, F3b2_sol
logical :: both
logical :: F3l_sol
real, parameter :: Pi = 3.14159265359
real, parameter :: ec = 1.60217733e-19
real, parameter :: eps0 = 8.854187817e-12
real, parameter :: kB = 1.380658e-23
real, external :: ran2, zbrent4, F3
real :: PiRatio
real, dimension(Nbs) :: Xb_old, Yb_old, Zb_old
real, dimension(Nbs) :: Xb_new, Yb_new, Zb_new
real, dimension(Nbs) :: Xb_tmp, Yb_tmp, Zb_tmp
real :: dphi0_max = 1.0 / 180.0 * Pi
real :: phi_span = 10.0 / 180.0 * Pi
real :: phi_w = 0.5 / 180.0 * Pi
real :: phi_inc = 0.01 / 180.0 * Pi
real :: a, b, phi0_m, phi1_m, phi1_o
real :: cth, sth, cphi, sphi
real :: F3a1_num, F3a2_num, F3b1_num, F3b2_num
real :: F3_lim
real :: dF31, dF32
real :: sboltz_m, sboltz_n, J_m, J_n
real :: x1, x2
real :: Random
real, dimension(1) :: ZERO2 = 0.0
real, dimension(3) :: q1, qp1, q2, r41, w, y, z
real, dimension(4) :: Xc, Yc, Zc
real, dimension(10) :: xb1, xb2
real, dimension(40) :: PROB, boltz
real, dimension(40) :: Ubend_m, Ubend_n, Ubt_m, Ubt_n
real, dimension(40) :: Utor_m, Utor_n
real, dimension(-1:4) :: l, th, th_m, phi
real, dimension(3,3) :: T0lab, T1lab, T0, T1, T2, MM
real, dimension(3,3) :: BB_m, BB_n
real, dimension(-2:4,3) :: r, u
real, dimension(-2:4,3,40) :: r_sol
real, dimension(-1:4,40) :: phi_sol
! The molecule must be unfolded from periodic boundary conditions.
Xb_new = Xb
Yb_new = Yb
Zb_new = Zb
call Outfold( Nbs, 0, BoxSize, Xb_new, Yb_new, Zb_new, &
ZERO2, ZERO2, ZERO2 )
Xb_old = Xb_new
Yb_old = Yb_new
Zb_old = Zb_new
success = .false.
dUintra = 0.0
ntry = 0
!do while ( .NOT. success )
ntry = ntry + 1
! Xb_new = Xb_old
! Yb_new = Yb_old
! Zb_new = Zb_old
do i = -2, 4
r(i,1) = Xb_new(mapf(i))
r(i,2) = Yb_new(mapf(i))
r(i,3) = Zb_new(mapf(i))
end do
do i = -1, 4
u(i,:) = r(i,:) - r(i-1,:)
l(i) = sqrt(dot_product(u(i,:),u(i,:)))
u(i,:) = u(i,:) / l(i)
end do
! Form matrix for the old Jacobian.
BB_m = 0.0
! BB_m(1:3,i) = ui x ( r4-ri )
do i = 1, 3
y(:) = r(4,:) - r(i,:)
BB_m(1,i) = u(i,2) * y(3) - u(i,3) * y(2)
BB_m(2,i) = u(i,3) * y(1) - u(i,1) * y(3)
BB_m(3,i) = u(i,1) * y(2) - u(i,2) * y(1)
end do
! Calculate the old bond angles.
do i = -1, 3
th_m(i) = dacos(dot_product(u(i,:),u(i+1,:)))
if( i >= 0 .AND. i <= 3 ) then
j = nint( th_m(i) * 180.0 / pi )
th_h(j) = th_h(j) + 1
end if
end do
! find phi0_m, phi1_m
! cos(phi(i)) = ( (u(i-1) x u(i)) / sin(th(i-1)) ) dot ( (u(i) x u(i+1)) / sin(th(i)) )
!do i = 0, 1
do i = 0, 3
sth = dsin(th_m(i-1))
z(1) = ( u(i-1,2) * u(i,3) - u(i-1,3) * u(i,2) ) / sth
z(2) = ( u(i-1,3) * u(i,1) - u(i-1,1) * u(i,3) ) / sth
z(3) = ( u(i-1,1) * u(i,2) - u(i-1,2) * u(i,1) ) / sth
sth = dsin(th_m(i))
y(1) = ( u(i,2) * u(i+1,3) - u(i,3) * u(i+1,2) ) / sth
y(2) = ( u(i,3) * u(i+1,1) - u(i,1) * u(i+1,3) ) / sth
y(3) = ( u(i,1) * u(i+1,2) - u(i,2) * u(i+1,1) ) / sth
a = -dot_product( z, y ) ! a = cos(phi)
w(1) = z(2) * y(3) - z(3) * y(2)
w(2) = z(3) * y(1) - z(1) * y(3)
w(3) = z(1) * y(2) - z(2) * y(1)
b = -dot_product( w(:), u(i,:) ) ! b = sin(phi)
if( a > 1.0 ) a = 1.0
if( a < -1.0 ) a = -1.0
if( a > 0 .AND. b > 0 ) then
phi(i) = dacos(a)
else if( a > 0 .AND. b < 0 ) then
phi(i) = -dacos(a)
else if( a < 0 .AND. b > 0 ) then
phi(i) = dacos(a)
else if( a < 0 .AND. b < 0 ) then
phi(i) = -dacos(a)
end if
j = nint( phi(i) * 180.0 / pi )
phi_h(j) = phi_h(j) + 1
end do
phi0_m = phi(0)
phi1_m = phi(1)
phi1_o = phi(1)
phi = 0.0
! Generate new bond angles.
th = th_m
do j = 0, 3
ibd = 0
do i = 1, Nbend
if( IBEND(1,i) == mapf(j-1) .AND. IBEND(2,i) == mapf(j) .AND. &
IBEND(3,i) == mapf(j+1) ) then
ibd = i
exit
else if( IBEND(3,i) == mapf(j-1) .AND. IBEND(2,i) == mapf(j) .AND. &
IBEND(1,i) == mapf(j+1) ) then
ibd = i
exit
end if
end do
i = -1
! i = 0
do while ( i == -1 .AND. ibd /= 0 )
a = Pi * ran2(Seed)
b = 0.5 * RBEND(1,ibd) * ( a - RBEND(2,ibd) ) ** 2.0
if( dsin(a) * exp( -beta * b ) > Ran2(Seed) ) then
i = 0
th(j) = Pi - a
end if
end do
end do
! Form T0lab = ( x0, y0, z0 )
! x0 = u0
! z0 = u0 x u-1 / sin(th(-1))
! y0 = z0 x x0 = z0 x u0
sth = dsin(th(-1))
z(1) = ( u(0,2) * u(-1,3) - u(0,3) * u(-1,2) ) / sth
z(2) = ( u(0,3) * u(-1,1) - u(0,1) * u(-1,3) ) / sth
z(3) = ( u(0,1) * u(-1,2) - u(0,2) * u(-1,1) ) / sth
y(1) = ( z(2) * u(0,3) - z(3) * u(0,2) )
y(2) = ( z(3) * u(0,1) - z(1) * u(0,3) )
y(3) = ( z(1) * u(0,2) - z(2) * u(0,1) )
T0lab(:,1) = u(0,:)
T0lab(1,2) = y(1)
T0lab(2,2) = y(2)
T0lab(3,2) = y(3)
T0lab(1,3) = z(1)
T0lab(2,3) = z(2)
T0lab(3,3) = z(3)
! select a phi0
phi(0) = ran2(Seed) * dphi0_max
if( ran2(Seed) > 0.5 ) phi(0) = -phi(0)
phi(0) = phi0_m + phi(0)
!phi(0) = phi0_m
! form T0, a 3x3 matrix.
! Ti = cos(th) sin(th) 0
! sin(th)cos(ph) -cos(th)cos(ph) sin(ph)
! sin(th)sin(ph) -cos(th)sin(ph) -cos(ph)
cth = dcos(th(0))
sth = dsin(th(0))
cphi = dcos(phi(0))
sphi = dsin(phi(0))
T0(1,1) = cth
T0(1,2) = sth
T0(1,3) = 0.0
T0(2,1) = sth * cphi
T0(2,2) = -cth * cphi
T0(2,3) = sphi
T0(3,1) = sth * sphi
T0(3,2) = -cth * sphi
T0(3,3) = -cphi
! Find r1
MM = matmul(T0lab, T0)
y = 0.0
y(1) = l(1)
z(:) = r(0,:) + matmul(MM(:,:), y(:))
r_sol(1,1,:) = z(1)
r_sol(1,2,:) = z(2)
r_sol(1,3,:) = z(3)
! Find the new u1
u(1,:) = z(:) - r(0,:)
l(1) = sqrt(dot_product(u(1,:),u(1,:)))
u(1,:) = u(1,:) / l(1)
! form T1lab
! x0 = u1
! z0 = u1 x u0 / sin(th0) ! says sin(th1) in paper !!
! y0 = z0 x x0 = z0 x u1
!sth = sin(th(0))
z(1) = ( u(1,2) * u(0,3) - u(1,3) * u(0,2) ) / sth
z(2) = ( u(1,3) * u(0,1) - u(1,1) * u(0,3) ) / sth
z(3) = ( u(1,1) * u(0,2) - u(1,2) * u(0,1) ) / sth
y(1) = ( z(2) * u(1,3) - z(3) * u(1,2) )
y(2) = ( z(3) * u(1,1) - z(1) * u(1,3) )
y(3) = ( z(1) * u(1,2) - z(2) * u(1,1) )
T1lab(:,1) = u(1,:)
T1lab(1,2) = y(1)
T1lab(2,2) = y(2)
T1lab(3,2) = y(3)
T1lab(1,3) = z(1)
T1lab(2,3) = z(2)
T1lab(3,3) = z(3)
r41(:) = matmul( transpose(T1lab(:,:)), r(4,:)-r(0,:) )
!u61(:) = matmul( transpose(T1lab(:,:)), u(6,:) )
cth = dcos(th(2))
sth = dsin(th(2))
q1(1) = l(3) + l(2)*cth
q1(2) = l(2)*sth
q1(3) = 0.0
qp1(1) = l(2) + l(3)*cth
qp1(2) = l(3)*sth
qp1(3) = 0.0
!cth = dcos(th(4))
!sth = dsin(th(4))
q2(1) = l(4)
q2(2) = 0.0
q2(3) = 0.0
!l1 = 0.0
!l1(1) = l(1)
! Search for solutions to F3.
! Start by looping over the 2 possible solution branches.
nsol = 0
!open(14, file='f5.dat')
do i = 1, 2
F3a1_sol = .false.
F3a2_sol = .false.
phi(1) = phi1_m - phi_span
do while ( phi(1) < phi1_m + phi_span + phi_inc )
nb = 0
gate = 0
F3a2_num = F3( i, l(1:4), th(1:3), phi(1:3), q1, q2, qp1, r41, &
gate, T1, T2, F3a2_sol )
if( F3a2_sol ) then
gate = 0
phi(1) = phi(1) + phi_inc
F3b2_num = F3( i, l(1:4), th(1:3), phi(1:3), &
q1, q2, qp1, r41, &
gate, T1, T2, F3b2_sol )
dF32 = ( F3b2_num - F3a2_num ) / phi_inc
end if
if( F3a1_sol .OR. F3a2_sol ) then
both = .false.
if( F3a1_sol .AND. F3a2_sol ) both = .true.
if( both .AND. F3a1_num * F3a2_num < 0.0 ) then
! Routine to make sure there is only one solution
! in the interval.
x1 = phi(1) - phi_inc - phi_w
x2 = phi(1) - phi_inc
nb = 10
gate = 0
a = phi(1)
call zbrak(F3,x1,x2,nblocks,xb1,xb2,nb, &
i, l(1:4), th(1:3), phi(1:3), &
q1, q2, qp1, r41, gate)
phi(1) = a
else if( both .AND. dF31 * dF32 < 0 .AND. f3a1_num * dF31 < 0 ) then
! Use a smaller range to test for solution.
x1 = phi(1) - phi_inc - phi_w
x2 = phi(1) - phi_inc
if( F3b2_sol ) x2 = x2 + phi_inc
nb = 10
gate = 0
a = phi(1)
call zbrak(F3,x1,x2,nblocks*2,xb1,xb2,nb, &
i, l(1:4), th(1:3), phi(1:3), &
q1, q2, qp1, r41, gate)
phi(1) = a
else if( .NOT. both .AND. F3a2_sol ) then
x1 = phi(1) - phi_inc
x2 = phi(1) - phi_w - phi_inc
if( x2 < phi1_m - phi_span ) then
phi(1) = phi(1) + phi_w
if( F3a2_sol ) phi(1) = phi(1) - phi_inc
F3a1_num = F3a2_num
F3a1_sol = F3a2_sol
dF31 = dF32
cycle
end if
a = phi(1)
F3l_sol = .false.
F3_lim = F3a1_num
j = 0
do while( .NOT. F3l_sol .AND. j < 20 )
gate = -nint(F3_lim)
if( j > 0 .AND. gate == 2 ) gate = 3
b = zbrent4(F3,x1,x2,1.0e-9, &
i, l(1:4), th(1:3), phi(1:3), &
q1, q2, qp1, r41, gate, &
T1, T2 )
if( nint( b ) == 1000 ) exit
if( gate == 3 ) gate = 2
gate = -gate
phi(1) = b
F3_lim = F3( i, l(1:4), th(1:3), phi(1:3), &
q1, q2, qp1, r41, &
gate, T1, T2, F3l_sol )
x2 = phi(1)
j = j + 1
if( j > 19 ) write(*,*) ' Caught in loop '
end do
if( nint( b ) == 1000 ) then
phi(1) = a
phi(1) = phi(1) + phi_w
if( F3a2_sol ) phi(1) = phi(1) - phi_inc
F3a1_num = F3a2_num
F3a1_sol = F3a2_sol
dF31 = dF32
cycle
end if
x2 = phi(1)
nb = 10
call zbrak(F3,x1,x2,nblocks,xb1,xb2,nb, &
i, l(1:4), th(1:3), phi(1:3), &
q1, q2, qp1, r41, gate)
phi(1) = a
else if ( .NOT. both .AND. F3a1_sol ) then
x1 = phi(1) - phi_w
x2 = phi(1)
if( x2 > phi1_m + phi_span + phi_inc ) then
phi(1) = phi(1) + phi_w
if( F3a2_sol ) phi(1) = phi(1) - phi_inc
F3a1_num = F3a2_num
F3a1_sol = F3a2_sol
dF31 = dF32
cycle
end if
a = phi(1)
F3l_sol = .false.
F3_lim = F3a2_num
j = 0
do while( .NOT. F3l_sol .AND. j < 20 )
gate = -nint(F3_lim)
if( j > 0 .AND. gate == 2 ) gate = 3
b = zbrent4(F3,x1,x2,1.0e-9, &
i, l(1:4), th(1:3), phi(1:3), &
q1, q2, qp1, r41, &
gate, T1, T2 )
if( nint( b ) == 1000 ) exit
if( gate == 3 ) gate = 2
gate = -gate
phi(1) = b
F3_lim = F3( i, l(1:4), th(1:3), phi(1:3), &
q1, q2, qp1, r41, &
gate, T1, T2, F3l_sol )
x2 = phi(1)
j = j + 1
if( j > 19 ) write(*,*) ' Caught in loop '
end do
if( nint( b ) == 1000 ) then
phi(1) = a
phi(1) = phi(1) + phi_w
if( F3a2_sol ) phi(1) = phi(1) - phi_inc
F3a1_num = F3a2_num
F3a1_sol = F3a2_sol
dF31 = dF32
cycle
end if
x2 = phi(1)
nb = 10
call zbrak(F3,x1,x2,nblocks,xb1,xb2,nb, &
i, l(1:4), th(1:3), phi(1:3), &
q1, q2, qp1, r41, gate)
phi(1) = a
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -