subroutine chainbuild C C Build the virtual polypeptide chain. Side-chain centroids are moveable. C As of 2/17/95. C implicit none 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. call chainbuild_cart return end C#ifdef DEBUG C if (perbox) then C first three CAs 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 none include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.LOCAL' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.INTERACT' integer i,j double precision e1(3),e2(3),e3(3) double precision be,be1,alfai double precision xp,yp,zp,cost2,sint2,cosomegi,sinomegi double precision dist,alpha,beta logical lprn,perbox,fail lprn=.false. c write (iout,*) "Calling chainbuild_extconf" 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 none integer i,j double precision cost,sint 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 none integer i,j double precision theti,phii,cost,sint,cosphi,sinphi 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 sc_coord_rebuild(i) C C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i). C implicit none include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.LOCAL' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.INTERACT' integer i,j,k double precision xx(3) double precision dsci,dsci_inv,alphi,omegi,cosalphi,sinalphi, & cosomegi,sinomegi,rp(3),theta2,cost2,sint2,rj double precision scalar double precision ref(3,3),scalp,sscalp,refnorm double precision alpha,beta 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) rp(1)= cosalphi rp(2)= sinalphi*cosomegi rp(3)=-sinalphi*sinomegi c Build the reference system do j=1,3 ref(j,1)=-dc_norm(j,i-1)+dc_norm(j,i) enddo refnorm=dsqrt(scalar(ref(1,1),ref(1,1))) do j=1,3 ref(j,1)=ref(j,1)/refnorm enddo scalp=scalar(ref(1,1),dc_norm(1,i)) sscalp=1.0d0/dsqrt(1.0d0-scalp*scalp) do j=1,3 ref(j,2)=(dc_norm(j,i)-scalp*ref(j,1))*sscalp enddo ref(1,3)= ref(2,1)*ref(3,2)-ref(3,1)*ref(2,2) ref(2,3)=-ref(1,1)*ref(3,2)+ref(3,1)*ref(1,2) ref(3,3)= ref(1,1)*ref(2,2)-ref(2,1)*ref(1,2) c do j=1,3 c write (iout,*) j,scalar(ref(1,j),ref(1,1)), c & scalar(ref(1,j),ref(1,2)),scalar(ref(1,j),ref(1,3)) c enddo c Bring the coordinates to the global reference system do j=1,3 dc_norm(j,nres+i)=0.0d0 do k=1,3 dc_norm(j,nres+i)=dc_norm(j,nres+i)+ref(j,k)*rp(k) enddo dc(j,nres+i)=dc_norm(j,nres+i)*dsci c(j,nres+i)=c(j,i)+dc(j,nres+i) enddo c write (iout,*) scalar(dc_norm(1,i+nres),dc_norm(1,i+nres)), c & dsqrt(scalar(dc(1,i+nres),dc(1,i+nres))) c Check the internal coordinates c c(:,2*nres+1)=ref(:,1)+c(:,i) c write (iout,*) "alpha",rad2deg*alphi, c & rad2deg*alpha(nres+i,i,2*nres+1) c write (iout,*) "omega",rad2deg*omegi, c & rad2deg*beta(nres+i,i,2*nres+1,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 none include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.LOCAL' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.INTERACT' integer i,j,k double precision xx(3) double precision dsci,dsci_inv,alphi,omegi,cosalphi,sinalphi, & cosomegi,sinomegi,xp,yp,zp,theta2,cost2,sint2,rj 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 #ifdef FIVEDIAG include 'COMMON.LANGEVIN.lang0.5diag' #else include 'COMMON.LANGEVIN.lang0' #endif #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 integer i,j 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