From 77b3def22adb62b65dbc315c01267bf4b77e6fad Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Thu, 9 Apr 2020 20:57:43 +0200 Subject: [PATCH] Adam's corrections --- source/cluster/wham/src-HCD/readpdb.F | 251 ++++++++---- .../unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos | 2 +- source/unres/src-HCD-5D/int_from_cart.F | 2 + source/unres/src-HCD-5D/readpdb-mult.F | 27 +- source/unres/src-HCD-5D/readpdb.F | 421 ++++++++++++-------- source/unres/src-HCD-5D/sc_move.F | 2 +- source/unres/src-HCD-5D/unres.F | 1 + source/wham/src-HCD/energy_p_new.F | 5 + source/wham/src-HCD/readpdb.F | 4 +- 9 files changed, 459 insertions(+), 256 deletions(-) diff --git a/source/cluster/wham/src-HCD/readpdb.F b/source/cluster/wham/src-HCD/readpdb.F index 98e538e..09e6e4e 100644 --- a/source/cluster/wham/src-HCD/readpdb.F +++ b/source/cluster/wham/src-HCD/readpdb.F @@ -19,7 +19,7 @@ C geometry. character*3 seq,res character*5 atom character*80 card - double precision sccor(3,20) + double precision sccor(3,50) integer rescode integer iterter(maxres) efree_temp=0.0d0 @@ -61,6 +61,7 @@ c write (2,*) "UNRES_PDB",unres_pdb itype(ires_old-1)=ntyp1 iterter(ires_old-1)=1 itype(ires_old)=ntyp1 + iterter(ires_old)=1 ishift1=ishift1+1 ibeg=2 ! write (iout,*) "Chain ended",ires,ishift,ires_old @@ -71,7 +72,7 @@ c write (2,*) "UNRES_PDB",unres_pdb else call sccenter(ires,iii,sccor) endif -c iii=0 + iii=0 endif ! Read free energy if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp @@ -167,45 +168,70 @@ c write (2,*) "iii",iii endif endif enddo - 10 continue -#ifdef DEBUG - write (iout,'(a,i5)') ' Number of residues found: ',ires -#endif - if (ires.eq.0) return -C Calculate the CM of the last side chain. - if (iii.gt.0) then - if (unres_pdb) then - do j=1,3 - dc(j,ires)=sccor(j,iii) - enddo - else - call sccenter(ires,iii,sccor) - endif - endif + 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) + + if (itype(i).eq.ntyp1) then + if (itype(i+1).eq.ntyp1) 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 +C if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the last dummy residue +C call refsys(i-3,i-2,i-1,e1,e2,e3,fail) +C if (fail) then +C e2(1)=0.0d0 +C e2(2)=1.0d0 +C e2(3)=0.0d0 +C endif !fail +C do j=1,3 +C c(j,i)=c(j,i-1)-1.9d0*e2(j) +C enddo +C else !unres_pdb + do j=1,3 + dcj=(c(j,i-2)-c(j,i-3))/2.0 + c(j,i)=c(j,i-1)+dcj + c(j,nres+i)=c(j,i) + enddo +C endif !unres_pdb + else !itype(i+1).eq.ntyp1 +C if (unres_pdb) then +C 2/15/2013 by Adam: corrected insertion of the first dummy residue +C call refsys(i+1,i+2,i+3,e1,e2,e3,fail) +C if (fail) then +C e2(1)=0.0d0 +C e2(2)=1.0d0 +C e2(3)=0.0d0 +C endif +C do j=1,3 +C c(j,i)=c(j,i+1)-1.9d0*e2(j) +C enddo +C else !unres_pdb + do j=1,3 + dcj=(c(j,i+3)-c(j,i+2))/2.0 + c(j,i)=c(j,i+1)-dcj + c(j,nres+i)=c(j,i) + enddo +C endif !unres_pdb + endif !itype(i+1).eq.ntyp1 + endif !itype.eq.ntyp1 + enddo +C Calculate the CM of the last side chain. + call sccenter(ires,iii,sccor) 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)-3.8d0*e2(j) - enddo - else do j=1,3 - dcj=c(j,nres-2)-c(j,nres-3) + dcj=(c(j,nres-2)-c(j,nres-3))/2.0 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 @@ -219,31 +245,12 @@ C 2/15/2013 by Adam: corrected insertion of the last dummy residue 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)-3.8d0*e2(j) - enddo - else do j=1,3 - dcj=c(j,4)-c(j,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 (lprn) then write (iout,'(/a)') @@ -321,7 +328,7 @@ c--------------------------------------------------------------------------- character*3 seq,res c character*5 atom character*80 card - dimension sccor(3,20) + dimension sccor(3,50) integer rescode logical lside,lprn if (lprn) then @@ -487,7 +494,7 @@ c--------------------------------------------------------------------------- implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' - dimension sccor(3,20) + dimension sccor(3,50) do j=1,3 sccmj=0.0D0 do i=1,nscat @@ -549,9 +556,11 @@ 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) - integer rescode - efree_temp=0.0d0 + double precision sccor(3,50) + integer rescode,iterter(maxres) + do i=1,maxres + iterter(i)=0 + enddo ibeg=1 ishift1=0 ishift=0 @@ -562,10 +571,27 @@ c write (2,*) "UNRES_PDB",unres_pdb lsecondary=.false. nhfrag=0 nbfrag=0 - do + do read (ipdbin,'(a80)',end=10) card -c write (iout,'(a)') card - if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10 + 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 @@ -579,9 +605,7 @@ 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. -c if (ibeg.eq.0) call sccenter(ires,iii,sccor) if (ibeg.eq.0) then -c write (iout,*) "Calculating sidechain center iii",iii if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) @@ -606,6 +630,13 @@ c write (iout,*) "BEG ires",ires 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) @@ -621,14 +652,14 @@ c & " ires_old",ires_old 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 +c if (card(27:27).eq."A" .or. card(27:27).eq."B") then c ishift1=ishift1+1 - endif +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 +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) @@ -653,13 +684,60 @@ c write (iout,*) "sidechain ",atom endif endif enddo - 10 continue -#ifdef DEBUG - write (iout,'(a,i5)') ' Number of residues found: ',ires -#endif - if (ires.eq.0) return + 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*e2(j) + 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 (iii.gt.0) then if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) @@ -667,19 +745,31 @@ C Calculate the CM of the last side chain. else call sccenter(ires,iii,sccor) endif - endif - nres=ires 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*e2(j) + enddo + else do j=1,3 - dcj=c(j,nres-2)-c(j,nres-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) @@ -701,11 +791,11 @@ 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)-3.8d0*e2(j) + c(j,1)=c(j,2)-1.9d0*e2(j) enddo else do j=1,3 - dcj=c(j,4)-c(j,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 @@ -751,16 +841,17 @@ C Calculate internal coordinates. 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 -c call chainbuild -C Copy the coordinates to reference coordinates - do i=1,2*nres + 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 - - ishift_pdb=ishift 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 ac97e59..db000f4 100644 --- a/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos +++ b/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos @@ -161,7 +161,7 @@ cartder.o : cartder.F readpdb.o : readpdb.F ${FC} ${FFLAGS2} ${CPPFLAGS} readpdb.F -readpdb.o : readpdb.F +readpdb-mult.o : readpdb-mult.F ${FC} ${FFLAGS2} ${CPPFLAGS} readpdb-mult.F sumsld.o : sumsld.f diff --git a/source/unres/src-HCD-5D/int_from_cart.F b/source/unres/src-HCD-5D/int_from_cart.F index 2ae703b..b36fae4 100644 --- a/source/unres/src-HCD-5D/int_from_cart.F +++ b/source/unres/src-HCD-5D/int_from_cart.F @@ -227,6 +227,8 @@ c--------------------------------------------------------------------------- integer i,j,ires,nscat double precision sccor(3,50) double precision sccmj + write (2,*) "sccenter",ires,nscat + call flush(2) do j=1,3 sccmj=0.0D0 do i=1,nscat diff --git a/source/unres/src-HCD-5D/readpdb-mult.F b/source/unres/src-HCD-5D/readpdb-mult.F index 9b99e64..cc9aec9 100644 --- a/source/unres/src-HCD-5D/readpdb-mult.F +++ b/source/unres/src-HCD-5D/readpdb-mult.F @@ -31,9 +31,10 @@ C geometry. lsecondary=.false. nhfrag=0 nbfrag=0 + iii=0 do read (ipdbin,'(a80)',end=10) card -! write (iout,'(a)') card +c write (iout,'(a)') card if (card(:5).eq.'HELIX') then nhfrag=nhfrag+1 lsecondary=.true. @@ -59,9 +60,10 @@ C geometry. itype(ires_old-1)=ntyp1 iterter(ires_old-1)=1 itype(ires_old)=ntyp1 + iterter(ires_old)=1 ishift1=ishift1+1 ibeg=2 -! write (iout,*) "Chain ended",ires,ishift,ires_old + write (iout,*) "Chain ended",ires,ishift,ires_old if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) @@ -69,7 +71,7 @@ C geometry. else call sccenter(ires,iii,sccor) endif -c iii=0 + iii=0 endif ! Read free energy if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp @@ -77,19 +79,20 @@ c iii=0 if (index(card(1:4),'ATOM').gt.0) then read (card(12:16),*) atom c write (2,'(a)') card -! write (iout,*) "! ",atom," !",ires +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 -! write (iout,*) "ires",ires,ires-ishift+ishift1, -! & " ires_old",ires_old -! write (iout,*) "ishift",ishift," ishift1",ishift1 -! write (iout,*) "IRES",ires-ishift+ishift1,ires_old +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 -! write (iout,*) "Calculating sidechain center iii",iii +c write (iout,*) "Calculating sidechain center iii",iii if (unres_pdb) then do j=1,3 dc(j,ires+nres)=sccor(j,iii) @@ -119,7 +122,7 @@ c write (iout,*) "BEG ires",ires ! Start a new chain ishift=-ires_old+ires-1 !!!!! ishift1=ishift1-1 !!!!! -! write (iout,*) "New chain started",ires,ishift,ishift1,"!" +c write (iout,*) "New chain started",ires,ishift,ishift1,"!" ires=ires-ishift+ishift1 ires_old=ires ibeg=0 @@ -136,11 +139,11 @@ c write (iout,*) "BEG ires",ires else ires=ires-ishift+ishift1 endif -! write (iout,*) "ires_old",ires_old," ires",ires +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 -! write (2,*) "ires",ires," res ",res!," ity"!,ity +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) diff --git a/source/unres/src-HCD-5D/readpdb.F b/source/unres/src-HCD-5D/readpdb.F index 9c8462d..1ffc884 100644 --- a/source/unres/src-HCD-5D/readpdb.F +++ b/source/unres/src-HCD-5D/readpdb.F @@ -1,7 +1,7 @@ subroutine readpdb C Read the PDB file and convert the peptide geometry into virtual-chain C geometry. - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.LOCAL' include 'COMMON.VAR' @@ -11,33 +11,27 @@ C geometry. include 'COMMON.GEO' include 'COMMON.NAMES' include 'COMMON.CONTROL' - include 'COMMON.DISTFIT' - include 'COMMON.SETUP' include 'COMMON.FRAG' - integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity, - & ishift_pdb - logical lprn /.false./,fail - double precision e1(3),e2(3),e3(3) - double precision dcj,efree_temp - character*3 seq,res - character*5 atom + include 'COMMON.SETUP' + include 'COMMON.SBRIDGE' + character*3 seq,atom,res character*80 card double precision sccor(3,50) - integer rescode - efree_temp=0.0d0 + double precision e1(3),e2(3),e3(3) + integer rescode,iterter(maxres),cou + logical fail + integer i,j,iii,ires,ires_old,ishift,ibeg + double precision dcj + bfac=0.0d0 + 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 -c write (iout,'(a)') card if (card(:5).eq.'HELIX') then nhfrag=nhfrag+1 lsecondary=.true. @@ -55,112 +49,145 @@ crc to be corrected !!! bfrag(4,nbfrag)=bfrag(2,nbfrag) crc---------------------------------------- endif - if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10 -c Read free energy - if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp + 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 + 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 + read (card(14:16),'(a3)') atom + if (atom.eq.'CA' .or. atom.eq.'CH3') then C Calculate the CM of the preceding residue. -c if (ibeg.eq.0) call sccenter(ires,iii,sccor) if (ibeg.eq.0) then -c write (iout,*) "Calculating sidechain center iii",iii if (unres_pdb) then do j=1,3 - dc(j,ires)=sccor(j,iii) + dc(j,ires+nres)=sccor(j,iii) enddo else - call sccenter(ires_old,iii,sccor) + call sccenter(ires,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 +c write (iout,'(a80)') card + read (card(23:26),*) ires + read (card(18:20),'(a3)') res + if (ibeg.eq.1) then 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 - ishift=ishift-(ires-ishift+ishift1-ires_old-1) - ires=ires-ishift+ishift1 - ires_old=ires + else if (ibeg.eq.2) then +c Start a new chain + ishift=-ires_old+ires-1 +c write (iout,*) "New chain started",ires,ishift + ibeg=0 endif - if (res.eq.'ACE' .or. res.eq.'NHE') then + ires=ires-ishift +c write (2,*) "ires",ires," ishift",ishift + if (res.eq.'ACE') 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 -c 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 -#ifdef DEBUG - write (iout,'(2i3,2x,a,3f8.3)') - & ires,itype(ires),res,(c(j,ires),j=1,3) -#endif - iii=iii+1 + read(card(61:66),*) bfac(ires) +c if(me.eq.king.or..not.out1file) +c & write (iout,'(2i3,2x,a,3f8.3)') +c & ires,itype(ires),res,(c(j,ires),j=1,3) + 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 + else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and. + & atom.ne.'N ' .and. atom.ne.'C ') then iii=iii+1 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) endif endif enddo - 10 continue -#ifdef DEBUG - write (iout,'(a,i5)') ' Number of residues found: ',ires -#endif - if (ires.eq.0) return + 10 if(me.eq.king.or..not.out1file) + & 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),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 + 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 + 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) + 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*e2(j) + 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 + write (iout,*) "After loop in readpbd" C Calculate the CM of the last side chain. - if (iii.gt.0) then if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) enddo - else + else call sccenter(ires,iii,sccor) endif - endif - nres=ires nsup=nres nstart_sup=1 if (itype(nres).ne.10) then @@ -175,11 +202,12 @@ 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)-3.8d0*e2(j) + 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) + 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 @@ -206,46 +234,27 @@ 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)-3.8d0*e2(j) + 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) + 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 (lprn) 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,i3,3f8.3,5x,3f8.3)') - & restyp(itype(ires)),ires,(c(j,ires),j=1,3), - & (c(j,ires+nres),j=1,3) - enddo - endif C Calculate internal coordinates. if(me.eq.king.or..not.out1file)then - write (iout,'(a)') - & "Backbone and SC coordinates as read from the PDB" do ires=1,nres write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3), & (c(j,nres+ires),j=1,3) enddo endif + call flush(iout) +c write(iout,*)"before int_from_cart nres",nres call int_from_cart(.true.,.false.) - call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) @@ -255,6 +264,8 @@ C Calculate internal coordinates. 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 @@ -264,15 +275,22 @@ C Calculate internal coordinates. 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,2*nres + 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 ' ', 6X,'X',11X,'Y',11X,'Z', + & 10X,'X',11X,'Y',11X,'Z') + 110 format (a,'(',i3,')',6f12.5) +cc enddiag do j=1,nbfrag do i=1,4 bfrag(i,j)=bfrag(i,j)-ishift @@ -284,14 +302,13 @@ C Copy the coordinates to reference coordinates hfrag(i,j)=hfrag(i,j)-ishift enddo enddo - ishift_pdb=ishift return end -c----------------------------------------------------------------------- +c--------------------------------------------------------------------------- subroutine readpdb_template(k) -C Read the PDB file with gaps for read_constr_homology with read2sigma +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) + implicit none include 'DIMENSIONS' include 'COMMON.LOCAL' include 'COMMON.VAR' @@ -301,20 +318,21 @@ C and convert the peptide geometry into virtual-chain geometry. include 'COMMON.GEO' include 'COMMON.NAMES' include 'COMMON.CONTROL' - include 'COMMON.DISTFIT' + include 'COMMON.FRAG' include 'COMMON.SETUP' - include 'COMMON.MD' - integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity, - & ishift_pdb + integer i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, + & ishift_pdb,ires_ca 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,50) - integer rescode - efree_temp=0.0d0 + double precision sccor(3,20) + integer rescode,iterter(maxres) + do i=1,maxres + iterter(i)=0 + enddo ibeg=1 ishift1=0 ishift=0 @@ -325,10 +343,27 @@ c write (2,*) "UNRES_PDB",unres_pdb lsecondary=.false. nhfrag=0 nbfrag=0 - do + do read (ipdbin,'(a80)',end=10) card -c write (iout,'(a)') card - if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10 + 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 @@ -342,9 +377,7 @@ 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. -c if (ibeg.eq.0) call sccenter(ires,iii,sccor) if (ibeg.eq.0) then -c write (iout,*) "Calculating sidechain center iii",iii if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) @@ -369,6 +402,13 @@ c write (iout,*) "BEG ires",ires 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) @@ -384,14 +424,14 @@ c & " ires_old",ires_old 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 +c if (card(27:27).eq."A" .or. card(27:27).eq."B") then c ishift1=ishift1+1 - endif +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 +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) @@ -416,13 +456,61 @@ c write (iout,*) "sidechain ",atom endif endif enddo - 10 continue -#ifdef DEBUG - write (iout,'(a,i5)') ' Number of residues found: ',ires -#endif - if (ires.eq.0) return + 10 if(me.eq.king.or..not.out1file) + & 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*e2(j) + 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 (iii.gt.0) then if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) @@ -430,19 +518,31 @@ C Calculate the CM of the last side chain. else call sccenter(ires,iii,sccor) endif - endif - nres=ires 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*e2(j) + enddo + else do j=1,3 - dcj=c(j,nres-2)-c(j,nres-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) @@ -464,18 +564,24 @@ 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)-3.8d0*e2(j) + c(j,1)=c(j,2)-1.9d0*e2(j) enddo else do j=1,3 - dcj=c(j,4)-c(j,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 (lprn) then + if (out_template_coord) then write (iout,'(/a)') & "Cartesian coordinates of the reference structure" write (iout,'(a,3(3x,a5),5x,3(3x,a5))') @@ -487,16 +593,7 @@ C Calculate internal coordinates. enddo endif C Calculate internal coordinates. - if(me.eq.king.or..not.out1file)then - write (iout,'(a)') - & "Backbone and SC coordinates as read from the PDB" - do ires=1,nres - write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') - & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3), - & (c(j,nres+ires),j=1,3) - enddo - endif - call int_from_cart(.true.,.false.) + call int_from_cart(.true.,.true.) call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) @@ -516,16 +613,18 @@ C Calculate internal coordinates. 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 -c call chainbuild -C Copy the coordinates to reference coordinates - do i=1,2*nres + 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 - - ishift_pdb=ishift return end + diff --git a/source/unres/src-HCD-5D/sc_move.F b/source/unres/src-HCD-5D/sc_move.F index 201e20f..7911b30 100644 --- a/source/unres/src-HCD-5D/sc_move.F +++ b/source/unres/src-HCD-5D/sc_move.F @@ -369,7 +369,7 @@ c parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) external funcgrad_restr1 external optsave #else - external func,gradient,fdum + external fdum external func_restr1,grad_restr1 logical not_done,change,reduce dimension iv(liv) diff --git a/source/unres/src-HCD-5D/unres.F b/source/unres/src-HCD-5D/unres.F index 08caade..be40e75 100644 --- a/source/unres/src-HCD-5D/unres.F +++ b/source/unres/src-HCD-5D/unres.F @@ -339,6 +339,7 @@ c print *,'Calling MINIMIZE.' potE=etot call cartoutx(0.0d0) endif + if (outpdb) call pdbout(etot,titel(:50),ipdb) if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) #ifdef LBFGS write (iout,'(a,a9)') 'LBFGS return code:',status diff --git a/source/wham/src-HCD/energy_p_new.F b/source/wham/src-HCD/energy_p_new.F index 835045b..bd69774 100644 --- a/source/wham/src-HCD/energy_p_new.F +++ b/source/wham/src-HCD/energy_p_new.F @@ -20,6 +20,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.SHIELD' include 'COMMON.CONTROL' include 'COMMON.TORCNSTR' + include 'COMMON.SAXS' double precision fact(6) c write(iout, '(a,i2)')'Calling etotal ipot=',ipot c call flush(iout) @@ -179,12 +180,16 @@ c write(iout,*) "TEST_ENE1 constr_homology=",constr_homology c write(iout,*) "TEST_ENE1 ehomology_constr=",ehomology_constr #ifdef DFA C BARTEK for dfa test! + edfadis=0.0d0 if (wdfa_dist.gt.0) call edfad(edfadis) c write(iout,*)'edfad is finished!', wdfa_dist,edfadis + edfator=0.0d0 if (wdfa_tor.gt.0) call edfat(edfator) c write(iout,*)'edfat is finished!', wdfa_tor,edfator + edfanei=0.0d0 if (wdfa_nei.gt.0) call edfan(edfanei) c write(iout,*)'edfan is finished!', wdfa_nei,edfanei + edfabet=0.0d0 if (wdfa_beta.gt.0) call edfab(edfabet) c write(iout,*)'edfab is finished!', wdfa_beta,edfabet #endif diff --git a/source/wham/src-HCD/readpdb.F b/source/wham/src-HCD/readpdb.F index 5993291..855d754 100644 --- a/source/wham/src-HCD/readpdb.F +++ b/source/wham/src-HCD/readpdb.F @@ -23,6 +23,7 @@ C geometry. logical lsecondary integer iterter(maxres) double precision efree_temp + iii=0 ibeg=1 ishift1=0 do @@ -53,6 +54,7 @@ C geometry. itype(ires_old-1)=ntyp1 iterter(ires_old-1)=1 itype(ires_old)=ntyp1 + iterter(ires_old)=1 ishift1=ishift1+1 ibeg=2 ! write (iout,*) "Chain ended",ires,ishift,ires_old @@ -63,7 +65,7 @@ C geometry. else call sccenter(ires,iii,sccor) endif -c iii=0 + iii=0 endif ! Read free energy if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp -- 1.7.9.5