From 915aca0187c1157b991cb5f13960e82e3007cca6 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Tue, 23 Jun 2020 10:24:09 +0200 Subject: [PATCH] pdbread-mult --- source/cluster/wham/src-HCD-5D/CMakeLists.txt | 4 +- .../wham/src-HCD-5D/Makefile-MPICH-ifort-okeanos | 2 +- source/cluster/wham/src-HCD-5D/readpdb-mult.F | 956 +++++++++++++++++++ .../unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos | 1 + source/unres/src-HCD-5D/parmread.F | 1 + source/unres/src-HCD-5D/readpdb-mult.F | 43 +- source/unres/src-HCD-5D/readrtns_CSA.F | 3 +- source/wham/src-HCD-5D/CMakeLists.txt | 4 +- .../wham/src-HCD-5D/Makefile_MPICH_ifort-okeanos | 11 +- source/wham/src-HCD-5D/readpdb-mult.F | 960 ++++++++++++++++++++ 10 files changed, 1965 insertions(+), 20 deletions(-) create mode 100644 source/cluster/wham/src-HCD-5D/readpdb-mult.F create mode 100644 source/wham/src-HCD-5D/readpdb-mult.F diff --git a/source/cluster/wham/src-HCD-5D/CMakeLists.txt b/source/cluster/wham/src-HCD-5D/CMakeLists.txt index 904cd4b..f7541e3 100644 --- a/source/cluster/wham/src-HCD-5D/CMakeLists.txt +++ b/source/cluster/wham/src-HCD-5D/CMakeLists.txt @@ -32,7 +32,7 @@ set(UNRES_CLUSTER_WHAM_M_SRC0 printmat.f probabl.F read_coords.F - readpdb.F + readpdb-mult.F readrtns.F rescode.f setup_var.f @@ -73,7 +73,7 @@ set(UNRES_CLUSTER_WHAM_M_PP_SRC rmscalc.F TMscore.F read_constr_homology.F - readpdb.F + readpdb-mult.F ) if(UNRES_DFA) diff --git a/source/cluster/wham/src-HCD-5D/Makefile-MPICH-ifort-okeanos b/source/cluster/wham/src-HCD-5D/Makefile-MPICH-ifort-okeanos index 7ce5b9c..c2d84d8 100644 --- a/source/cluster/wham/src-HCD-5D/Makefile-MPICH-ifort-okeanos +++ b/source/cluster/wham/src-HCD-5D/Makefile-MPICH-ifort-okeanos @@ -19,7 +19,7 @@ LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a object = main_clust.o arcos.o cartprint.o chainbuild.o convert.o initialize_p.o \ matmult.o readrtns.o pinorm.o rescode.o intcor.o timing.o misc.o \ - geomout.o readpdb.o read_coords.o parmread.o probabl.o fitsq.o hc.o \ + geomout.o readpdb-mult.o read_coords.o parmread.o probabl.o fitsq.o hc.o \ track.o wrtclust.o srtclust.o noyes.o contact.o printmat.o \ int_from_cart1.o energy_p_new.o boxshift.o icant.o proc_proc.o \ work_partition.o setup_var.o read_ref_str.o gnmr1.o permut.o \ diff --git a/source/cluster/wham/src-HCD-5D/readpdb-mult.F b/source/cluster/wham/src-HCD-5D/readpdb-mult.F new file mode 100644 index 0000000..c6139ee --- /dev/null +++ b/source/cluster/wham/src-HCD-5D/readpdb-mult.F @@ -0,0 +1,956 @@ + subroutine readpdb +C Read the PDB file and convert the peptide geometry into virtual-chain +C geometry. + implicit none + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + include 'COMMON.SBRIDGE' + include 'COMMON.FRAG' + character*3 seq,atom,res + character*80 card + double precision e1(3),e2(3),e3(3) + double precision sccor(3,50) + integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old + double precision dcj + integer rescode,kkk,lll,icha,cou,kupa,iprzes + logical lsecondary,sccalc,fail,zero + integer iterter(maxres) + double precision efree_temp + iii=0 + ibeg=1 + ishift1=0 + sccalc=.false. + bfac=0.0d0 + do i=1,maxres + iterter(i)=0 + enddo + ibeg=1 + ishift1=0 + lsecondary=.false. + nhfrag=0 + nbfrag=0 + iii=0 + sccalc=.false. + do + read (ipdbin,'(a80)',end=10) card +c write (iout,'(a)') card +c call flush(iout) + if (card(:5).eq.'HELIX') then + nhfrag=nhfrag+1 + lsecondary=.true. + read(card(22:25),*) hfrag(1,nhfrag) + read(card(34:37),*) hfrag(2,nhfrag) + endif + if (card(:5).eq.'SHEET') then + nbfrag=nbfrag+1 + lsecondary=.true. + read(card(24:26),*) bfrag(1,nbfrag) + read(card(35:37),*) bfrag(2,nbfrag) +!rc---------------------------------------- +!rc to be corrected !!! + bfrag(3,nbfrag)=bfrag(1,nbfrag) + bfrag(4,nbfrag)=bfrag(2,nbfrag) +!rc---------------------------------------- + endif + if (card(:3).eq.'END') then + goto 10 + else if (card(:3).eq.'TER') then +! End current chain + ires_old=ires+2 + itype(ires_old-1)=ntyp1 + iterter(ires_old-1)=1 + itype(ires_old)=ntyp1 + iterter(ires_old)=1 +c ishift1=ishift1+1 + ibeg=2 + write (iout,*) "Chain ended",ires,ishift,ires_old,ibeg + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + iii=0 + sccalc=.true. + endif +! Read free energy +c if (index(card,"FREE ENERGY").gt.0) then +c ifree=index(card,"FREE ENERGY")+12 +c read(card(ifree:),*,err=1115,end=1115) efree_temp +c 1115 continue +c endif +! Fish out the ATOM cards. + if (index(card(1:4),'ATOM').gt.0) then + sccalc=.false. + read (card(12:16),*) atom +c write (2,'(a)') card +c write (iout,*) "ibeg",ibeg +c write (iout,*) "! ",atom," !",ires +! if (atom.eq.'CA' .or. atom.eq.'CH3') then + read (card(23:26),*) ires + read (card(18:20),'(a3)') res +c write (iout,*) "ires",ires,ires-ishift+ishift1, +c & " ires_old",ires_old +c write (iout,*) "ishift",ishift," ishift1",ishift1 +c write (iout,*) "IRES",ires-ishift+ishift1,ires_old + if (ires-ishift+ishift1.ne.ires_old) then +! Calculate the CM of the preceding residue. +! if (ibeg.eq.0) call sccenter(ires,iii,sccor) + if (ibeg.eq.0) then +c write (iout,*) "Calculating sidechain center iii",iii +c write (iout,*) "ires",ires + if (unres_pdb) then +c write (iout,'(i5,3f10.5)') ires,(sccor(j,iii),j=1,3) + do j=1,3 + dc(j,ires_old)=sccor(j,iii) + enddo + else + call sccenter(ires_old,iii,sccor) + endif + iii=0 + sccalc=.true. + endif +! Start new residue. +c write (iout,*) "ibeg",ibeg + if (res.eq.'Cl-' .or. res.eq.'Na+') then + ires=ires_old + cycle + else if (ibeg.eq.1) then +c write (iout,*) "BEG ires",ires + ishift=ires-1 + if (res.ne.'GLY' .and. res.ne. 'ACE') then + ishift=ishift-1 + itype(1)=ntyp1 + endif + ires=ires-ishift+ishift1 + ires_old=ires +! write (iout,*) "ishift",ishift," ires",ires,& +! " ires_old",ires_old + ibeg=0 + else if (ibeg.eq.2) then +! Start a new chain + ishift=-ires_old+ires-1 !!!!! +c ishift1=ishift1-1 !!!!! +c write (iout,*) "New chain started",ires,ires_old,ishift, +c & ishift1 + ires=ires-ishift+ishift1 + write (iout,*) "New chain started ires",ires + ires_old=ires +c ires=ires_old+1 + ibeg=0 + else + ishift=ishift-(ires-ishift+ishift1-ires_old-1) + ires=ires-ishift+ishift1 + ires_old=ires + endif + if (res.eq.'ACE' .or. res.eq.'NHE') then + itype(ires)=10 + else + itype(ires)=rescode(ires,res,0) + endif + else + ires=ires-ishift+ishift1 + endif +c write (iout,*) "ires_old",ires_old," ires",ires + if (card(27:27).eq."A" .or. card(27:27).eq."B") then +! ishift1=ishift1+1 + endif +c write (2,*) "ires",ires," res ",res!," ity"!,ity + if (atom.eq.'CA' .or. atom.eq.'CH3' .or. + & res.eq.'NHE'.and.atom(:2).eq.'HN') then + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) +c write (iout,*) "backbone ",atom +c write (iout,*) ires,res,(c(j,ires),j=1,3) +#ifdef DEBUG + write (iout,'(i6,i3,2x,a,3f8.3)') + & ires,itype(ires),res,(c(j,ires),j=1,3) +#endif + iii=iii+1 + do j=1,3 + sccor(j,iii)=c(j,ires) + enddo +c write (2,*) card(23:27),ires,itype(ires),iii + else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. + & atom.ne.'N' .and. atom.ne.'C' .and. + & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. + & atom.ne.'OXT' .and. atom(:2).ne.'3H') then +! write (iout,*) "sidechain ",atom + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) +c write (2,*) "iii",iii + endif + endif + enddo + 10 write (iout,'(a,i5)') ' Nres: ',ires +c write (iout,*) "iii",iii +C Calculate dummy residue coordinates inside the "chain" of a multichain +C system + nres=ires +c write (iout,*) "dc" +c do i=1,nres +c write (iout,'(i5,3f10.5)') i,(dc(j,i),j=1,3) +c enddo + do i=2,nres-1 +c write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i) + if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then + if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then +C 16/01/2014 by Adasko: Adding to dummy atoms in the chain +C first is connected prevous chain (itype(i+1).eq.ntyp1)=true +C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the last dummy residue +c print *,i,'tu dochodze' + call refsys(i-3,i-2,i-1,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif !fail +c print *,i,'a tu?' + do j=1,3 + c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + enddo + else !unres_pdb + do j=1,3 + dcj=(c(j,i-2)-c(j,i-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i-1)+dcj + c(j,nres+i)=c(j,i) + dC(j,i)=c(j,i) + enddo + endif !unres_pdb + else !itype(i+1).eq.ntyp1 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(i+1,i+2,i+3,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + enddo + else !unres_pdb + do j=1,3 + dcj=(c(j,i+3)-c(j,i+2))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i+1)-dcj + c(j,nres+i)=c(j,i) + dC(j,i)=c(j,i) + enddo + endif !unres_pdb + endif !itype(i+1).eq.ntyp1 + endif !itype.eq.ntyp1 + enddo + write (iout,*) "After loop in readpbd" +C Calculate the CM of the last side chain. + if (.not. sccalc) then + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else +c write (iout,*) "Calling sccenter iii",iii + call sccenter(ires,iii,sccor) + endif + endif + nsup=nres + nstart_sup=1 + if (itype(nres).ne.10) then + nres=nres+1 + itype(nres)=ntyp1 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the last dummy residue + call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + enddo + else + do j=1,3 + dcj=(c(j,nres-2)-c(j,nres-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,nres)=c(j,nres-1)+dcj + c(j,2*nres)=c(j,nres) + enddo + endif + endif + do i=2,nres-1 + do j=1,3 + c(j,i+nres)=dc(j,i) + enddo + enddo + do j=1,3 + c(j,nres+1)=c(j,1) + c(j,2*nres)=c(j,nres) + enddo + if (itype(1).eq.ntyp1) then + nsup=nsup-1 + nstart_sup=2 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(2,3,4,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0) + enddo + else + do j=1,3 + dcj=(c(j,4)-c(j,3))/2.0 + c(j,1)=c(j,2)-dcj + c(j,nres+1)=c(j,1) + enddo + endif + endif +C Calculate internal coordinates. + write (iout,'(/a)') + & "Cartesian coordinates of the reference structure" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(c(j,ires),j=1,3), + & (c(j,ires+nres),j=1,3) + zero=.false. + enddo + do ires=1,nres + zero=zero.or.itype(ires).eq.0 + enddo + if (zero) then + write (iout,'(2a)') "Gaps in PDB coordinates detected;", + & " look for ZERO in the control output above." + write (iout,'(2a)') "Repair the PDB file using MODELLER", + & " or other softwared and resubmit." + call flush(iout) + stop + endif +c write(iout,*)"before int_from_cart nres",nres + call int_from_cart(.true.,.false.) + do i=1,nres + thetaref(i)=theta(i) + phiref(i)=phi(i) + enddo + dc(:,0)=c(:,1) + do i=1,nres-1 + do j=1,3 + dc(j,i)=c(j,i+1)-c(j,i) + dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) + enddo +c write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3), +c & vbld_inv(i+1) + enddo + do i=2,nres-1 + do j=1,3 + dc(j,i+nres)=c(j,i+nres)-c(j,i) + dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) + enddo +c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3), +c & vbld_inv(i+nres) + enddo + call sc_loc_geom(.false.) + call int_from_cart1(.false.) +c call chainbuild +C Copy the coordinates to reference coordinates + do i=1,nres + do j=1,3 + cref(j,i)=c(j,i) + cref(j,i+nres)=c(j,i+nres) + enddo + enddo + 100 format (//' alpha-carbon coordinates ', + & ' centroid coordinates'/ + 1 ' ', 7X,'X',11X,'Y',11X,'Z', + & 10X,'X',11X,'Y',11X,'Z') + 110 format (a,'(',i4,')',6f12.5) +cc enddiag + do j=1,nbfrag + do i=1,4 + bfrag(i,j)=bfrag(i,j)-ishift + enddo + enddo + + do j=1,nhfrag + do i=1,2 + hfrag(i,j)=hfrag(i,j)-ishift + enddo + enddo + return + end +c--------------------------------------------------------------------------- + subroutine readpdb_template(k) +C Read the PDB file for read_constr_homology with read2sigma +C and convert the peptide geometry into virtual-chain geometry. + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + include 'COMMON.CONTROL' + include 'COMMON.SETUP' + integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity + logical lprn /.false./,fail + double precision e1(3),e2(3),e3(3) + double precision dcj,efree_temp + character*3 seq,res + character*5 atom + character*80 card + double precision sccor(3,20) + integer rescode,iterter(maxres) + logical zero + do i=1,maxres + iterter(i)=0 + enddo + ibeg=1 + ishift1=0 + ishift=0 +c write (2,*) "UNRES_PDB",unres_pdb + ires=0 + ires_old=0 + iii=0 + lsecondary=.false. + nhfrag=0 + nbfrag=0 + do + read (ipdbin,'(a80)',end=10) card + if (card(:3).eq.'END') then + goto 10 + else if (card(:3).eq.'TER') then +C End current chain + ires_old=ires+2 + itype(ires_old-1)=ntyp1 + iterter(ires_old-1)=1 + itype(ires_old)=ntyp1 + iterter(ires_old)=1 + ibeg=2 +c write (iout,*) "Chain ended",ires,ishift,ires_old + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + endif +C Fish out the ATOM cards. + if (index(card(1:4),'ATOM').gt.0) then + read (card(12:16),*) atom +c write (iout,*) "! ",atom," !",ires +c if (atom.eq.'CA' .or. atom.eq.'CH3') then + read (card(23:26),*) ires + read (card(18:20),'(a3)') res +c write (iout,*) "ires",ires,ires-ishift+ishift1, +c & " ires_old",ires_old +c write (iout,*) "ishift",ishift," ishift1",ishift1 +c write (iout,*) "IRES",ires-ishift+ishift1,ires_old + if (ires-ishift+ishift1.ne.ires_old) then +C Calculate the CM of the preceding residue. + if (ibeg.eq.0) then + if (unres_pdb) then + do j=1,3 + dc(j,ires_old)=sccor(j,iii) + enddo + else + call sccenter(ires_old,iii,sccor) + endif + iii=0 + endif +C Start new residue. + if (res.eq.'Cl-' .or. res.eq.'Na+') then + ires=ires_old + cycle + else if (ibeg.eq.1) then +c write (iout,*) "BEG ires",ires + ishift=ires-1 + if (res.ne.'GLY' .and. res.ne. 'ACE') then + ishift=ishift-1 + itype(1)=ntyp1 + endif + ires=ires-ishift+ishift1 + ires_old=ires +c write (iout,*) "ishift",ishift," ires",ires, +c & " ires_old",ires_old +c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift + ibeg=0 + else if (ibeg.eq.2) then +c Start a new chain + ishift=-ires_old+ires-1 + ires=ires_old+1 +c write (iout,*) "New chain started",ires,ishift + ibeg=0 + else + ishift=ishift-(ires-ishift+ishift1-ires_old-1) + ires=ires-ishift+ishift1 + ires_old=ires + endif + if (res.eq.'ACE' .or. res.eq.'NHE') then + itype(ires)=10 + else + itype(ires)=rescode(ires,res,0) + endif + else + ires=ires-ishift+ishift1 + endif +c write (iout,*) "ires_old",ires_old," ires",ires +c if (card(27:27).eq."A" .or. card(27:27).eq."B") then +c ishift1=ishift1+1 +c endif +c write (2,*) "ires",ires," res ",res," ity",ity + if (atom.eq.'CA' .or. atom.eq.'CH3' .or. + & res.eq.'NHE'.and.atom(:2).eq.'HN') then + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) +c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3) +#ifdef DEBUG + write (iout,'(2i3,2x,a,3f8.3)') + & ires,itype(ires),res,(c(j,ires),j=1,3) +#endif + iii=iii+1 + do j=1,3 + sccor(j,iii)=c(j,ires) + enddo + if (ishift.ne.0) then + ires_ca=ires+ishift-ishift1 + else + ires_ca=ires + endif +c write (*,*) card(23:27),ires,itype(ires) + else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. + & atom.ne.'N' .and. atom.ne.'C' .and. + & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. + & atom.ne.'OXT' .and. atom(:2).ne.'3H') then +c write (iout,*) "sidechain ",atom + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif + endif + enddo + 10 write (iout,'(a,i5)') ' Nres: ',ires +C Calculate dummy residue coordinates inside the "chain" of a multichain +C system + nres=ires + do i=2,nres-1 +c write (iout,*) i,itype(i),itype(i+1) + if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then + if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then +C 16/01/2014 by Adasko: Adding to dummy atoms in the chain +C first is connected prevous chain (itype(i+1).eq.ntyp1)=true +C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the last dummy residue + call refsys(i-3,i-2,i-1,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif !fail + do j=1,3 + c(j,i)=c(j,i-1)-1.9d0*e2(j) + enddo + else !unres_pdb + do j=1,3 + dcj=(c(j,i-2)-c(j,i-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i-1)+dcj + c(j,nres+i)=c(j,i) + enddo + endif !unres_pdb + else !itype(i+1).eq.ntyp1 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(i+1,i+2,i+3,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + enddo + else !unres_pdb + do j=1,3 + dcj=(c(j,i+3)-c(j,i+2))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i+1)-dcj + c(j,nres+i)=c(j,i) + enddo + endif !unres_pdb + endif !itype(i+1).eq.ntyp1 + endif !itype.eq.ntyp1 + enddo +C Calculate the CM of the last side chain. + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + nsup=nres + nstart_sup=1 + if (itype(nres).ne.10) then + nres=nres+1 + itype(nres)=ntyp1 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the last dummy residue + call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + enddo + else + do j=1,3 + dcj=(c(j,nres-2)-c(j,nres-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,nres)=c(j,nres-1)+dcj + c(j,2*nres)=c(j,nres) + enddo + endif + endif + do i=2,nres-1 + do j=1,3 + c(j,i+nres)=dc(j,i) + enddo + enddo + do j=1,3 + c(j,nres+1)=c(j,1) + c(j,2*nres)=c(j,nres) + enddo + if (itype(1).eq.ntyp1) then + nsup=nsup-1 + nstart_sup=2 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(2,3,4,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0) + enddo + else + do j=1,3 + dcj=(c(j,4)-c(j,3))/2.0 + c(j,1)=c(j,2)-dcj + c(j,nres+1)=c(j,1) + enddo + endif + endif +C Copy the coordinates to reference coordinates +c do i=1,2*nres +c do j=1,3 +c cref(j,i)=c(j,i) +c enddo +c enddo +C Calculate internal coordinates. + if (out_template_coord) then + write (iout,'(/a)') + & "Cartesian coordinates of the reference structure" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(c(j,ires),j=1,3), + & (c(j,ires+nres),j=1,3) + enddo + endif + zero=.false. + do ires=1,nres + zero=zero.or.itype(ires).eq.0 + enddo + if (zero) then + write (iout,'(2a)') "Gaps in PDB coordinates detected;", + & " look for ZERO in the control output above." + write (iout,'(2a)') "Repair the PDB file using MODELLER", + & " or other softwared and resubmit." + call flush(iout) + stop + endif +C Calculate internal coordinates. + call int_from_cart(.true.,out_template_coord) + call sc_loc_geom(.false.) + do i=1,nres + thetaref(i)=theta(i) + phiref(i)=phi(i) + enddo + dc(:,0)=c(:,1) + do i=1,nres-1 + do j=1,3 + dc(j,i)=c(j,i+1)-c(j,i) + dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) + enddo + enddo + do i=2,nres-1 + do j=1,3 + dc(j,i+nres)=c(j,i+nres)-c(j,i) + dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) + enddo +c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3), +c & vbld_inv(i+nres) + enddo + do i=1,nres + do j=1,3 + cref(j,i)=c(j,i) + cref(j,i+nres)=c(j,i+nres) + enddo + enddo + do i=1,2*nres + do j=1,3 + chomo(j,i,k)=c(j,i) + enddo + enddo + + return + end +c--------------------------------------------------------------------------- + subroutine int_from_cart(lside,lprn) + implicit none + include 'DIMENSIONS' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + character*3 seq,atom,res + character*80 card + double precision sccor(3,50) + integer rescode + double precision dist,alpha,beta,di + integer i,j,iti + logical lside,lprn + if (lprn) then + write (iout,'(/a)') + & 'Internal coordinates calculated from crystal structure.' + if (lside) then + write (iout,'(8a)') ' Res ',' dvb',' Theta', + & ' Phi',' Dsc_id',' Dsc',' Alpha', + & ' Omega' + else + write (iout,'(4a)') ' Res ',' dvb',' Theta', + & ' Phi' + endif + endif + do i=2,nres + iti=itype(i) +c write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1) + if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and. + & (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then + write (iout,'(a,i4)') 'Bad Cartesians for residue',i +c stop + endif + vbld(i)=dist(i-1,i) + vbld_inv(i)=1.0d0/vbld(i) + theta(i+1)=alpha(i-1,i,i+1) + if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1) + enddo +c if (itype(1).eq.ntyp1) then +c do j=1,3 +c c(j,1)=c(j,2)+(c(j,3)-c(j,4)) +c enddo +c endif +c if (itype(nres).eq.ntyp1) then +c do j=1,3 +c c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3)) +c enddo +c endif + if (lside) then + do i=2,nres-1 + do j=1,3 + c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1)) + enddo + iti=itype(i) + di=dist(i,nres+i) + vbld(i+nres)=di + if (itype(i).ne.10) then + vbld_inv(i+nres)=1.0d0/di + else + vbld_inv(i+nres)=0.0d0 + endif + if (iti.ne.10) then + alph(i)=alpha(nres+i,i,maxres2) + omeg(i)=beta(nres+i,i,maxres2,i+1) + endif + if (iti.ne.10) then + alph(i)=alpha(nres+i,i,maxres2) + omeg(i)=beta(nres+i,i,maxres2,i+1) + endif + if (lprn) + & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1), + & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di, + & rad2deg*alph(i),rad2deg*omeg(i) + enddo + else if (lprn) then + do i=2,nres + iti=itype(i) + write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1), + & rad2deg*theta(i),rad2deg*phi(i) + enddo + endif + return + end +c--------------------------------------------------------------------------- + subroutine sccenter(ires,nscat,sccor) + implicit none + include 'DIMENSIONS' + include 'COMMON.CHAIN' + integer ires,nscat,i,j + double precision sccor(3,50),sccmj + do j=1,3 + sccmj=0.0D0 + do i=1,nscat + sccmj=sccmj+sccor(j,i) + enddo + dc(j,ires)=sccmj/nscat + enddo + return + end +c--------------------------------------------------------------------------- + subroutine sc_loc_geom(lprn) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + include 'COMMON.CONTROL' + include 'COMMON.SETUP' + double precision x_prime(3),y_prime(3),z_prime(3) + logical lprn + do i=1,nres-1 + do j=1,3 + dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i)) + enddo + enddo + do i=2,nres-1 + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + do j=1,3 + dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i)) + enddo + else + do j=1,3 + dc_norm(j,i+nres)=0.0d0 + enddo + endif + enddo + do i=2,nres-1 + costtab(i+1) =dcos(theta(i+1)) + sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1)) + cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1))) + sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1))) + cosfac2=0.5d0/(1.0d0+costtab(i+1)) + cosfac=dsqrt(cosfac2) + sinfac2=0.5d0/(1.0d0-costtab(i+1)) + sinfac=dsqrt(sinfac2) + it=itype(i) + if (it.ne.10 .and. itype(i).ne.ntyp1) then +c +C Compute the axes of tghe local cartesian coordinates system; store in +c x_prime, y_prime and z_prime +c + do j=1,3 + x_prime(j) = 0.00 + y_prime(j) = 0.00 + z_prime(j) = 0.00 + enddo + do j = 1,3 + x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac + y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac + enddo +c write (iout,*) "x_prime",(x_prime(j),j=1,3) +c write (iout,*) "y_prime",(y_prime(j),j=1,3) + call vecpr(x_prime,y_prime,z_prime) +c write (iout,*) "z_prime",(z_prime(j),j=1,3) +c +C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i), +C to local coordinate system. Store in xx, yy, zz. +c + xx=0.0d0 + yy=0.0d0 + zz=0.0d0 + do j = 1,3 + xx = xx + x_prime(j)*dc_norm(j,i+nres) + yy = yy + y_prime(j)*dc_norm(j,i+nres) + zz = zz + z_prime(j)*dc_norm(j,i+nres) + enddo + + xxref(i)=xx + yyref(i)=yy + zzref(i)=zz + else + xxref(i)=0.0d0 + yyref(i)=0.0d0 + zzref(i)=0.0d0 + endif + enddo + if (lprn) then + write (iout,*) "xxref,yyref,zzref" + do i=2,nres + iti=itype(i) + write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i), + & zzref(i) + enddo + endif + return + end +c--------------------------------------------------------------------------- + subroutine bond_regular + implicit none + include 'DIMENSIONS' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.CHAIN' + integer i,i1,i2 + do i=1,nres-1 + vbld(i+1)=vbl + vbld_inv(i+1)=vblinv + vbld(i+1+nres)=dsc(iabs(itype(i+1))) + vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1))) +c print *,vbld(i+1),vbld(i+1+nres) + enddo +c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain + do i=1,nchain + i1=chain_border(1,i) + i2=chain_border(2,i) + if (i1.gt.1) then + vbld(i1)=vbld(i1)/2 + vbld_inv(i1)=vbld_inv(i1)*2 + endif + if (i2.lt.nres) then + vbld(i2+1)=vbld(i2+1)/2 + vbld_inv(i2+1)=vbld_inv(i2+1)*2 + endif + enddo + return + end + diff --git a/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos b/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos index 853e319..3fcb971 100644 --- a/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos +++ b/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos @@ -9,6 +9,7 @@ OPTE = -c -O3 -ipo -mcmodel=medium -shared-intel -dynamic OPT2 = -O2 -ip -mcmodel=medium -shared-intel -dynamic OPT0 = -g -O0 -mcmodel=medium -shared-intel -dynamic OPT1 = -g -CA -CB -mcmodel=medium -shared-intel -dynamic +#OPT = -g -CA -CB -mcmodel=medium -shared-intel -dynamic FFLAGS = -c ${OPT} FFLAGSE = -c ${OPTE} diff --git a/source/unres/src-HCD-5D/parmread.F b/source/unres/src-HCD-5D/parmread.F index 7550fd5..721d05b 100644 --- a/source/unres/src-HCD-5D/parmread.F +++ b/source/unres/src-HCD-5D/parmread.F @@ -462,6 +462,7 @@ C here will be the apropriate recalibrating for D-aminoacid c write (2,*) "Start reading THETA_PDB",ithep_pdb do i=1,ntyp c write (2,*) 'i=',i + call flush(iout) read (ithep_pdb,*,err=111,end=111) & a0thet(i),(athet(j,i,1,1),j=1,2), & (bthet(j,i,1,1),j=1,2) diff --git a/source/unres/src-HCD-5D/readpdb-mult.F b/source/unres/src-HCD-5D/readpdb-mult.F index 41fe7f6..76cb6b6 100644 --- a/source/unres/src-HCD-5D/readpdb-mult.F +++ b/source/unres/src-HCD-5D/readpdb-mult.F @@ -64,9 +64,9 @@ c call flush(iout) iterter(ires_old-1)=1 itype(ires_old)=ntyp1 iterter(ires_old)=1 - ishift1=ishift1+1 +c ishift1=ishift1+1 ibeg=2 - write (iout,*) "Chain ended",ires,ishift,ires_old + write (iout,*) "Chain ended",ires,ishift,ires_old,ibeg if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) @@ -95,8 +95,8 @@ c write (iout,*) "! ",atom," !",ires read (card(18:20),'(a3)') res c write (iout,*) "ires",ires,ires-ishift+ishift1, c & " ires_old",ires_old -c write (iout,*) "ishift",ishift," ishift1",ishift1 -c write (iout,*) "IRES",ires-ishift+ishift1,ires_old +c write (iout,*) "ishift",ishift," ishift1",ishift1 +c write (iout,*) "IRES",ires-ishift+ishift1,ires_old if (ires-ishift+ishift1.ne.ires_old) then ! Calculate the CM of the preceding residue. ! if (ibeg.eq.0) call sccenter(ires,iii,sccor) @@ -115,6 +115,7 @@ c write (iout,'(i5,3f10.5)') ires,(sccor(j,iii),j=1,3) sccalc=.true. endif ! Start new residue. +c write (iout,*) "ibeg",ibeg if (res.eq.'Cl-' .or. res.eq.'Na+') then ires=ires_old cycle @@ -133,10 +134,13 @@ c write (iout,*) "BEG ires",ires else if (ibeg.eq.2) then ! Start a new chain ishift=-ires_old+ires-1 !!!!! - ishift1=ishift1-1 !!!!! -c write (iout,*) "New chain started",ires,ishift,ishift1,"!" +c ishift1=ishift1-1 !!!!! +c write (iout,*) "New chain started",ires,ires_old,ishift, +c & ishift1 ires=ires-ishift+ishift1 + write (iout,*) "New chain started ires",ires ires_old=ires +c ires=ires_old+1 ibeg=0 else ishift=ishift-(ires-ishift+ishift1-ires_old-1) @@ -159,7 +163,8 @@ c write (2,*) "ires",ires," res ",res!," ity"!,ity if (atom.eq.'CA' .or. atom.eq.'CH3' .or. & res.eq.'NHE'.and.atom(:2).eq.'HN') then read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) -! write (iout,*) "backbone ",atom +c write (iout,*) "backbone ",atom +c write (iout,*) ires,res,(c(j,ires),j=1,3) #ifdef DEBUG write (iout,'(i6,i3,2x,a,3f8.3)') & ires,itype(ires),res,(c(j,ires),j=1,3) @@ -229,7 +234,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,i)=c(j,i+1)-1.9d0*e2(j) + c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) enddo else !unres_pdb do j=1,3 @@ -412,8 +417,9 @@ C and convert the peptide geometry into virtual-chain geometry. character*3 seq,res character*5 atom character*80 card - double precision sccor(3,20) + double precision sccor(3,50) integer rescode,iterter(maxres) + logical zero do i=1,maxres iterter(i)=0 enddo @@ -581,7 +587,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,i)=c(j,i+1)-1.9d0*e2(j) + c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) enddo else !unres_pdb do j=1,3 @@ -616,7 +622,7 @@ C 2/15/2013 by Adam: corrected insertion of the last dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,nres)=c(j,nres-1)-1.9d0*e2(j) + c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) enddo else do j=1,3 @@ -648,7 +654,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,1)=c(j,2)-1.9d0*e2(j) + c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0) enddo else do j=1,3 @@ -676,6 +682,18 @@ C Calculate internal coordinates. & (c(j,ires+nres),j=1,3) enddo endif + zero=.false. + do ires=1,nres + zero=zero.or.itype(ires).eq.0 + enddo + if (zero) then + write (iout,'(2a)') "Gaps in PDB coordinates detected;", + & " look for ZERO in the control output above." + write (iout,'(2a)') "Repair the PDB file using MODELLER", + & " or other softwared and resubmit." + call flush(iout) + stop + endif C Calculate internal coordinates. call int_from_cart(.true.,out_template_coord) call sc_loc_geom(.false.) @@ -683,6 +701,7 @@ C Calculate internal coordinates. thetaref(i)=theta(i) phiref(i)=phi(i) enddo + dc(:,0)=c(:,1) do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index e5f0b41..eeaf74c 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -1184,7 +1184,8 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb 335 continue unres_pdb=.false. nres_temp=nres - call readpdb +c call readpdb + call readpdb_template(nmodel_start+1) close(ipdbin) if (nres.ge.nres_temp) then nmodel_start=nmodel_start+1 diff --git a/source/wham/src-HCD-5D/CMakeLists.txt b/source/wham/src-HCD-5D/CMakeLists.txt index d7bdfcd..42d52e9 100644 --- a/source/wham/src-HCD-5D/CMakeLists.txt +++ b/source/wham/src-HCD-5D/CMakeLists.txt @@ -48,7 +48,7 @@ set(UNRES_WHAM_M_SRC0 PMFprocess.F oligomer.F readrtns_compar.F - readpdb.F + readpdb-mult.F fitsq.f contact.f elecont.f @@ -90,7 +90,7 @@ set(UNRES_WHAM_M_PP_SRC PMFprocess.F read_constr_homology.F read_dist_constr.F - readpdb.F + readpdb-mult.F read_ref_str.F readrtns_compar.F readrtns.F diff --git a/source/wham/src-HCD-5D/Makefile_MPICH_ifort-okeanos b/source/wham/src-HCD-5D/Makefile_MPICH_ifort-okeanos index b04295c..c5e3133 100644 --- a/source/wham/src-HCD-5D/Makefile_MPICH_ifort-okeanos +++ b/source/wham/src-HCD-5D/Makefile_MPICH_ifort-okeanos @@ -1,10 +1,12 @@ BIN = ~/bin FC = ftn -OPT = -mcmodel=medium -shared-intel -O3 -dynamic +OPT = -mcmodel=medium -shared-intel -O2 -dynamic +OPTE = -mcmodel=medium -shared-intel -O3 -dynamic #OPT = -O3 -intel-static -mcmodel=medium #OPT = -O3 -ip -w #OPT = -g -CB -mcmodel=medium -shared-intel -dynamic FFLAGS = ${OPT} -c -I. -I./include_unres -I$(INSTALL_DIR)/include +FFLAGSE = ${OPTE} -c -I. -I./include_unres -I$(INSTALL_DIR)/include LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a .f.o: @@ -59,7 +61,7 @@ objects = \ objects_compar = \ readrtns_compar.o \ - readpdb.o fitsq.o contact.o \ + readpdb-mult.o fitsq.o contact.o \ elecont.o contfunc.o cont_frag.o conf_compar.o match_contact.o \ angnorm.o odlodc.o promienie.o qwolynes.o read_ref_str.o \ rmscalc.o secondary.o proc_cont.o define_pairs.o mysort.o @@ -156,6 +158,11 @@ NEWCORR5D_DFA: ${objects} ${objects_compar} dfa.o xdrf/libxdrf.a xdrf/libxdrf.a: cd xdrf && make +energy_p_new.o: energy_p_new.F + ${FC} ${FFLAGSE} ${CPPFLAGS} energy_p_new.F + +wham_calc1.o: wham_calc1.F + ${FC} ${FFLAGSE} ${CPPFLAGS} wham_calc1.F clean: /bin/rm -f *.o && /bin/rm -f compinfo && cd xdrf && make clean diff --git a/source/wham/src-HCD-5D/readpdb-mult.F b/source/wham/src-HCD-5D/readpdb-mult.F new file mode 100644 index 0000000..37b15c1 --- /dev/null +++ b/source/wham/src-HCD-5D/readpdb-mult.F @@ -0,0 +1,960 @@ + subroutine readpdb +C Read the PDB file and convert the peptide geometry into virtual-chain +C geometry. + implicit none + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.CONTROL' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + include 'COMMON.SBRIDGE' + include 'COMMON.FRAG' + character*3 seq,atom,res + character*80 card + double precision e1(3),e2(3),e3(3) + double precision sccor(3,50) + integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old + double precision dcj + integer rescode,kkk,lll,icha,cou,kupa,iprzes + logical lsecondary,sccalc,fail,zero + integer iterter(maxres) + double precision efree_temp + iii=0 + ibeg=1 + ishift1=0 + sccalc=.false. + bfac=0.0d0 + do i=1,maxres + iterter(i)=0 + enddo + ibeg=1 + ishift1=0 + lsecondary=.false. + nhfrag=0 + nbfrag=0 + iii=0 + sccalc=.false. + do + read (ipdbin,'(a80)',end=10) card +c write (iout,'(a)') card +c call flush(iout) + if (card(:5).eq.'HELIX') then + nhfrag=nhfrag+1 + lsecondary=.true. + read(card(22:25),*) hfrag(1,nhfrag) + read(card(34:37),*) hfrag(2,nhfrag) + endif + if (card(:5).eq.'SHEET') then + nbfrag=nbfrag+1 + lsecondary=.true. + read(card(24:26),*) bfrag(1,nbfrag) + read(card(35:37),*) bfrag(2,nbfrag) +!rc---------------------------------------- +!rc to be corrected !!! + bfrag(3,nbfrag)=bfrag(1,nbfrag) + bfrag(4,nbfrag)=bfrag(2,nbfrag) +!rc---------------------------------------- + endif + if (card(:3).eq.'END') then + goto 10 + else if (card(:3).eq.'TER') then +! End current chain + ires_old=ires+2 + itype(ires_old-1)=ntyp1 + iterter(ires_old-1)=1 + itype(ires_old)=ntyp1 + iterter(ires_old)=1 +c ishift1=ishift1+1 + ibeg=2 + write (iout,*) "Chain ended",ires,ishift,ires_old,ibeg + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + iii=0 + sccalc=.true. + endif +! Read free energy +c if (index(card,"FREE ENERGY").gt.0) then +c ifree=index(card,"FREE ENERGY")+12 +c read(card(ifree:),*,err=1115,end=1115) efree_temp +c 1115 continue +c endif +! Fish out the ATOM cards. + if (index(card(1:4),'ATOM').gt.0) then + sccalc=.false. + read (card(12:16),*) atom +c write (2,'(a)') card +c write (iout,*) "ibeg",ibeg +c write (iout,*) "! ",atom," !",ires +! if (atom.eq.'CA' .or. atom.eq.'CH3') then + read (card(23:26),*) ires + read (card(18:20),'(a3)') res +c write (iout,*) "ires",ires,ires-ishift+ishift1, +c & " ires_old",ires_old +c write (iout,*) "ishift",ishift," ishift1",ishift1 +c write (iout,*) "IRES",ires-ishift+ishift1,ires_old + if (ires-ishift+ishift1.ne.ires_old) then +! Calculate the CM of the preceding residue. +! if (ibeg.eq.0) call sccenter(ires,iii,sccor) + if (ibeg.eq.0) then +c write (iout,*) "Calculating sidechain center iii",iii +c write (iout,*) "ires",ires + if (unres_pdb) then +c write (iout,'(i5,3f10.5)') ires,(sccor(j,iii),j=1,3) + do j=1,3 + dc(j,ires_old)=sccor(j,iii) + enddo + else + call sccenter(ires_old,iii,sccor) + endif + iii=0 + sccalc=.true. + endif +! Start new residue. +c write (iout,*) "ibeg",ibeg + if (res.eq.'Cl-' .or. res.eq.'Na+') then + ires=ires_old + cycle + else if (ibeg.eq.1) then +c write (iout,*) "BEG ires",ires + ishift=ires-1 + if (res.ne.'GLY' .and. res.ne. 'ACE') then + ishift=ishift-1 + itype(1)=ntyp1 + endif + ires=ires-ishift+ishift1 + ires_old=ires +! write (iout,*) "ishift",ishift," ires",ires,& +! " ires_old",ires_old + ibeg=0 + else if (ibeg.eq.2) then +! Start a new chain + ishift=-ires_old+ires-1 !!!!! +c ishift1=ishift1-1 !!!!! +c write (iout,*) "New chain started",ires,ires_old,ishift, +c & ishift1 + ires=ires-ishift+ishift1 + write (iout,*) "New chain started ires",ires + ires_old=ires +c ires=ires_old+1 + ibeg=0 + else + ishift=ishift-(ires-ishift+ishift1-ires_old-1) + ires=ires-ishift+ishift1 + ires_old=ires + endif + if (res.eq.'ACE' .or. res.eq.'NHE') then + itype(ires)=10 + else + itype(ires)=rescode(ires,res,0) + endif + else + ires=ires-ishift+ishift1 + endif +c write (iout,*) "ires_old",ires_old," ires",ires + if (card(27:27).eq."A" .or. card(27:27).eq."B") then +! ishift1=ishift1+1 + endif +c write (2,*) "ires",ires," res ",res!," ity"!,ity + if (atom.eq.'CA' .or. atom.eq.'CH3' .or. + & res.eq.'NHE'.and.atom(:2).eq.'HN') then + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) +c write (iout,*) "backbone ",atom +c write (iout,*) ires,res,(c(j,ires),j=1,3) +#ifdef DEBUG + write (iout,'(i6,i3,2x,a,3f8.3)') + & ires,itype(ires),res,(c(j,ires),j=1,3) +#endif + iii=iii+1 + do j=1,3 + sccor(j,iii)=c(j,ires) + enddo +c write (2,*) card(23:27),ires,itype(ires),iii + else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. + & atom.ne.'N' .and. atom.ne.'C' .and. + & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. + & atom.ne.'OXT' .and. atom(:2).ne.'3H') then +! write (iout,*) "sidechain ",atom + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) +c write (2,*) "iii",iii + endif + endif + enddo + 10 write (iout,'(a,i5)') ' Nres: ',ires +c write (iout,*) "iii",iii +C Calculate dummy residue coordinates inside the "chain" of a multichain +C system + nres=ires +c write (iout,*) "dc" +c do i=1,nres +c write (iout,'(i5,3f10.5)') i,(dc(j,i),j=1,3) +c enddo + do i=2,nres-1 +c write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i) + if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then + if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then +C 16/01/2014 by Adasko: Adding to dummy atoms in the chain +C first is connected prevous chain (itype(i+1).eq.ntyp1)=true +C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the last dummy residue +c print *,i,'tu dochodze' + call refsys(i-3,i-2,i-1,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif !fail +c print *,i,'a tu?' + do j=1,3 + c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + enddo + else !unres_pdb + do j=1,3 + dcj=(c(j,i-2)-c(j,i-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i-1)+dcj + c(j,nres+i)=c(j,i) + dC(j,i)=c(j,i) + enddo + endif !unres_pdb + else !itype(i+1).eq.ntyp1 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(i+1,i+2,i+3,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + enddo + else !unres_pdb + do j=1,3 + dcj=(c(j,i+3)-c(j,i+2))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i+1)-dcj + c(j,nres+i)=c(j,i) + dC(j,i)=c(j,i) + enddo + endif !unres_pdb + endif !itype(i+1).eq.ntyp1 + endif !itype.eq.ntyp1 + enddo + write (iout,*) "After loop in readpbd" +C Calculate the CM of the last side chain. + if (.not. sccalc) then + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else +c write (iout,*) "Calling sccenter iii",iii + call sccenter(ires,iii,sccor) + endif + endif + nsup=nres + nstart_sup=1 + if (itype(nres).ne.10) then + nres=nres+1 + itype(nres)=ntyp1 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the last dummy residue + call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + enddo + else + do j=1,3 + dcj=(c(j,nres-2)-c(j,nres-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,nres)=c(j,nres-1)+dcj + c(j,2*nres)=c(j,nres) + enddo + endif + endif + do i=2,nres-1 + do j=1,3 + c(j,i+nres)=dc(j,i) + enddo + enddo + do j=1,3 + c(j,nres+1)=c(j,1) + c(j,2*nres)=c(j,nres) + enddo + if (itype(1).eq.ntyp1) then + nsup=nsup-1 + nstart_sup=2 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(2,3,4,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0) + enddo + else + do j=1,3 + dcj=(c(j,4)-c(j,3))/2.0 + c(j,1)=c(j,2)-dcj + c(j,nres+1)=c(j,1) + enddo + endif + endif +C Calculate internal coordinates. + write (iout,'(/a)') + & "Cartesian coordinates of the reference structure" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(c(j,ires),j=1,3), + & (c(j,ires+nres),j=1,3) + zero=.false. + enddo + do ires=1,nres + zero=zero.or.itype(ires).eq.0 + enddo + if (zero) then + write (iout,'(2a)') "Gaps in PDB coordinates detected;", + & " look for ZERO in the control output above." + write (iout,'(2a)') "Repair the PDB file using MODELLER", + & " or other softwared and resubmit." + call flush(iout) + stop + endif +c write(iout,*)"before int_from_cart nres",nres + call int_from_cart(.true.,.false.) + do i=1,nres + thetaref(i)=theta(i) + phiref(i)=phi(i) + enddo + dc(:,0)=c(:,1) + do i=1,nres-1 + do j=1,3 + dc(j,i)=c(j,i+1)-c(j,i) + dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) + enddo +c write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3), +c & vbld_inv(i+1) + enddo + do i=2,nres-1 + do j=1,3 + dc(j,i+nres)=c(j,i+nres)-c(j,i) + dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) + enddo +c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3), +c & vbld_inv(i+nres) + enddo + call sc_loc_geom(.false.) + call int_from_cart1(.false.) +c call chainbuild +C Copy the coordinates to reference coordinates + do i=1,nres + do j=1,3 + cref(j,i)=c(j,i) + cref(j,i+nres)=c(j,i+nres) + enddo + enddo + 100 format (//' alpha-carbon coordinates ', + & ' centroid coordinates'/ + 1 ' ', 7X,'X',11X,'Y',11X,'Z', + & 10X,'X',11X,'Y',11X,'Z') + 110 format (a,'(',i4,')',6f12.5) +cc enddiag + do j=1,nbfrag + do i=1,4 + bfrag(i,j)=bfrag(i,j)-ishift + enddo + enddo + + do j=1,nhfrag + do i=1,2 + hfrag(i,j)=hfrag(i,j)-ishift + enddo + enddo + return + end +c--------------------------------------------------------------------------- + subroutine readpdb_template(k) +C Read the PDB file for read_constr_homology with read2sigma +C and convert the peptide geometry into virtual-chain geometry. + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + include 'COMMON.CONTROL' + include 'COMMON.SETUP' + integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity + logical lprn /.false./,fail + double precision e1(3),e2(3),e3(3) + double precision dcj,efree_temp + character*3 seq,res + character*5 atom + character*80 card + double precision sccor(3,20) + integer rescode,iterter(maxres) + logical zero + do i=1,maxres + iterter(i)=0 + enddo + ibeg=1 + ishift1=0 + ishift=0 +c write (2,*) "UNRES_PDB",unres_pdb + ires=0 + ires_old=0 + iii=0 + lsecondary=.false. + nhfrag=0 + nbfrag=0 + do + read (ipdbin,'(a80)',end=10) card + if (card(:3).eq.'END') then + goto 10 + else if (card(:3).eq.'TER') then +C End current chain + ires_old=ires+2 + itype(ires_old-1)=ntyp1 + iterter(ires_old-1)=1 + itype(ires_old)=ntyp1 + iterter(ires_old)=1 + ibeg=2 +c write (iout,*) "Chain ended",ires,ishift,ires_old + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + endif +C Fish out the ATOM cards. + if (index(card(1:4),'ATOM').gt.0) then + read (card(12:16),*) atom +c write (iout,*) "! ",atom," !",ires +c if (atom.eq.'CA' .or. atom.eq.'CH3') then + read (card(23:26),*) ires + read (card(18:20),'(a3)') res +c write (iout,*) "ires",ires,ires-ishift+ishift1, +c & " ires_old",ires_old +c write (iout,*) "ishift",ishift," ishift1",ishift1 +c write (iout,*) "IRES",ires-ishift+ishift1,ires_old + if (ires-ishift+ishift1.ne.ires_old) then +C Calculate the CM of the preceding residue. + if (ibeg.eq.0) then + if (unres_pdb) then + do j=1,3 + dc(j,ires_old)=sccor(j,iii) + enddo + else + call sccenter(ires_old,iii,sccor) + endif + iii=0 + endif +C Start new residue. + if (res.eq.'Cl-' .or. res.eq.'Na+') then + ires=ires_old + cycle + else if (ibeg.eq.1) then +c write (iout,*) "BEG ires",ires + ishift=ires-1 + if (res.ne.'GLY' .and. res.ne. 'ACE') then + ishift=ishift-1 + itype(1)=ntyp1 + endif + ires=ires-ishift+ishift1 + ires_old=ires +c write (iout,*) "ishift",ishift," ires",ires, +c & " ires_old",ires_old +c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift + ibeg=0 + else if (ibeg.eq.2) then +c Start a new chain + ishift=-ires_old+ires-1 + ires=ires_old+1 +c write (iout,*) "New chain started",ires,ishift + ibeg=0 + else + ishift=ishift-(ires-ishift+ishift1-ires_old-1) + ires=ires-ishift+ishift1 + ires_old=ires + endif + if (res.eq.'ACE' .or. res.eq.'NHE') then + itype(ires)=10 + else + itype(ires)=rescode(ires,res,0) + endif + else + ires=ires-ishift+ishift1 + endif +c write (iout,*) "ires_old",ires_old," ires",ires +c if (card(27:27).eq."A" .or. card(27:27).eq."B") then +c ishift1=ishift1+1 +c endif +c write (2,*) "ires",ires," res ",res," ity",ity + if (atom.eq.'CA' .or. atom.eq.'CH3' .or. + & res.eq.'NHE'.and.atom(:2).eq.'HN') then + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) +c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3) +#ifdef DEBUG + write (iout,'(2i3,2x,a,3f8.3)') + & ires,itype(ires),res,(c(j,ires),j=1,3) +#endif + iii=iii+1 + do j=1,3 + sccor(j,iii)=c(j,ires) + enddo + if (ishift.ne.0) then + ires_ca=ires+ishift-ishift1 + else + ires_ca=ires + endif +c write (*,*) card(23:27),ires,itype(ires) + else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. + & atom.ne.'N' .and. atom.ne.'C' .and. + & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. + & atom.ne.'OXT' .and. atom(:2).ne.'3H') then +c write (iout,*) "sidechain ",atom + iii=iii+1 + read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) + endif + endif + enddo + 10 write (iout,'(a,i5)') ' Nres: ',ires +C Calculate dummy residue coordinates inside the "chain" of a multichain +C system + nres=ires + do i=2,nres-1 +c write (iout,*) i,itype(i),itype(i+1) + if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then + if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then +C 16/01/2014 by Adasko: Adding to dummy atoms in the chain +C first is connected prevous chain (itype(i+1).eq.ntyp1)=true +C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the last dummy residue + call refsys(i-3,i-2,i-1,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif !fail + do j=1,3 + c(j,i)=c(j,i-1)-1.9d0*e2(j) + enddo + else !unres_pdb + do j=1,3 + dcj=(c(j,i-2)-c(j,i-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i-1)+dcj + c(j,nres+i)=c(j,i) + enddo + endif !unres_pdb + else !itype(i+1).eq.ntyp1 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(i+1,i+2,i+3,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + enddo + else !unres_pdb + do j=1,3 + dcj=(c(j,i+3)-c(j,i+2))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,i)=c(j,i+1)-dcj + c(j,nres+i)=c(j,i) + enddo + endif !unres_pdb + endif !itype(i+1).eq.ntyp1 + endif !itype.eq.ntyp1 + enddo +C Calculate the CM of the last side chain. + if (unres_pdb) then + do j=1,3 + dc(j,ires)=sccor(j,iii) + enddo + else + call sccenter(ires,iii,sccor) + endif + nsup=nres + nstart_sup=1 + if (itype(nres).ne.10) then + nres=nres+1 + itype(nres)=ntyp1 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the last dummy residue + call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) + enddo + else + do j=1,3 + dcj=(c(j,nres-2)-c(j,nres-3))/2.0 + if (dcj.eq.0) dcj=1.23591524223 + c(j,nres)=c(j,nres-1)+dcj + c(j,2*nres)=c(j,nres) + enddo + endif + endif + do i=2,nres-1 + do j=1,3 + c(j,i+nres)=dc(j,i) + enddo + enddo + do j=1,3 + c(j,nres+1)=c(j,1) + c(j,2*nres)=c(j,nres) + enddo + if (itype(1).eq.ntyp1) then + nsup=nsup-1 + nstart_sup=2 + if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue + call refsys(2,3,4,e1,e2,e3,fail) + if (fail) then + e2(1)=0.0d0 + e2(2)=1.0d0 + e2(3)=0.0d0 + endif + do j=1,3 + c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0) + enddo + else + do j=1,3 + dcj=(c(j,4)-c(j,3))/2.0 + c(j,1)=c(j,2)-dcj + c(j,nres+1)=c(j,1) + enddo + endif + endif +C Copy the coordinates to reference coordinates +c do i=1,2*nres +c do j=1,3 +c cref(j,i)=c(j,i) +c enddo +c enddo +C Calculate internal coordinates. + if (out_template_coord) then + write (iout,'(/a)') + & "Cartesian coordinates of the reference structure" + write (iout,'(a,3(3x,a5),5x,3(3x,a5))') + & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + do ires=1,nres + write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)') + & restyp(itype(ires)),ires,(c(j,ires),j=1,3), + & (c(j,ires+nres),j=1,3) + enddo + endif + zero=.false. + do ires=1,nres + zero=zero.or.itype(ires).eq.0 + enddo + if (zero) then + write (iout,'(2a)') "Gaps in PDB coordinates detected;", + & " look for ZERO in the control output above." + write (iout,'(2a)') "Repair the PDB file using MODELLER", + & " or other softwared and resubmit." + call flush(iout) + stop + endif +C Calculate internal coordinates. + call int_from_cart(.true.,out_template_coord) + call sc_loc_geom(.false.) + do i=1,nres + thetaref(i)=theta(i) + phiref(i)=phi(i) + enddo + dc(:,0)=c(:,1) + do i=1,nres-1 + do j=1,3 + dc(j,i)=c(j,i+1)-c(j,i) + dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) + enddo + enddo + do i=2,nres-1 + do j=1,3 + dc(j,i+nres)=c(j,i+nres)-c(j,i) + dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) + enddo +c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3), +c & vbld_inv(i+nres) + enddo + do i=1,nres + do j=1,3 + cref(j,i)=c(j,i) + cref(j,i+nres)=c(j,i+nres) + enddo + enddo + do i=1,2*nres + do j=1,3 + chomo(j,i,k)=c(j,i) + enddo + enddo + + return + end +c--------------------------------------------------------------------------- + subroutine int_from_cart(lside,lprn) + implicit none + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + character*3 seq,atom,res + character*80 card + double precision sccor(3,50) + integer rescode + double precision dist,alpha,beta,di + integer i,j,iti + logical lside,lprn + if (lprn) then + write (iout,'(/a)') + & 'Internal coordinates calculated from crystal structure.' + if (lside) then + write (iout,'(8a)') ' Res ',' dvb',' Theta', + & ' Phi',' Dsc_id',' Dsc',' Alpha', + & ' Omega' + else + write (iout,'(4a)') ' Res ',' dvb',' Theta', + & ' Phi' + endif + endif + do i=2,nres + iti=itype(i) +c write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1) + if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and. + & (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then + write (iout,'(a,i4)') 'Bad Cartesians for residue',i +c stop + endif + vbld(i)=dist(i-1,i) + vbld_inv(i)=1.0d0/vbld(i) + theta(i+1)=alpha(i-1,i,i+1) + if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1) + enddo +c if (itype(1).eq.ntyp1) then +c do j=1,3 +c c(j,1)=c(j,2)+(c(j,3)-c(j,4)) +c enddo +c endif +c if (itype(nres).eq.ntyp1) then +c do j=1,3 +c c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3)) +c enddo +c endif + if (lside) then + do i=2,nres-1 + do j=1,3 + c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1)) + enddo + iti=itype(i) + di=dist(i,nres+i) + vbld(i+nres)=di + if (itype(i).ne.10) then + vbld_inv(i+nres)=1.0d0/di + else + vbld_inv(i+nres)=0.0d0 + endif + if (iti.ne.10) then + alph(i)=alpha(nres+i,i,maxres2) + omeg(i)=beta(nres+i,i,maxres2,i+1) + endif + if (iti.ne.10) then + alph(i)=alpha(nres+i,i,maxres2) + omeg(i)=beta(nres+i,i,maxres2,i+1) + endif + if (lprn) + & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1), + & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di, + & rad2deg*alph(i),rad2deg*omeg(i) + enddo + else if (lprn) then + do i=2,nres + iti=itype(i) + write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1), + & rad2deg*theta(i),rad2deg*phi(i) + enddo + endif + return + end +c--------------------------------------------------------------------------- + subroutine sccenter(ires,nscat,sccor) + implicit none + include 'DIMENSIONS' + include 'COMMON.CHAIN' + integer ires,nscat,i,j + double precision sccor(3,50),sccmj + do j=1,3 + sccmj=0.0D0 + do i=1,nscat + sccmj=sccmj+sccor(j,i) + enddo + dc(j,ires)=sccmj/nscat + enddo + return + end +c--------------------------------------------------------------------------- + subroutine sc_loc_geom(lprn) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.NAMES' + include 'COMMON.CONTROL' + include 'COMMON.SETUP' + double precision x_prime(3),y_prime(3),z_prime(3) + logical lprn + do i=1,nres-1 + do j=1,3 + dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i)) + enddo + enddo + do i=2,nres-1 + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + do j=1,3 + dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i)) + enddo + else + do j=1,3 + dc_norm(j,i+nres)=0.0d0 + enddo + endif + enddo + do i=2,nres-1 + costtab(i+1) =dcos(theta(i+1)) + sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1)) + cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1))) + sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1))) + cosfac2=0.5d0/(1.0d0+costtab(i+1)) + cosfac=dsqrt(cosfac2) + sinfac2=0.5d0/(1.0d0-costtab(i+1)) + sinfac=dsqrt(sinfac2) + it=itype(i) + if (it.ne.10 .and. itype(i).ne.ntyp1) then +c +C Compute the axes of tghe local cartesian coordinates system; store in +c x_prime, y_prime and z_prime +c + do j=1,3 + x_prime(j) = 0.00 + y_prime(j) = 0.00 + z_prime(j) = 0.00 + enddo + do j = 1,3 + x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac + y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac + enddo +c write (iout,*) "x_prime",(x_prime(j),j=1,3) +c write (iout,*) "y_prime",(y_prime(j),j=1,3) + call vecpr(x_prime,y_prime,z_prime) +c write (iout,*) "z_prime",(z_prime(j),j=1,3) +c +C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i), +C to local coordinate system. Store in xx, yy, zz. +c + xx=0.0d0 + yy=0.0d0 + zz=0.0d0 + do j = 1,3 + xx = xx + x_prime(j)*dc_norm(j,i+nres) + yy = yy + y_prime(j)*dc_norm(j,i+nres) + zz = zz + z_prime(j)*dc_norm(j,i+nres) + enddo + + xxref(i)=xx + yyref(i)=yy + zzref(i)=zz + else + xxref(i)=0.0d0 + yyref(i)=0.0d0 + zzref(i)=0.0d0 + endif + enddo + if (lprn) then + write (iout,*) "xxref,yyref,zzref" + do i=2,nres + iti=itype(i) + write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i), + & zzref(i) + enddo + endif + return + end +c--------------------------------------------------------------------------- + subroutine bond_regular + implicit none + include 'DIMENSIONS' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.CHAIN' + integer i,i1,i2 + do i=1,nres-1 + vbld(i+1)=vbl + vbld_inv(i+1)=vblinv + vbld(i+1+nres)=dsc(iabs(itype(i+1))) + vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1))) +c print *,vbld(i+1),vbld(i+1+nres) + enddo +c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain + do i=1,nchain + i1=chain_border(1,i) + i2=chain_border(2,i) + if (i1.gt.1) then + vbld(i1)=vbld(i1)/2 + vbld_inv(i1)=vbld_inv(i1)*2 + endif + if (i2.lt.nres) then + vbld(i2+1)=vbld(i2+1)/2 + vbld_inv(i2+1)=vbld_inv(i2+1)*2 + endif + enddo + return + end + -- 1.7.9.5