include 'COMMON.CHAIN'
include 'COMMON.HEADER'
include 'COMMON.FFIELD'
- include 'COMMON.FREE'
include 'COMMON.INTERACT'
include "COMMON.SPLITELE"
include 'COMMON.SHIELD'
+ include 'COMMON.FREE'
+ include 'COMMON.LANGEVIN'
character*320 controlcard,ucase
#ifdef MPL
include 'COMMON.INFO'
min_var=(index(controlcard,'MINVAR').gt.0)
plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
- call readi(controlcard,'NCUT',ncut,1)
+ call readi(controlcard,'NCUT',ncut,0)
+ call readi(controlcard,'NCLUST',nclust,5)
call readi(controlcard,'SYM',symetr,1)
write (iout,*) 'sym', symetr
call readi(controlcard,'NSTART',nstart,0)
lgrp=(index(controlcard,'LGRP').gt.0)
caonly=(index(controlcard,'CA_ONLY').gt.0)
print_dist=(index(controlcard,'PRINT_DIST').gt.0)
- call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
+ if (ncut.gt.0)
+ & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
call readi(controlcard,'IOPT',iopt,2)
lside = index(controlcard,"SIDE").gt.0
efree = index(controlcard,"EFREE").gt.0
print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
print_homology_models=
& index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
+ call readi(controlcard,'NSAXS',nsaxs,0)
+ call readi(controlcard,'SAXS_MODE',saxs_mode,0)
+ call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
+ call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
+ write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
+ & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
if (min_var) iopt=1
return
end
C Read weights of the subsequent energy terms.
call card_concat(weightcard)
write(iout,*) weightcard
-C call reada(weightcard,'WSC',wsc,1.0d0)
- write(iout,*) wsc
+ call reada(weightcard,'WSC',wsc,1.0d0)
call reada(weightcard,'WLONG',wsc,wsc)
call reada(weightcard,'WSCP',wscp,1.0d0)
call reada(weightcard,'WELEC',welec,1.0D0)
call reada(weightcard,'WTORD',wtor_d,1.0D0)
call reada(weightcard,'WANG',wang,1.0D0)
call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+ call reada(weightcard,'WSAXS',wsaxs,0.0D0)
call reada(weightcard,'SCAL14',scal14,0.4D0)
call reada(weightcard,'SCALSCP',scalscp,1.0d0)
call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
if (nstart.lt.nnt) nstart=nnt
if (nend.gt.nct .or. nend.eq.0) nend=nct
write (iout,*) "nstart",nstart," nend",nend
+ write (iout,*) "calling read_saxs_consrtr",nsaxs
+ if (nsaxs.gt.0) call read_saxs_constr
nres0=nres
if (constr_homology.gt.0) then
call read_constr_homology(print_homology_restraints)
return
end
+c-------------------------------------------------------------------------------
+ subroutine read_saxs_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'
+ double precision cm(3)
+c read(inp,*) nsaxs
+ write (iout,*) "Calling read_saxs nsaxs",nsaxs
+ call flush(iout)
+ if (saxs_mode.eq.0) then
+c SAXS distance distribution
+ do i=1,nsaxs
+ read(inp,*) distsaxs(i),Psaxs(i)
+ enddo
+ Cnorm = 0.0d0
+ do i=1,nsaxs
+ Cnorm = Cnorm + Psaxs(i)
+ enddo
+ write (iout,*) "Cnorm",Cnorm
+ do i=1,nsaxs
+ Psaxs(i)=Psaxs(i)/Cnorm
+ enddo
+ write (iout,*) "Normalized distance distribution from SAXS"
+ do i=1,nsaxs
+ write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
+ enddo
+ Wsaxs0=0.0d0
+ do i=1,nsaxs
+ Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
+ enddo
+ write (iout,*) "Wsaxs0",Wsaxs0
+ else
+c SAXS "spheres".
+ do i=1,nsaxs
+ read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
+ enddo
+ do j=1,3
+ cm(j)=0.0d0
+ enddo
+ do i=1,nsaxs
+ do j=1,3
+ cm(j)=cm(j)+Csaxs(j,i)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=cm(j)/nsaxs
+ enddo
+ do i=1,nsaxs
+ do j=1,3
+ Csaxs(j,i)=Csaxs(j,i)-cm(j)
+ enddo
+ enddo
+ write (iout,*) "SAXS sphere coordinates"
+ do i=1,nsaxs
+ write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
+ enddo
+ endif
+ return
+ end
c====-------------------------------------------------------------------
subroutine read_constr_homology(lprn)
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.INTERACT'
+ include 'COMMON.NAMES'
include 'COMMON.HOMRESTR'
c
c For new homol impl
integer ilen
external ilen
logical lprn
- logical unres_pdb
+ logical unres_pdb,liiflag
c
c FP - Nov. 2014 Temporary specifications for new vars
c
- double precision rescore_tmp,x12,y12,z12,rescore2_tmp
+ double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
+ & rescore3_tmp
double precision, dimension (max_template,maxres) :: rescore
double precision, dimension (max_template,maxres) :: rescore2
+ double precision, dimension (max_template,maxres) :: rescore3
character*24 tpl_k_rescore
c -----------------------------------------------------------------
c Reading multiple PDB ref structures and calculation of retraints
c
c New
c
- lim_theta=0
- lim_xx=0
c
c Reading HM global scores (prob not required)
c
tpl_k_rescore="template"//kic2//".sco"
unres_pdb=.false.
- call readpdb
+ call readpdb_template(k)
do i=1,2*nres
do j=1,3
crefjlee(j,i)=c(j,i)
c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
open (ientin,file=tpl_k_rescore,status='old')
if (nnt.gt.1) rescore(k,1)=0.0d0
- do irec=nnt,maxdim ! loop for reading res sim
+ do irec=nnt,nct ! loop for reading res sim
if (read2sigma) then
read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
- & idomain_tmp
+ & rescore3_tmp,idomain_tmp
i_tmp=i_tmp+nnt-1
idomain(k,i_tmp)=idomain_tmp
rescore(k,i_tmp)=rescore_tmp
rescore2(k,i_tmp)=rescore2_tmp
+ rescore3(k,i_tmp)=rescore3_tmp
+ write(iout,'(a7,i5,3f10.5,i5)') "rescore",
+ & i_tmp,rescore2_tmp,rescore_tmp,
+ & rescore3_tmp,idomain_tmp
else
idomain(k,irec)=1
read (ientin,*,end=1401) rescore_tmp
c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
c Instead of res sim other local measure of b/b str reliability possible
- sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
+ if (sigma_dih(k,i).ne.0)
+ & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
enddo
lim_dih=nct-nnt-2
sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
& rescore(k,i-2))/3.0
c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
- sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+ if (sigma_theta(k,i).ne.0)
+ & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
c rescore(k,i-2) ! right expression ?
c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
enddo
endif
- lim_theta=nct-nnt-1
if (waga_d.gt.0.0d0) then
c open (ientin,file=tpl_k_sigma_d,status='old')
c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
- sigma_d(k,i)=rescore(k,i) ! right expression ?
- sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+ sigma_d(k,i)=rescore3(k,i) ! right expression ?
+ if (sigma_d(k,i).ne.0)
+ & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
c read (ientin,*) sigma_d(k,i) ! 1st variant
- if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
enddo
- lim_xx=nct-nnt+1
endif
enddo
c
c
if (waga_dist.ne.0.0d0) then
ii=0
+ liiflag=.true.
do i=nnt,nct-2
do j=i+2,nct
ii=ii+1
- if (ii_in_use(ii).eq.0) then
- do ki=ii,lim_odl-1
- ires_homo(ki)=ires_homo(ki+1)
- jres_homo(ki)=jres_homo(ki+1)
- ii_in_use(ki)=ii_in_use(ki+1)
+ if (ii_in_use(ii).eq.0.and.liiflag) then
+ liiflag=.false.
+ iistart=ii
+ endif
+ if (ii_in_use(ii).ne.0.and..not.liiflag.or.
+ & .not.liiflag.and.ii.eq.lim_odl) then
+ if (ii.eq.lim_odl) then
+ iishift=ii-iistart+1
+ else
+ iishift=ii-iistart
+ endif
+ liiflag=.true.
+ do ki=iistart,lim_odl-iishift
+ ires_homo(ki)=ires_homo(ki+iishift)
+ jres_homo(ki)=jres_homo(ki+iishift)
+ ii_in_use(ki)=ii_in_use(ki+iishift)
do k=1,constr_homology
- odl(k,ki)=odl(k,ki+1)
- sigma_odl(k,ki)=sigma_odl(k,ki+1)
- l_homo(k,ki)=l_homo(k,ki+1)
+ odl(k,ki)=odl(k,ki+iishift)
+ sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
+ l_homo(k,ki)=l_homo(k,ki+iishift)
enddo
enddo
- ii=ii-1
- lim_odl=lim_odl-1
+ ii=ii-iishift
+ lim_odl=lim_odl-iishift
endif
enddo
enddo
& ki=1,constr_homology)
enddo
write (iout,*) "Dihedral angle restraints from templates"
- do i=nnt+3,lim_dih
- write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+ do i=nnt+3,nct
+ write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+ & (rad2deg*dih(ki,i),
& rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
enddo
write (iout,*) "Virtual-bond angle restraints from templates"
- do i=nnt+2,lim_theta
- write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+ do i=nnt+2,nct
+ write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
+ & (rad2deg*thetatpl(ki,i),
& rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
enddo
write (iout,*) "SC restraints from templates"
- do i=nnt,lim_xx
+ do i=nnt,nct
write(iout,'(i5,100(4f8.2,4x))') i,
& (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
& 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
projects / unres.git / blobdiff