unres

diff --git a/source/unres/src-HCD-5D/COMMON.INTERACT b/source/unres/src-HCD-5D/COMMON.INTERACT
index 9b023e5..2a55393 100644 (file)
--- a/source/unres/src-HCD-5D/COMMON.INTERACT
+++ b/source/unres/src-HCD-5D/COMMON.INTERACT
@@ -1,6 +1,6 @@
       double precision aa,bb,augm,aad,bad,app,bpp,ale6,ael3,ael6,
      &aa_lip,bb_lip,aa_aq,bb_aq,sc_aa_tube_par,sc_bb_tube_par,
-     & pep_aa_tube,pep_bb_tube
+     & pep_aa_tube,pep_bb_tube,alpha_GB,alpha_GB1
 
       integer expon,expon2
       integer nnt,nct,nint_gr,istart,iend,itype,itel,itypro,
@@ -11,7 +11,7 @@
       common /interact/aa_aq(ntyp,ntyp),bb_aq(ntyp,ntyp),
      & aa_lip(ntyp,ntyp),bb_lip(ntyp,ntyp),
      & sc_aa_tube_par(ntyp),sc_bb_tube_par(ntyp),
-     & pep_aa_tube,pep_bb_tube,
+     & pep_aa_tube,pep_bb_tube,alpha_GB,alpha_GB1,
      & augm(ntyp,ntyp),
      & aad(ntyp,2),bad(ntyp,2),app(2,2),bpp(2,2),ael6(2,2),ael3(2,2),
      & expon,expon2,nnt,nct,nint_gr(maxres),istart(maxres,maxint_gr),

diff --git a/source/unres/src-HCD-5D/elecont.f b/source/unres/src-HCD-5D/elecont.f
index 7c024ea..79d6a7e 100644 (file)
--- a/source/unres/src-HCD-5D/elecont.f
+++ b/source/unres/src-HCD-5D/elecont.f
@@ -428,7 +428,7 @@ cd            write (iout,*) i1,j1,not_done
             iii1=max0(ii1-1,1)
             do ij=iii1,i1
              nsec(ij)=nsec(ij)+1
-             if (nsec(ij).le.2) then
+             if (nsec(ij).le.2 .and. nsec(ij).gt.0) then
               isec(ij,nsec(ij))=nbeta
              endif
             enddo

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 2d94dc0..9ec8107 100644 (file)
--- a/source/unres/src-HCD-5D/energy_p_new_barrier.F
+++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F
@@ -2013,6 +2013,11 @@ C
      & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
       double precision boxshift
+      double precision x0,y0,r012,rij12,facx0,
+     &  facx02,afacx0,bfacx0,abfacx0,Afac,BBfac,Afacsig,Bfacsig
+c      alpha_GB=0.5d0
+c      alpha_GB=0.01d0
+c      alpha_GB1=1.0d0+1.0d0/alpha_GB 
       evdw=0.0D0
 ccccc      energy_dec=.false.
 C      print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -2173,67 +2178,150 @@ c     &       " sig",sig," sig0ij",sig0ij
 c for diagnostics; uncomment
 c            rij_shift=1.2*sig0ij
 C I hate to put IF's in the loops, but here don't have another choice!!!!
-              if (rij_shift.le.0.0D0) then
-                evdw=1.0D20
+c              if (rij_shift.le.0.0D0) then
+              x0=alpha_GB*(sig-sig0ij)
+              if (energy_dec) write (iout,*) i,j," x0",x0
+              if (rij_shift.le.x0) then
+c                sig=2.0d0*sig0ij
+                sigder=-sig*sigsq
+c                sigder=0.0d0
+                fac=rij**expon
+                rij12=fac*fac
+c                rij12=1.0d0
+                x0=alpha_GB*(sig-sig0ij)
+                facx0=1.0d0/x0**expon
+                facx02=facx0*facx0
+                r012=((1.0d0+alpha_GB)*(sig-sig0ij))**(2*expon)
+                afacx0=aa*facx02
+                bfacx0=bb*facx0
+                abfacx0=afacx0+0.5d0*bfacx0
+                Afac=alpha_GB1*abfacx0
+                Afacsig=0.5d0*alpha_GB1*bfacx0/(sig-sig0ij)
+                BBfac=Afac-(afacx0+bfacx0)
+c                BBfac=0.0d0
+                Bfacsig=(-alpha_GB1*(abfacx0+afacx0)+
+     &              (afacx0+afacx0+bfacx0))/(sig-sig0ij)
+c                Bfacsig=0.0d0
+                Afac=Afac*r012
+                Afacsig=Afacsig*r012
+c                Afac=1.0d0
+c                Afacsig=0.0d0
+c    w(x)=4*eps*((1.0+1.0/alpha_GB)*(y0**12-0.5*y0**6)*(r0/x)**12-(1+1/alpha)*(y0**12-0.5*y0**6)+y0**12-y0**6)
+c                eps1=1.0d0
+c                eps2rt=1.0d0
+c                eps3rt=1.0d0
+                e1 = eps1*eps2rt*eps3rt*Afac*rij12
+                e2 = -eps1*eps2rt*eps3rt*BBfac
+                evdwij = e1+e2
+                eps2der=evdwij*eps3rt
+                eps3der=evdwij*eps2rt
+c                eps2der=0.0d0
+c                eps3der=0.0d0
+c                eps1_om12=0.0d0
+                evdwij=evdwij*eps2rt*eps3rt
+c                Afacsig=0.d0
+c                Bfacsig=0.0d0
+                if (lprn) then
+                  write (iout,*) "aa",aa," bb",bb," sig0ij",sig0ij
+                  sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+                  epsi=bb**2/aa
+                  write (iout,'(2(a3,i3,2x),18(0pf9.5))')
+     &             restyp(itypi),i,restyp(itypj),j,
+     &             epsi,sigm,chi1,chi2,chip1,chip2,
+     &             eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+     &             eps1*eps2rt**2*eps3rt**2,om1,om2,om12,
+     &             1.0D0/rij,rij_shift,
+     &             evdwij
+                endif
+                if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') 
+     &          'RE r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
+                evdw=evdw+evdwij*sss
+C Calculate gradient components.
+                e1=e1*eps2rt*eps3rt
+                sigder=-expon*eps1*eps2rt*eps2rt*eps3rt*eps3rt
+     &            *(Afacsig*rij12-Bfacsig)*sigder
+                fac=-2.0d0*expon*e1*rij*rij
+c              print '(2i4,6f8.4)',i,j,sss,sssgrad*
+c     &        evdwij,fac,sigma(itypi,itypj),expon
+                fac=fac+evdwij*sssgrad/sss*rij
+c                fac=0.0d0
+c                write (iout,*) "sigder",sigder," fac",fac," e1",e1,
+c     &              " e2",e2," sss",sss," sssgrad",sssgrad,"esp123",
+c     &              eps1*eps2rt**2*eps3rt**2
+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
+C              gg_lipi(3)=0.0d0
+C              gg_lipj(3)=0.0d0
+                gg(1)=xj*fac
+                gg(2)=yj*fac
+                gg(3)=zj*fac
+c                evdw=1.0D20
 cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 cd     &        restyp(itypi),i,restyp(itypj),j,
 cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
-                return
-              endif
-              sigder=-sig*sigsq
+c                return
+              else
+                rij_shift=1.0D0/rij_shift 
+                sigder=-sig*sigsq
 c---------------------------------------------------------------
-              rij_shift=1.0D0/rij_shift 
-              fac=rij_shift**expon
+                fac=rij_shift**expon
 C here to start with
 C            if (c(i,3).gt.
-              faclip=fac
-              e1=fac*fac*aa
-              e2=fac*bb
-              evdwij=eps1*eps2rt*eps3rt*(e1+e2)
-              eps2der=evdwij*eps3rt
-              eps3der=evdwij*eps2rt
+                faclip=fac
+                e1=fac*fac*aa
+                e2=fac*bb
+                evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+                eps2der=evdwij*eps3rt
+                eps3der=evdwij*eps2rt
 C       write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij,
 C     &((sslipi+sslipj)/2.0d0+
 C     &(2.0d0-sslipi-sslipj)/2.0d0)
 c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
-              evdwij=evdwij*eps2rt*eps3rt
-              evdw=evdw+evdwij*sss
-              if (lprn) then
-                sigm=dabs(aa/bb)**(1.0D0/6.0D0)
-                epsi=bb**2/aa
-                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-     &           restyp(itypi),i,restyp(itypj),j,
-     &           epsi,sigm,chi1,chi2,chip1,chip2,
-     &           eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
-     &           om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
-     &           evdwij
-              endif
+                evdwij=evdwij*eps2rt*eps3rt
+                evdw=evdw+evdwij*sss
+                if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') 
+     &          'GB r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
+                if (lprn) then
+                  write (iout,*) "aa",aa," bb",bb," sig0ij",sig0ij
+                  sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+                  epsi=bb**2/aa
+                  write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+     &             restyp(itypi),i,restyp(itypj),j,
+     &             epsi,sigm,chi1,chi2,chip1,chip2,
+     &             eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+     &             om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+     &             evdwij
+                endif
 
-              if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)') 
-     &          'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
 
 C Calculate gradient components.
-              e1=e1*eps1*eps2rt**2*eps3rt**2
-              fac=-expon*(e1+evdwij)*rij_shift
-              sigder=fac*sigder
-              fac=rij*fac
-c            print '(2i4,6f8.4)',i,j,sss,sssgrad*
-c     &      evdwij,fac,sigma(itypi,itypj),expon
-              fac=fac+evdwij*sssgrad/sss*rij
-c            fac=0.0d0
+                e1=e1*eps1*eps2rt**2*eps3rt**2
+                fac=-expon*(e1+evdwij)*rij_shift
+                sigder=fac*sigder
+                fac=rij*fac
+c              print '(2i4,6f8.4)',i,j,sss,sssgrad*
+c     &        evdwij,fac,sigma(itypi,itypj),expon
+                fac=fac+evdwij*sssgrad/sss*rij
+c              fac=0.0d0
 C Calculate the radial part of the gradient
-              gg_lipi(3)=eps1*(eps2rt*eps2rt)
+                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
-C            gg_lipi(3)=0.0d0
-C            gg_lipj(3)=0.0d0
-              gg(1)=xj*fac
-              gg(2)=yj*fac
-              gg(3)=zj*fac
+                gg_lipj(3)=ssgradlipj*gg_lipi(3)
+                gg_lipi(3)=gg_lipi(3)*ssgradlipi
+C              gg_lipi(3)=0.0d0
+C              gg_lipj(3)=0.0d0
+                gg(1)=xj*fac
+                gg(2)=yj*fac
+                gg(3)=zj*fac
+              endif
 C Calculate angular part of the gradient.
 c            call sc_grad_scale(sss)
               call sc_grad

diff --git a/source/unres/src-HCD-5D/geomout.F b/source/unres/src-HCD-5D/geomout.F
index 3dcde10..553dbf6 100644 (file)
--- a/source/unres/src-HCD-5D/geomout.F
+++ b/source/unres/src-HCD-5D/geomout.F
@@ -145,7 +145,7 @@ cmodel      write (iunit,'(a5,i6)') 'MODEL',1
       endif
       do i=1,nss
        if (dyn_ss) then
-        write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
+        write (iunit,30) ica(iss(idssb(i)))+1,ica(iss(jdssb(i)))+1
        else
         write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
        endif
@@ -153,7 +153,7 @@ cmodel      write (iunit,'(a5,i6)') 'MODEL',1
       write (iunit,'(a6)') 'ENDMDL'     
   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
-  30  FORMAT ('CONECT',8I7)
+  30  FORMAT ('CONECT',8I5)
       return
       end
 c------------------------------------------------------------------------------
@@ -335,6 +335,7 @@ c-----------------------------------------------------------------
 #else
       parameter (me=0)
 #endif
+      include 'COMMON.CONTROL'
       include 'COMMON.CHAIN'
       include 'COMMON.INTERACT'
       include 'COMMON.NAMES'
@@ -359,8 +360,13 @@ c-----------------------------------------------------------------
       call xdrfint_(ixdrf, nss, iret) 
       do j=1,nss
        if (dyn_ss) then
-        call xdrfint_(ixdrf, idssb(j)+nres, iret)
-        call xdrfint_(ixdrf, jdssb(j)+nres, iret)
+        if (modecalc.eq.14) then
+          call xdrfint_(ixdrf, idssb(j), iret)
+          call xdrfint_(ixdrf, jdssb(j), iret)
+        else
+          call xdrfint_(ixdrf, iss(idssb(j))+nres, iret)
+          call xdrfint_(ixdrf, iss(jdssb(j))+nres, iret)
+        endif
        else
         call xdrfint_(ixdrf, ihpb(j), iret)
         call xdrfint_(ixdrf, jhpb(j), iret)
@@ -391,8 +397,13 @@ c     & " nss",nss
       call xdrfint(ixdrf, nss, iret) 
       do j=1,nss
        if (dyn_ss) then
-        call xdrfint(ixdrf, idssb(j)+nres, iret)
-        call xdrfint(ixdrf, jdssb(j)+nres, iret)
+        if (modecalc.eq.14) then
+          call xdrfint(ixdrf, idssb(j), iret)
+          call xdrfint(ixdrf, jdssb(j), iret)
+        else
+          call xdrfint(ixdrf, iss(idssb(j))+nres, iret)
+          call xdrfint(ixdrf, iss(jdssb(j))+nres, iret)
+        endif
        else
         call xdrfint(ixdrf, ihpb(j), iret)
         call xdrfint(ixdrf, jhpb(j), iret)

diff --git a/source/unres/src-HCD-5D/readpdb-mult.F b/source/unres/src-HCD-5D/readpdb-mult.F
index 12d011e..e603267 100644 (file)
--- a/source/unres/src-HCD-5D/readpdb-mult.F
+++ b/source/unres/src-HCD-5D/readpdb-mult.F
@@ -28,6 +28,7 @@ C geometry.
          iterter(i)=0
       enddo
       ibeg=1
+      ishift=0
       ishift1=0
       lsecondary=.false.
       nhfrag=0

diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F
index 58b86cc..aeafdc3 100644 (file)
--- a/source/unres/src-HCD-5D/readrtns_CSA.F
+++ b/source/unres/src-HCD-5D/readrtns_CSA.F
@@ -196,6 +196,12 @@ c      print *,'AFMlog',AFMlog,selfguide,"KUPA"
       call readi(controlcard,'TUBEMOD',tubelog,0)
 c      write (iout,*) TUBElog,"TUBEMODE"
       call readi(controlcard,'IPRINT',iprint,0)
+c 6/22/2021 AL: alpha_GB: parameter to switch between the GB SC-SC
+c    interaction potential and the all-repulsive potential with singularity
+c    at zero site-site distance
+      call reada(controlcard,'ALPHA_GB',alpha_GB,1.0d-2)
+      alpha_GB1 = 1.0d0+1.0d0/alpha_GB
+      write (iout,*) "alpha_GB",alpha_GB," alpha_GB1",alpha_GB1
 C SHIELD keyword sets if the shielding effect of side-chains is used
 C 0 denots no shielding is used all peptide are equally despite the 
 C solvent accesible area
@@ -3066,7 +3072,8 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
       character*2 kic2
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
-      integer ki,i,ii,j,k,l,ii_in_use(maxdim_cont),i_tmp,idomain_tmp,
+      integer ki,i,ii,j,k,l,ii_in_use(maxchain*maxdim),i_tmp,
+     & idomain_tmp,
      & irec,ik,iistart,nres_temp
       integer ilen
       external ilen

diff --git a/source/unres/src-HCD-5D/unres.F b/source/unres/src-HCD-5D/unres.F
index 2e0ebaf..4c60573 100644 (file)
--- a/source/unres/src-HCD-5D/unres.F
+++ b/source/unres/src-HCD-5D/unres.F
@@ -294,10 +294,10 @@ c      print *,"after etotal"
       call enerprint(energy(0))
       if (secondary_str) then
       call hairpin(.true.,nharp,iharp)
-c        print *,'after hairpin'
+      print *,'after hairpin'
       call secondary2(.true.)
       endif
-c        print *,'after secondary'
+      print *,'after secondary'
       if (minim) then
 crc overlap test
         if (indpdb.ne.0 .and. .not.dccart) then