use io_clust
!#define CLUSTER
use io_units
- use io_base, only: permut
+ use io_base, only: permut, read_dist_constr,gen_dist_constr2
use geometry_data, only: nres,theta,phi,alph,omeg,&
c,cref
use energy_data, only: nnt,nct
! write (iout,*) "Before permut"
! write (iout,*) "symetr", symetr
! call flush(iout)
- call permut(symetr)
+! call permut(symetr)
! write (iout,*) "after permut"
! call flush(iout)
print *,'MAIN: nnt=',nnt,' nct=',nct
write(iout,*) "after daread_ccoords"
ind1=0
DO I=1,NCON_work-1
+! print *,"Tu1",i
if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
do k=1,2*nres
do l=1,3
enddo
enddo
DO J=I+1,NCON_work
+! print *,"tu2",j
IND=IOFFSET(NCON_work,I,J)
#ifdef MPI
if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
real(kind=8) :: przes(3),obrot(3,3)
real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2)
integer :: i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
- integer :: iaperm,ibezperm,run
+ integer :: iaperm,ibezperm,run,nequiv(50)
real(kind=8) :: rms,rmsmina
-! write (iout,*) "tu dochodze"
+! print *, "tu dochodze",symetr
rmsmina=10d10
nperm=1
do i=1,symetr
nperm=i*nperm
enddo
! write (iout,*) "nperm",nperm
- call permut(symetr)
+ if (symetr.ne.1) call permut(nequiv(i),nperm,tabperm)
+! call permut(symetr)
+! print *,"after permut"
! write (iout,*) "tabperm", tabperm(1,1)
+ if (symetr.ne.1) then
do kkk=1,nperm
if (lside) then
ii=0
if (non_conv) print *,non_conv,icon,jcon
if (rmsmina.gt.rms) rmsmina=rms
enddo
+ else
+! print *,"in else difconf"
+ nstart=2
+ kkk=1
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=allcart(j,i,jcon)
+ enddo
+ enddo
+ call fitsq(rms,c(1,nstart),cref(1,nstart,kkk),nend-nstart+1,&
+ przes,&
+ obrot,non_conv)
+ rmsmina=rms
+ endif
difconf=dsqrt(rmsmina)
return
end function difconf
!c if (constr_homology) tole=dmax1(tole,1.5d0)
write (iout,*) "with_homology_constr ",with_dihed_constr,&
" CONSTR_HOMOLOGY",constr_homology
+ write(iout,*) "CONSTR_DIST=",constr_dist
read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
indpdb,constr_dist,raw_psipred, with_theta_constr
use geometry, only: chainbuild,alloc_geo_arrays
use energy, only: alloc_ener_arrays
- use io_base, only: rescode
+ use io_base, only: rescode,read_dist_constr,gen_dist_constr2
use control, only: setup_var,init_int_table
use conform_compar, only: contact
! implicit none
!C enddo
endif ! ntheta_constr.gt.0
endif! with_theta_constr
+ if (constr_dist.gt.0) then
+ call read_dist_constr
+
+ link_start=1
+ link_end=nhpb
+ clink_start=1
+ clink_end=cnhpb
+ endif
+
+
+
if (constr_homology.gt.0) then
!c write (iout,*) "About to call read_constr_homology"
!c call flush(iout)
10 stop "Error in fragment file"
end subroutine read_klapaucjusz
-!-----------------------------------------------------------------------------
+
+
+! subroutine read_dist_constr
+! use MPI_data
+!! use control
+! use geometry, only: dist
+! use geometry_data
+! use control_data
+! use energy_data
+!! implicit real*8 (a-h,o-z)
+!! include 'DIMENSIONS'
+!#ifdef MPI
+! include 'mpif.h'
+!#endif
+!! include 'COMMON.SETUP'
+!! include 'COMMON.CONTROL'
+!! include 'COMMON.CHAIN'
+!! include 'COMMON.IOUNITS'
+!! include 'COMMON.SBRIDGE'
+! integer,dimension(2,100) :: ifrag_,ipair_
+! real(kind=8),dimension(100) :: wfrag_,wpair_
+! character(len=640) :: controlcard
+!
+!!el local variables
+! integer :: i,k,j,ii,jj,itemp,mnumkk,mnumjj,k1,j1
+! integer :: nfrag_,npair_,ndist_
+! real(kind=8) :: dist_cut,ddjk
+!
+!! write (iout,*) "Calling read_dist_constr"
+!! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+!! call flush(iout)
+! if(.not.allocated(cihpb)) allocate(cihpb(maxdim))
+! if(.not.allocated(cjhpb)) allocate(cjhpb(maxdim))
+! if(.not.allocated(cdhpb)) allocate(cdhpb(maxdim))
+! if(.not.allocated(cdhpb1)) allocate(cdhpb1(maxdim))
+! if(.not.allocated(cforcon)) allocate(cforcon(maxdim))
+! if(.not.allocated(cfordepth)) allocate(cfordepth(maxdim))
+! if(.not.allocated(cibecarb)) allocate(cibecarb(maxdim))
+! if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
+! call gen_dist_constr2
+! go to 1712
+! endif
+! cnhpb=0
+! call card_concat(controlcard,.true.)
+! 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
+!! if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
+!! if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
+!! if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
+!! if(.not.allocated(forcon)) allocate(forcon(maxdim))
+!
+! call flush(iout)
+! do i=1,nfrag_
+! 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,nres
+! j1=j
+! k1=k
+! if (j.gt.nres) j1=j-nres
+! if (k.gt.nres) k1=k-nres
+! mnumkk=molnum(k1)
+! mnumjj=molnum(j1)
+!
+! if ((itype(k1,mnumkk).gt.ntyp_molec(mnumkk)).or.&
+! (itype(j1,mnumjj).gt.ntyp_molec(mnumjj))) cycle
+! write (iout,*) "j",j," k",k,itype(k1,mnumkk),itype(j1,mnumjj)
+! ddjk=dist(j,k)
+! if (constr_dist.eq.1) then
+! cnhpb=cnhpb+1
+! cihpb(cnhpb)=j
+! cjhpb(cnhpb)=k
+! cdhpb(cnhpb)=ddjk
+! cforcon(cnhpb)=wfrag_(i)
+! else if (constr_dist.eq.2) then
+! if (ddjk.le.dist_cut) then
+! print *,"tu",nhpb,i,j,k,maxdim
+! cnhpb=cnhpb+1
+! cihpb(cnhpb)=j
+! cjhpb(cnhpb)=k
+! cdhpb(cnhpb)=ddjk
+! cforcon(cnhpb)=wfrag_(i)
+! endif
+! else
+! cnhpb=cnhpb+1
+! cihpb(cnhpb)=j
+! cjhpb(cnhpb)=k
+! cdhpb(cnhpb)=ddjk
+! cforcon(cnhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+! endif
+!#ifdef MPI
+! if (.not.out1file .or. me.eq.king) &
+! write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
+! cnhpb,cihpb(cnhpb),cjhpb(cnhpb),cdhpb(cnhpb),cforcon(cnhpb)
+!#else
+! write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
+! cnhpb,cihpb(cnhpb),cjhpb(cnhpb),cdhpb(cnhpb),cforcon(cnhpb)
+!#endif
+! 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)
+! cnhpb=cnhpb+1
+! cihpb(cnhpb)=j
+! cjhpb(cnhpb)=k
+! cforcon(cnhpb)=wpair_(i)
+! cdhpb(cnhpb)=dist(j,k)
+!#ifdef MPI
+! if (.not.out1file .or. me.eq.king) &
+! write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+! nhpb,cihpb(nhpb),cjhpb(nhpb),cdhpb(nhpb),cforcon(nhpb)
+!#else
+! write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+! nhpb,cihpb(nhpb),cjhpb(nhpb),cdhpb(nhpb),cforcon(nhpb)
+!#endif
+! enddo
+! enddo
+! endif
+! enddo
+! do i=1,ndist_
+!! read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+!! if (forcon(nhpb+1).gt.0.0d0) then
+!! nhpb=nhpb+1
+!! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+! if (constr_dist.eq.11) then
+! read (inp,*) cihpb(cnhpb+1),cjhpb(cnhpb+1),cdhpb(i),cdhpb1(i), &
+! cibecarb(i),cforcon(cnhpb+1),cfordepth(cnhpb+1)
+!!EL fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
+! cfordepth(cnhpb+1)=cfordepth(cnhpb+1)**(0.25d0)
+! cforcon(cnhpb+1)=cforcon(cnhpb+1)**(0.25d0)
+! else
+! print *,"in else"
+! read (inp,*) cihpb(cnhpb+1),cjhpb(cnhpb+1),cdhpb(i),cdhpb1(i), &
+! cibecarb(i),cforcon(cnhpb+1)
+! print *,"afterread",cihpb(cnhpb+1)
+! endif
+! if (cforcon(cnhpb+1).gt.0.0d0) then
+! cnhpb=cnhpb+1
+! if (cibecarb(i).gt.0) then
+! cihpb(cnhpb)=cihpb(cnhpb)+nres
+! cjhpb(cnhpb)=cjhpb(cnhpb)+nres
+! endif
+! if (cdhpb(cnhpb).eq.0.0d0) &
+! cdhpb(cnhpb)=dist(ihpb(cnhpb),jhpb(cnhpb))
+! endif
+!
+!#ifdef MPI
+! if (.not.out1file .or. me.eq.king) &
+! write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+! cnhpb,cihpb(cnhpb),cjhpb(cnhpb),cdhpb(cnhpb),cforcon(cnhpb)
+!#else
+! write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+! cnhpb,cihpb(cnhpb),cjhpb(cnhpb),cdhpb(cnhpb),cforcon(cnhpb)
+!#endif
+! enddo
+! 1712 continue
+! call flush(iout)
+! return
+! end subroutine read_dist_constr
+!
+!! 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
+!! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+! subroutine gen_dist_constr2
+! use MPI_data
+!! use control
+! use geometry, only: dist
+! use geometry_data
+! use control_data
+! use energy_data
+! integer :: i,j
+! real(kind=8) :: distance
+! if (constr_dist.eq.11) then
+! do i=nstart_sup,nstart_sup+nsup-1
+! do j=i+2,nstart_sup+nsup-1
+! distance=dist(i,j)
+! if (distance.le.15.0) then
+! cnhpb=cnhpb+1
+! cihpb(nhpb)=i+nstart_seq-nstart_sup
+! cjhpb(nhpb)=j+nstart_seq-nstart_sup
+! cforcon(nhpb)=sqrt(0.04*distance)
+! cfordepth(nhpb)=sqrt(40.0/distance)
+! cdhpb(nhpb)=distance-0.1d0
+! cdhpb1(nhpb)=distance+0.1d0
+!
+!#ifdef MPI
+! if (.not.out1file .or. me.eq.king) &
+! write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
+! cnhpb,cihpb(nhpb),cjhpb(nhpb),cdhpb(nhpb),cforcon(nhpb)
+!#else
+! write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
+! cnhpb,cihpb(nhpb),cjhpb(nhpb),cdhpb(nhpb),cforcon(nhpb)
+!#endif
+! endif
+! enddo
+! enddo
+! else
+! do i=nstart_sup,nstart_sup+nsup-1
+! do j=i+2,nstart_sup+nsup-1
+! cnhpb=cnhpb+1
+! cihpb(nhpb)=i+nstart_seq-nstart_sup
+! cjhpb(nhpb)=j+nstart_seq-nstart_sup
+! cforcon(nhpb)=weidis
+! cdhpb(nhpb)=dist(i,j)
+! enddo
+! enddo
+! endif
+! return
+! end subroutine gen_dist_constr2
+!
+!!-----------------------------------------------------------------------------
end module io_clust
nperm=i*nperm
enddo
c write (iout,*) "nperm",nperm
- call permut(symetr)
+! call permut(symetr)
+ call permut(nequiv(i),nperm,tabperm)
c write (iout,*) "tabperm", tabperm(1,1)
do kkk=1,nperm
if (lside) then
enddo
call int_from_cart1(.false.)
! call etotal(energia(0),fT)
- call etotal(energia)
+ call etotal(energia(0))
totfree(i)=energia(0)
! write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
! write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
etors=enetb(13,i)
etors_d=enetb(14,i)
ehpb=enetb(15,i)
+! print *,"ehpb",ehpb
! estr=enetb(18,i)
estr=enetb(17,i)
! esccor=enetb(19,i)
set (CMAKE_Fortran_FLAGS_RELEASE " ")
set (CMAKE_Fortran_FLAGS_DEBUG "-O0 -g -traceback")
# set(FFLAGS0 "-fpp -c -CB -g -ip " )
- set(FFLAGS0 "-O3 -ip -fpp -mcmodel=large" )
+ set(FFLAGS0 "-CB -g -ip -fpp -mcmodel=large" )
# set(FFLAGS0 "-O0 -CB -CA -g" )
set(FFLAGS1 "-fpp -c -O " )
set(FFLAGS2 "-fpp -c -g -CA -CB ")
integer :: i,j,k,iti,mnum,term
real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
mag1,mag2,v(3)
-!#ifdef DEBUG
+#ifdef DEBUG
write (iout,*) "Velocities, kietic"
do i=0,nres
write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
(d_t(j,i+nres),j=1,3)
enddo
-!#endif
+#endif
KEt_p=0.0d0
KEt_sc=0.0d0
! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
do i=nnt,term
mnum=molnum(i)
if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
- write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
+! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
if (mnum.gt.4) then
do j=1,3
v(j)=incr(j)+0.5d0*d_t(j,i)
incr(j)=incr(j)+d_t(j,i)
enddo
enddo
- write(iout,*) 'KEt_p', KEt_p
+! write(iout,*) 'KEt_p', KEt_p
! The translational part for the side chain virtual bond
! Only now we can initialize incr with zeros. It must be equal
! to the velocities of the first Calpha.
v(j)=incr(j)+d_t(j,nres+i)
enddo
endif
- write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
+! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
enddo
enddo
! goto 111
- write(iout,*) 'KEt_sc', KEt_sc
+! write(iout,*) 'KEt_sc', KEt_sc
! The part due to stretching and rotation of the peptide groups
KEr_p=0.0D0
do i=nnt,nct-1
do j=1,3
incr(j)=d_t(j,i)
enddo
- write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3),KEr_p,Ip(mnum)
+! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3),KEr_p,Ip(mnum)
KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
+incr(3)*incr(3))
enddo
! goto 111
- write(iout,*) 'KEr_p', KEr_p
+! write(iout,*) 'KEr_p', KEr_p
! The rotational part of the side chain virtual bond
KEr_sc=0.0D0
do i=nnt,nct
do j=1,3
incr(j)=d_t(j,nres+i)
enddo
- write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
+! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
incr(3)*incr(3))
endif
enddo
! The total kinetic energy
111 continue
- write(iout,*) 'KEr_sc', KEr_sc
+! write(iout,*) 'KEr_sc', KEr_sc
KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
- write (iout,*) "KE_total",KE_total
+! write (iout,*) "KE_total",KE_total
return
end subroutine kinetic
!-----------------------------------------------------------------------------
innt=iposd_chain(ichain)
! innt_org=
innt_org=chain_border(1,ichain)
+! print *,"befor molnum0",innt_org
if ((molnum(innt_org).eq.5).or.(molnum(innt_org).eq.4)) go to 137
if(.not.allocated(ghalf)) print *,"COCO"
if(.not.allocated(Ghalf)) allocate(Ghalf(1300*(1300+1)/2))
write(iout,*) i,ghalf(i)
enddo
#endif
- mnum=molnum(innt+1)
+! print *,"befor molnum0.5",innt_org+1
+ mnum=molnum(innt_org+1)
call gldiag(1300*2,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
write(iout,*) "after internal",ierr,iwork
#ifdef DEBUG
do i=1,n
do k=1,3
ii=ii+1
+! print *, "before molnum",innt_org+1
mnum=molnum(innt_org+1)
if (molnum(innt_org+1).ge.4) geigen(i)=3.0/msc(itype(innt_org+i-1,mnum),mnum)
!should not it be other way around
do j=1,3
ind=ind+1
d_t(j,i)=d_t_work(ind)
- enddo
+ enddo !j
+! print *, "before molnum2",i
mnum=molnum(i)
+ print *, "after molnum2",i
if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum).and.mnum.le.2) then
do j=1,3
ind=ind+1
d_t(j,i+nres)=d_t_work(ind)
- enddo
+ enddo !j
endif
- enddo
+ enddo !i
enddo
+ print *,"before large1"
if (large) then
write (iout,*)
write (iout,*) "Random velocities in the Calpha,SC space"
restyp(itype(i,mnum),mnum),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
enddo
endif
+ print *,"before kineticCASC"
call kinetic_CASC(Ek1)
+ print *,"after kineticCASC"
+
!
! Transform the velocities to virtual-bond space
!
enddo
ibeg=inct+1
do i=innt,inct
+! print *,"before molnum3",i
mnum=molnum(i)
if (iabs(itype(i,1).eq.10).or.mnum.ge.3) then
!c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
enddo
- print *,"after msc",Im
+! print *,"after msc",Im
do i=nnt,nct-1
mnum=molnum(i)
mnum1=molnum(i+1)
Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*&
vbld(i+1)*vbld(i+1)*0.25d0
enddo
- print *,"after Ip",Im
+! print *,"after Ip",Im
do i=nnt,nct
mnum=molnum(i)
mnum1=molnum(i+1)
dc_norm(2,inres))*vbld(inres)*vbld(inres)
Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*&
dc_norm(3,inres))*vbld(inres)*vbld(inres)
- print *,i,inres,vbld(inres),iti,dc_norm(1,inres),Im(1,1)
+! print *,i,inres,vbld(inres),iti,dc_norm(1,inres),Im(1,1)
endif
enddo
- print *,"after ISC",Im
- print *,"before angnom",L
- print *,"before angnom",Im
+! print *,"after ISC",Im
+! print *,"before angnom",L
+! print *,"before angnom",Im
call angmom(cm,L)
Im(2,1)=Im(1,2)
Im(3,1)=Im(1,3)
enddo
enddo
!c Finding the eigenvectors and eignvalues of the inertia tensor
- print *,"before djacob",Imcp,eigvec,eigval
+! print *,"before djacob",Imcp,eigvec,eigval
call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
- print *, "after djacob"
+! print *, "after djacob"
do i=1,3
if (dabs(eigval(i)).gt.1.0d-15) then
Id(i,i)=1.0d0/eigval(i)
"WLT "," "," ","WTUBE " ,&
"WVDWPPNUCL","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",&
"WELSB ","WBOND_NUCL","WANG_NUCL ","WSBLOC ","WTOR_NUCL ",&
- "WTORD_NUCL","WCORR_NUCL","WCORR3_NUC","WNULL ","WNULL ",&
+ "WTORD_NUCL","WCORR_NUCL","WCORR3NUCL","WNULL ","WNULL ",&
"WCATPROT ","WCATCAT ","WNULL ","WNULL ","WNULL ",&
"WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO ","WCATNUCL ",&
"H_CONS ","WNULL ","WNULL ","WNULL ","WNULL ",&
! print *,"before",ees,evdw1,ecorr
! write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
if (nres_molec(2).gt.0) then
+! write(iout,*) "welpsb2,",welpsb
call ebond_nucl(estr_nucl)
call ebend_nucl(ebe_nucl)
call etor_nucl(etors_nucl)
wtor=weights(13)*fact(1)
wtor_d=weights(14)*fact(2)
wsccor=weights(21)*fact(1)
- welpsb=weights(28)*fact(1)
+ welpsb=weights(29)*fact(1)
wcorr_nucl= weights(37)*fact(1)
wcorr3_nucl=weights(38)*fact(2)
wtor_nucl= weights(35)*fact(1)
!c Restraints from contact prediction
dd=dist(ii,jj)
if (constr_dist.eq.11) then
- ehpb=ehpb+fordepth(i)**4.0d0 &
+ ehpb=ehpb+cfordepth(i)**4.0d0 &
*rlornmr1(dd,cdhpb(i),cdhpb1(i),cforcon(i))
- fac=fordepth(i)**4.0d0 &
+ fac=cfordepth(i)**4.0d0 &
*rlornmr1prim(dd,cdhpb(i),cdhpb1(i),cforcon(i))/dd
if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, &
ehpb,cfordepth(i),dd
else
if (cdhpb1(i).gt.0.0d0) then
ehpb=ehpb+2*cforcon(i)*gnmr1(dd,cdhpb(i),cdhpb1(i))
- fac=forcon(i)*gnmr1prim(dd,cdhpb(i),cdhpb1(i))/dd
+ fac=cforcon(i)*gnmr1prim(dd,cdhpb(i),cdhpb1(i))/dd
!c write (iout,*) "beta nmr",
!c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
else
if (constr_dist.eq.11) then
ehpb=ehpb+cfordepth(i)**4.0d0 &
- *rlornmr1(dd,dhpb(i),cdhpb1(i),cforcon(i))
- fac=fordepth(i)**4.0d0 &
+ *rlornmr1(dd,cdhpb(i),cdhpb1(i),cforcon(i))
+ fac=cfordepth(i)**4.0d0 &
*rlornmr1prim(dd,cdhpb(i),cdhpb1(i),cforcon(i))/dd
if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, &
- ehpb,fordepth(i),dd
+ ehpb,cfordepth(i),dd
else
if (cdhpb1(i).gt.0.0d0) then
ehpb=ehpb+2*cforcon(i)*gnmr1(dd,cdhpb(i),cdhpb1(i))
- fac=forcon(i)*gnmr1prim(dd,cdhpb(i),cdhpb1(i))/dd
+ fac=cforcon(i)*gnmr1prim(dd,cdhpb(i),cdhpb1(i))/dd
!c write (iout,*) "alph nmr",
!c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
else
waga=cforcon(i)
!C Calculate the contribution to energy.
ehpb=ehpb+waga*rdis*rdis
- ! if (energy_dec)
-! write (iout,'(a6,2i5,5f10.3)') "edisip",ii,jj, &
-! ehpb,dd,cdhpb(i),waga,rdis
+ if (energy_dec) &
+ write (iout,'(a6,2i5,5f10.3)') "edisip",ii,jj, &
+ ehpb,dd,cdhpb(i),waga,rdis
!c write (iout,*) "alpha reg",dd,waga*rdis*rdis
!C
!c AL 3/30/2022 handle the cases of an isolated-residue chain
if (i.eq.nnt .and. itype(i+1,1).eq.ntyp1) cycle
if (i.eq.nct .and. itype(i-1,1).eq.ntyp1) cycle
+ if (itype(i-1,1).eq.ntyp1 .and. itype(i+1,1).eq.ntyp1) cycle
+! if (i.eq.nct .and. itype(i-1,1).eq.ntyp1) cycle
+
! costtab(i+1) =dcos(theta(i+1))
if (it.eq.10) goto 1
#ifdef SC_END
*cossc1)*gradene
gsclocx(:,i)=gsclocx(:,i)+vbld_inv(i+nres)*&
(dC_norm(:,i-1)-dC_norm(:,i+nres)*cossc1)*gradene
-#ifdef ENERGY_DEC
+!#ifdef ENERGY_DEC
if (energy_dec) write (2,'(2hC ,a3,i6,2(a,f10.5))')&
restyp(iti,1),i," angle",rad2deg*dacos(cossc1)," escloc",sumene
-#endif
+!#endif
else if (i.eq.nnt .or. itype(i-1,1).eq.ntyp1) then
!c AL 3/30/2022 handle a sidechain of a loose N-end
cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
*cossc)*gradene
gsclocx(:,i)=gsclocx(:,i)+vbld_inv(i+nres)*&
(dC_norm(:,i)-dC_norm(:,i+nres)*cossc)*gradene
-#ifdef ENERGY_DEC
- if (energy_dec) write (2,'(2hN ,a3,i6,2(a,f10.5))')
- & restyp(iti),i," angle",rad2deg*dacos(cossc)," escloc",sumene
-#endif
+!#ifdef ENERGY_DEC
+ if (energy_dec) write (2,'(2hN ,a3,i6,2(a,f10.5))')&
+ restyp(iti,1),i," angle",rad2deg*dacos(cossc)," escloc",sumene
+!#endif
else
#endif
!
!(maxres)
allocate(jcont_hb(maxconts,nres))
#endif
+ allocate(jcont_hb(maxconts,nres))
allocate(num_cont_hb(nres))
!(maxconts,maxres)
! common /rotat/
! print *,"tu?",i,istart_nucl(i,iint),iend_nucl(i,iint)
do j=istart_nucl(i,iint),iend_nucl(i,iint)
ind=ind+1
-! print *,"JESTEM"
+! print *,"JESTEM",i,j
itypj=itype(j,2)
if (itypj.eq.ntyp1_molec(2)) cycle
dscj_inv=vbld_inv(j+nres)
!C Calculate angular part of the gradient.
call sc_grad_nucl
call eelsbij(eelij,num_conti2)
- if (energy_dec .and. &
- (j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2)) &
- write (istat,'(e14.5)') evdwij
+! if (energy_dec .and. &
+! (j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2)) &
+! write (istat,'(e14.5)') evdwij
eelsb=eelsb+eelij
enddo ! j
enddo ! iint
!C 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+facfac-fac1
- if (energy_dec) then
- if(j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2) &
- write (istat,'(2a1,i4,1x,2a1,i4,4f10.5,3e12.5,$)') &
- sugartyp(istype(i)),restyp(itypi,2),i,sugartyp(istype(j)),&
- restyp(itypj,2),j,1.0d0/rij,cosa,cosb,cosg,fac*r3ij, &
- (4.0D0+facfac-fac1)*r6ij,(2.0d0-2.0d0*facfac+fac1)*r6ij
- write (iout,'(a6,2i5,e15.3)') 'ees',i,j,eesij
- endif
+! if (energy_dec) then
+! if(j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2) &
+! write (istat,'(2a1,i4,1x,2a1,i4,4f10.5,3e12.5,$)') &
+! sugartyp(istype(i)),restyp(itypi,2),i,sugartyp(istype(j)),&
+! restyp(itypj,2),j,1.0d0/rij,cosa,cosb,cosg,fac*r3ij, &
+! (4.0D0+facfac-fac1)*r6ij,(2.0d0-2.0d0*facfac+fac1)*r6ij
+! write (iout,'(a6,2i5,e15.3)') 'ees',i,j,eesij
+! endif
!C
!C Calculate contributions to the Cartesian gradient.
gelsbc(k,j)=gelsbc(k,j)+ggg(k)
gelsbc(k,i)=gelsbc(k,i)-ggg(k)
enddo
-! IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and.
- IF ( j.gt.i+1 .and.&
- num_conti.le.maxcont) THEN
+ IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and. &
+ ( j.gt.i+1 .and. num_conti2.le.maxcont)) THEN
!C
!C Calculate the contact function. The ith column of the array JCONT will
!C contain the numbers of atoms that make contacts with the atom I (of numbers
num_conti=num_conti+1
num_conti2=num_conti2+1
- if (num_conti.gt.maxconts) then
+ if (num_conti2.gt.maxconts) then
write (iout,*) 'WARNING - max. # of contacts exceeded;',&
' will skip next contacts for this conf.',maxconts
else
- jcont_hb(num_conti,i)=j
+ jcont_hb(num_conti2,i)=j
!c write (iout,*) "num_conti",num_conti,
!c & " jcont_hb",jcont_hb(num_conti,i)
!C Calculate contact energies
ecosbm=ecosb1-ecosb2
ecosgm=ecosg1-ecosg2
!C End diagnostics
- facont_hb(num_conti,i)=fcont
+ facont_hb(num_conti2,i)=fcont
fprimcont=fprimcont/rij
do k=1,3
gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
gggm(2)=gggm(2)+ees0mijp*yj
gggm(3)=gggm(3)+ees0mijp*zj
!C Derivatives due to the contact function
- gacont_hbr(1,num_conti,i)=fprimcont*xj
- gacont_hbr(2,num_conti,i)=fprimcont*yj
- gacont_hbr(3,num_conti,i)=fprimcont*zj
+ gacont_hbr(1,num_conti2,i)=fprimcont*xj
+ gacont_hbr(2,num_conti2,i)=fprimcont*yj
+ gacont_hbr(3,num_conti2,i)=fprimcont*zj
do k=1,3
!c
!c Gradient of the correlation terms
!c
- gacontp_hb1(k,num_conti,i)= &
+ gacontp_hb1(k,num_conti2,i)= &
(ecosap*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres)) &
+ ecosbp*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
- gacontp_hb2(k,num_conti,i)= &
+ gacontp_hb2(k,num_conti2,i)= &
(ecosap*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres)) &
+ ecosgp*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
- gacontp_hb3(k,num_conti,i)=gggp(k)
- gacontm_hb1(k,num_conti,i)= &
+ gacontp_hb3(k,num_conti2,i)=gggp(k)
+ gacontm_hb1(k,num_conti2,i)= &
(ecosam*(dc_norm(k,j+nres)-cosa*dc_norm(k,i+nres)) &
+ ecosbm*(erij(k)-cosb*dc_norm(k,i+nres)))*vbld_inv(i+nres)
- gacontm_hb2(k,num_conti,i)= &
+ gacontm_hb2(k,num_conti2,i)= &
(ecosam*(dc_norm(k,i+nres)-cosa*dc_norm(k,j+nres))&
+ ecosgm*(erij(k)-cosg*dc_norm(k,j+nres)))*vbld_inv(j+nres)
- gacontm_hb3(k,num_conti,i)=gggm(k)
+ gacontm_hb3(k,num_conti2,i)=gggm(k)
enddo
endif
endif
! double precision x(maxvar)
character(len=256) :: pdbfile
character(len=800) :: weightcard
- character(len=80) :: weightcard_t!,ucase
+ character(len=80) :: weightcard_t,chartest!,ucase
! integer,dimension(:),allocatable :: itype_pdb !(maxres)
! common /pizda/ itype_pdb
logical :: fail !seq_comp,
'Scaling factor of 1,4 SC-p interactions:',scal14
write (iout,'(a,f8.3)') &
'General scaling factor of SC-p interactions:',scalscp
+ write(iout,*) "welpsb",welpsb
endif
r0_corr=cutoff_corr-delt_corr
do i=1,ntyp
enddo
endif
endif
- print *,"CZY TU DOCHODZE"
+ print *,"CZY TU DOCHODZE",nucleic
if (indpdb.eq.0) then
nres_molec(:)=0
allocate(sequence(maxres,5))
if (nucleic) then
! Read sequence if not taken from the pdb file.
molec=2
+ print *,"jest before read nres"
read (inp,*) nres_molec(molec)
print *,'nres=',nres_molec(molec)
! allocate(sequence(maxres,5))
enddo
do i=2,nres-1
alph(i)=110d0*deg2rad
- if (molnum(i).eq.2) alph(i)=30d0*deg2rad
+ if (molnum(i).eq.2) alph(i)=60d0*deg2rad
enddo
do i=2,nres-1
omeg(i)=-120d0*deg2rad
- if (molnum(i).eq.2) omeg(i)=60d0*deg2rad
+ if (molnum(i).eq.2) omeg(i)=30d0*deg2rad
if (itype(i,1).le.0) omeg(i)=-omeg(i)
enddo
call chainbuild
read (inp,*) cihpb(cnhpb+1),cjhpb(cnhpb+1),cdhpb(i),cdhpb1(i), &
cibecarb(i),cforcon(cnhpb+1),cfordepth(cnhpb+1)
!EL fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
- fordepth(nhpb+1)=cfordepth(cnhpb+1)**(0.25d0)
- forcon(nhpb+1)=cforcon(cnhpb+1)**(0.25d0)
+ cfordepth(cnhpb+1)=cfordepth(cnhpb+1)**(0.25d0)
+ cforcon(cnhpb+1)=cforcon(cnhpb+1)**(0.25d0)
else
-!C print *,"in else"
+ print *,"in else"
read (inp,*) cihpb(cnhpb+1),cjhpb(cnhpb+1),cdhpb(i),cdhpb1(i), &
cibecarb(i),cforcon(cnhpb+1)
+ print *,"afterread",cihpb(cnhpb+1)
endif
- if (cforcon(nhpb+1).gt.0.0d0) then
+ if (cforcon(cnhpb+1).gt.0.0d0) then
cnhpb=cnhpb+1
- if (ibecarb(i).gt.0) then
- cihpb(i)=ihpb(i)+nres
- cjhpb(i)=jhpb(i)+nres
+ if (cibecarb(i).gt.0) then
+ cihpb(cnhpb)=cihpb(cnhpb)+nres
+ cjhpb(cnhpb)=cjhpb(cnhpb)+nres
endif
- if (cdhpb(nhpb).eq.0.0d0) &
+ if (cdhpb(cnhpb).eq.0.0d0) &
cdhpb(cnhpb)=dist(ihpb(cnhpb),jhpb(cnhpb))
endif
do i=nnt,nct
! write(iout,*), "coord",c(1,i),c(2,i),c(3,i)
iti=itype(i,molnum(i))
- print *,i,molnum(i)
+! print *,i,molnum(i)
if (molnum(i+1).eq.0) then
iti1=ntyp1_molec(molnum(i))
else
COS45 = .70710678
S45SQ = 0.50
C45SQ = 0.50
- print *,"in math",N,NMAX
+! print *,"in math",N,NMAX
! UNIT EIGENVECTOR MATRIX
DO 70 I = 1,N
DO 7 J = I,N
1000 FORMAT (/1X,'NONCONVERGENT JACOBI. LARGEST OFF-DIAGONAL ELE',&
'MENT = ',1PE14.7)
! ARRANGEMENT OF EIGENVALUES IN ASCENDING ORDER
- print *,"before first AJJ?"
+! print *,"before first AJJ?"
300 DO 14 I=1,N
14 AJJ(I)=A(I,I)
LT=N+1
17 AIIMIN=AJJ(I)
IT=I
16 CONTINUE
- print *,IT,"IT"
+! print *,IT,"IT"
IN=L
AII(IN)=AIIMIN
AJJ(IT)=1.0E+30
!c set default parameters for the line search
!c
if (stpmax .eq. 0.0d0) stpmax = 5.0d0
+ if (stpmax .eq. 0.0d0) stpmax = 2.0d0
stpmin = 1.0d-16
cappa = 0.9d0
slpmax = 10000.0d0
blank = ' '
if (stpmin .eq. 0.0d0) stpmin = 1.0d-16
if (stpmax .eq. 0.0d0) stpmax = 2.0d0
+ if (stpmax .eq. 0.0d0) stpmax = 1.5d0
if (cappa .eq. 0.0d0) cappa = 0.1d0
if (slpmax .eq. 0.0d0) slpmax = 10000.0d0
if (angmax .eq. 0.0d0) angmax = 180.0d0
integer :: errmsg_count,maxerrmsg_count=100000
!el real(kind=8) :: rmsnat,gyrate
!el external rmsnat,gyrate
- real(kind=8) :: tole=0.0d0
+ real(kind=8) :: tole=0.5d0
integer i,itj,ii,iii,j,k,l,licz
integer ir,ib,ipar,iparm
integer iscor,islice
wtor_d,wsccor,wbond
#endif
call restore_parm(iparm)
+ print *,"after restore parm"
+ print *,beta_h(ib,ipar)
call rescale_weights(1.0d0/(beta_h(ib,ipar)*1.987D-3))
#ifdef DEBUG
write (iout,*) "before etot w=",1.0d0/(beta_h(ib,ipar)*1.987D-3)
wtor_d,wsccor,wbond
#endif
! call etotal(energia(0),fT)
+ print *,"before etotal"
call etotal(energia(0))
+ print *,"after etotal"
!write(iout,*)"check c and dc after etotal",1.0d0/(0.001987*beta_h(ib,ipar))
!do k=1,2*nres+2
!write(iout,*)k,"c=",(c(l,k),l=1,3)
!
use energy_data
use geometry_data, only:nres,deg2rad,c,dc,nres_molec,crefjlee,cref
- use control_data, only:iscode,dyn_ss,pdbref,indpdb
+ use control_data, only:iscode,pdbref,indpdb
use io_base, only:rescode
use control, only:setup_var,init_int_table,hpb_partition
use geometry, only:alloc_geo_arrays
allocate(ww(max_eneW))
do i=1,n_eneW
key = wname(i)(:ilen(wname(i)))
+
call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
+ write(iout,*) "NAMES ",key(:ilen(key)),ww(i),i
enddo
write (iout,*) "iparm",iparm," myparm",myparm
wcatcat=ww(42)
wcatprot=ww(41)
wcorr3_nucl=ww(38)
+ write(iout,*) "CZY TU BLAD", ww(38)
wcorr_nucl=ww(37)
wtor_d_nucl=ww(36)
wtor_nucl=ww(35)
use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
scelemode,TUBEmode,tor_mode,energy_dec,r_cut_ang,r_cut_mart,&
- rlamb_mart
+ rlamb_mart,constr_dist
use energy_data, only:distchainmax,tubeR0,tubecenter,dyn_ss,constr_homology
use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
use work_part
!
use io_units
- use control_data, only:indpdb
+ use control_data, only:indpdb,constr_dist
#ifdef MPI
use mpi_data
! use mpi_
logical :: punch_dist,print_rms,caonly,verbose,merge_helices,&
bxfile,cxfile,histfile,entfile,zscfile,rmsrgymap,&
with_dihed_constr,check_conf,histout
- integer :: icomparfunc,pdbint,ensembles,constr_dist,oldion,shield_mode
+ integer :: icomparfunc,pdbint,ensembles,oldion,shield_mode
!---------------------------------------------------------------------------
! COMMON.OBCINKA
! common /obcinka/
projects / unres4.git / commitdiff