From 25149eedfd23dd08ea519ea9a3c9f3cc0279805f Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Sat, 6 Jun 2020 17:14:46 +0200 Subject: [PATCH] Adam's cluster & unres corrections --- source/cluster/wham/src-HCD-5D/energy_p_new.F | 2 +- source/unres/src-HCD-5D/boxshift.f | 101 ---- source/unres/src-HCD-5D/energy_p_new_barrier.F | 6 + source/unres/src-HCD-5D/readpdb.F | 631 ------------------------ source/unres/src-HCD-5D/readrtns_CSA.F | 10 +- 5 files changed, 15 insertions(+), 735 deletions(-) delete mode 100644 source/unres/src-HCD-5D/boxshift.f delete mode 100644 source/unres/src-HCD-5D/readpdb.F diff --git a/source/cluster/wham/src-HCD-5D/energy_p_new.F b/source/cluster/wham/src-HCD-5D/energy_p_new.F index 119bad6..db2e043 100644 --- a/source/cluster/wham/src-HCD-5D/energy_p_new.F +++ b/source/cluster/wham/src-HCD-5D/energy_p_new.F @@ -2899,7 +2899,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eelloc',i,j,eel_loc_ij c if (eel_loc_ij.ne.0) diff --git a/source/unres/src-HCD-5D/boxshift.f b/source/unres/src-HCD-5D/boxshift.f deleted file mode 100644 index 29d3406..0000000 --- a/source/unres/src-HCD-5D/boxshift.f +++ /dev/null @@ -1,101 +0,0 @@ - -c------------------------------------------------------------------------ - double precision function boxshift(x,boxsize) - implicit none - double precision x,boxsize - double precision xtemp - xtemp=dmod(x,boxsize) - if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then - boxshift=xtemp-boxsize - else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then - boxshift=xtemp+boxsize - else - boxshift=xtemp - endif - return - end -c-------------------------------------------------------------------------- - subroutine closest_img(xi,yi,zi,xj,yj,zj) - include 'DIMENSIONS' - include 'COMMON.CHAIN' - integer xshift,yshift,zshift,subchap - double precision dist_init,xj_safe,yj_safe,zj_safe, - & xj_temp,yj_temp,zj_temp,dist_temp - xj_safe=xj - yj_safe=yj - zj_safe=zj - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - subchap=0 - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 - xj=xj_safe+xshift*boxxsize - yj=yj_safe+yshift*boxysize - zj=zj_safe+zshift*boxzsize - dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - if(dist_temp.lt.dist_init) then - dist_init=dist_temp - xj_temp=xj - yj_temp=yj - zj_temp=zj - subchap=1 - endif - enddo - enddo - enddo - if (subchap.eq.1) then - xj=xj_temp-xi - yj=yj_temp-yi - zj=zj_temp-zi - else - xj=xj_safe-xi - yj=yj_safe-yi - zj=zj_safe-zi - endif - return - end -c-------------------------------------------------------------------------- - subroutine to_box(xi,yi,zi) - implicit none - include 'DIMENSIONS' - include 'COMMON.CHAIN' - double precision xi,yi,zi - xi=dmod(xi,boxxsize) - if (xi.lt.0.0d0) xi=xi+boxxsize - yi=dmod(yi,boxysize) - if (yi.lt.0.0d0) yi=yi+boxysize - zi=dmod(zi,boxzsize) - if (zi.lt.0.0d0) zi=zi+boxzsize - return - end -c-------------------------------------------------------------------------- - subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi) - implicit none - include 'DIMENSIONS' - include 'COMMON.CHAIN' - double precision xi,yi,zi,sslipi,ssgradlipi - double precision fracinbuf - double precision sscalelip,sscagradlip - - if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then -C the energy transfer exist - if (zi.lt.buflipbot) then -C what fraction I am in - fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick) -C lipbufthick is thickenes of lipid buffore - sslipi=sscalelip(fracinbuf) - ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) - sslipi=sscalelip(fracinbuf) - ssgradlipi=sscagradlip(fracinbuf)/lipbufthick - else - sslipi=1.0d0 - ssgradlipi=0.0 - endif - else - sslipi=0.0d0 - ssgradlipi=0.0 - endif - return - end diff --git a/source/unres/src-HCD-5D/energy_p_new_barrier.F b/source/unres/src-HCD-5D/energy_p_new_barrier.F index 3f5429d..2c701ca 100644 --- a/source/unres/src-HCD-5D/energy_p_new_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F @@ -2308,6 +2308,12 @@ C Calculate gradient components. fac=rij*fac-2*expon*rrij*e_augm fac=fac+(evdwij+e_augm)*sssgrad/sss*rij C Calculate the radial part of the gradient + gg_lipi(3)=eps1*(eps2rt*eps2rt) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + gg_lipj(3)=ssgradlipj*gg_lipi(3) + gg_lipi(3)=gg_lipi(3)*ssgradlipi gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac diff --git a/source/unres/src-HCD-5D/readpdb.F b/source/unres/src-HCD-5D/readpdb.F deleted file mode 100644 index c56b1df..0000000 --- a/source/unres/src-HCD-5D/readpdb.F +++ /dev/null @@ -1,631 +0,0 @@ - subroutine readpdb -C Read the PDB file and convert the peptide geometry into virtual-chain -C geometry. - 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' - include 'COMMON.CONTROL' - include 'COMMON.FRAG' - include 'COMMON.SETUP' - include 'COMMON.SBRIDGE' - character*3 seq,atom,res - character*80 card - double precision sccor(3,50) - 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 - lsecondary=.false. - nhfrag=0 - nbfrag=0 - ires=0 - do - read (ipdbin,'(a80)',end=10) card - 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) -crc---------------------------------------- -crc to be corrected !!! - bfrag(3,nbfrag)=bfrag(1,nbfrag) - bfrag(4,nbfrag)=bfrag(2,nbfrag) -crc---------------------------------------- - endif - 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(14:16),'(a3)') atom - if (atom.eq.'CA' .or. atom.eq.'CH3') 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+nres)=sccor(j,iii) - enddo - else - call sccenter(ires,iii,sccor) - endif - endif -C Start new residue. -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 -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 -c write (iout,*) "New chain started",ires,ishift - ibeg=0 - endif - 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 - read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) - 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 - 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 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 (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 Calculate internal coordinates. - if(me.eq.king.or..not.out1file)then - 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.) - do i=1,nres - thetaref(i)=theta(i) - phiref(i)=phi(i) - enddo - 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 ' ', 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 - 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 none - 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.FRAG' - include 'COMMON.SETUP' - 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,20) - integer rescode,iterter(maxres) - 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)=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 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 (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*e2(j) - 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*e2(j) - 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 -C Calculate internal coordinates. - call int_from_cart(.true.,.true.) - call sc_loc_geom(.false.) - do i=1,nres - thetaref(i)=theta(i) - phiref(i)=phi(i) - enddo - 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 - diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index d76b29e..f120ec7 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -1215,6 +1215,8 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb C endif +c write (iout,*) "iranconf",iranconf," extconf",extconf, +c & " start_from_models",start_from_model if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and. & modecalc.ne.10) then @@ -3053,7 +3055,8 @@ c double precision, dimension (max_template,maxres) :: rescore2 double precision, dimension (max_template,maxres) :: rescore3 double precision distal - character*24 pdbfile,tpl_k_rescore + character*24 tpl_k_rescore + character*256 pdbfile c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints c not using pre-computed ones stored in files model_ki_{dist,angle} @@ -3132,6 +3135,7 @@ c do k=1,constr_homology read(inp,'(a)') pdbfile + pdbfiles_chomo(k)=pdbfile if(me.eq.king .or. .not. out1file) & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file', & pdbfile(:ilen(pdbfile)) @@ -3636,7 +3640,8 @@ c double precision rescore_tmp,x12,y12,z12,rescore2_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 - character*24 pdbfile,tpl_k_rescore + character*24 tpl_k_rescore + character*256 pdbfile c c For new homol impl @@ -3655,6 +3660,7 @@ c Read pdb files read(ientin,'(a)') pdbfile write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', & pdbfile(:ilen(pdbfile)) + pdbfiles_chomo(k)=pdbfile open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a,5x,a)') 'Error opening PDB file', -- 1.7.9.5