& constr_homology,homol_nset,iset
real*8 waga_homology
real*8 waga_dist, waga_angle,waga_theta, waga_d, dist_cut
+ & ,dist2_cut
logical refstr,pdbref,punch_dist,print_dist,caonly,lside,
& lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx,
& with_dihed_constr,out1file,print_homology_restraints,
& print_contact_map,print_homology_models
+ & ,read2sigma,l_homo
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,constr_dist,
& with_dihed_constr, constr_homology,homol_nset,out1file,
- & print_contact_map
+ & print_contact_map,read2sigma,l_homo(max_template,maxdim)
common /cntrlr/ waga_homology(10),
- & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,iset,
+ & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,dist2_cut,
+ & iset,
& print_homology_restraints,print_homology_models
& WSHET(MAXRES,MAXRES), EDFABET,
& CK(4),SCK(4),S1(4),S2(4)
c & ,DFAEXP(15001),
+
+ DATA CK/1.0D0,1.58740105197D0,2.08008382305D0,2.51984209979D0/
+ DATA SCK/1.0D0,1.25992104989D0,1.44224957031D0,1.58740105197D0/
+ DATA S1/3.75D0,5.75D0,7.75D0,9.75D0/
+ DATA S2/4.25D0,6.25D0,8.25D0,10.25D0/
parameter (maxprocs=16)
C Max. number of AA residues
integer maxres,maxres2
- parameter (maxres=600)
+ parameter (maxres=800)
C Appr. max. number of interaction sites
parameter (maxres2=2*maxres)
C Max. number of variables
- block data dfadata
- include 'DIMENSIONS'
- include 'COMMON.INTERACT'
- include 'COMMON.DFA'
- DATA CK/1.0D0,1.58740105197D0,2.08008382305D0,2.51984209979D0/
- DATA SCK/1.0D0,1.25992104989D0,1.44224957031D0,1.58740105197D0/
- DATA S1/3.75D0,5.75D0,7.75D0,9.75D0/
- DATA S2/4.25D0,6.25D0,8.25D0,10.25D0/
- end
-c------------------------------------------------------------------------
-
subroutine init_dfa_vars
include 'DIMENSIONS'
c write(*,*) 'case1:',tmp_n
- t1dx=t1dx+0.0d0
- t1dy=t1dy+0.0d0
- t1dz=t1dz+0.0d0
+cc t1dx=t1dx+0.0d0
+cc t1dy=t1dy+0.0d0
+cc t1dz=t1dz+0.0d0
t2dx(j)=0.0d0
t2dy(j)=0.0d0
t2dz(j)=0.0d0
tmp_n=tmp_n+1.0d0
c write(*,*) 'case4:',tmp_n
- t1dx=t1dx+0.0d0
- t1dy=t1dy+0.0d0
- t1dz=t1dz+0.0d0
+cc t1dx=t1dx+0.0d0
+cc t1dy=t1dy+0.0d0
+cc t1dz=t1dz+0.0d0
t2dx(j)=0.0d0
t2dy(j)=0.0d0
t2dz(j)=0.0d0
double precision enebet
enebet=0.0d0
- bx=0.0d0;by=0.0d0;bz=0.0d0
- shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
+c bx=0.0d0;by=0.0d0;bz=0.0d0
+c shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
gdfab=0.0d0
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc******************************************************************************
c real*8 dfaexp(15001)
real*8 uum, uup
real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
+ real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
+ common /shetf/ shetfx,shetfy,shetfz
+
common /sheca/ bx,by,bz
common /phys1/ inb,nmax,iselect
common /kyori2/ dis
common /sinsin/ sth
real*8 r_pair_mat(maxca,maxca)
+ real*8 e_gcont,fprim_gcont,de_gcont
ci integer istrand(maxca,maxca)
ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
ci common /shetest/ istrand,istrand_p,istrand_m
do i=1,inb-7
do j=i+4,inb-3
+
+ if (dis(i,j).lt.dfa_cutoff) then
+ call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
+
ip=i+1
ipp=i+2
jp=j+1
vbet(i,j)=uup+uum
vbetp=vbetp+uup
vbetm=vbetm+uum
- vbeta=vbeta+vbet(i,j)
+ vbeta=vbeta+vbet(i,j)*e_gcont
+
+
+ if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then
+c gradient correction from gcont
+ de_gcont=vbet(i,j)*fprim_gcont/dis(i,j)
+ shetfx(i)=shetfx(i) + de_gcont*rx(i,j)
+ shetfy(i)=shetfy(i) + de_gcont*ry(i,j)
+ shetfz(i)=shetfz(i) + de_gcont*rz(i,j)
+
+ shetfx(j)=shetfx(j) - de_gcont*rx(i,j)
+ shetfy(j)=shetfy(j) - de_gcont*ry(i,j)
+ shetfz(j)=shetfz(j) - de_gcont*rz(i,j)
+
+c energy correction from gcont
+ vbet(i,j)=vbet(i,j)*e_gcont
+ vbetap(i,j)=vbetap(i,j)*e_gcont
+ vbetap1(i,j)=vbetap1(i,j)*e_gcont
+ vbetap2(i,j)=vbetap2(i,j)*e_gcont
+ vbetam(i,j)=vbetam(i,j)*e_gcont
+ vbetam1(i,j)=vbetam1(i,j)*e_gcont
+ vbetam2(i,j)=vbetam2(i,j)*e_gcont
+ endif
+
ci elseif(istrand(i,j).eq.0)then
ci vbet(i,j)=0
c write(*,*) 'vbetm:',vbetm
c write(*,*) 'vbet(i,j):',vbet(i,j)
c stop
-
+
+ else
+ vbetap(i,j)=0
+ vbetap1(i,j)=0
+ vbetap2(i,j)=0
+ vbetam(i,j)=0
+ vbetam1(i,j)=0
+ vbetam2(i,j)=0
+ vbet(i,j)=0
+ endif
enddo
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbet(maxca,maxca)
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
+ real*8 e_gcont,fprim_gcont
C--------------------------------------------------------------------------------
do i=4,inb-4
im3=i-3
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
+ real*8 e_gcont,fprim_gcont
c********************************************************************************
do i=3,inb-5
imm=i-2
im=i-1
do j=i+2,inb-3
+
+ if (dis(imm,j).lt.dfa_cutoff) then
+ call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
jp=j+1
jpp=j+2
! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
! shefz(i,5)=shefz(i,5)
! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
-
+
+ endif
ci endif
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
+ real*8 e_gcont,fprim_gcont
C********************************************************************************
do i=2,inb-6
ip=i+1
im=i-1
do j=i+3,inb-3
+
+ if (dis(im,j).lt.dfa_cutoff) then
+ call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
jp=j+1
jpp=j+2
! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
! shefz(i,6)=shefz(i,6)
! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
-
+
+ endif
ci endif
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8
+ real*8 e_gcont,fprim_gcont
C********************************************************************************
do j=7,inb-1
jm=j-1
jmm=j-2
do i=1,j-6
+
+ if (dis(i,jmm).lt.dfa_cutoff) then
+ call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
ip=i+1
ipp=i+2
! shefz(j,11)=shefz(j,11)
! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
+ endif
ci endif
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
+ real*8 e_gcont,fprim_gcont
!c*************************************************************************c
do j=6,inb-2
jp=j+1
jm=j-1
do i=1,j-5
+
+ if (dis(i,jm).lt.dfa_cutoff) then
+ call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
ip=i+1
ipp=i+2
shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
$ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
+ endif
+
ci endif
ENDDO
dij=dist(i,j)
c write (iout,*) "dij(",i,j,") =",dij
do k=1,constr_homology
+ if(.not.l_homo(k,ii)) cycle
distance(k)=odl(k,ii)-dij
c write (iout,*) "distance(",k,") =",distance(k)
c
endif
enddo
- min_odl=minval(distancek)
+c min_odl=minval(distancek)
+ do kk=1,constr_homology
+ if(l_homo(kk,ii)) then
+ min_odl=distancek(kk)
+ exit
+ endif
+ enddo
+ do kk=1,constr_homology
+ if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl)
+ & min_odl=distancek(kk)
+ enddo
c write (iout,* )"min_odl",min_odl
#ifdef DEBUG
write (iout,*) "ij dij",i,j,dij
c Nie wiem po co to liczycie jeszcze raz!
c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/
c & (2*(sigma_odl(i,j,k))**2))
+ if(.not.l_homo(k,ii)) cycle
if (waga_dist.ge.0.0d0) then
c
c For Gaussian-type Urestr
c & *waga_dist)+min_odl
c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist
c
+ if(.not.l_homo(k,ii)) cycle
if (waga_dist.ge.0.0d0) then
c For Gaussian-type Urestr
c
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+ if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
esccor=0.0D0
do i=itau_start,itau_end
esccor_ii=0.0D0
+ if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
phii=phi(i)
character*2 kic2
character*24 model_ki_dist, model_ki_angle
character*500 controlcard
- integer ki, i, j, k, l
+ 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
c
c FP - Nov. 2014 Temporary specifications for new vars
c
- double precision rescore_tmp,x12,y12,z12
+ 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 Reading multiple PDB ref structures and calculation of retraints
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)
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)
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 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
write (iout,*) "CONSTR_HOMOLOGY",constr_homology
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)') 'Error opening PDB file.'
+ 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
+ & pdbfile(:ilen(pdbfile))
stop
34 continue
c print *,'Begin reading pdb data'
write(kic2,'(bz,i2.2)') k
tpl_k_rescore="template"//kic2//".sco"
-c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
-c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
-c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
-c tpl_k_sigma_d="template"//kic2//".sigma_d"
unres_pdb=.false.
call readpdb
c
c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
open (ientin,file=tpl_k_rescore,status='old')
- do irec=1,maxdim ! loop for reading res sim
- if (irec.eq.1) then
- rescore(k,irec)=0.0d0
- goto 1301
- endif
- read (ientin,*,end=1401) rescore_tmp
+ if (nnt.gt.1) rescore(k,1)=0.0d0
+ do irec=nnt,maxdim ! loop for reading res sim
+ if (read2sigma) then
+ read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_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
+ 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
+ rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
- 1301 continue
+ endif
enddo
1401 continue
- close (ientin)
-c open (ientin,file=tpl_k_sigma_odl,status='old')
-c do irec=1,maxdim ! loop for reading sigma_odl
-c read (ientin,*,end=1401) i, j,
-c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
-c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
-c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
-c enddo
-c 1401 continue
-c close (ientin)
+ close (ientin)
if (waga_dist.ne.0.0d0) then
ii=0
- do i = nnt,nct-2 ! right? without parallel.
- do j=i+2,nct ! right?
-c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
-c do j=i+2,nres ! ibid
-c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
-c do j=i+2,nct ! ibid
+ 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
-c
-c Attempt to replace dist(i,j) by its definition in ...
-c
- 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)
odl(k,ii)=distal
-c
-c odl(k,ii)=dist(i,j)
-c write (iout,*) "dist(",i,j,") =",dist(i,j)
-c write (iout,*) "distal = ",distal
-c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
-c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
-c & "rescore(",k,j,") =",rescore(k,j)
-c
-c Calculation of sigma from res sim
-c
-c if (odl(k,ii).le.6.0d0) then
-c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
-c Other functional forms possible depending on odl(k,ii), eg.
-c
- if (odl(k,ii).le.dist_cut) then
- sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
-c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
- else
+ 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))* ! sigma ~ rescore ~ error
+ 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))* ! sigma ~ rescore ~ error
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
& dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
#endif
-c Following expr replaced by a positive exp argument
-c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
-c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
-
-c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
-c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
+ 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
-c
- sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
-c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
-c
-c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
-c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
enddo
-c read (ientin,*) sigma_odl(k,ii) ! 1st variant
enddo
-c lim_odl=ii
-c if (constr_homology.gt.0) call homology_partition
+ lim_odl=ii
endif
c
c Theta, dihedral and SC retraints
c enddo
c1402 continue
c close (ientin)
- do i = nnt+3,nct ! right? without parallel.
-c do i=1,nres ! alternative for bounds acc to readpdb?
-c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
-c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
+ 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 & "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 ! right expression ?
-c
+ & 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 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))
c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
- if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
-c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
enddo
+ lim_dih=nct-nnt-2
endif
if (waga_theta.gt.0.0d0) then
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-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 ! right expression ?
+ & 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))
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)
- if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
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 & sigma_d(k,i+nnt-1)
c enddo
c1404 continue
- close (ientin)
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) goto 1 ! right?
+ 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 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?
- 1 continue
enddo
+ lim_xx=nct-nnt+1
endif
- close(ientin)
enddo
- if (waga_dist.ne.0.0d0) lim_odl=ii
+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
+ 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)
+ 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)
+ enddo
+ enddo
+ ii=ii-1
+ lim_odl=lim_odl-1
+ endif
+ enddo
+ enddo
+ endif
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,
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,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
- & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
+ 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,lim_dih
- write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+ write (iout,'(i5,100(2f8.2,4x))') 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,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+ write (iout,'(i5,100(2f8.2,4x))') 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
- write(iout,'(i5,10(4f8.2,4x))') i,
+ 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
c -----------------------------------------------------------------
return
end
-
-
+c----------------------------------------------------------------------
& ensembles,constr_dist,constr_homology,homol_nset,
& iset,ihset
real*8 waga_homology
- real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut
+ real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut,
+ & dist2_cut
logical refstr,pdbref,punch_dist,print_rms,caonly,verbose,
& merge_helices,bxfile,cxfile,histfile,entfile,zscfile,
- & rmsrgymap,with_dihed_constr,check_conf,histout,out1file
+ & rmsrgymap,with_dihed_constr,check_conf,histout,out1file,
+ & read2sigma,l_homo
common /cntrl/ iscode,indpdb,refstr,pdbref,outpdb,outmol2,
& punch_dist,print_rms,caonly,verbose,icomparfunc,pdbint,
& merge_helices,bxfile,cxfile,histfile,entfile,zscfile,rmsrgymap,
& ensembles,with_dihed_constr,check_conf,histout,constr_dist,
- & constr_homology,out1file,homol_nset
+ & constr_homology,out1file,homol_nset,read2sigma
common /homol/ waga_homology(maxR),
- & waga_dist,waga_angle,waga_theta,waga_d,dist_cut,iset,ihset
+ & waga_dist,waga_angle,waga_theta,waga_d,dist_cut,dist2_cut,
+ & iset,ihset,l_homo(max_template,maxdim)
& WSHET(MAXRES,MAXRES), EDFABET,
& CK(4),SCK(4),S1(4),S2(4)
c & ,DFAEXP(15001),
+
+ DATA CK/1.0D0,1.58740105197D0,2.08008382305D0,2.51984209979D0/
+ DATA SCK/1.0D0,1.25992104989D0,1.44224957031D0,1.58740105197D0/
+ DATA S1/3.75D0,5.75D0,7.75D0,9.75D0/
+ DATA S2/4.25D0,6.25D0,8.25D0,10.25D0/
c parameter (max_cg_procs=maxprocs)
C Max. number of AA residues
integer maxres
- parameter (maxres=600)
+ parameter (maxres=800)
c parameter (maxres=400)
C Appr. max. number of interaction sites
integer maxres2
- block data dfadata
- include 'DIMENSIONS'
- include 'COMMON.INTERACT'
- include 'COMMON.DFA'
- DATA CK/1.0D0,1.58740105197D0,2.08008382305D0,2.51984209979D0/
- DATA SCK/1.0D0,1.25992104989D0,1.44224957031D0,1.58740105197D0/
- DATA S1/3.75D0,5.75D0,7.75D0,9.75D0/
- DATA S2/4.25D0,6.25D0,8.25D0,10.25D0/
- end
-c------------------------------------------------------------------------
-
subroutine init_dfa_vars
include 'DIMENSIONS'
c write(*,*) 'case1:',tmp_n
- t1dx=t1dx+0.0d0
- t1dy=t1dy+0.0d0
- t1dz=t1dz+0.0d0
+cc t1dx=t1dx+0.0d0
+cc t1dy=t1dy+0.0d0
+cc t1dz=t1dz+0.0d0
t2dx(j)=0.0d0
t2dy(j)=0.0d0
t2dz(j)=0.0d0
tmp_n=tmp_n+1.0d0
c write(*,*) 'case4:',tmp_n
- t1dx=t1dx+0.0d0
- t1dy=t1dy+0.0d0
- t1dz=t1dz+0.0d0
+cc t1dx=t1dx+0.0d0
+cc t1dy=t1dy+0.0d0
+cc t1dz=t1dz+0.0d0
t2dx(j)=0.0d0
t2dy(j)=0.0d0
t2dz(j)=0.0d0
double precision enebet
enebet=0.0d0
- bx=0.0d0;by=0.0d0;bz=0.0d0
- shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
+c bx=0.0d0;by=0.0d0;bz=0.0d0
+c shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
gdfab=0.0d0
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc******************************************************************************
c real*8 dfaexp(15001)
real*8 uum, uup
real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
+ real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
+ common /shetf/ shetfx,shetfy,shetfz
+
common /sheca/ bx,by,bz
common /phys1/ inb,nmax,iselect
common /kyori2/ dis
common /sinsin/ sth
real*8 r_pair_mat(maxca,maxca)
+ real*8 e_gcont,fprim_gcont,de_gcont
ci integer istrand(maxca,maxca)
ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
ci common /shetest/ istrand,istrand_p,istrand_m
do i=1,inb-7
do j=i+4,inb-3
+
+ if (dis(i,j).lt.dfa_cutoff) then
+ call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
+
ip=i+1
ipp=i+2
jp=j+1
vbet(i,j)=uup+uum
vbetp=vbetp+uup
vbetm=vbetm+uum
- vbeta=vbeta+vbet(i,j)
+ vbeta=vbeta+vbet(i,j)*e_gcont
+
+
+ if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then
+c gradient correction from gcont
+ de_gcont=vbet(i,j)*fprim_gcont/dis(i,j)
+ shetfx(i)=shetfx(i) + de_gcont*rx(i,j)
+ shetfy(i)=shetfy(i) + de_gcont*ry(i,j)
+ shetfz(i)=shetfz(i) + de_gcont*rz(i,j)
+
+ shetfx(j)=shetfx(j) - de_gcont*rx(i,j)
+ shetfy(j)=shetfy(j) - de_gcont*ry(i,j)
+ shetfz(j)=shetfz(j) - de_gcont*rz(i,j)
+
+c energy correction from gcont
+ vbet(i,j)=vbet(i,j)*e_gcont
+ vbetap(i,j)=vbetap(i,j)*e_gcont
+ vbetap1(i,j)=vbetap1(i,j)*e_gcont
+ vbetap2(i,j)=vbetap2(i,j)*e_gcont
+ vbetam(i,j)=vbetam(i,j)*e_gcont
+ vbetam1(i,j)=vbetam1(i,j)*e_gcont
+ vbetam2(i,j)=vbetam2(i,j)*e_gcont
+ endif
+
ci elseif(istrand(i,j).eq.0)then
ci vbet(i,j)=0
c write(*,*) 'vbetm:',vbetm
c write(*,*) 'vbet(i,j):',vbet(i,j)
c stop
-
+
+ else
+ vbetap(i,j)=0
+ vbetap1(i,j)=0
+ vbetap2(i,j)=0
+ vbetam(i,j)=0
+ vbetam1(i,j)=0
+ vbetam2(i,j)=0
+ vbet(i,j)=0
+ endif
enddo
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbet(maxca,maxca)
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
+ real*8 e_gcont,fprim_gcont
C--------------------------------------------------------------------------------
do i=4,inb-4
im3=i-3
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
+ real*8 e_gcont,fprim_gcont
c********************************************************************************
do i=3,inb-5
imm=i-2
im=i-1
do j=i+2,inb-3
+
+ if (dis(imm,j).lt.dfa_cutoff) then
+ call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
jp=j+1
jpp=j+2
! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
! shefz(i,5)=shefz(i,5)
! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
-
+
+ endif
ci endif
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
+ real*8 e_gcont,fprim_gcont
C********************************************************************************
do i=2,inb-6
ip=i+1
im=i-1
do j=i+3,inb-3
+
+ if (dis(im,j).lt.dfa_cutoff) then
+ call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
jp=j+1
jpp=j+2
! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
! shefz(i,6)=shefz(i,6)
! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
-
+
+ endif
ci endif
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8
+ real*8 e_gcont,fprim_gcont
C********************************************************************************
do j=7,inb-1
jm=j-1
jmm=j-2
do i=1,j-6
+
+ if (dis(i,jmm).lt.dfa_cutoff) then
+ call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
ip=i+1
ipp=i+2
! shefz(j,11)=shefz(j,11)
! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
+ endif
ci endif
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
+ real*8 e_gcont,fprim_gcont
!c*************************************************************************c
do j=6,inb-2
jp=j+1
jm=j-1
do i=1,j-5
+
+ if (dis(i,jm).lt.dfa_cutoff) then
+ call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
ip=i+1
ipp=i+2
shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
$ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
+ endif
+
ci endif
ENDDO
dij=dist(i,j)
c write (iout,*) "dij(",i,j,") =",dij
do k=1,constr_homology
+ if(.not.l_homo(k,ii)) cycle
distance(k)=odl(k,ii)-dij
c write (iout,*) "distance(",k,") =",distance(k)
c
endif
enddo
- min_odl=minval(distancek)
+c min_odl=minval(distancek)
+ do kk=1,constr_homology
+ if(l_homo(kk,ii)) then
+ min_odl=distancek(kk)
+ exit
+ endif
+ enddo
+ do kk=1,constr_homology
+ if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl)
+ & min_odl=distancek(kk)
+ enddo
c write (iout,* )"min_odl",min_odl
#ifdef DEBUG
write (iout,*) "ij dij",i,j,dij
c Nie wiem po co to liczycie jeszcze raz!
c odleg3=-waga_dist(iset)*((distance(i,j,k)**2)/
c & (2*(sigma_odl(i,j,k))**2))
+ if(.not.l_homo(k,ii)) cycle
if (waga_dist.ge.0.0d0) then
c
c For Gaussian-type Urestr
c & *waga_dist)+min_odl
c sgodl=-godl(k)*distance(k)*sigma_odl(k,ii)*waga_dist
c
+ if(.not.l_homo(k,ii)) cycle
if (waga_dist.ge.0.0d0) then
c For Gaussian-type Urestr
c
double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
& escloc,ehomology_constr,
& ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
- & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt
+ & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt,
+ & edfadis,edfator,edfanei,edfabet
integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
double precision qfree,sumprob,eini,efree,rmsdev
character*80 bxname
esccor=enetb(19,i,iparm)
edihcnstr=enetb(20,i,iparm)
ehomology_constr=enetb(22,i,iparm)
+ edfadis=enetb(23,i,iparm)
+ edfator=enetb(24,i,iparm)
+ edfanei=enetb(25,i,iparm)
+ edfabet=enetb(26,i,iparm)
if (homol_nset.gt.1)
& ehomology_constr=waga_homology(homol_nset)*ehomology_constr
#ifdef SPLITELE
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
& +wbond*estr+ehomology_constr
+ & +wdfa_dist*edfadis
+ & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
#else
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
& +wbond*estr+ehomology_constr
+ & +wdfa_dist*edfadis
+ & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
#endif
#ifdef MPI
Fdimless_(i)=
character*2 kic2
character*24 model_ki_dist, model_ki_angle
character*500 controlcard
- integer ki, i, j, k, l
+ 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
c
c FP - Nov. 2014 Temporary specifications for new vars
c
- double precision rescore_tmp,x12,y12,z12
+ 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 Reading multiple PDB ref structures and calculation of retraints
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 readi(controlcard,"HOMOL_NSET",homol_nset,1)
+ call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
+ call readi(controlcard,"HOMOL_SET",homol_nset,1)
+ read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
call readi(controlcard,"IHSET",ihset,1)
if (homol_nset.gt.1)then
call card_concat(controlcard,.true.)
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 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
do k=1,constr_homology
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)') 'Error opening PDB file.'
+ 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
+ & pdbfile(:ilen(pdbfile))
stop
34 continue
c print *,'Begin reading pdb data'
write(kic2,'(bz,i2.2)') k
tpl_k_rescore="template"//kic2//".sco"
-c tpl_k_sigma_odl="template"//kic2//".sigma_odl"
-c tpl_k_sigma_dih="template"//kic2//".sigma_dih"
-c tpl_k_sigma_theta="template"//kic2//".sigma_theta"
-c tpl_k_sigma_d="template"//kic2//".sigma_d"
unres_pdb=.false.
call readpdb
c
c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
open (ientin,file=tpl_k_rescore,status='old')
- do irec=1,maxdim ! loop for reading res sim
- if (irec.eq.1) then
- rescore(k,irec)=0.0d0
- goto 1301
- endif
- read (ientin,*,end=1401) rescore_tmp
+ if (nnt.gt.1) rescore(k,1)=0.0d0
+ do irec=nnt,maxdim ! loop for reading res sim
+ if (read2sigma) then
+ read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_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
+ 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)
- 1301 continue
+ 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)
-c open (ientin,file=tpl_k_sigma_odl,status='old')
-c do irec=1,maxdim ! loop for reading sigma_odl
-c read (ientin,*,end=1401) i, j,
-c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
-c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
-c & sigma_odl_temp(i+nnt-1,j+nnt-1,k)
-c enddo
-c 1401 continue
-c close (ientin)
+ close (ientin)
if (waga_dist.ne.0.0d0) then
ii=0
- do i = nnt,nct-2 ! right? without parallel.
- do j=i+2,nct ! right?
-c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
-c do j=i+2,nres ! ibid
-c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
-c do j=i+2,nct ! ibid
+ 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
-c
-c Attempt to replace dist(i,j) by its definition in ...
-c
- 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)
odl(k,ii)=distal
-c
-c odl(k,ii)=dist(i,j)
-c write (iout,*) "dist(",i,j,") =",dist(i,j)
-c write (iout,*) "distal = ",distal
-c write (iout,*) "odl(",k,ii,") =",odl(k,ii)
-c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
-c & "rescore(",k,j,") =",rescore(k,j)
-c
-c Calculation of sigma from res sim
-c
-c if (odl(k,ii).le.6.0d0) then
-c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
-c Other functional forms possible depending on odl(k,ii), eg.
-c
- if (odl(k,ii).le.dist_cut) then
- sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
-c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
-c write (iout,*) "c sigma_odl",k,ii,sigma_odl(k,ii)
- else
+ 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))* ! sigma ~ rescore ~ error
+ 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))* ! sigma ~ rescore ~ error
+ sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
& dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
-c write (iout,*) "d sigma_odl",k,ii,sigma_odl(k,ii),
-c & odl(k,ii),dist_cut
#endif
-c Following expr replaced by a positive exp argument
-c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
-c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
-
-c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
-c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
+ 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
-c
- sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
-c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
-c
-c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
-c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity
enddo
-c read (ientin,*) sigma_odl(k,ii) ! 1st variant
enddo
-c lim_odl=ii
-c if (constr_homology.gt.0) call homology_partition
+ lim_odl=ii
endif
c
c Theta, dihedral and SC retraints
c enddo
c1402 continue
c close (ientin)
- do i = nnt+3,nct ! right? without parallel.
-c do i=1,nres ! alternative for bounds acc to readpdb?
-c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
-c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
+ 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)
sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
& rescore(k,i-2)+rescore(k,i-3))/4.0
-c
+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 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))
c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
- if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
-c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
enddo
+ lim_dih=nct-nnt-2
endif
if (waga_theta.gt.0.0d0) then
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 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
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)
- if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
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 & sigma_d(k,i+nnt-1)
c enddo
c1404 continue
- close (ientin)
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) goto 1 ! right?
+ 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 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?
- 1 continue
enddo
+ lim_xx=nct-nnt+1
endif
- close(ientin)
enddo
- if (waga_dist.ne.0.0d0) lim_odl=ii
+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
+ 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)
+ 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)
+ enddo
+ enddo
+ ii=ii-1
+ lim_odl=lim_odl-1
+ endif
+ enddo
+ enddo
+ endif
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,
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,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
- & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
+ 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,lim_dih
- write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+ write (iout,'(i5,100(2f8.2,4x))') 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,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
+ write (iout,'(i5,100(2f8.2,4x))') 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
- write(iout,'(i5,10(4f8.2,4x))') i,
+ 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
-#ifdef AIX
- call flush_(iout)
-#else
- call flush(iout)
-#endif
c -----------------------------------------------------------------
return
end
projects / unres.git / commitdiff