do i=1,nres
itype(i)=rescode(i,sequence(i),iscode)
enddo
+ if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
+ write (iout,*)
+ & "Glycine is the first full residue, initial dummy deleted"
+ do i=1,nres
+ itype(i)=itype(i+1)
+ enddo
+ nres=nres-1
+ endif
+ if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
+ write (iout,*)
+ & "Glycine is the last full residue, terminal dummy deleted"
+ nres=nres-1
+ endif
+
C Assign initial virtual bond lengths
do i=2,nres
vbld(i)=vbl
call readi(controlcard,"NPAIR",npair_,0)
call readi(controlcard,"NDIST",ndist_,0)
call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+
call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
endif
enddo
do i=1,ndist_
+ if (constr_dist.eq.11) then
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+ & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
+ fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
+ else
read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
& ibecarb(i),forcon(nhpb+1)
+ endif
if (forcon(nhpb+1).gt.0.0d0) then
nhpb=nhpb+1
if (ibecarb(i).gt.0) then
if (.not.out1file .or. me.eq.king) then
#endif
do i=1,nhpb
+ if (constr_dist.eq.11) then
+ write (iout,'(a,3i5,2f8.2,i2,2f10.1)') "+dist.constr11 ",
+ & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i),
+ & fordepth(i)
+ else
write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
& i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
+ endif
enddo
call flush(iout)
#ifdef MPI
include 'COMMON.MD'
include 'COMMON.GEO'
include 'COMMON.INTERACT'
+ include 'COMMON.NAMES'
c
c For new homol impl
c
character*2 kic2
character*24 model_ki_dist, model_ki_angle
character*500 controlcard
- integer ki, i, j, k, l, ii_in_use(maxdim)
+ integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
logical lprn /.true./
+ integer ilen
+ external ilen
+ logical liiflag
c
c FP - Nov. 2014 Temporary specifications for new vars
c
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)
+ start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
+ if(.not.read2sigma.and.start_from_model) then
+ write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
+ start_from_model=.false.
+ endif
+ if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON'
+ if(start_from_model .and. rest) then
+ write(iout,*) 'START_FROM_MODELS is OFF'
+ write(iout,*) 'remove restart keyword from input'
+ endif
if (homol_nset.gt.1)then
call card_concat(controlcard)
read(controlcard,*) (waga_homology(i),i=1,homol_nset)
lim_odl=0
lim_dih=0
c
-c New
-c
- lim_theta=0
- lim_xx=0
+ write(iout,*) 'nnt=',nnt,'nct=',nct
c
+ do i = nnt,nct
+ do k=1,constr_homology
+ idomain(k,i)=0
+ enddo
+ enddo
+
ii=0
do i = nnt,nct-2
do j=i+2,nct
do k=1,constr_homology
+ read(inp,'(a)') 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
tpl_k_rescore="template"//kic2//".sco"
+ unres_pdb=.false.
+ if (read2sigma) then
+ call readpdb_template(k)
+ else
+ call readpdb
+ 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(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')
- do irec=1,maxdim ! loop for reading res sim
- if (irec.eq.1) then
- rescore(k,irec)=0.0d0
- goto 1301
- endif
+ 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) rescore2_tmp,rescore_tmp,
- & idomain(k,irec)
- rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0)
- rescore2(k,irec)=0.5d0*(rescore2_tmp+0.5d0)
+ 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
+ write(iout,'(a7,i5,2f10.5,i5)') "rescore",
+ & i_tmp,rescore2_tmp,rescore_tmp,
+ & idomain_tmp
else
idomain(k,irec)=1
read (ientin,*,end=1401) rescore_tmp
rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
endif
- 1301 continue
enddo
1401 continue
- close (ientin)
-
- read(inp,'(a)') 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))
- open(ipdbin,file=pdbfile,status='old',err=33)
- goto 34
- 33 write (iout,'(a)') 'Error opening PDB file.'
- stop
- 34 continue
-c print *,'Begin reading pdb data'
-
- unres_pdb=.false.
- if (read2sigma) then
- call readpdb_template(k)
- else
- call readpdb
- endif
-c
-C Copy the coordinates from reference coordinates (?)
- do i=1,2*nres
- do j=1,3
- c(j,i)=cref(j,i)
-c write (iout,*) "c(",j,i,") =",c(j,i)
- enddo
- enddo
+ 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 (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+ & .and. distal.le.dist2_cut ) then
- 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
-
- if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then
-
ii=ii+1
ii_in_use(ii)=1
l_homo(k,ii)=.true.
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 (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) ! 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
- endif
- 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)
-
-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
+ 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.
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
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) ! right expression ?
-c
+ 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
- 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
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) ! right expression ?
- sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
+ 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
- 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
-
do i = nnt,nct ! right? without parallel.
c do i=2,nres-1 ! alternative for bounds acc to readpdb?
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))
+ 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
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(2f8.2,1x,l1,4x))')
+ 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),
+ 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,10(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
- write(iout,'(i5,10(4f8.2,4x))') i,
+ 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