double precision betaT
- integer iscode,indpdb,outpdb,outmol2,iopt,nstart,nend,symetr
+ integer iscode,indpdb,outpdb,outmol2,iopt,nstart,nend,symetr,
+ & constr_dist
logical refstr,pdbref,punch_dist,print_dist,caonly,lside,
- & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx
+ & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx,
+ & with_dihed_constr
common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2,
& punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int,
- & from_cart,from_bx,from_cx,efree,iopt,nstart,nend,symetr
+ & from_cart,from_bx,from_cx,efree,iopt,nstart,nend,symetr,
+ & constr_dist
integer ns,nss,nfree,iss
common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,
& ns,nss,nfree,iss(maxss)
- double precision dhpb,dhpb1,forcon
+ double precision dhpb,dhpb1,forcon,fordepth
integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb
common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim),
+ & fordepth(maxdim),
& ihpb(maxdim),jhpb(maxdim),nhpb
double precision weidis
common /restraints/ weidis
endif
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
- if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
- & iabs(itype(jjj)).eq.1) then
+C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+C & iabs(itype(jjj)).eq.1) then
+C call ssbond_ene(iii,jjj,eij)
+C ehpb=ehpb+2*eij
+C else
+ if (.not.dyn_ss .and. i.le.nss) then
+ if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+ & iabs(itype(jjj)).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
- else
+ endif !ii.gt.neres
+ else if (ii.gt.nres .and. jj.gt.nres) then
+c Restraints from contact prediction
+ dd=dist(ii,jj)
+ if (constr_dist.eq.11) then
+C ehpb=ehpb+fordepth(i)**4.0d0
+C & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ ehpb=ehpb+fordepth(i)**4.0d0
+ & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ fac=fordepth(i)**4.0d0
+ & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
+C & ehpb,fordepth(i),dd
+C write(iout,*) ehpb,"atu?"
+C ehpb,"tu?"
+C fac=fordepth(i)**4.0d0
+C & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+ else !constr_dist.eq.11
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "beta nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else !dhpb(i).gt.0.00
+
C Calculate the distance between the two points and its difference from the
C target distance.
dd=dist(ii,jj)
C Evaluate gradient.
C
fac=waga*rdis/dd
+ endif !dhpb(i).gt.0
+ endif
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
do j=1,3
ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
enddo
endif
+ else !ii.gt.nres
+C write(iout,*) "before"
+ dd=dist(ii,jj)
+C write(iout,*) "after",dd
+ if (constr_dist.eq.11) then
+ ehpb=ehpb+fordepth(i)**4.0d0
+ & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ fac=fordepth(i)**4.0d0
+ & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+C ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i))
+C fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd
+C print *,ehpb,"tu?"
+C write(iout,*) ehpb,"btu?",
+C & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i)
+C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
+C & ehpb,fordepth(i),dd
+ else
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "alph nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ rdis=dd-dhpb(i)
+C Get the force constant corresponding to this distance.
+ waga=forcon(i)
+C Calculate the contribution to energy.
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "alpha reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+ fac=waga*rdis/dd
+ endif
+ endif
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
+cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
+C If this is a SC-SC distance, we need to calculate the contributions to the
+C Cartesian gradient in the SC vectors (ghpbx).
+ if (iii.lt.ii) then
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
+ endif
do j=iii,jjj-1
do k=1,3
ghpbc(k,j)=ghpbc(k,j)+ggg(k)
enddo
endif
enddo
- ehpb=0.5D0*ehpb
+ if (constr_dist.ne.11) ehpb=0.5D0*ehpb
return
end
C--------------------------------------------------------------------------
pdbref=(index(controlcard,'PDBREF').gt.0)
iscode=index(controlcard,'ONE_LETTER')
tree=(index(controlcard,'MAKE_TREE').gt.0)
+ with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
+ call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+ write (iout,*) "with_dihed_constr ",with_dihed_constr,
+ & " CONSTR_DIST",constr_dist
+ call flush(iout)
min_var=(index(controlcard,'MINVAR').gt.0)
plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
include 'COMMON.CONTROL'
include 'COMMON.CONTACTS'
include 'COMMON.TIME1'
+ include 'COMMON.TORCNSTR'
#ifdef MPL
include 'COMMON.INFO'
#endif
print *,'Call Read_Bridge.'
call read_bridge
+C this fragment reads diheadral constrains
+ if (with_dihed_constr) then
+
+ read (inp,*) ndih_constr
+ if (ndih_constr.gt.0) then
+ read (inp,*) ftors
+ write (iout,*) 'FTORS',ftors
+C ftors is the force constant for torsional quartic constrains
+ read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
+ write (iout,*)
+ & 'There are',ndih_constr,' constraints on phi angles.'
+ do i=1,ndih_constr
+ write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
+ enddo
+ do i=1,ndih_constr
+ phi0(i)=deg2rad*phi0(i)
+ drange(i)=deg2rad*drange(i)
+ enddo
+ endif ! endif ndif_constr.gt.0
+ endif ! with_dihed_constr
nnt=1
nct=nres
print *,'NNT=',NNT,' NCT=',NCT
dyn_ss_mask(iss(i))=.true.
enddo
endif
+c Read distance restraints
+ if (constr_dist.gt.0) then
+ call read_dist_constr
+ call hpb_partition
+ endif
return
end
c-----------------------------------------------------------------------------
read (rekord(iread:),*) wartosc
return
end
+C----------------------------------------------------------------------
+ subroutine multreadi(rekord,lancuch,tablica,dim,default)
+ implicit none
+ integer dim,i
+ integer tablica(dim),default
+ character*(*) rekord,lancuch
+ character*80 aux
+ integer ilen,iread
+ external ilen
+ do i=1,dim
+ tablica(i)=default
+ enddo
+ iread=index(rekord,lancuch(:ilen(lancuch))//"=")
+ if (iread.eq.0) return
+ iread=iread+ilen(lancuch)+1
+ read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
+ 10 return
+ end
+
c----------------------------------------------------------------------------
subroutine card_concat(card)
include 'DIMENSIONS'
#endif
return
end
+ subroutine read_dist_constr
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ integer ifrag_(2,100),ipair_(2,100)
+ double precision wfrag_(100),wpair_(100)
+ character*500 controlcard
+ logical lprn /.true./
+ write (iout,*) "Calling read_dist_constr"
+C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+C call flush(iout)
+ write(iout,*) "TU sie wywalam?"
+ call card_concat(controlcard)
+ write (iout,*) controlcard
+ call flush(iout)
+ call readi(controlcard,"NFRAG",nfrag_,0)
+ call readi(controlcard,"NPAIR",npair_,0)
+ call readi(controlcard,"NDIST",ndist_,0)
+ call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+ call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
+ call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
+ call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
+ call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
+ write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+ write (iout,*) "IFRAG"
+ do i=1,nfrag_
+ write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ enddo
+ write (iout,*) "IPAIR"
+ do i=1,npair_
+ write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+ enddo
+ call flush(iout)
+ do i=1,nfrag_
+ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
+ if (ifrag_(2,i).gt.nstart_sup+nsup-1)
+ & ifrag_(2,i)=nstart_sup+nsup-1
+c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ call flush(iout)
+ if (wfrag_(i).gt.0.0d0) then
+ do j=ifrag_(1,i),ifrag_(2,i)-1
+ do k=j+1,ifrag_(2,i)
+ write (iout,*) "j",j," k",k
+ ddjk=dist(j,k)
+ if (constr_dist.eq.1) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ else if (constr_dist.eq.2) then
+ if (ddjk.le.dist_cut) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ endif
+ else
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ endif
+ if (lprn)
+ & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,npair_
+ if (wpair_(i).gt.0.0d0) then
+ ii = ipair_(1,i)
+ jj = ipair_(2,i)
+ if (ii.gt.jj) then
+ itemp=ii
+ ii=jj
+ jj=itemp
+ endif
+ do j=ifrag_(1,ii),ifrag_(2,ii)
+ do k=ifrag_(1,jj),ifrag_(2,jj)
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ forcon(nhpb)=wpair_(i)
+ dhpb(nhpb)=dist(j,k)
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,ndist_
+ if (constr_dist.eq.11) then
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+ & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
+ fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
+C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ else
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+ endif
+ if (forcon(nhpb+1).gt.0.0d0) then
+ nhpb=nhpb+1
+ if (ibecarb(i).gt.0) then
+ ihpb(i)=ihpb(i)+nres
+ jhpb(i)=jhpb(i)+nres
+ endif
+ if (dhpb(nhpb).eq.0.0d0)
+ & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ endif
+C endif
+ enddo
+ call hpb_partition
+ call flush(iout)
+ return
+ end
else
difi=0.0
endif
-cd write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
-cd & rad2deg*phi0(i), rad2deg*drange(i),
-cd & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
+ if (energy_dec) then
+ write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc",
+ & i,itori,rad2deg*phii,
+ & rad2deg*phi0(i), rad2deg*drange(i),
+ & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
+ endif
enddo
cd write (iout,*) 'edihcnstr',edihcnstr
return
edihcnstr=edihcnstr+0.25d0*ftors*difi**4
gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
endif
-! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
-! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
+ write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih",
+ & i,itori,rad2deg*phii,
+ & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
enddo
! write (iout,*) 'edihcnstr',edihcnstr
return
else
difi=0.0d0
endif
+ write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih",
+ & i,itori,rad2deg*phii,
+ & rad2deg*difi,0.25d0*ftors*difi**4
c write (iout,'(2i5,4f10.5,e15.5)') i,itori,phii,phi0(i),difi,
c & drange(i),edihi
! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
projects / unres.git / commitdiff