Merge branch 'devel' of mmka:unres into devel
[unres.git] / source / cluster / wham / src / chainbuild.f
1       subroutine chainbuild
2
3 C Build the virtual polypeptide chain. Side-chain centroids are moveable.
4 C As of 2/17/95.
5 C
6       implicit real*8 (a-h,o-z)
7       include 'DIMENSIONS'
8       include 'COMMON.CHAIN'
9       include 'COMMON.LOCAL'
10       include 'COMMON.GEO'
11       include 'COMMON.VAR'
12       include 'COMMON.IOUNITS'
13       include 'COMMON.NAMES'
14       include 'COMMON.INTERACT'
15       logical lprn
16 C Set lprn=.true. for debugging
17       lprn = .false.
18 C
19 C Define the origin and orientation of the coordinate system and locate the
20 C first three CA's and SC(2).
21 C
22       call orig_frame
23 *
24 * Build the alpha-carbon chain.
25 *
26       do i=4,nres
27         call locate_next_res(i)
28       enddo     
29 C
30 C First and last SC must coincide with the corresponding CA.
31 C
32       do j=1,3
33         dc(j,nres+1)=0.0D0
34         dc_norm(j,nres+1)=0.0D0
35         dc(j,nres+nres)=0.0D0
36         dc_norm(j,nres+nres)=0.0D0
37         c(j,nres+1)=c(j,1)
38         c(j,nres+nres)=c(j,nres)
39       enddo
40 *
41 * Temporary diagnosis
42 *
43       if (lprn) then
44
45       call cartprint
46       write (iout,'(/a)') 'Recalculated internal coordinates'
47       do i=2,nres-1
48         do j=1,3
49           c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
50         enddo
51         be=0.0D0
52         if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
53         be1=rad2deg*beta(nres+i,i,maxres2,i+1)
54         alfai=0.0D0
55         if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
56         write (iout,1212) restyp(itype(i)),i,dist(i-1,i),
57      &  alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,maxres2),be1
58       enddo   
59  1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
60
61       endif
62
63       return
64       end
65 c-------------------------------------------------------------------------
66       subroutine orig_frame
67 C
68 C Define the origin and orientation of the coordinate system and locate 
69 C the first three atoms.
70 C
71       implicit real*8 (a-h,o-z)
72       include 'DIMENSIONS'
73       include 'COMMON.CHAIN'
74       include 'COMMON.LOCAL'
75       include 'COMMON.GEO'
76       include 'COMMON.VAR'
77       cost=dcos(theta(3))
78       sint=dsin(theta(3))
79       t(1,1,1)=-cost
80       t(1,2,1)=-sint 
81       t(1,3,1)= 0.0D0
82       t(2,1,1)=-sint
83       t(2,2,1)= cost
84       t(2,3,1)= 0.0D0
85       t(3,1,1)= 0.0D0
86       t(3,2,1)= 0.0D0
87       t(3,3,1)= 1.0D0
88       r(1,1,1)= 1.0D0
89       r(1,2,1)= 0.0D0
90       r(1,3,1)= 0.0D0
91       r(2,1,1)= 0.0D0
92       r(2,2,1)= 1.0D0
93       r(2,3,1)= 0.0D0
94       r(3,1,1)= 0.0D0
95       r(3,2,1)= 0.0D0
96       r(3,3,1)= 1.0D0
97       do i=1,3
98         do j=1,3
99           rt(i,j,1)=t(i,j,1)
100         enddo
101       enddo
102       do i=1,3
103         do j=1,3
104           prod(i,j,1)=0.0D0
105           prod(i,j,2)=t(i,j,1)
106         enddo
107         prod(i,i,1)=1.0D0
108       enddo   
109       c(1,1)=0.0D0
110       c(2,1)=0.0D0
111       c(3,1)=0.0D0
112       c(1,2)=vbl
113       c(2,2)=0.0D0
114       c(3,2)=0.0D0
115       dc(1,1)=vbl
116       dc(2,1)=0.0D0
117       dc(3,1)=0.0D0
118       dc_norm(1,1)=1.0D0
119       dc_norm(2,1)=0.0D0
120       dc_norm(3,1)=0.0D0
121       do j=1,3
122         dc_norm(j,2)=prod(j,1,2)
123         dc(j,2)=vbl*prod(j,1,2)
124         c(j,3)=c(j,2)+dc(j,2)
125       enddo
126       call locate_side_chain(2)
127       return
128       end
129 c-----------------------------------------------------------------------------
130       subroutine locate_next_res(i)
131 C
132 C Locate CA(i) and SC(i-1)
133 C
134       implicit real*8 (a-h,o-z)
135       include 'DIMENSIONS'
136       include 'COMMON.CHAIN'
137       include 'COMMON.LOCAL'
138       include 'COMMON.GEO'
139       include 'COMMON.VAR'
140       include 'COMMON.IOUNITS'
141       include 'COMMON.NAMES'
142       include 'COMMON.INTERACT'
143 C
144 C Define the rotation matrices corresponding to CA(i)
145 C
146       theti=theta(i)
147       phii=phi(i)
148       cost=dcos(theti)
149       sint=dsin(theti)
150       cosphi=dcos(phii)
151       sinphi=dsin(phii)
152 * Define the matrices of the rotation about the virtual-bond valence angles
153 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
154 * program), R(i,j,k), and, the cumulative matrices of rotation RT
155       t(1,1,i-2)=-cost
156       t(1,2,i-2)=-sint 
157       t(1,3,i-2)= 0.0D0
158       t(2,1,i-2)=-sint
159       t(2,2,i-2)= cost
160       t(2,3,i-2)= 0.0D0
161       t(3,1,i-2)= 0.0D0
162       t(3,2,i-2)= 0.0D0
163       t(3,3,i-2)= 1.0D0
164       r(1,1,i-2)= 1.0D0
165       r(1,2,i-2)= 0.0D0
166       r(1,3,i-2)= 0.0D0
167       r(2,1,i-2)= 0.0D0
168       r(2,2,i-2)=-cosphi
169       r(2,3,i-2)= sinphi
170       r(3,1,i-2)= 0.0D0
171       r(3,2,i-2)= sinphi
172       r(3,3,i-2)= cosphi
173       rt(1,1,i-2)=-cost
174       rt(1,2,i-2)=-sint
175       rt(1,3,i-2)=0.0D0
176       rt(2,1,i-2)=sint*cosphi
177       rt(2,2,i-2)=-cost*cosphi
178       rt(2,3,i-2)=sinphi
179       rt(3,1,i-2)=-sint*sinphi
180       rt(3,2,i-2)=cost*sinphi
181       rt(3,3,i-2)=cosphi
182       call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
183       do j=1,3
184         dc_norm(j,i-1)=prod(j,1,i-1)
185         dc(j,i-1)=vbl*prod(j,1,i-1)
186         c(j,i)=c(j,i-1)+dc(j,i-1)
187       enddo
188 cd    print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
189
190 C Now calculate the coordinates of SC(i-1)
191 C
192       call locate_side_chain(i-1)
193       return
194       end
195 c-----------------------------------------------------------------------------
196       subroutine locate_side_chain(i)
197
198 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
199 C
200       implicit real*8 (a-h,o-z)
201       include 'DIMENSIONS'
202       include 'COMMON.CHAIN'
203       include 'COMMON.LOCAL'
204       include 'COMMON.GEO'
205       include 'COMMON.VAR'
206       include 'COMMON.IOUNITS'
207       include 'COMMON.NAMES'
208       include 'COMMON.INTERACT'
209       dimension xx(3)
210
211       dsci=dsc(itype(i))
212       dsci_inv=dsc_inv(itype(i))
213       alphi=alph(i)
214       omegi=omeg(i)
215       cosalphi=dcos(alphi)
216       sinalphi=dsin(alphi)
217       cosomegi=dcos(omegi)
218       sinomegi=dsin(omegi) 
219       xp= dsci*cosalphi
220       yp= dsci*sinalphi*cosomegi
221       zp=-dsci*sinalphi*sinomegi
222 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
223 * X-axis aligned with the vector DC(*,i)
224       theta2=pi-0.5D0*theta(i+1)
225       cost2=dcos(theta2)
226       sint2=dsin(theta2)
227       xx(1)= xp*cost2+yp*sint2
228       xx(2)=-xp*sint2+yp*cost2
229       xx(3)= zp
230 cd    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
231 cd   &   xp,yp,zp,(xx(k),k=1,3)
232       do j=1,3
233         xloc(j,i)=xx(j)
234       enddo
235 * Bring the SC vectors to the common coordinate system.
236       xx(1)=xloc(1,i)
237       xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
238       xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
239       do j=1,3
240         xrot(j,i)=xx(j)
241       enddo
242       do j=1,3
243         rj=0.0D0
244         do k=1,3
245           rj=rj+prod(j,k,i-1)*xx(k)
246         enddo
247         dc(j,nres+i)=rj
248         dc_norm(j,nres+i)=rj*dsci_inv
249         c(j,nres+i)=c(j,i)+rj
250       enddo
251       return
252       end