3 C Build the virtual polypeptide chain. Side-chain centroids are moveable.
6 implicit real*8 (a-h,o-z)
12 include 'COMMON.IOUNITS'
13 include 'COMMON.NAMES'
14 include 'COMMON.INTERACT'
15 double precision e1(3),e2(3),e3(3)
16 logical lprn,perbox,fail
17 C Set lprn=.true. for debugging
21 print *, 'enter chainbuild'
29 print *,'before refsys'
30 call refsys(2,3,4,e1,e2,e3,fail)
31 print *,'after refsys'
40 print *,'dc',dc(1,0),dc(2,0),dc(3,0)
58 veclen=veclen+(c(i,2)-c(i,1))**2
82 call locate_side_chain(2)
86 if (theti.ne.theti) theti=100.0
88 if (phii.ne.phii) phii=180.0
97 * Define the matrices of the rotation about the virtual-bond valence angles
98 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
99 * program), R(i,j,k), and, the cumulative matrices of rotation RT
121 rt(2,1,i-2)=sint*cosphi
122 rt(2,2,i-2)=-cost*cosphi
124 rt(3,1,i-2)=-sint*sinphi
125 rt(3,2,i-2)=cost*sinphi
127 call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
129 dc_norm(j,i-1)=prod(j,1,i-1)
130 dc(j,i-1)=vbld(i)*prod(j,1,i-1)
132 call locate_side_chain(i-1)
136 C Define the origin and orientation of the coordinate system and locate the
137 C first three CA's and SC(2).
141 * Build the alpha-carbon chain.
144 call locate_next_res(i)
147 C First and last SC must coincide with the corresponding CA.
151 dc_norm(j,nres+1)=0.0D0
152 dc(j,nres+nres)=0.0D0
153 dc_norm(j,nres+nres)=0.0D0
155 c(j,nres+nres)=c(j,nres)
158 * Temporary diagnosis
163 write (iout,'(/a)') 'Recalculated internal coordinates'
166 c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
169 if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
170 be1=rad2deg*beta(nres+i,i,maxres2,i+1)
172 if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
173 write (iout,1212) restyp(itype(i)),i,dist(i-1,i),
174 & alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,maxres2),be1
176 1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
183 c-------------------------------------------------------------------------
184 subroutine orig_frame
186 C Define the origin and orientation of the coordinate system and locate
187 C the first three atoms.
189 implicit real*8 (a-h,o-z)
191 include 'COMMON.CHAIN'
192 include 'COMMON.LOCAL'
246 dc_norm(j,2)=prod(j,1,2)
247 dc(j,2)=vbld(3)*prod(j,1,2)
248 c(j,3)=c(j,2)+dc(j,2)
250 call locate_side_chain(2)
253 c-----------------------------------------------------------------------------
254 subroutine locate_next_res(i)
256 C Locate CA(i) and SC(i-1)
258 implicit real*8 (a-h,o-z)
260 include 'COMMON.CHAIN'
261 include 'COMMON.LOCAL'
264 include 'COMMON.IOUNITS'
265 include 'COMMON.NAMES'
266 include 'COMMON.INTERACT'
268 C Define the rotation matrices corresponding to CA(i)
272 if (theti.ne.theti) theti=100.0
274 if (phii.ne.phii) phii=180.0
283 * Define the matrices of the rotation about the virtual-bond valence angles
284 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
285 * program), R(i,j,k), and, the cumulative matrices of rotation RT
307 rt(2,1,i-2)=sint*cosphi
308 rt(2,2,i-2)=-cost*cosphi
310 rt(3,1,i-2)=-sint*sinphi
311 rt(3,2,i-2)=cost*sinphi
313 call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
315 dc_norm(j,i-1)=prod(j,1,i-1)
316 dc(j,i-1)=vbld(i)*prod(j,1,i-1)
317 c(j,i)=c(j,i-1)+dc(j,i-1)
319 cd print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
321 C Now calculate the coordinates of SC(i-1)
323 call locate_side_chain(i-1)
326 c-----------------------------------------------------------------------------
327 subroutine locate_side_chain(i)
329 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
331 implicit real*8 (a-h,o-z)
333 include 'COMMON.CHAIN'
334 include 'COMMON.LOCAL'
337 include 'COMMON.IOUNITS'
338 include 'COMMON.NAMES'
339 include 'COMMON.INTERACT'
343 c dsci_inv=dsc_inv(itype(i))
345 dsci_inv=vbld_inv(i+nres)
349 if (alphi.ne.alphi) alphi=100.0
350 if (omegi.ne.omegi) omegi=-100.0
360 yp= dsci*sinalphi*cosomegi
361 zp=-dsci*sinalphi*sinomegi
362 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
363 * X-axis aligned with the vector DC(*,i)
364 theta2=pi-0.5D0*theta(i+1)
367 xx(1)= xp*cost2+yp*sint2
368 xx(2)=-xp*sint2+yp*cost2
370 cd print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
371 cd & xp,yp,zp,(xx(k),k=1,3)
375 * Bring the SC vectors to the common coordinate system.
377 xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
378 xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
385 rj=rj+prod(j,k,i-1)*xx(k)
388 dc_norm(j,nres+i)=rj*dsci_inv
389 c(j,nres+i)=c(j,i)+rj
393 c------------------------------------------
399 include 'COMMON.CONTROL'
403 include 'COMMON.LANGEVIN'
405 include 'COMMON.LANGEVIN.lang0'
407 include 'COMMON.CHAIN'
408 include 'COMMON.DERIV'
410 include 'COMMON.LOCAL'
411 include 'COMMON.INTERACT'
412 include 'COMMON.IOUNITS'
413 include 'COMMON.NAMES'
414 include 'COMMON.TIME1'
415 include 'COMMON.REMD'
416 include 'COMMON.SETUP'
417 include 'COMMON.MUCA'
418 include 'COMMON.HAIRPIN'
422 C write(*,*) 'initial', i,j,c(j,i)
425 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
427 if (allareout.eq.1) then
428 ireturnval=int(c(j,i)/boxxsize)
429 if (c(j,i).le.0) ireturnval=ireturnval-1
430 do k=chain_beg,chain_end
431 c(j,k)=c(j,k)-ireturnval*boxxsize
432 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
438 if (int(c(j,i)/boxxsize).eq.0) allareout=0
441 if (allareout.eq.1) then
442 ireturnval=int(c(j,i)/boxxsize)
443 if (c(j,i).le.0) ireturnval=ireturnval-1
445 c(j,k)=c(j,k)-ireturnval*boxxsize
446 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
451 C write(*,*) 'befor no jump', i,j,c(j,i)
455 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
456 difference=abs(c(j,i-1)-c(j,i))
457 C print *,'diff', difference
458 if (difference.gt.boxxsize/2.0) then
459 if (c(j,i-1).gt.c(j,i)) then
468 c(j,i)=c(j,i)+nojumpval*boxxsize
469 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
473 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
474 difference=abs(c(j,i-1)-c(j,i))
475 if (difference.gt.boxxsize/2.0) then
476 if (c(j,i-1).gt.c(j,i)) then
485 c(j,i)=c(j,i)+nojumpval*boxxsize
486 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
490 C write(*,*) 'after no jump', i,j,c(j,i)
497 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
499 if (allareout.eq.1) then
500 ireturnval=int(c(j,i)/boxysize)
501 if (c(j,i).le.0) ireturnval=ireturnval-1
502 do k=chain_beg,chain_end
503 c(j,k)=c(j,k)-ireturnval*boxysize
504 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
510 if (int(c(j,i)/boxysize).eq.0) allareout=0
513 if (allareout.eq.1) then
514 ireturnval=int(c(j,i)/boxysize)
515 if (c(j,i).le.0) ireturnval=ireturnval-1
517 c(j,k)=c(j,k)-ireturnval*boxysize
518 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
523 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
524 difference=abs(c(j,i-1)-c(j,i))
525 if (difference.gt.boxysize/2.0) then
526 if (c(j,i-1).gt.c(j,i)) then
535 c(j,i)=c(j,i)+nojumpval*boxysize
536 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
540 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
541 difference=abs(c(j,i-1)-c(j,i))
542 if (difference.gt.boxysize/2.0) then
543 if (c(j,i-1).gt.c(j,i)) then
552 c(j,i)=c(j,i)+nojumpval*boxysize
553 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
558 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
560 if (allareout.eq.1) then
561 ireturnval=int(c(j,i)/boxysize)
562 if (c(j,i).le.0) ireturnval=ireturnval-1
563 do k=chain_beg,chain_end
564 c(j,k)=c(j,k)-ireturnval*boxzsize
565 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
571 if (int(c(j,i)/boxzsize).eq.0) allareout=0
574 if (allareout.eq.1) then
575 ireturnval=int(c(j,i)/boxzsize)
576 if (c(j,i).le.0) ireturnval=ireturnval-1
578 c(j,k)=c(j,k)-ireturnval*boxzsize
579 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
584 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
585 difference=abs(c(j,i-1)-c(j,i))
586 if (difference.gt.(boxzsize/2.0)) then
587 if (c(j,i-1).gt.c(j,i)) then
596 c(j,i)=c(j,i)+nojumpval*boxzsize
597 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
601 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
602 difference=abs(c(j,i-1)-c(j,i))
603 if (difference.gt.boxzsize/2.0) then
604 if (c(j,i-1).gt.c(j,i)) then
613 c(j,i)=c(j,i)+nojumpval*boxzsize
614 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize