subroutine chainbuild C C Build the virtual polypeptide chain. Side-chain centroids are moveable. C As of 2/17/95. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.LOCAL' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.INTERACT' double precision e1(3),e2(3),e3(3) logical lprn,perbox,fail C Set lprn=.true. for debugging lprn = .false. perbox=.false. fail=.false. print *, 'enter chainbuild' call chainbuild_cart return end C#ifdef DEBUG C if (perbox) then C first three CA's and SC(2). C subroutine chainbuild_extconf C C Build the virtual polypeptide chain. Side-chain centroids are moveable. C As of 2/17/95. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.LOCAL' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.INTERACT' double precision e1(3),e2(3),e3(3) logical lprn,perbox,fail call orig_frame * * Build the alpha-carbon chain. * do i=4,nres call locate_next_res(i) enddo C C First and last SC must coincide with the corresponding CA. C do j=1,3 dc(j,nres+1)=0.0D0 dc_norm(j,nres+1)=0.0D0 dc(j,nres+nres)=0.0D0 dc_norm(j,nres+nres)=0.0D0 c(j,nres+1)=c(j,1) c(j,nres+nres)=c(j,nres) enddo * * Temporary diagnosis * if (lprn) then call cartprint write (iout,'(/a)') 'Recalculated internal coordinates' do i=2,nres-1 do j=1,3 c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1)) enddo be=0.0D0 if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i) be1=rad2deg*beta(nres+i,i,maxres2,i+1) alfai=0.0D0 if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i) write (iout,1212) restyp(itype(i)),i,dist(i-1,i), & alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,maxres2),be1 enddo 1212 format (a3,'(',i3,')',2(f10.5,2f10.2)) C endif endif return end C#endif c------------------------------------------------------------------------- subroutine orig_frame C C Define the origin and orientation of the coordinate system and locate C the first three atoms. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.LOCAL' include 'COMMON.GEO' include 'COMMON.VAR' cost=dcos(theta(3)) sint=dsin(theta(3)) t(1,1,1)=-cost t(1,2,1)=-sint t(1,3,1)= 0.0D0 t(2,1,1)=-sint t(2,2,1)= cost t(2,3,1)= 0.0D0 t(3,1,1)= 0.0D0 t(3,2,1)= 0.0D0 t(3,3,1)= 1.0D0 r(1,1,1)= 1.0D0 r(1,2,1)= 0.0D0 r(1,3,1)= 0.0D0 r(2,1,1)= 0.0D0 r(2,2,1)= 1.0D0 r(2,3,1)= 0.0D0 r(3,1,1)= 0.0D0 r(3,2,1)= 0.0D0 r(3,3,1)= 1.0D0 do i=1,3 do j=1,3 rt(i,j,1)=t(i,j,1) enddo enddo do i=1,3 do j=1,3 prod(i,j,1)=0.0D0 prod(i,j,2)=t(i,j,1) enddo prod(i,i,1)=1.0D0 enddo c(1,1)=0.0D0 c(2,1)=0.0D0 c(3,1)=0.0D0 c(1,2)=vbld(2) c(2,2)=0.0D0 c(3,2)=0.0D0 dc(1,0)=0.0d0 dc(2,0)=0.0D0 dc(3,0)=0.0D0 dc(1,1)=vbld(2) dc(2,1)=0.0D0 dc(3,1)=0.0D0 dc_norm(1,0)=0.0D0 dc_norm(2,0)=0.0D0 dc_norm(3,0)=0.0D0 dc_norm(1,1)=1.0D0 dc_norm(2,1)=0.0D0 dc_norm(3,1)=0.0D0 do j=1,3 dc_norm(j,2)=prod(j,1,2) dc(j,2)=vbld(3)*prod(j,1,2) c(j,3)=c(j,2)+dc(j,2) enddo call locate_side_chain(2) return end c----------------------------------------------------------------------------- subroutine locate_next_res(i) C C Locate CA(i) and SC(i-1) C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.LOCAL' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.INTERACT' C C Define the rotation matrices corresponding to CA(i) C #ifdef OSF theti=theta(i) if (theti.ne.theti) theti=100.0 phii=phi(i) if (phii.ne.phii) phii=180.0 #else theti=theta(i) phii=phi(i) #endif cost=dcos(theti) sint=dsin(theti) cosphi=dcos(phii) sinphi=dsin(phii) * Define the matrices of the rotation about the virtual-bond valence angles * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this * program), R(i,j,k), and, the cumulative matrices of rotation RT t(1,1,i-2)=-cost t(1,2,i-2)=-sint t(1,3,i-2)= 0.0D0 t(2,1,i-2)=-sint t(2,2,i-2)= cost t(2,3,i-2)= 0.0D0 t(3,1,i-2)= 0.0D0 t(3,2,i-2)= 0.0D0 t(3,3,i-2)= 1.0D0 r(1,1,i-2)= 1.0D0 r(1,2,i-2)= 0.0D0 r(1,3,i-2)= 0.0D0 r(2,1,i-2)= 0.0D0 r(2,2,i-2)=-cosphi r(2,3,i-2)= sinphi r(3,1,i-2)= 0.0D0 r(3,2,i-2)= sinphi r(3,3,i-2)= cosphi rt(1,1,i-2)=-cost rt(1,2,i-2)=-sint rt(1,3,i-2)=0.0D0 rt(2,1,i-2)=sint*cosphi rt(2,2,i-2)=-cost*cosphi rt(2,3,i-2)=sinphi rt(3,1,i-2)=-sint*sinphi rt(3,2,i-2)=cost*sinphi rt(3,3,i-2)=cosphi call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1)) do j=1,3 dc_norm(j,i-1)=prod(j,1,i-1) dc(j,i-1)=vbld(i)*prod(j,1,i-1) c(j,i)=c(j,i-1)+dc(j,i-1) enddo cd print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3) C C Now calculate the coordinates of SC(i-1) C call locate_side_chain(i-1) return end c----------------------------------------------------------------------------- subroutine locate_side_chain(i) C C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i). C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.LOCAL' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.INTERACT' dimension xx(3) c dsci=dsc(itype(i)) c dsci_inv=dsc_inv(itype(i)) dsci=vbld(i+nres) dsci_inv=vbld_inv(i+nres) #ifdef OSF alphi=alph(i) omegi=omeg(i) if (alphi.ne.alphi) alphi=100.0 if (omegi.ne.omegi) omegi=-100.0 #else alphi=alph(i) omegi=omeg(i) #endif cosalphi=dcos(alphi) sinalphi=dsin(alphi) cosomegi=dcos(omegi) sinomegi=dsin(omegi) xp= dsci*cosalphi yp= dsci*sinalphi*cosomegi zp=-dsci*sinalphi*sinomegi * Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its * X-axis aligned with the vector DC(*,i) theta2=pi-0.5D0*theta(i+1) cost2=dcos(theta2) sint2=dsin(theta2) xx(1)= xp*cost2+yp*sint2 xx(2)=-xp*sint2+yp*cost2 xx(3)= zp cd print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i, cd & xp,yp,zp,(xx(k),k=1,3) do j=1,3 xloc(j,i)=xx(j) enddo * Bring the SC vectors to the common coordinate system. xx(1)=xloc(1,i) xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1) xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1) do j=1,3 xrot(j,i)=xx(j) enddo do j=1,3 rj=0.0D0 do k=1,3 rj=rj+prod(j,k,i-1)*xx(k) enddo dc(j,nres+i)=rj dc_norm(j,nres+i)=rj*dsci_inv c(j,nres+i)=c(j,i)+rj enddo return end c------------------------------------------ subroutine returnbox include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' #ifndef LANG0 include 'COMMON.LANGEVIN' #else include 'COMMON.LANGEVIN.lang0' #endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' include 'COMMON.REMD' include 'COMMON.SETUP' include 'COMMON.MUCA' include 'COMMON.HAIRPIN' C change suggested by Ana - begin integer allareout C change suggested by Ana - end j=1 chain_beg=1 C do i=1,nres C write(*,*) 'initial', i,j,c(j,i) C enddo C change suggested by Ana - begin allareout=1 C change suggested by Ana -end do i=1,nres-1 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then chain_end=i if (allareout.eq.1) then ireturnval=int(c(j,i)/boxxsize) if (c(j,i).le.0) ireturnval=ireturnval-1 do k=chain_beg,chain_end c(j,k)=c(j,k)-ireturnval*boxxsize c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize enddo C Suggested by Ana if (chain_beg.eq.1) & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize C Suggested by Ana -end endif chain_beg=i+1 allareout=1 else if (int(c(j,i)/boxxsize).eq.0) allareout=0 endif enddo if (allareout.eq.1) then ireturnval=int(c(j,i)/boxxsize) if (c(j,i).le.0) ireturnval=ireturnval-1 do k=chain_beg,nres c(j,k)=c(j,k)-ireturnval*boxxsize c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize enddo endif C NO JUMP C do i=1,nres C write(*,*) 'befor no jump', i,j,c(j,i) C enddo nojumpval=0 do i=2,nres if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then difference=abs(c(j,i-1)-c(j,i)) C print *,'diff', difference if (difference.gt.boxxsize/2.0) then if (c(j,i-1).gt.c(j,i)) then nojumpval=1 else nojumpval=-1 endif else nojumpval=0 endif endif c(j,i)=c(j,i)+nojumpval*boxxsize c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize enddo nojumpval=0 do i=2,nres if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then difference=abs(c(j,i-1)-c(j,i)) if (difference.gt.boxxsize/2.0) then if (c(j,i-1).gt.c(j,i)) then nojumpval=1 else nojumpval=-1 endif else nojumpval=0 endif endif c(j,i)=c(j,i)+nojumpval*boxxsize c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize enddo C do i=1,nres C write(*,*) 'after no jump', i,j,c(j,i) C enddo C NOW Y dimension C suggesed by Ana begins allareout=1 C suggested by Ana ends j=2 chain_beg=1 do i=1,nres-1 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then chain_end=i if (allareout.eq.1) then ireturnval=int(c(j,i)/boxysize) if (c(j,i).le.0) ireturnval=ireturnval-1 do k=chain_beg,chain_end c(j,k)=c(j,k)-ireturnval*boxysize c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize enddo C Suggested by Ana if (chain_beg.eq.1) & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize C Suggested by Ana -end endif chain_beg=i+1 allareout=1 else if (int(c(j,i)/boxysize).eq.0) allareout=0 endif enddo if (allareout.eq.1) then ireturnval=int(c(j,i)/boxysize) if (c(j,i).le.0) ireturnval=ireturnval-1 do k=chain_beg,nres c(j,k)=c(j,k)-ireturnval*boxysize c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize enddo endif nojumpval=0 do i=2,nres if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then difference=abs(c(j,i-1)-c(j,i)) if (difference.gt.boxysize/2.0) then if (c(j,i-1).gt.c(j,i)) then nojumpval=1 else nojumpval=-1 endif else nojumpval=0 endif endif c(j,i)=c(j,i)+nojumpval*boxysize c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize enddo nojumpval=0 do i=2,nres if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then difference=abs(c(j,i-1)-c(j,i)) if (difference.gt.boxysize/2.0) then if (c(j,i-1).gt.c(j,i)) then nojumpval=1 else nojumpval=-1 endif else nojumpval=0 endif endif c(j,i)=c(j,i)+nojumpval*boxysize c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize enddo C Now Z dimension C Suggested by Ana -begins allareout=1 C Suggested by Ana -ends j=3 chain_beg=1 do i=1,nres-1 if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then chain_end=i if (allareout.eq.1) then ireturnval=int(c(j,i)/boxysize) if (c(j,i).le.0) ireturnval=ireturnval-1 do k=chain_beg,chain_end c(j,k)=c(j,k)-ireturnval*boxzsize c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize enddo C Suggested by Ana if (chain_beg.eq.1) & dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize C Suggested by Ana -end endif chain_beg=i+1 allareout=1 else if (int(c(j,i)/boxzsize).eq.0) allareout=0 endif enddo if (allareout.eq.1) then ireturnval=int(c(j,i)/boxzsize) if (c(j,i).le.0) ireturnval=ireturnval-1 do k=chain_beg,nres c(j,k)=c(j,k)-ireturnval*boxzsize c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize enddo endif nojumpval=0 do i=2,nres if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then difference=abs(c(j,i-1)-c(j,i)) if (difference.gt.(boxzsize/2.0)) then if (c(j,i-1).gt.c(j,i)) then nojumpval=1 else nojumpval=-1 endif else nojumpval=0 endif endif c(j,i)=c(j,i)+nojumpval*boxzsize c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize enddo nojumpval=0 do i=2,nres if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then difference=abs(c(j,i-1)-c(j,i)) if (difference.gt.boxzsize/2.0) then if (c(j,i-1).gt.c(j,i)) then nojumpval=1 else nojumpval=-1 endif else nojumpval=0 endif endif c(j,i)=c(j,i)+nojumpval*boxzsize c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize enddo return end