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