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'
#endif
- integer i
-
+ integer i,i1,i2,it1,it2
+ double precision pi
read (INP,'(a80)') titel
call card_concat(controlcard)
call readi(controlcard,'RESCALE',rescale_mode,2)
call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
write (iout,*) "DISTCHAINMAX",distchainmax
+C Reading the dimensions of box in x,y,z coordinates
+ call reada(controlcard,'BOXX',boxxsize,100.0d0)
+ call reada(controlcard,'BOXY',boxysize,100.0d0)
+ call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+c Cutoff range for interactions
+ call reada(controlcard,"R_CUT",r_cut,15.0d0)
+ call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+ call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+ call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+ if (lipthick.gt.0.0d0) then
+ bordliptop=(boxzsize+lipthick)/2.0
+ bordlipbot=bordliptop-lipthick
+C endif
+ if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
+ & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+ buflipbot=bordlipbot+lipbufthick
+ bufliptop=bordliptop-lipbufthick
+ if ((lipbufthick*2.0d0).gt.lipthick)
+ &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+ endif
+ write(iout,*) "bordliptop=",bordliptop
+ write(iout,*) "bordlipbot=",bordlipbot
+ write(iout,*) "bufliptop=",bufliptop
+ write(iout,*) "buflipbot=",buflipbot
+C Shielding mode
+ call readi(controlcard,'SHIELD',shield_mode,0)
+ write (iout,*) "SHIELD MODE",shield_mode
+ if (shield_mode.gt.0) then
+ pi=3.141592d0
+C VSolvSphere the volume of solving sphere
+C print *,pi,"pi"
+C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
+C there will be no distinction between proline peptide group and normal peptide
+C group in case of shielding parameters
+ VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+ VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+ write (iout,*) VSolvSphere,VSolvSphere_div
+C long axis of side chain
+ do i=1,ntyp
+ long_r_sidechain(i)=vbldsc0(1,i)
+ short_r_sidechain(i)=sigma0(i)
+ enddo
+ buff_shield=1.0d0
+ endif
call readi(controlcard,'PDBOUT',outpdb,0)
call readi(controlcard,'MOL2OUT',outmol2,0)
refstr=(index(controlcard,'REFSTR').gt.0)
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
+ with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
+ write (iout,*) "with_theta_constr ",with_theta_constr
+ 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)
- 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
write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
lprint_cart=index(controlcard,"PRINT_CART") .gt.0
lprint_int=index(controlcard,"PRINT_INT") .gt.0
+ call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
+ write (iout,*) "with_homology_constr ",with_dihed_constr,
+ & " CONSTR_HOMOLOGY",constr_homology
+ print_homology_restraints=
+ & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0
+ print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0
+ print_homology_models=
+ & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0
+ read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").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
include 'COMMON.CONTROL'
include 'COMMON.CONTACTS'
include 'COMMON.TIME1'
+ include 'COMMON.TORCNSTR'
+ include 'COMMON.SHIELD'
#ifdef MPL
include 'COMMON.INFO'
#endif
double precision x(maxvar)
integer itype_pdb(maxres)
logical seq_comp
- integer i,j
+ integer i,j,kkk,i1,i2,it1,it2
C
C Body
C
C Read weights of the subsequent energy terms.
call card_concat(weightcard)
+ write(iout,*) weightcard
call reada(weightcard,'WSC',wsc,1.0d0)
call reada(weightcard,'WLONG',wsc,wsc)
call reada(weightcard,'WSCP',wscp,1.0d0)
call reada(weightcard,'WTURN4',wturn4,1.0D0)
call reada(weightcard,'WTURN6',wturn6,1.0D0)
call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
+ call reada(weightcard,'WSCCOR',wsccor,1.0D0)
call reada(weightcard,'WBOND',wbond,1.0D0)
call reada(weightcard,'WTOR',wtor,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)
call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
if (index(weightcard,'SOFT').gt.0) ipot=6
+ call reada(weightcard,"D0CM",d0cm,3.78d0)
+ call reada(weightcard,"AKCM",akcm,15.1d0)
+ call reada(weightcard,"AKTH",akth,11.0d0)
+ call reada(weightcard,"AKCT",akct,12.0d0)
+ call reada(weightcard,"V1SS",v1ss,-1.08d0)
+ call reada(weightcard,"V2SS",v2ss,7.61d0)
+ call reada(weightcard,"V3SS",v3ss,13.7d0)
+ call reada(weightcard,"EBR",ebr,-5.50D0)
+ call reada(weightcard,'WSHIELD',wshield,1.0d0)
+ write(iout,*) 'WSHIELD',wshield
+ call reada(weightcard,'WLT',wliptran,0.0D0)
+ call reada(weightcard,"ATRISS",atriss,0.301D0)
+ call reada(weightcard,"BTRISS",btriss,0.021D0)
+ call reada(weightcard,"CTRISS",ctriss,1.001D0)
+ call reada(weightcard,"DTRISS",dtriss,1.001D0)
+ write (iout,*) "ATRISS=", atriss
+ write (iout,*) "BTRISS=", btriss
+ write (iout,*) "CTRISS=", ctriss
+ write (iout,*) "DTRISS=", dtriss
+ dyn_ss=(index(weightcard,'DYN_SS').gt.0)
+ do i=1,maxres
+ dyn_ss_mask(i)=.false.
+ enddo
+ do i=1,maxres-1
+ do j=i+1,maxres
+ dyn_ssbond_ij(i,j)=1.0d300
+ enddo
+ enddo
+ call reada(weightcard,"HT",Ht,0.0D0)
+ if (dyn_ss) then
+ ss_depth=ebr/wsc-0.25*eps(1,1)
+ Ht=Ht/wsc-0.25*eps(1,1)
+ akcm=akcm*wstrain/wsc
+ akth=akth*wstrain/wsc
+ akct=akct*wstrain/wsc
+ v1ss=v1ss*wstrain/wsc
+ v2ss=v2ss*wstrain/wsc
+ v3ss=v3ss*wstrain/wsc
+ else
+ ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+ endif
+ write (iout,'(/a)') "Disulfide bridge parameters:"
+ write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
+ write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
+ write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
+ write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
+ & ' v3ss:',v3ss
+
C 12/1/95 Added weight for the multi-body term WCORR
call reada(weightcard,'WCORRH',wcorr,1.0D0)
if (wcorr4.gt.0.0d0) wcorr=wcorr4
weights(16)=wvdwpp
weights(17)=wbond
weights(18)=scal14
+ weights(19)=wsccor
write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
& wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
- & wturn4,wturn6
+ & wturn4,wturn6,wsccor
10 format (/'Energy-term weights (unscaled):'//
& 'WSCC= ',f10.6,' (SC-SC)'/
& 'WSCP= ',f10.6,' (SC-p)'/
& 'WCORR6= ',f10.6,' (multi-body 6th order)'/
& 'WTURN3= ',f10.6,' (turns, 3rd order)'/
& 'WTURN4= ',f10.6,' (turns, 4th order)'/
- & 'WTURN6= ',f10.6,' (turns, 6th order)')
+ & 'WTURN6= ',f10.6,' (turns, 6th order)'/
+ & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
+
if (wcorr4.gt.0.0d0) then
write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
& 'between contact pairs of peptide groups'
do i=1,nres
#ifdef PROCOR
- if (itype(i).eq.21 .or. itype(i+1).eq.21) then
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
#else
- if (itype(i).eq.21) then
+ if (itype(i).eq.ntyp1) then
#endif
itel(i)=0
#ifdef PROCOR
- else if (itype(i+1).ne.20) then
+ else if (iabs(itype(i+1)).ne.20) then
#else
- else if (itype(i).ne.20) then
+ else if (iabs(itype(i)).ne.20) then
#endif
itel(i)=1
else
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
+C read (inp,*) ftors
+C write (iout,*) 'FTORS',ftors
+C ftors is the force constant for torsional quartic constrains
+ read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
+ & i=1,ndih_constr)
+ write (iout,*)
+ & 'There are',ndih_constr,' constraints on phi angles.'
+ do i=1,ndih_constr
+ write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
+ & ftors(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
+ if (with_theta_constr) then
+C with_theta_constr is keyword allowing for occurance of theta constrains
+ read (inp,*) ntheta_constr
+C ntheta_constr is the number of theta constrains
+ if (ntheta_constr.gt.0) then
+C read (inp,*) ftors
+ read (inp,*) (itheta_constr(i),theta_constr0(i),
+ & theta_drange(i),for_thet_constr(i),
+ & i=1,ntheta_constr)
+C the above code reads from 1 to ntheta_constr
+C itheta_constr(i) residue i for which is theta_constr
+C theta_constr0 the global minimum value
+C theta_drange is range for which there is no energy penalty
+C for_thet_constr is the force constant for quartic energy penalty
+C E=k*x**4
+C if(me.eq.king.or..not.out1file)then
+ write (iout,*)
+ & 'There are',ntheta_constr,' constraints on phi angles.'
+ do i=1,ntheta_constr
+ write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
+ & theta_drange(i),
+ & for_thet_constr(i)
+ enddo
+C endif
+ do i=1,ntheta_constr
+ theta_constr0(i)=deg2rad*theta_constr0(i)
+ theta_drange(i)=deg2rad*theta_drange(i)
+ enddo
+C if(me.eq.king.or..not.out1file)
+C & write (iout,*) 'FTORS',ftors
+C do i=1,ntheta_constr
+C ii = itheta_constr(i)
+C thetabound(1,ii) = phi0(i)-drange(i)
+C thetabound(2,ii) = phi0(i)+drange(i)
+C enddo
+ endif ! ntheta_constr.gt.0
+ endif! with_theta_constr
+
nnt=1
nct=nres
print *,'NNT=',NNT,' NCT=',NCT
- if (itype(1).eq.21) nnt=2
- if (itype(nres).eq.21) nct=nct-1
+ if (itype(1).eq.ntyp1) nnt=2
+ if (itype(nres).eq.ntyp1) nct=nct-1
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)
+ endif
+
c if (pdbref) then
c read(inp,'(a)') pdbfile
c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
nstart_sup=nnt
nstart_seq=nnt
nsup=nct-nnt+1
+ kkk=1
do i=1,2*nres
do j=1,3
- cref(j,i)=c(j,i)
+ cref(j,i,kkk)=c(j,i)
enddo
enddo
endif
call contact(.true.,ncont_ref,icont_ref)
endif
+ if (ns.gt.0) then
+C write (iout,'(/a,i3,a)')
+C & 'The chain contains',ns,' disulfide-bridging cysteines.'
+ write (iout,'(20i4)') (iss(i),i=1,ns)
+ if (dyn_ss) then
+ write(iout,*)"Running with dynamic disulfide-bond formation"
+ else
+ write (iout,'(/a/)') 'Pre-formed links are:'
+ do i=1,nss
+ i1=ihpb(i)-nres
+ i2=jhpb(i)-nres
+ it1=itype(i1)
+ it2=itype(i2)
+ write (iout,'(2a,i3,3a,i3,a,3f10.3)')
+ & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
+ & ebr,forcon(i)
+ enddo
+ write (iout,'(a)')
+ endif
+ endif
+ if (ns.gt.0.and.dyn_ss) then
+ do i=nss+1,nhpb
+ ihpb(i-nss)=ihpb(i)
+ jhpb(i-nss)=jhpb(i)
+ forcon(i-nss)=forcon(i)
+ dhpb(i-nss)=dhpb(i)
+ enddo
+ nhpb=nhpb-nss
+ nss=0
+ call hpb_partition
+ do i=1,ns
+ 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-----------------------------------------------------------------------------
do i=1,ns
if (itype(iss(i)).ne.1) then
write (iout,'(2a,i3,a)')
- & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
+ & 'Do you REALLY think that the residue ',
+ & restyp(itype(iss(i))),i,
& ' can form a disulfide bridge?!!!'
write (*,'(2a,i3,a)')
- & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
+ & 'Do you REALLY think that the residue ',
+ & restyp(itype(iss(i))),i,
& ' can form a disulfide bridge?!!!'
#ifdef MPL
call mp_stopall(error_msg)
enddo
write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
20 continue
- dhpb(i)=dbr
- forcon(i)=fbr
+C dhpb(i)=dbr
+C forcon(i)=fbr
enddo
do i=1,nss
ihpb(i)=ihpb(i)+nres
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'
open (isidep1,file=sidepname,status="old")
call getenv('SCCORPAR',sccorname)
open (isccor,file=sccorname,status="old")
+ call getenv('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old')
#ifndef OLDSCP
C
C 8/9/01 In the newest version SCp interaction constants are read from a file
#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./
+ logical normalize,next
+ integer restr_type
+ double precision xlink(4,0:4) /
+c a b c sigma
+ & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
+ & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
+ & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
+ & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
+ & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
+ write (iout,*) "Calling read_dist_constr"
+c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+c call flush(iout)
+ next=.true.
+
+ DO WHILE (next)
+
+ call card_concat(controlcard)
+ next = index(controlcard,"NEXT").gt.0
+ call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
+ write (iout,*) "restr_type",restr_type
+ call readi(controlcard,"NFRAG",nfrag_,0)
+ 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)
+ if (restr_type.eq.10)
+ & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
+ 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)
+ normalize = index(controlcard,"NORMALIZE").gt.0
+ write (iout,*) "WBOLTZD",wboltzd
+c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+c write (iout,*) "IFRAG"
+c do i=1,nfrag_
+c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+c enddo
+c write (iout,*) "IPAIR"
+c do i=1,npair_
+c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+c enddo
+ if (nfrag_.gt.0) then
+ nres0=nres
+ read(inp,'(a)') pdbfile
+ write (iout,*)
+ & "Distance restraints will be constructed from structure ",pdbfile
+ open(ipdbin,file=pdbfile,status='old',err=11)
+ call readpdb(.true.)
+ nres=nres0
+ close(ipdbin)
+ endif
+ 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)
+c call flush(iout)
+ if (wfrag_(i).eq.0.0d0) cycle
+ do j=ifrag_(1,i),ifrag_(2,i)-1
+ do k=j+1,ifrag_(2,i)
+c write (iout,*) "j",j," k",k
+ ddjk=dist(j,k)
+ if (restr_type.eq.1) then
+ nhpb=nhpb+1
+ irestr_type(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
+ irestr_type(nhpb)=1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ endif
+ else if (restr_type.eq.3) then
+ nhpb=nhpb+1
+ irestr_type(nhpb)=1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ endif
+ write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ enddo
+ enddo
+ enddo
+ do i=1,npair_
+ if (wpair_(i).eq.0.0d0) cycle
+ 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)
+ if (restr_type.eq.1) then
+ nhpb=nhpb+1
+ irestr_type(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
+ irestr_type(nhpb)=1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ endif
+ else if (restr_type.eq.3) then
+ nhpb=nhpb+1
+ irestr_type(nhpb)=1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ endif
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ enddo
+ enddo
+ enddo
+
+c print *,ndist_
+ write (iout,*) "Distance restraints as read from input"
+ do i=1,ndist_
+ if (restr_type.eq.11) then
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
+ & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
+c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
+ if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
+ nhpb=nhpb+1
+ irestr_type(nhpb)=11
+ write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+ & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
+ if (ibecarb(nhpb).gt.0) then
+ ihpb(nhpb)=ihpb(nhpb)+nres
+ jhpb(nhpb)=jhpb(nhpb)+nres
+ endif
+ else if (constr_dist.eq.10) then
+c Cross-lonk Markov-like potential
+ call card_concat(controlcard)
+ call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
+ call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
+ ibecarb(nhpb+1)=0
+ if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
+ if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
+ if (index(controlcard,"ZL").gt.0) then
+ link_type=1
+ else if (index(controlcard,"ADH").gt.0) then
+ link_type=2
+ else if (index(controlcard,"PDH").gt.0) then
+ link_type=3
+ else if (index(controlcard,"DSS").gt.0) then
+ link_type=4
+ else
+ link_type=0
+ endif
+ call reada(controlcard,"AXLINK",dhpb(nhpb+1),
+ & xlink(1,link_type))
+ call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
+ & xlink(2,link_type))
+ call reada(controlcard,"CXLINK",fordepth(nhpb+1),
+ & xlink(3,link_type))
+ call reada(controlcard,"SIGMA",forcon(nhpb+1),
+ & xlink(4,link_type))
+ call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
+c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
+c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
+ if (forcon(nhpb+1).le.0.0d0 .or.
+ & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
+ nhpb=nhpb+1
+ irestr_type(nhpb)=10
+ if (ibecarb(nhpb).gt.0) then
+ ihpb(nhpb)=ihpb(nhpb)+nres
+ jhpb(nhpb)=jhpb(nhpb)+nres
+ endif
+ write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+ & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
+ & irestr_type(nhpb)
+ else
+C print *,"in else"
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
+ & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
+ if (forcon(nhpb+1).gt.0.0d0) then
+ nhpb=nhpb+1
+ if (dhpb1(nhpb).eq.0.0d0) then
+ irestr_type(nhpb)=1
+ else
+ irestr_type(nhpb)=2
+ endif
+ if (ibecarb(nhpb).gt.0) then
+ ihpb(nhpb)=ihpb(nhpb)+nres
+ jhpb(nhpb)=jhpb(nhpb)+nres
+ endif
+ if (dhpb(nhpb).eq.0.0d0)
+ & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ endif
+ write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
+ endif
+C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+C if (forcon(nhpb+1).gt.0.0d0) then
+C nhpb=nhpb+1
+C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ enddo
+
+ ENDDO ! next
+
+ fordepthmax=0.0d0
+ if (normalize) then
+ do i=nss+1,nhpb
+ if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
+ & fordepthmax=fordepth(i)
+ enddo
+ do i=nss+1,nhpb
+ if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
+ enddo
+ endif
+ if (nhpb.gt.nss) then
+ write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
+ & "The following",nhpb-nss,
+ & " distance restraints have been imposed:",
+ & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
+ & " score"," type"
+ do i=nss+1,nhpb
+ write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
+ & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
+ & irestr_type(i)
+ enddo
+ endif
+ call hpb_partition
+ call flush(iout)
+ return
+ 11 write (iout,*)"read_dist_restr: error reading reference structure"
+ stop
+ 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 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.INTERACT'
+ include 'COMMON.NAMES'
+ include 'COMMON.HOMRESTR'
+c
+c For new homol impl
+c
+ include 'COMMON.VAR'
+c include 'include_unres/COMMON.VAR'
+c
+
+c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
+c & dist_cut
+c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
+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, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
+ integer idomain(max_template,maxres)
+ integer ilen
+ external ilen
+ logical lprn
+ logical unres_pdb,liiflag
+c
+c FP - Nov. 2014 Temporary specifications for new vars
+c
+ 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 not using pre-computed ones stored in files model_ki_{dist,angle}
+c FP (Nov., 2014)
+c -----------------------------------------------------------------
+c
+c
+c Alternative: reading from input
+#ifdef DEBUG
+ write (iout,*) "BEGIN READ HOMOLOGY INFO"
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+#endif
+ call card_concat(controlcard)
+ call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
+ call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
+ call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
+ call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
+ call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
+ call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
+ dist1cut=(index(controlcard,'DIST1CUT').gt.0)
+ call readi(controlcard,"HOMOL_NSET",homol_nset,1)
+ read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
+ if (homol_nset.gt.1)then
+ call readi(controlcard,"ISET",iset,1)
+ call card_concat(controlcard)
+ read(controlcard,*) (waga_homology(i),i=1,homol_nset)
+ else
+ iset=1
+ waga_homology(1)=1.0
+ endif
+c
+#ifdef DEBUG
+ write(iout,*) "read_constr_homology iset",iset
+ write(iout,*) "waga_homology(",iset,")",waga_homology(iset)
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+#endif
+cd write (iout,*) "nnt",nnt," nct",nct
+cd call flush(iout)
+
+
+ lim_odl=0
+ lim_dih=0
+c
+c New
+c
+c
+c Reading HM global scores (prob not required)
+c
+ do i = nnt,nct
+ do k=1,constr_homology
+ idomain(k,i)=0
+ enddo
+ enddo
+c open (4,file="HMscore")
+c do k=1,constr_homology
+c read (4,*,end=521) hmscore_tmp
+c hmscore(k)=hmscore_tmp ! Another transformation can be used
+c write(*,*) "Model", k, ":", hmscore(k)
+c enddo
+c521 continue
+
+ ii=0
+ do i = nnt,nct-2
+ do j=i+2,nct
+ ii=ii+1
+ ii_in_use(ii)=0
+ enddo
+ enddo
+c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+
+ if (read_homol_frag) then
+ call read_klapaucjusz
+ else
+ write (iout,*) "CONSTR_HOMOLOGY",constr_homology
+ do k=1,constr_homology
+
+ read(inp,'(a)') pdbfile
+c write (iout,*) "k ",k," pdbfile ",pdbfile
+c Next stament causes error upon compilation (?)
+c if(me.eq.king.or. .not. out1file)
+c write (iout,'(2a)') 'PDB data will be read from file ',
+c & pdbfile(:ilen(pdbfile))
+ write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
+ & pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=33)
+ goto 34
+ 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
+ & pdbfile(:ilen(pdbfile))
+ stop
+ 34 continue
+c print *,'Begin reading pdb data'
+c
+c Files containing res sim or local scores (former containing sigmas)
+c
+
+ write(kic2,'(bz,i2.2)') k
+
+ tpl_k_rescore="template"//kic2//".sco"
+
+ unres_pdb=.false.
+ call readpdb_template(k)
+ do i=1,2*nres
+ do j=1,3
+ crefjlee(j,i)=c(j,i)
+ enddo
+ enddo
+#ifdef DEBUG
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
+ & (crefjlee(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "READ HOMOLOGY INFO"
+ write (iout,*) "read_constr_homology x: after reading pdb file"
+ write (iout,*) "waga_homology(",iset,")",waga_homology(iset)
+ write (iout,*) "waga_dist",waga_dist
+ write (iout,*) "waga_angle",waga_angle
+ write (iout,*) "waga_theta",waga_theta
+ write (iout,*) "waga_d",waga_d
+ write (iout,*) "dist_cut",dist_cut
+#endif
+#ifdef AIX
+ call flush_(iout)
+#else
+ call flush(iout)
+#endif
+
+c
+c Distance restraints
+c
+c ... --> odl(k,ii)
+C Copy the coordinates from reference coordinates (?)
+ do i=1,2*nres
+ do j=1,3
+c c(j,i)=cref(j,i)
+c write (iout,*) "c(",j,i,") =",c(j,i)
+ enddo
+ enddo
+c
+c From read_dist_constr (commented out 25/11/2014 <-> res sim)
+c
+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,nct ! loop for reading res sim
+ if (read2sigma) then
+ read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_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 rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
+ rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
+c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
+ endif
+ enddo
+ 1401 continue
+ close (ientin)
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ do i = nnt,nct-2
+ do j=i+2,nct
+
+ x12=c(1,i)-c(1,j)
+ y12=c(2,i)-c(2,j)
+ z12=c(3,i)-c(3,j)
+ distal=dsqrt(x12*x12+y12*y12+z12*z12)
+c write (iout,*) k,i,j,distal,dist2_cut
+ if (dist1cut .and. k.gt.1) then
+ ii=ii+1
+ if (l_homo(1,ii)) then
+ ii_in_use(ii)=1
+ l_homo(k,ii)=.true.
+ ires_homo(ii)=i
+ jres_homo(ii)=j
+ odl(k,ii)=distal
+ sigma_odl(k,ii)=sigma_odl(1,ii)
+ else
+ l_homo(k,ii)=.false.
+ endif
+ else
+ if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+ & .and. distal.le.dist2_cut ) then
+
+ ii=ii+1
+ ii_in_use(ii)=1
+ l_homo(k,ii)=.true.
+
+c write (iout,*) "k",k
+c write (iout,*) "i",i," j",j," constr_homology",
+c & constr_homology
+ ires_homo(ii)=i
+ jres_homo(ii)=j
+ odl(k,ii)=distal
+ if (read2sigma) then
+ sigma_odl(k,ii)=0
+ do ik=i,j
+ sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
+ enddo
+ sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
+ if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
+ & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+ else
+ if (odl(k,ii).le.dist_cut) then
+ sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
+ else
+#ifdef OLDSIGMA
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
+#else
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+#endif
+ endif
+ endif
+ sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
+ else
+ ii=ii+1
+ l_homo(k,ii)=.false.
+ endif
+ endif
+ enddo
+ enddo
+ lim_odl=ii
+ endif
+c
+c Theta, dihedral and SC retraints
+c
+ if (waga_angle.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_dih,status='old')
+c do irec=1,maxres-3 ! loop for reading sigma_dih
+c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
+c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
+c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_dih(k,i+nnt-1)
+c enddo
+c1402 continue
+c close (ientin)
+ do i = nnt+3,nct
+ if (idomain(k,i).eq.0) then
+ sigma_dih(k,i)=0.0
+ cycle
+ endif
+ dih(k,i)=phiref(i) ! right?
+c read (ientin,*) sigma_dih(k,i) ! original variant
+c write (iout,*) "dih(",k,i,") =",dih(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c & "rescore(",k,i-1,") =",rescore(k,i-1),
+c & "rescore(",k,i-2,") =",rescore(k,i-2),
+c & "rescore(",k,i-3,") =",rescore(k,i-3)
+
+ sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2)+rescore(k,i-3))/4.0
+c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
+c write (iout,*) "Raw sigmas for dihedral angle restraints"
+c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
+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
+ 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
+ endif
+
+ if (waga_theta.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_theta,status='old')
+c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
+c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
+c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_theta(k,i+nnt-1)
+c enddo
+c1403 continue
+c close (ientin)
+
+ do i = nnt+2,nct ! right? without parallel.
+c do i = i=1,nres ! alternative for bounds acc to readpdb?
+c do i=ithet_start,ithet_end ! with FG parallel.
+ if (idomain(k,i).eq.0) then
+ sigma_theta(k,i)=0.0
+ cycle
+ endif
+ thetatpl(k,i)=thetaref(i)
+c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
+c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+c & "rescore(",k,i-1,") =",rescore(k,i-1),
+c & "rescore(",k,i-2,") =",rescore(k,i-2)
+c read (ientin,*) sigma_theta(k,i) ! 1st variant
+ 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
+ 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
+
+ if (waga_d.gt.0.0d0) then
+c open (ientin,file=tpl_k_sigma_d,status='old')
+c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
+c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
+c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
+c & sigma_d(k,i+nnt-1)
+c enddo
+c1404 continue
+
+ do i = nnt,nct ! right? without parallel.
+c do i=2,nres-1 ! alternative for bounds acc to readpdb?
+c do i=loc_start,loc_end ! with FG parallel.
+ if (itype(i).eq.10) cycle
+ if (idomain(k,i).eq.0 ) then
+ sigma_d(k,i)=0.0
+ cycle
+ endif
+ xxtpl(k,i)=xxref(i)
+ yytpl(k,i)=yyref(i)
+ zztpl(k,i)=zzref(i)
+c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
+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)=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
+ enddo
+ endif
+ enddo
+c
+c remove distance restraints not used in any model from the list
+c shift data in all arrays
+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.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+iishift)
+ sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
+ l_homo(k,ki)=l_homo(k,ki+iishift)
+ enddo
+ enddo
+ ii=ii-iishift
+ lim_odl=lim_odl-iishift
+ endif
+ enddo
+ enddo
+ endif
+
+ endif ! .not. klapaucjusz
+
+ if (constr_homology.gt.0) call homology_partition
+ if (constr_homology.gt.0) call init_int_table
+cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
+cd & "lim_xx=",lim_xx
+c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+c
+c Print restraints
+c
+ if (.not.lprn) return
+cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write (iout,*) "Distance restraints from templates"
+ do ii=1,lim_odl
+ write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
+ & ii,ires_homo(ii),jres_homo(ii),
+ & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
+ & ki=1,constr_homology)
+ enddo
+ write (iout,*) "Dihedral angle restraints from templates"
+ 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,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,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)
+ enddo
+ endif
+c -----------------------------------------------------------------
+ return
+ end
+c----------------------------------------------------------------------
+ subroutine read_klapaucjusz
+
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.INTERACT'
+ include 'COMMON.NAMES'
+ include 'COMMON.HOMRESTR'
+ character*256 fragfile
+ integer ninclust(maxclust),inclust(max_template,maxclust),
+ & nresclust(maxclust),iresclust(maxres,maxclust)
+
+ character*2 kic2
+ character*24 model_ki_dist, model_ki_angle
+ character*500 controlcard
+ integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
+ integer idomain(max_template,maxres)
+ logical lprn /.true./
+ integer ilen
+ external ilen
+ logical unres_pdb,liiflag
+c
+c
+ double precision rescore_tmp,x12,y12,z12,rescore2_tmp
+ double precision, dimension (max_template,maxres) :: rescore
+ double precision, dimension (max_template,maxres) :: rescore2
+ character*24 tpl_k_rescore
+
+c
+c For new homol impl
+c
+ include 'COMMON.VAR'
+c
+ double precision chomo(3,maxres2+2,max_template)
+ call getenv("FRAGFILE",fragfile)
+ open(ientin,file=fragfile,status="old",err=10)
+ read(ientin,*) constr_homology,nclust
+ l_homo = .false.
+ sigma_theta=0.0
+ sigma_d=0.0
+ sigma_dih=0.0
+c Read pdb files
+ do k=1,constr_homology
+ read(ientin,'(a)') pdbfile
+ write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
+ & pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=33)
+ goto 34
+ 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
+ & pdbfile(:ilen(pdbfile))
+ stop
+ 34 continue
+ unres_pdb=.false.
+ call readpdb_template(k)
+ do i=1,2*nres
+ do j=1,3
+ chomo(j,i,k)=c(j,i)
+ enddo
+ enddo
+ do i=1,nres
+ rescore(k,i)=0.2d0
+ rescore2(k,i)=1.0d0
+ enddo
+ enddo
+c Read clusters
+ do i=1,nclust
+ read(ientin,*) ninclust(i),nresclust(i)
+ read(ientin,*) (inclust(k,i),k=1,ninclust(i))
+ read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
+ enddo
+c
+c Loop over clusters
+c
+ do l=1,nclust
+ do ll = 1,ninclust(l)
+
+ k = inclust(ll,l)
+ do i=1,nres
+ idomain(k,i)=0
+ enddo
+ do i=1,nresclust(l)
+ if (nnt.gt.1) then
+ idomain(k,iresclust(i,l)+1) = 1
+ else
+ idomain(k,iresclust(i,l)) = 1
+ endif
+ enddo
+c
+c Distance restraints
+c
+c ... --> odl(k,ii)
+C Copy the coordinates from reference coordinates (?)
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=chomo(j,i,k)
+c write (iout,*) "c(",j,i,") =",c(j,i)
+ enddo
+ enddo
+ call int_from_cart1(.false.)
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ do i = nnt,nct-2
+ do j=i+2,nct
+
+ x12=c(1,i)-c(1,j)
+ y12=c(2,i)-c(2,j)
+ z12=c(3,i)-c(3,j)
+ distal=dsqrt(x12*x12+y12*y12+z12*z12)
+c write (iout,*) k,i,j,distal,dist2_cut
+
+ if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+ & .and. distal.le.dist2_cut ) then
+
+ ii=ii+1
+ ii_in_use(ii)=1
+ l_homo(k,ii)=.true.
+
+c write (iout,*) "k",k
+c write (iout,*) "i",i," j",j," constr_homology",
+c & constr_homology
+ ires_homo(ii)=i
+ jres_homo(ii)=j
+ odl(k,ii)=distal
+ if (read2sigma) then
+ sigma_odl(k,ii)=0
+ do ik=i,j
+ sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
+ enddo
+ sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
+ if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
+ & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+ else
+ if (odl(k,ii).le.dist_cut) then
+ sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
+ else
+#ifdef OLDSIGMA
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
+#else
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
+ & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+#endif
+ endif
+ endif
+ sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
+ else
+ ii=ii+1
+c l_homo(k,ii)=.false.
+ endif
+ enddo
+ enddo
+ lim_odl=ii
+ endif
+c
+c Theta, dihedral and SC retraints
+c
+ if (waga_angle.gt.0.0d0) then
+ do i = nnt+3,nct
+ if (idomain(k,i).eq.0) then
+c sigma_dih(k,i)=0.0
+ cycle
+ endif
+ dih(k,i)=phiref(i)
+ sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2)+rescore(k,i-3))/4.0
+c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
+c & " sigma_dihed",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))
+ enddo
+ lim_dih=nct-nnt-2
+ endif
+
+ if (waga_theta.gt.0.0d0) then
+ do i = nnt+2,nct
+ if (idomain(k,i).eq.0) then
+c sigma_theta(k,i)=0.0
+ cycle
+ endif
+ thetatpl(k,i)=thetaref(i)
+ sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
+ & rescore(k,i-2))/3.0
+ if (sigma_theta(k,i).ne.0)
+ & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+ enddo
+ endif
+
+ if (waga_d.gt.0.0d0) then
+ do i = nnt,nct
+ if (itype(i).eq.10) cycle
+ if (idomain(k,i).eq.0 ) then
+c sigma_d(k,i)=0.0
+ cycle
+ endif
+ xxtpl(k,i)=xxref(i)
+ yytpl(k,i)=yyref(i)
+ zztpl(k,i)=zzref(i)
+ sigma_d(k,i)=rescore(k,i)
+ if (sigma_d(k,i).ne.0)
+ & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
+ if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
+ enddo
+ endif
+ enddo ! l
+ enddo ! ll
+c
+c remove distance restraints not used in any model from the list
+c shift data in all arrays
+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.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+iishift)
+ sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
+ l_homo(k,ki)=l_homo(k,ki+iishift)
+ enddo
+ enddo
+ ii=ii-iishift
+ lim_odl=lim_odl-iishift
+ endif
+ enddo
+ enddo
+ endif
+
+ return
+ 10 stop "Error infragment file"
+ end
+c----------------------------------------------------------------------
projects / unres.git / blobdiff