📄 conrot.f90
字号:
phi_sol(3,nsol) = phi(3)
! determine r_sol's
MM = matmul( matmul( T0lab, T0 ), T1 )
y = 0.0
y(1) = l(2)
r_sol(2,:,nsol) = r_sol(1,:,nsol) + matmul( MM(:,:), y(:) )
MM = matmul( MM, T2 )
y = 0.0
y(1) = l(3)
r_sol(3,:,nsol) = r_sol(2,:,nsol) + matmul( MM(:,:), y(:) )
end do
end if
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
end do
end do
isol = 0
PROB = 0.0
PROB(1:nsol) = abs( phi_sol(1,1:nsol) - phi1_o )
if( minval(PROB(1:nsol)) < 5.0e-7 ) isol = minloc( PROB(1:nsol) )
sol_m = isol(1)
if( sol_m == 0 ) then
! write(*,*) ' Original solution not found in Conrot ! '
! write(*,"(1x,A, F10.7)") ' phi1_old = ', phi1_o
! write(*,"(1x,A, 30(2x,F10.7))") ' phi1_new = ', phi_sol(1,1:nsol)
! open(23, file='nosol.dat')
! read(23,*) j
! close(23)
! open(23, file='nosol.dat')
! write(23,*) j + 1
! close(23)
! cycle
return
end if
! Find the old torsion energy with the exact values of the old coordinates.
! If the values from resolving the solution are used, small differences result,
! which accumulate over time.
Ubend_m = 0.0
Utor_m = 0.0
do j = 1, nsol
if( j == sol_m ) then
do i = 1, 3
Xb_new(mapf(i)) = Xb_old(mapf(i))
Yb_new(mapf(i)) = Yb_old(mapf(i))
Zb_new(mapf(i)) = Zb_old(mapf(i))
end do
else
do i = 1, 3
Xb_new(mapf(i)) = r_sol(i,1,j)
Yb_new(mapf(i)) = r_sol(i,2,j)
Zb_new(mapf(i)) = r_sol(i,3,j)
end do
end if
do i = 1, Nbend
ibd = 0
if( IBEND(1,i) == mapf(1) .OR. IBEND(1,i) == mapf(2) .OR. &
IBEND(1,i) == mapf(3) ) then
ibd = i
else if( IBEND(2,i) == mapf(1) .OR. IBEND(2,i) == mapf(2) .OR. &
IBEND(2,i) == mapf(3) ) then
ibd = i
else if( IBEND(3,i) == mapf(1) .OR. IBEND(3,i) == mapf(2) .OR. &
IBEND(3,i) == mapf(3) ) then
ibd = i
end if
if( ibd == i ) then
Xc(1:3) = Xb_new( IBEND(1:3,i) )
Yc(1:3) = Yb_new( IBEND(1:3,i) )
Zc(1:3) = Zb_new( IBEND(1:3,i) )
call IntraBend( Xc, Yc, Zc, RBEND(1,i), RBEND(2,i), b )
Ubend_m(j) = Ubend_m(j) + b
if( IBEND(1,i) == mapf(3) .AND. IBEND(2,i) == mapf(-2) .AND. &
IBEND(3,i) == mapf(-1) ) then
Ubt_m(j) = Ubt_m(j) + b
else if( IBEND(3,i) == mapf(3) .AND. IBEND(2,i) == mapf(-2) .AND. &
IBEND(1,i) == mapf(-1) ) then
Ubt_m(j) = Ubt_m(j) + b
end if
end if
end do
do i = 1, Ntor
ibd = 0
if( ITOR(1,i) == mapf(1) .OR. ITOR(1,i) == mapf(2) .OR. &
ITOR(1,i) == mapf(3) ) then
ibd = i
else if( ITOR(2,i) == mapf(1) .OR. ITOR(2,i) == mapf(2) .OR. &
ITOR(2,i) == mapf(3) ) then
ibd = i
else if( ITOR(3,i) == mapf(1) .OR. ITOR(3,i) == mapf(2) .OR. &
ITOR(3,i) == mapf(3) ) then
ibd = i
else if( ITOR(4,i) == mapf(1) .OR. ITOR(4,i) == mapf(2) .OR. &
ITOR(4,i) == mapf(3) ) then
ibd = i
end if
if( ibd == i ) then
Xc(1:4) = Xb_new( ITOR(1:4,i) )
Yc(1:4) = Yb_new( ITOR(1:4,i) )
Zc(1:4) = Zb_new( ITOR(1:4,i) )
call IntraTorsion( Xc, Yc, Zc, 2, RTOR(1:4,i), b )
Utor_m(j) = Utor_m(j) + b
end if
end do
end do
boltz(:) = exp( -beta * ( Ubt_m(:) + Utor_m(:) ) )
! boltz(:) = exp( -beta * Utor_m(:) )
sboltz_m = sum(boltz(1:nsol))
!Utor_m(sol_m) = a
!return
! Calculate the Jacobians.
call ludcmp( BB_m, 3, 3, indx, J_m )
do i = 1, 3
J_m = J_m * BB_m(i,i)
end do
J_m = abs( 1.0 / J_m )
call ludcmp( BB_n, 3, 3, indx, J_n )
do i = 1, 3
J_n = J_n * BB_n(i,i)
end do
J_n = abs( 1.0 / J_n )
PiRatio = log( sboltz_n / sboltz_m * J_n / J_m )
if( Piratio > log(ran2(Seed)) ) then
! accept move
success = .true.
dUintra = Ubend_n(sol_n) - Ubend_m(sol_m) + Utor_n(sol_n) - Utor_m(sol_m)
! write(*,*) ' Success', dUintra
Xb = Xb_tmp
Yb = Yb_tmp
Zb = Zb_tmp
end if
! Xb_new = Xb_old
! Yb_new = Yb_old
! Zb_new = Zb_old
!end do
return
end subroutine conrot4
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! function to calculate F3
real function F3( branch, l, th, phi, q1, q2, qp1, r41, &
gate, T1, T2, sol )
implicit none
integer :: branch, gate
real, dimension(4) :: l
real, dimension(3) :: th, phi
real, dimension(3) :: q1, q2, qp1, r41, y
real, dimension(3,3) :: T1, T2
logical :: sol
! local stuff
real :: pi, c1, h, a, b, cth, sth, cphi, sphi
real, dimension(3) :: t, q3
! intial settings
sol = .false.
pi = dacos(-1.0)
cth = dcos(th(1))
sth = dsin(th(1))
cphi = dcos(phi(1))
sphi = dsin(phi(1))
T1(1,1) = cth
T1(1,2) = sth
T1(1,3) = 0.0
T1(2,1) = sth * cphi
T1(2,2) = -cth * cphi
T1(2,3) = sphi
T1(3,1) = sth * sphi
T1(3,2) = -cth * sphi
T1(3,3) = -cphi
y = 0.0
y(1) = l(1)
t = matmul( transpose(T1), r41 - y )
! determine c1 to find if a solution exists for phi(2)
c1 = ( dot_product(q1,q1) - dot_product(q2,q2) + dot_product(t,t) - 2*t(1)*qp1(1) ) / &
( 2*qp1(2)*sqrt( dot_product(t,t) - t(1)*t(1) ) )
select case ( gate )
case ( 0, -2, +2 )
if( abs(c1) > 1.0 ) then
F3 = -1.0
return
end if
case ( -1, +3 )
c1 = sign(1.0,c1)
case ( +1 )
F3 = abs(c1) - 1.0
sol = .true.
return
end select
sol = .true.
h = 0.0
if( t(3) < 0.0 ) h = pi
if( branch == 1 ) then
phi(2) = dasin(c1) - datan(t(2)/t(3)) - h
else
phi(2) = pi - dasin(c1) - datan(t(2)/t(3)) - h
end if
if( phi(2) > pi ) phi(2) = phi(2) - 2*pi
if( phi(2) < -pi ) phi(2) = phi(2) + 2*pi
cth = dcos(th(2))
sth = dsin(th(2))
cphi = dcos(phi(2))
sphi = dsin(phi(2))
T2(1,1) = cth
T2(1,2) = sth
T2(1,3) = 0.0
T2(2,1) = sth * cphi
T2(2,2) = -cth * cphi
T2(2,3) = sphi
T2(3,1) = sth * sphi
T2(3,2) = -cth * sphi
T2(3,3) = -cphi
q3 = matmul( transpose(T2), t ) - q1
cth = dcos(th(3))
sth = dsin(th(3))
a = q2(3)
b = sth*q2(1) - cth*q2(2)
sphi = ( a*q3(2) + b*q3(3) ) / ( a*a + b*b )
cphi = ( q3(2) -a*sphi ) / b
if( cphi > 1.0 ) cphi = 1.0
if( cphi < -1.0 ) cphi = -1.0
if( cphi > 0 .AND. sphi > 0 ) then
phi(3) = dacos(cphi)
else if( cphi > 0 .AND. sphi < 0 ) then
phi(3) = -dacos(cphi)
else if( cphi < 0 .AND. sphi > 0 ) then
phi(3) = dacos(cphi)
else if( cphi < 0 .AND. sphi < 0 ) then
phi(3) = -dacos(cphi)
end if
y = 0.0
y(1) = l(2)
y = t - y
y = matmul( transpose(T2), y )
y(1) = y(1) - l(3)
F3 = y(1) - cth*l(4)
end function F3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE zbrak(fx,x1,x2,n,xb1,xb2,nb, &
branch, l, th, phi, &
q1, q2, qp1, r41, gate)
INTEGER n,nb
REAL x1,x2,xb1(nb),xb2(nb),fx
EXTERNAL fx
INTEGER i,nbb
REAL dx,fc,fp,x
! additions
integer :: branch, gate, gate_tmp
real, dimension(4) :: l
real, dimension(3) :: th, phi
real, dimension(3) :: q1, q2, qp1, r41
real, dimension(3,3) :: T1, T2
logical :: sol_fp, sol_fc
nbb=0
x=x1
dx=(x2-x1)/n
gate_tmp = 0
phi(1) = x
fp=fx(branch, l, th, phi, q1, q2, qp1, r41, &
gate_tmp, T1, T2, sol_fp)
do 11 i=1,n
x=x+dx
if( i == n .AND. gate /= 0 ) then
x = x2
gate_tmp = gate
end if
phi(1) = x
fc=fx(branch, l, th, phi, q1, q2, qp1, r41, &
gate_tmp, T1, T2, sol_fc)
! if( .NOT. sol_fc ) then
!
! write(*,*) ' No solution in zbrak, i= ', i
! exit
!
! end if
if( i == n ) gate_tmp = 0
if( fc*fp < 0.0 .AND. sol_fp .AND. sol_fc ) then
if(i==n) gate_tmp = gate
nbb=nbb+1
xb1(nbb)=x-dx
xb2(nbb)=x
if(nbb.eq.nb)goto 1
endif
fp=fc
sol_fp = sol_fc
11 continue
1 continue
nb=nbb
gate = gate_tmp
return
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real function zbrent4(func,x1,x2,tol, &
branch, l, th, phi, q1, q2, qp1, &
r41, gate, T1, T2 )
INTEGER ITMAX
REAL tol,x1,x2,func,EPS
EXTERNAL func
PARAMETER (ITMAX=100,EPS=3.e-8)
INTEGER iter
REAL a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm
! additions
integer :: branch, gate, gate_tmp
real, dimension(4) :: l
real, dimension(3) :: th, phi
real, dimension(3) :: q1, q2, qp1, r41
real, dimension(3,3) :: T1, T2
logical :: sol
a=x1
b=x2
gate_tmp = 0
if( gate > 0 ) gate_tmp = gate
if( gate_tmp == 3 ) then
gate_tmp = 2
end if
! fa=func(a) ! gate = 0 for x1
phi(1) = a
fa=func(branch, l, th, phi, q1, q2, qp1, r41, &
gate_tmp, T1, T2, sol)
! fb=func(b) ! if a "wall" exists, it exists at x2 ... gate can be -1, -2
phi(1) = b
fb=func(branch, l, th, phi, q1, q2, qp1, r41, &
gate, T1, T2, sol)
! if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.))pause
! *'root must be bracketed for zbrent'
c=b
fc=fb
do 12 iter=1,ITMAX
if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then
c=a
fc=fa
d=b-a
e=d
endif
if(abs(fc).lt.abs(fb)) then
a=b
b=c
c=a
fa=fb
fb=fc
fc=fa
endif
tol1=2.*EPS*abs(b)+0.5*tol
xm=.5*(c-b)
if(abs(xm).le.tol1 .or. fb.eq.0.)then
zbrent4=b
! if( abs(fb) > 1.0e-6 ) write(*,*) ' fb= ', fb
return
endif
if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
s=fb/fa
if(a.eq.c) then
p=2.*xm*s
q=1.-s
else
q=fa/fc
r=fb/fc
p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
q=(q-1.)*(r-1.)*(s-1.)
endif
if(p.gt.0.) q=-q
p=abs(p)
if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
e=d
d=p/q
else
d=xm
e=d
endif
else
d=xm
e=d
endif
a=b
fa=fb
if(abs(d) .gt. tol1) then
b=b+d
else
b=b+sign(tol1,xm)
endif
! fb=func(b)
phi(1) = b
fb=func(branch, l, th, phi, q1, q2, qp1, r41, &
gate_tmp, T1, T2, sol)
if( .NOT. sol ) then
write(*,*) ' Solution not found in zbrent'
zbrent4=1000.0
return
end if
12 continue
write(*,*) ' zbrent exceeding maximum iterations'
stop
END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE ludcmp(a,n,np,indx,d)
INTEGER n,np,indx(n),NMAX
REAL d,a(np,np),TINY
PARAMETER (NMAX=500,TINY=1.0e-20)
INTEGER i,imax,j,k
REAL aamax,dum,sum,vv(NMAX)
d=1.
do 22 i=1,n
aamax=0.
do 21 j=1,n
if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
21 continue
if (aamax.eq.0.) pause 'singular matrix in ludcmp'
vv(i)=1./aamax
22 continue
do 29 j=1,n
do 24 i=1,j-1
sum=a(i,j)
do 23 k=1,i-1
sum=sum-a(i,k)*a(k,j)
23 continue
a(i,j)=sum
24 continue
aamax=0.
do 26 i=j,n
sum=a(i,j)
do 25 k=1,j-1
sum=sum-a(i,k)*a(k,j)
25 continue
a(i,j)=sum
dum=vv(i)*abs(sum)
if (dum.ge.aamax) then
imax=i
aamax=dum
endif
26 continue
if (j.ne.imax)then
do 27 k=1,n
dum=a(imax,k)
a(imax,k)=a(j,k)
a(j,k)=dum
27 continue
d=-d
vv(imax)=vv(j)
endif
indx(j)=imax
if(a(j,j).eq.0.)a(j,j)=TINY
if(j.ne.n)then
dum=1./a(j,j)
do 28 i=j+1,n
a(i,j)=a(i,j)*dum
28 continue
endif
29 continue
return
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -