Adam's cluster & unres corrections

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 (file)
--- 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 (file)
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 (file)
--- 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 (file)
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 (file)
--- 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',