Fixed the addition of dummy chain ends to structures read from UNRES PDB files
[unres.git] / source / unres / src_MD / bond_move.f
1       subroutine bond_move(nbond,nstart,psi,lprint,error)
2 C Move NBOND fragment starting from the CA(nstart) by angle PSI.
3       implicit real*8 (a-h,o-z)
4       include 'DIMENSIONS'
5       integer nbond,nstart
6       double precision psi
7       logical fail,error,lprint
8       include 'COMMON.GEO'
9       include 'COMMON.CHAIN'
10       include 'COMMON.VAR'
11       include 'COMMON.IOUNITS'
12       include 'COMMON.MCM'
13       dimension x(3),e(3,3),e1(3),e2(3),e3(3),rot(3,3),trans(3,3)
14       error=.false.
15       nend=nstart+nbond
16       if (print_mc.gt.2) then
17       write (iout,*) 'nstart=',nstart,' nend=',nend,' nbond=',nbond
18       write (iout,*) 'psi=',psi
19       write (iout,'(a)') 'Original coordinates of the fragment'
20       do i=nstart,nend
21         write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
22       enddo
23       endif
24       if (nstart.lt.1 .or. nend .gt.nres .or. nbond.lt.2 .or. 
25      & nbond.ge.nres-1) then
26         write (iout,'(a)') 'Bad data in BOND_MOVE.'
27         error=.true.
28         return
29       endif
30 C Generate the reference system.
31       i2=nend
32       i3=nstart
33       i4=nstart+1
34       call refsys(i2,i3,i4,e1,e2,e3,error) 
35 C Return, if couldn't define the reference system.
36       if (error) return
37 C Compute the transformation matrix.
38       cospsi=dcos(psi)
39       sinpsi=dsin(psi)
40       rot(1,1)=1.0D0
41       rot(1,2)=0.0D0
42       rot(1,3)=0.0D0
43       rot(2,1)=0.0D0
44       rot(2,2)=cospsi
45       rot(2,3)=-sinpsi
46       rot(3,1)=0.0D0
47       rot(3,2)=sinpsi
48       rot(3,3)=cospsi
49       do i=1,3
50         e(1,i)=e1(i)
51         e(2,i)=e2(i)
52         e(3,i)=e3(i)
53       enddo
54
55       if (print_mc.gt.2) then
56       write (iout,'(a)') 'Reference system and matrix r:'
57       do i=1,3
58         write(iout,'(i5,2(3f10.5,5x))')i,(e(i,j),j=1,3),(rot(i,j),j=1,3)
59       enddo
60       endif
61
62       call matmult(rot,e,trans)
63       do i=1,3
64         do j=1,3
65           e(i,1)=e1(i)
66           e(i,2)=e2(i)
67           e(i,3)=e3(i)
68         enddo
69       enddo
70       call matmult(e,trans,trans)
71
72       if (lprint) then
73       write (iout,'(a)') 'The trans matrix:'
74       do i=1,3
75         write (iout,'(i5,3f10.5)') i,(trans(i,j),j=1,3)
76       enddo
77       endif
78
79       do i=nstart,nend
80         do j=1,3
81           rij=c(j,nstart)
82           do k=1,3
83             rij=rij+trans(j,k)*(c(k,i)-c(k,nstart))
84           enddo
85           x(j)=rij
86         enddo
87         do j=1,3
88           c(j,i)=x(j)
89         enddo
90       enddo
91
92       if (lprint) then
93       write (iout,'(a)') 'Rotated coordinates of the fragment'
94       do i=nstart,nend
95         write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
96       enddo
97       endif
98
99 c     call int_from_cart(.false.,lprint)
100       if (nstart.gt.1) then
101         theta(nstart+1)=alpha(nstart-1,nstart,nstart+1)
102         phi(nstart+2)=beta(nstart-1,nstart,nstart+1,nstart+2)
103         if (nstart.gt.2) phi(nstart+1)=
104      &                beta(nstart-2,nstart-1,nstart,nstart+1)
105       endif
106       if (nend.lt.nres) then
107         theta(nend+1)=alpha(nend-1,nend,nend+1)
108         phi(nend+1)=beta(nend-2,nend-1,nend,nend+1)
109         if (nend.lt.nres-1) phi(nend+2)=
110      &                beta(nend-1,nend,nend+1,nend+2)
111       endif
112       if (print_mc.gt.2) then
113       write (iout,'(/a,i3,a,i3,a/)') 
114      & 'Moved internal coordinates of the ',nstart,'-',nend,
115      & ' fragment:'
116       do i=nstart+1,nstart+2
117         write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
118       enddo
119       do i=nend+1,nend+2
120         write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
121       enddo
122       endif
123       return
124       end