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