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'
25 print *,'before refsys'
26 call refsys(2,3,4,e1,e2,e3,fail)
27 print *,'after refsys'
36 print *,'dc',dc(1,0),dc(2,0),dc(3,0)
54 veclen=veclen+(c(i,2)-c(i,1))**2
78 call locate_side_chain(2)
82 if (theti.ne.theti) theti=100.0
84 if (phii.ne.phii) phii=180.0
93 * Define the matrices of the rotation about the virtual-bond valence angles
94 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
95 * program), R(i,j,k), and, the cumulative matrices of rotation RT
117 rt(2,1,i-2)=sint*cosphi
118 rt(2,2,i-2)=-cost*cosphi
120 rt(3,1,i-2)=-sint*sinphi
121 rt(3,2,i-2)=cost*sinphi
123 call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
125 dc_norm(j,i-1)=prod(j,1,i-1)
126 dc(j,i-1)=vbld(i)*prod(j,1,i-1)
128 call locate_side_chain(i-1)
132 C Define the origin and orientation of the coordinate system and locate the
133 C first three CA's and SC(2).
137 * Build the alpha-carbon chain.
140 call locate_next_res(i)
143 C First and last SC must coincide with the corresponding CA.
147 dc_norm(j,nres+1)=0.0D0
148 dc(j,nres+nres)=0.0D0
149 dc_norm(j,nres+nres)=0.0D0
151 c(j,nres+nres)=c(j,nres)
154 * Temporary diagnosis
159 write (iout,'(/a)') 'Recalculated internal coordinates'
162 c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
165 if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
166 be1=rad2deg*beta(nres+i,i,maxres2,i+1)
168 if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
169 write (iout,1212) restyp(itype(i)),i,dist(i-1,i),
170 & alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,maxres2),be1
172 1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
178 c-------------------------------------------------------------------------
179 subroutine orig_frame
181 C Define the origin and orientation of the coordinate system and locate
182 C the first three atoms.
184 implicit real*8 (a-h,o-z)
186 include 'COMMON.CHAIN'
187 include 'COMMON.LOCAL'
241 dc_norm(j,2)=prod(j,1,2)
242 dc(j,2)=vbld(3)*prod(j,1,2)
243 c(j,3)=c(j,2)+dc(j,2)
245 call locate_side_chain(2)
248 c-----------------------------------------------------------------------------
249 subroutine locate_next_res(i)
251 C Locate CA(i) and SC(i-1)
253 implicit real*8 (a-h,o-z)
255 include 'COMMON.CHAIN'
256 include 'COMMON.LOCAL'
259 include 'COMMON.IOUNITS'
260 include 'COMMON.NAMES'
261 include 'COMMON.INTERACT'
263 C Define the rotation matrices corresponding to CA(i)
267 if (theti.ne.theti) theti=100.0
269 if (phii.ne.phii) phii=180.0
278 * Define the matrices of the rotation about the virtual-bond valence angles
279 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
280 * program), R(i,j,k), and, the cumulative matrices of rotation RT
302 rt(2,1,i-2)=sint*cosphi
303 rt(2,2,i-2)=-cost*cosphi
305 rt(3,1,i-2)=-sint*sinphi
306 rt(3,2,i-2)=cost*sinphi
308 call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
310 dc_norm(j,i-1)=prod(j,1,i-1)
311 dc(j,i-1)=vbld(i)*prod(j,1,i-1)
312 c(j,i)=c(j,i-1)+dc(j,i-1)
314 cd print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
316 C Now calculate the coordinates of SC(i-1)
318 call locate_side_chain(i-1)
321 c-----------------------------------------------------------------------------
322 subroutine locate_side_chain(i)
324 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
326 implicit real*8 (a-h,o-z)
328 include 'COMMON.CHAIN'
329 include 'COMMON.LOCAL'
332 include 'COMMON.IOUNITS'
333 include 'COMMON.NAMES'
334 include 'COMMON.INTERACT'
338 c dsci_inv=dsc_inv(itype(i))
340 dsci_inv=vbld_inv(i+nres)
344 if (alphi.ne.alphi) alphi=100.0
345 if (omegi.ne.omegi) omegi=-100.0
355 yp= dsci*sinalphi*cosomegi
356 zp=-dsci*sinalphi*sinomegi
357 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
358 * X-axis aligned with the vector DC(*,i)
359 theta2=pi-0.5D0*theta(i+1)
362 xx(1)= xp*cost2+yp*sint2
363 xx(2)=-xp*sint2+yp*cost2
365 cd print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
366 cd & xp,yp,zp,(xx(k),k=1,3)
370 * Bring the SC vectors to the common coordinate system.
372 xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
373 xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
380 rj=rj+prod(j,k,i-1)*xx(k)
383 dc_norm(j,nres+i)=rj*dsci_inv
384 c(j,nres+i)=c(j,i)+rj
388 c------------------------------------------
394 include 'COMMON.CONTROL'
398 include 'COMMON.LANGEVIN'
400 include 'COMMON.LANGEVIN.lang0'
402 include 'COMMON.CHAIN'
403 include 'COMMON.DERIV'
405 include 'COMMON.LOCAL'
406 include 'COMMON.INTERACT'
407 include 'COMMON.IOUNITS'
408 include 'COMMON.NAMES'
409 include 'COMMON.TIME1'
410 include 'COMMON.REMD'
411 include 'COMMON.SETUP'
412 include 'COMMON.MUCA'
413 include 'COMMON.HAIRPIN'
417 C write(*,*) 'initial', i,j,c(j,i)
420 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
422 if (allareout.eq.1) then
423 ireturnval=int(c(j,i)/boxxsize)
424 if (c(j,i).le.0) ireturnval=ireturnval-1
425 do k=chain_beg,chain_end
426 c(j,k)=c(j,k)-ireturnval*boxxsize
427 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
433 if (int(c(j,i)/boxxsize).eq.0) allareout=0
436 if (allareout.eq.1) then
437 ireturnval=int(c(j,i)/boxxsize)
438 if (c(j,i).le.0) ireturnval=ireturnval-1
440 c(j,k)=c(j,k)-ireturnval*boxxsize
441 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
446 C write(*,*) 'befor no jump', i,j,c(j,i)
450 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
451 difference=abs(c(j,i-1)-c(j,i))
452 C print *,'diff', difference
453 if (difference.gt.boxxsize/2.0) then
454 if (c(j,i-1).gt.c(j,i)) then
463 c(j,i)=c(j,i)+nojumpval*boxxsize
464 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
468 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
469 difference=abs(c(j,i-1)-c(j,i))
470 if (difference.gt.boxxsize/2.0) then
471 if (c(j,i-1).gt.c(j,i)) then
480 c(j,i)=c(j,i)+nojumpval*boxxsize
481 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
485 C write(*,*) 'after no jump', i,j,c(j,i)
492 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
494 if (allareout.eq.1) then
495 ireturnval=int(c(j,i)/boxysize)
496 if (c(j,i).le.0) ireturnval=ireturnval-1
497 do k=chain_beg,chain_end
498 c(j,k)=c(j,k)-ireturnval*boxysize
499 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
505 if (int(c(j,i)/boxysize).eq.0) allareout=0
508 if (allareout.eq.1) then
509 ireturnval=int(c(j,i)/boxysize)
510 if (c(j,i).le.0) ireturnval=ireturnval-1
512 c(j,k)=c(j,k)-ireturnval*boxysize
513 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
518 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
519 difference=abs(c(j,i-1)-c(j,i))
520 if (difference.gt.boxysize/2.0) then
521 if (c(j,i-1).gt.c(j,i)) then
530 c(j,i)=c(j,i)+nojumpval*boxysize
531 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
535 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
536 difference=abs(c(j,i-1)-c(j,i))
537 if (difference.gt.boxysize/2.0) then
538 if (c(j,i-1).gt.c(j,i)) then
547 c(j,i)=c(j,i)+nojumpval*boxysize
548 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
553 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
555 if (allareout.eq.1) then
556 ireturnval=int(c(j,i)/boxysize)
557 if (c(j,i).le.0) ireturnval=ireturnval-1
558 do k=chain_beg,chain_end
559 c(j,k)=c(j,k)-ireturnval*boxzsize
560 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
566 if (int(c(j,i)/boxzsize).eq.0) allareout=0
569 if (allareout.eq.1) then
570 ireturnval=int(c(j,i)/boxzsize)
571 if (c(j,i).le.0) ireturnval=ireturnval-1
573 c(j,k)=c(j,k)-ireturnval*boxzsize
574 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
579 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
580 difference=abs(c(j,i-1)-c(j,i))
581 if (difference.gt.(boxzsize/2.0)) then
582 if (c(j,i-1).gt.c(j,i)) then
591 c(j,i)=c(j,i)+nojumpval*boxzsize
592 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
596 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
597 difference=abs(c(j,i-1)-c(j,i))
598 if (difference.gt.boxzsize/2.0) then
599 if (c(j,i-1).gt.c(j,i)) then
608 c(j,i)=c(j,i)+nojumpval*boxzsize
609 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize