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
26 C first three CAs and SC(2).
28 subroutine chainbuild_extconf
30 C Build the virtual polypeptide chain. Side-chain centroids are moveable.
33 implicit real*8 (a-h,o-z)
35 include 'COMMON.CHAIN'
36 include 'COMMON.LOCAL'
39 include 'COMMON.IOUNITS'
40 include 'COMMON.NAMES'
41 include 'COMMON.INTERACT'
42 double precision e1(3),e2(3),e3(3)
43 logical lprn /.false./,perbox,fail
45 c write (iout,*) "Calling chainbuild_extconf"
48 * Build the alpha-carbon chain.
51 call locate_next_res(i)
54 C First and last SC must coincide with the corresponding CA.
58 dc_norm(j,nres+1)=0.0D0
60 dc_norm(j,nres+nres)=0.0D0
62 c(j,nres+nres)=c(j,nres)
70 write (iout,*) 'dc_norm'
72 write (iout,'(a3,1x,i3,3f10.5,5x,3f10.5)')
73 & restyp(itype(i)),i,(dc_norm(j,i),j=1,3),
74 & (dc_norm(j,i+nres),j=1,3)
76 write (iout,'(/a)') 'Recalculated internal coordinates'
79 c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
82 if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
83 be1=rad2deg*beta(nres+i,i,maxres2,i+1)
85 if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
86 write (iout,1212) restyp(itype(i)),i,dist(i-1,i),
87 & alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,maxres2),be1
89 1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
96 c-------------------------------------------------------------------------
99 C Define the origin and orientation of the coordinate system and locate
100 C the first three atoms.
102 implicit real*8 (a-h,o-z)
104 include 'COMMON.CHAIN'
105 include 'COMMON.LOCAL'
159 dc_norm(j,2)=prod(j,1,2)
160 dc(j,2)=vbld(3)*prod(j,1,2)
161 c(j,3)=c(j,2)+dc(j,2)
163 call locate_side_chain(2)
166 c-----------------------------------------------------------------------------
167 subroutine locate_prev_res(i)
169 C Locate CA(i) and SC(i-1)
171 implicit real*8 (a-h,o-z)
173 include 'COMMON.CHAIN'
174 include 'COMMON.LOCAL'
177 include 'COMMON.IOUNITS'
178 include 'COMMON.NAMES'
179 include 'COMMON.INTERACT'
181 C Define the rotation matrices corresponding to CA(i)
185 if (theti.ne.theti) theti=100.0
187 if (phii.ne.phii) phii=180.0
196 * Define the matrices of the rotation about the virtual-bond valence angles
197 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
198 * program), R(i,j,k), and, the cumulative matrices of rotation RT
220 rt(2,1,i+2)=sint*cosphi
221 rt(2,2,i+2)=-cost*cosphi
223 rt(3,1,i+2)=-sint*sinphi
224 rt(3,2,i+2)=cost*sinphi
226 call matmult(prod(1,1,i+2),rt(1,1,i+2),prod(1,1,i+1))
228 dc_norm(j,i)=-prod(j,1,i+1)
229 dc(j,i)=-vbld(i)*prod(j,1,i+1)
230 c(j,i)=c(j,i+1)+dc(j,i)
232 cd print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i+1),j=1,3),(c(j,i),j=1,3)
234 C Now calculate the coordinates of SC(i+1)
236 call locate_side_chain(i+1)
239 c-----------------------------------------------------------------------------
240 subroutine locate_next_res(i)
242 C Locate CA(i) and SC(i-1)
244 implicit real*8 (a-h,o-z)
246 include 'COMMON.CHAIN'
247 include 'COMMON.LOCAL'
250 include 'COMMON.IOUNITS'
251 include 'COMMON.NAMES'
252 include 'COMMON.INTERACT'
254 C Define the rotation matrices corresponding to CA(i)
258 if (theti.ne.theti) theti=100.0
260 if (phii.ne.phii) phii=180.0
269 * Define the matrices of the rotation about the virtual-bond valence angles
270 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
271 * program), R(i,j,k), and, the cumulative matrices of rotation RT
293 rt(2,1,i-2)=sint*cosphi
294 rt(2,2,i-2)=-cost*cosphi
296 rt(3,1,i-2)=-sint*sinphi
297 rt(3,2,i-2)=cost*sinphi
299 call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
301 dc_norm(j,i-1)=prod(j,1,i-1)
302 dc(j,i-1)=vbld(i)*prod(j,1,i-1)
303 c(j,i)=c(j,i-1)+dc(j,i-1)
305 cd print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
307 C Now calculate the coordinates of SC(i-1)
309 call locate_side_chain(i-1)
312 c-----------------------------------------------------------------------------
313 subroutine locate_side_chain(i)
315 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
317 implicit real*8 (a-h,o-z)
319 include 'COMMON.CHAIN'
320 include 'COMMON.LOCAL'
323 include 'COMMON.IOUNITS'
324 include 'COMMON.NAMES'
325 include 'COMMON.INTERACT'
329 c dsci_inv=dsc_inv(itype(i))
331 dsci_inv=vbld_inv(i+nres)
335 if (alphi.ne.alphi) alphi=100.0
336 if (omegi.ne.omegi) omegi=-100.0
346 yp= dsci*sinalphi*cosomegi
347 zp=-dsci*sinalphi*sinomegi
348 * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
349 * X-axis aligned with the vector DC(*,i)
350 theta2=pi-0.5D0*theta(i+1)
353 xx(1)= xp*cost2+yp*sint2
354 xx(2)=-xp*sint2+yp*cost2
356 cd print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
357 cd & xp,yp,zp,(xx(k),k=1,3)
361 * Bring the SC vectors to the common coordinate system.
363 xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
364 xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
371 rj=rj+prod(j,k,i-1)*xx(k)
374 dc_norm(j,nres+i)=rj*dsci_inv
375 c(j,nres+i)=c(j,i)+rj
379 c------------------------------------------
385 include 'COMMON.CONTROL'
389 include 'COMMON.LANGEVIN'
391 include 'COMMON.LANGEVIN.lang0'
393 include 'COMMON.CHAIN'
394 include 'COMMON.DERIV'
396 include 'COMMON.LOCAL'
397 include 'COMMON.INTERACT'
398 include 'COMMON.IOUNITS'
399 include 'COMMON.NAMES'
400 include 'COMMON.TIME1'
401 include 'COMMON.REMD'
402 include 'COMMON.SETUP'
403 include 'COMMON.MUCA'
404 include 'COMMON.HAIRPIN'
405 C change suggested by Ana - begin
407 C change suggested by Ana - end
411 C write(*,*) 'initial', i,j,c(j,i)
413 C change suggested by Ana - begin
415 C change suggested by Ana -end
417 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
419 if (allareout.eq.1) then
420 ireturnval=int(c(j,i)/boxxsize)
421 if (c(j,i).le.0) ireturnval=ireturnval-1
422 do k=chain_beg,chain_end
423 c(j,k)=c(j,k)-ireturnval*boxxsize
424 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
428 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
429 C Suggested by Ana -end
434 if (int(c(j,i)/boxxsize).eq.0) allareout=0
437 if (allareout.eq.1) then
438 ireturnval=int(c(j,i)/boxxsize)
439 if (c(j,i).le.0) ireturnval=ireturnval-1
441 c(j,k)=c(j,k)-ireturnval*boxxsize
442 c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
447 C write(*,*) 'befor no jump', i,j,c(j,i)
451 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
452 difference=abs(c(j,i-1)-c(j,i))
453 C print *,'diff', difference
454 if (difference.gt.boxxsize/2.0) then
455 if (c(j,i-1).gt.c(j,i)) then
464 c(j,i)=c(j,i)+nojumpval*boxxsize
465 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
469 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
470 difference=abs(c(j,i-1)-c(j,i))
471 if (difference.gt.boxxsize/2.0) then
472 if (c(j,i-1).gt.c(j,i)) then
481 c(j,i)=c(j,i)+nojumpval*boxxsize
482 c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
486 C write(*,*) 'after no jump', i,j,c(j,i)
490 C suggesed by Ana begins
492 C suggested by Ana ends
496 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
498 if (allareout.eq.1) then
499 ireturnval=int(c(j,i)/boxysize)
500 if (c(j,i).le.0) ireturnval=ireturnval-1
501 do k=chain_beg,chain_end
502 c(j,k)=c(j,k)-ireturnval*boxysize
503 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
507 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
508 C Suggested by Ana -end
513 if (int(c(j,i)/boxysize).eq.0) allareout=0
516 if (allareout.eq.1) then
517 ireturnval=int(c(j,i)/boxysize)
518 if (c(j,i).le.0) ireturnval=ireturnval-1
520 c(j,k)=c(j,k)-ireturnval*boxysize
521 c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
526 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
527 difference=abs(c(j,i-1)-c(j,i))
528 if (difference.gt.boxysize/2.0) then
529 if (c(j,i-1).gt.c(j,i)) then
538 c(j,i)=c(j,i)+nojumpval*boxysize
539 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
543 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
544 difference=abs(c(j,i-1)-c(j,i))
545 if (difference.gt.boxysize/2.0) then
546 if (c(j,i-1).gt.c(j,i)) then
555 c(j,i)=c(j,i)+nojumpval*boxysize
556 c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
559 C Suggested by Ana -begins
561 C Suggested by Ana -ends
565 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
567 if (allareout.eq.1) then
568 ireturnval=int(c(j,i)/boxysize)
569 if (c(j,i).le.0) ireturnval=ireturnval-1
570 do k=chain_beg,chain_end
571 c(j,k)=c(j,k)-ireturnval*boxzsize
572 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
576 & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
577 C Suggested by Ana -end
582 if (int(c(j,i)/boxzsize).eq.0) allareout=0
585 if (allareout.eq.1) then
586 ireturnval=int(c(j,i)/boxzsize)
587 if (c(j,i).le.0) ireturnval=ireturnval-1
589 c(j,k)=c(j,k)-ireturnval*boxzsize
590 c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
595 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
596 difference=abs(c(j,i-1)-c(j,i))
597 if (difference.gt.(boxzsize/2.0)) then
598 if (c(j,i-1).gt.c(j,i)) then
607 c(j,i)=c(j,i)+nojumpval*boxzsize
608 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
612 if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
613 difference=abs(c(j,i-1)-c(j,i))
614 if (difference.gt.boxzsize/2.0) then
615 if (c(j,i-1).gt.c(j,i)) then
624 c(j,i)=c(j,i)+nojumpval*boxzsize
625 c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize