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'
414 C change suggested by Ana - begin
416 C change suggested by Ana - end
420 C write(*,*) 'initial', i,j,c(j,i)
422 C change suggested by Ana - begin
424 C change suggested by Ana -end
426 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
428 if (allareout.eq.1) then
429 ireturnval=int(c(j,i)/boxxsize)
430 if (c(j,i).le.0) ireturnval=ireturnval-1
431 do k=chain_beg,chain_end
432 c(j,k)=c(j,k)-ireturnval*boxxsize
433 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
437 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
438 C Suggested by Ana -end
443 if (int(c(j,i)/boxxsize).eq.0) allareout=0
446 if (allareout.eq.1) then
447 ireturnval=int(c(j,i)/boxxsize)
448 if (c(j,i).le.0) ireturnval=ireturnval-1
450 c(j,k)=c(j,k)-ireturnval*boxxsize
451 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
456 C write(*,*) 'befor no jump', i,j,c(j,i)
460 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
461 difference=abs(c(j,i-1)-c(j,i))
462 C print *,'diff', difference
463 if (difference.gt.boxxsize/2.0) then
464 if (c(j,i-1).gt.c(j,i)) then
473 c(j,i)=c(j,i)+nojumpval*boxxsize
474 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
478 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
479 difference=abs(c(j,i-1)-c(j,i))
480 if (difference.gt.boxxsize/2.0) then
481 if (c(j,i-1).gt.c(j,i)) then
490 c(j,i)=c(j,i)+nojumpval*boxxsize
491 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
495 C write(*,*) 'after no jump', i,j,c(j,i)
499 C suggesed by Ana begins
501 C suggested by Ana ends
505 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
507 if (allareout.eq.1) then
508 ireturnval=int(c(j,i)/boxysize)
509 if (c(j,i).le.0) ireturnval=ireturnval-1
510 do k=chain_beg,chain_end
511 c(j,k)=c(j,k)-ireturnval*boxysize
512 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
516 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
517 C Suggested by Ana -end
522 if (int(c(j,i)/boxysize).eq.0) allareout=0
525 if (allareout.eq.1) then
526 ireturnval=int(c(j,i)/boxysize)
527 if (c(j,i).le.0) ireturnval=ireturnval-1
529 c(j,k)=c(j,k)-ireturnval*boxysize
530 c(j,k+nres)=c(j,k+nres)-ireturnval*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
552 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
553 difference=abs(c(j,i-1)-c(j,i))
554 if (difference.gt.boxysize/2.0) then
555 if (c(j,i-1).gt.c(j,i)) then
564 c(j,i)=c(j,i)+nojumpval*boxysize
565 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
568 C Suggested by Ana -begins
570 C Suggested by Ana -ends
574 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
576 if (allareout.eq.1) then
577 ireturnval=int(c(j,i)/boxysize)
578 if (c(j,i).le.0) ireturnval=ireturnval-1
579 do k=chain_beg,chain_end
580 c(j,k)=c(j,k)-ireturnval*boxzsize
581 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
585 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
586 C Suggested by Ana -end
591 if (int(c(j,i)/boxzsize).eq.0) allareout=0
594 if (allareout.eq.1) then
595 ireturnval=int(c(j,i)/boxzsize)
596 if (c(j,i).le.0) ireturnval=ireturnval-1
598 c(j,k)=c(j,k)-ireturnval*boxzsize
599 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
604 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
605 difference=abs(c(j,i-1)-c(j,i))
606 if (difference.gt.(boxzsize/2.0)) then
607 if (c(j,i-1).gt.c(j,i)) then
616 c(j,i)=c(j,i)+nojumpval*boxzsize
617 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
621 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
622 difference=abs(c(j,i-1)-c(j,i))
623 if (difference.gt.boxzsize/2.0) then
624 if (c(j,i-1).gt.c(j,i)) then
633 c(j,i)=c(j,i)+nojumpval*boxzsize
634 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize