From: czarek Date: Thu, 17 Feb 2022 09:09:54 +0000 (+0100) Subject: unres X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=refs%2Fheads%2Fadam unres --- diff --git a/source/unres/src-HCD-5D/COMMON.INTERACT b/source/unres/src-HCD-5D/COMMON.INTERACT index 9b023e5..2a55393 100644 --- 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 --- 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 --- 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 --- 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 --- 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 --- 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 --- 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