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'
419 C change suggested by Ana - begin
421 C change suggested by Ana - end
425 C write(*,*) 'initial', i,j,c(j,i)
427 C change suggested by Ana - begin
429 C change suggested by Ana -end
431 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
433 if (allareout.eq.1) then
434 ireturnval=int(c(j,i)/boxxsize)
435 if (c(j,i).le.0) ireturnval=ireturnval-1
436 do k=chain_beg,chain_end
437 c(j,k)=c(j,k)-ireturnval*boxxsize
438 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
442 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
443 C Suggested by Ana -end
448 if (int(c(j,i)/boxxsize).eq.0) allareout=0
451 if (allareout.eq.1) then
452 ireturnval=int(c(j,i)/boxxsize)
453 if (c(j,i).le.0) ireturnval=ireturnval-1
455 c(j,k)=c(j,k)-ireturnval*boxxsize
456 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
461 C write(*,*) 'befor no jump', i,j,c(j,i)
465 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
466 difference=abs(c(j,i-1)-c(j,i))
467 C print *,'diff', difference
468 if (difference.gt.boxxsize/2.0) then
469 if (c(j,i-1).gt.c(j,i)) then
478 c(j,i)=c(j,i)+nojumpval*boxxsize
479 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
483 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
484 difference=abs(c(j,i-1)-c(j,i))
485 if (difference.gt.boxxsize/2.0) then
486 if (c(j,i-1).gt.c(j,i)) then
495 c(j,i)=c(j,i)+nojumpval*boxxsize
496 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
500 C write(*,*) 'after no jump', i,j,c(j,i)
504 C suggesed by Ana begins
506 C suggested by Ana ends
510 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
512 if (allareout.eq.1) then
513 ireturnval=int(c(j,i)/boxysize)
514 if (c(j,i).le.0) ireturnval=ireturnval-1
515 do k=chain_beg,chain_end
516 c(j,k)=c(j,k)-ireturnval*boxysize
517 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
521 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
522 C Suggested by Ana -end
527 if (int(c(j,i)/boxysize).eq.0) allareout=0
530 if (allareout.eq.1) then
531 ireturnval=int(c(j,i)/boxysize)
532 if (c(j,i).le.0) ireturnval=ireturnval-1
534 c(j,k)=c(j,k)-ireturnval*boxysize
535 c(j,k+nres)=c(j,k+nres)-ireturnval*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
557 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
558 difference=abs(c(j,i-1)-c(j,i))
559 if (difference.gt.boxysize/2.0) then
560 if (c(j,i-1).gt.c(j,i)) then
569 c(j,i)=c(j,i)+nojumpval*boxysize
570 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
573 C Suggested by Ana -begins
575 C Suggested by Ana -ends
579 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
581 if (allareout.eq.1) then
582 ireturnval=int(c(j,i)/boxysize)
583 if (c(j,i).le.0) ireturnval=ireturnval-1
584 do k=chain_beg,chain_end
585 c(j,k)=c(j,k)-ireturnval*boxzsize
586 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
590 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
591 C Suggested by Ana -end
596 if (int(c(j,i)/boxzsize).eq.0) allareout=0
599 if (allareout.eq.1) then
600 ireturnval=int(c(j,i)/boxzsize)
601 if (c(j,i).le.0) ireturnval=ireturnval-1
603 c(j,k)=c(j,k)-ireturnval*boxzsize
604 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
609 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
610 difference=abs(c(j,i-1)-c(j,i))
611 if (difference.gt.(boxzsize/2.0)) then
612 if (c(j,i-1).gt.c(j,i)) then
621 c(j,i)=c(j,i)+nojumpval*boxzsize
622 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
626 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
627 difference=abs(c(j,i-1)-c(j,i))
628 if (difference.gt.boxzsize/2.0) then
629 if (c(j,i-1).gt.c(j,i)) then
638 c(j,i)=c(j,i)+nojumpval*boxzsize
639 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize