fail=.true.
ii=0
do while (fail .and. ii .le. maxsi)
- call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail)
+ call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail,molnum(i))
ii = ii+1
enddo
endif
if (refstr) then
call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')&
- itime,totT,EK,potE,totE,&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
potEcomp(23),me
format1="a133"
else
!C print *,'A CHUJ',potEcomp(23)
write (line1,'(i10,f15.2,7f12.3,i5,$)') &
- itime,totT,EK,potE,totE,&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
kinetic_T,t_bath,gyrate(),&
potEcomp(23),me
format1="a114"
call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, &
f9.3,i5,$)') &
- itime,totT,EK,potE,totE,&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
distance,potEcomp(23),me
format1="a133"
else
!C print *,'A CHUJ',potEcomp(23)
write (line1,'(i10,f15.2,8f12.3,i5,$)')&
- itime,totT,EK,potE,totE, &
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51), &
kinetic_T,t_bath,gyrate(),&
distance,potEcomp(23),me
format1="a114"
endif
+ else if (velnanoconst.ne.0 ) then
+ if (refstr) then
+ call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
+ write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, &
+ f9.3,i5,$)') &
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
+ rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
+ vecsim,vectrue,me
+ format1="a133"
+!C print *,"CHUJOWO"
+ else
+!C print *,'A CHUJ',potEcomp(23)
+ write (line1,'(i10,f15.2,8f12.3,i5,$)')&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51), &
+ kinetic_T,t_bath,gyrate(),&
+ vecsim,vectrue,me
+ format1="a114"
+ endif
else
if (refstr) then
call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') &
- itime,totT,EK,potE,totE,&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
format1="a133"
else
write (line1,'(i10,f15.2,7f12.3,i5,$)') &
- itime,totT,EK,potE,totE,&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
amax,kinetic_T,t_bath,gyrate(),me
format1="a114"
endif
use MPI_data
use compare, only:seq_comp,contact
use control
+ implicit none
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
#ifdef MPI
! integer ilen
!el external ilen
!el local varables
- integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2
+ integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2,mnum
real(kind=8),dimension(3,maxres2+2) :: c_alloc
real(kind=8),dimension(3,0:maxres2) :: dc_alloc
integer :: iti,nsi,maxsi,itrial,itmp
real(kind=8) :: wlong,scalscp,co,ssscale,phihel,phibet,sigmahel,&
- sigmabet,sumv
+ sigmabet,sumv,nres_temp
allocate(weights(n_ene))
!-----------------------------
allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
+! print *,"KUR..",wtor_nucl
call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
- call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,0.0D0)
+ call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.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,'WBOND_NUCL',wbond_nucl,0.0D0)
call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
+ call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
+ call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0)
call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
call reada(weightcard,'WSCPHO',wscpho,0.0d0)
call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
weights(17)=wbond
weights(18)=scal14
weights(21)=wsccor
+ weights(26)=wvdwpp_nucl
+ weights(27)=welpp
+ weights(28)=wvdwpsb
+ weights(29)=welpsb
+ weights(30)=wvdwsb
+ weights(31)=welsb
+ weights(32)=wbond_nucl
+ weights(33)=wang_nucl
+ weights(34)=wsbloc
+ weights(35)=wtor_nucl
+ weights(36)=wtor_d_nucl
+ weights(37)=wcorr_nucl
+ weights(38)=wcorr3_nucl
+ weights(41)=wcatcat
+ weights(42)=wcatprot
+ weights(46)=wscbase
+ weights(47)=wpepbase
+ weights(48)=wscpho
+ weights(49)=wpeppho
+ weights(50)=wcatnucl
+ weights(56)=wcat_tran
if(me.eq.king.or..not.out1file) &
write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
34 continue
! print *,'Begin reading pdb data'
call readpdb
+ if (.not.allocated(crefjlee)) allocate (crefjlee(3,2*nres+2))
+ 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
+#endif
+
! call int_from_cart1(.true.)
! print *,'Finished reading pdb data'
if(me.eq.king.or..not.out1file) &
- write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
+ write (iout,'(a,i7,a,i7)')'nsup=',nsup,&
' nstart_sup=',nstart_sup !,"ergwergewrgae"
!el if(.not.allocated(itype_pdb))
allocate(itype_pdb(nres))
nsi=0
fail=.true.
do while (fail.and.nsi.le.maxsi)
- call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
+ call gen_side(iti,theta(i+1),alph(i),omeg(i),fail,molnum(i))
nsi=nsi+1
enddo
if(fail) write(iout,*)'Adding sidechain failed for res ',&
enddo
endif
endif
-
+ print *,"CZY TU DOCHODZE"
if (indpdb.eq.0) then
nres_molec(:)=0
allocate(sequence(maxres,5))
! Read sequence if not taken from the pdb file.
molec=2
read (inp,*) nres_molec(molec)
-! print *,'nres=',nres
+ print *,'nres=',nres_molec(molec)
! allocate(sequence(maxres,5))
! if (iscode.gt.0) then
read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
+ print *,"KUR**"
+ print *,(sequence(i,molec),i=1,nres_molec(molec))
! Convert sequence to numeric code
do i=1,nres_molec(molec)
itmp=itmp+1
istype(itmp)=sugarcode(sequence(i,molec)(1:1),i)
- itype(itmp,molec)=rescode(i,sequence(i,molec)(2:2),iscode,molec)
+ sequence(i,molec)=sequence(i,molec)(1:2)
+ itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
+ write(iout,*) "NUCLE=", itype(itmp,molec)
enddo
endif
endif
if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
do i=2,nres
+ mnum=molnum(i)
+ if (molnum(i).eq.1) then
vbld(i)=vbl
vbld_inv(i)=vblinv
+
+ else
+ write(iout,*) "typy",itype(i,mnum),ntyp1_molec(mnum),i
+ vbld(i)=6.0
+ vbld_inv(i)=1.0/6.0
+ if ((itype(i,mnum).eq.ntyp1_molec(mnum)).or.&
+ (itype(i-1,mnum).eq.ntyp1_molec(mnum))) then
+ vbld(i)=1.9
+ vbld_inv(i)=1.0/1.9
+ endif
+ endif
enddo
do i=2,nres-1
- print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i
+ mnum=molnum(i)
+ if (molnum(i).eq.1) then
+! print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i
vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
+ else
+ write(iout,*) "typy2",itype(i,mnum),ntyp1_molec(mnum),i
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
+ vbld_inv(i+nres)=1.0
+ vbld(i+nres)=0.0
+ else
+ vbld(i+nres)=vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
+ vbld_inv(i+nres)=1.0/vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
+ endif
+ endif
! write (iout,*) "i",i," itype",itype(i,1),
! & " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
enddo
endif
call read_bridge
!--------------------------------
- print *,"tu dochodze"
+! print *,"tu dochodze"
! znamy nres oraz nss można zaalokowac potrzebne tablice
call alloc_geo_arrays
+ print *,"after alloc_geo_arrays"
call alloc_ener_arrays
+ print *,"after alloc_ener_arrays"
!--------------------------------
! 8/13/98 Set limits to generating the dihedral angles
do i=1,nres
enddo
else if (ndih_constr.lt.0) then
raw_psipred=.true.
+ allocate(idih_constr(nres))
allocate(secprob(3,nres))
allocate(vpsipred(3,nres))
allocate(sdihed(2,nres))
endif
#endif
nct=nres
+ allocate(ireschain(nres))
+ ireschain=0
+ write(iout,*),"before seq2chains",ireschain
+ call seq2chains
+ write(iout,*) "after seq2chains",nchain
+ allocate ( chain_border1(2,nchain))
+ chain_border1(1,1)=1
+ chain_border1(2,1)=chain_border(2,1)+1
+ do i=2,nchain-1
+ chain_border1(1,i)=chain_border(1,i)-1
+ chain_border1(2,i)=chain_border(2,i)+1
+ enddo
+ if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
+ chain_border1(2,nchain)=nres
+ write(iout,*) "nres",nres," nchain",nchain
+ do i=1,nchain
+ write(iout,*)"chain",i,chain_length(i),chain_border(1,i),&
+ chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
+ enddo
+ allocate(tabpermchain(50,5040))
+ call chain_symmetry(npermchain,tabpermchain)
print *,'NNT=',NNT,' NCT=',NCT
if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
if (pdbref) then
if(me.eq.king.or..not.out1file) &
- write (iout,'(a,i3)') 'nsup=',nsup
+ write (iout,'(a,i7)') 'nsup=',nsup
nstart_seq=nnt
if (nsup.le.(nct-nnt+1)) then
do i=0,nct-nnt+1-nsup
'Error - sequences to be superposed do not match.'
endif
111 continue
- if (nsup.eq.0) nsup=nct-nnt
+ if (nsup.eq.0) nsup=nct-nnt+1
if (nstart_sup.eq.0) nstart_sup=nnt
if (nstart_seq.eq.0) nstart_seq=nnt
if(me.eq.king.or..not.out1file) &
endif
if(me.eq.king.or..not.out1file)then
write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
- write (iout,*) 'IZ_SC=',iz_sc
+ write (iout,*) 'IZ_SC=',iz_sc, 'NCT=',nct
endif
!----------------------
call init_int_table
restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
enddo
endif
+ if (constr_homology.gt.0) then
+! write (iout,*) "Calling read_constr_homology"
+! call flush(iout)
+ call read_constr_homology
+ if (indpdb.gt.0 .or. pdbref) then
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=crefjlee(j,i)
+ cref(j,i,1)=crefjlee(j,i)
+ enddo
+ enddo
+ endif
+#define DEBUG
+#ifdef DEBUG
+ write (iout,*) "sc_loc_geom: Array C"
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),&
+ (c(j,i+nres),j=1,3)
+ enddo
+ write (iout,*) "Array Cref"
+ do i=1,nres
+ write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),&
+ (cref(j,i+nres,1),j=1,3)
+ enddo
+#endif
+ call int_from_cart1(.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+ enddo
+ else
+ homol_nset=0
+ if (start_from_model) then
+ nmodel_start=0
+ do
+ read(inp,'(a)',end=332,err=332) pdbfile
+ if (me.eq.king .or. .not. out1file)&
+ write (iout,'(a,5x,a)') 'Opening PDB file',&
+ pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=336)
+ goto 335
+ 336 write (iout,'(a,5x,a)') 'Error opening PDB file',&
+ pdbfile(:ilen(pdbfile))
+ call flush(iout)
+ stop
+ 335 continue
+ unres_pdb=.false.
+ nres_temp=nres
+! call readpdb
+ call readpdb_template(nmodel_start+1)
+ close(ipdbin)
+ if (nres.ge.nres_temp) then
+ nmodel_start=nmodel_start+1
+ pdbfiles_chomo(nmodel_start)=pdbfile
+ do i=1,2*nres
+ do j=1,3
+ chomo(j,i,nmodel_start)=c(j,i)
+ enddo
+ enddo
+ else
+ if (me.eq.king .or. .not. out1file) &
+ write (iout,'(a,2i7,1x,a)') &
+ "Different number of residues",nres_temp,nres, &
+ " model skipped."
+ endif
+ nres=nres_temp
+ enddo
+ 332 continue
+ if (nmodel_start.eq.0) then
+ if (me.eq.king .or. .not. out1file) &
+ write (iout,'(a)') &
+ "No valid starting model found START_FROM_MODELS is OFF"
+ start_from_model=.false.
+ endif
+ write (iout,*) "nmodel_start",nmodel_start
+ endif
+ endif
+
endif
if (constr_dist.gt.0) call read_dist_constr
write (iout,*) "After read_dist_constr nhpb",nhpb
if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
+ if (velnanoconst.ne.0) call read_afmnano
call hpb_partition
if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
write (iout,'(a)') 'Extended chain initial geometry.'
do i=3,nres
theta(i)=90d0*deg2rad
+ if (molnum(i).eq.2) theta(i)=160d0*deg2rad
enddo
do i=4,nres
phi(i)=180d0*deg2rad
+ if (molnum(i).eq.2) phi(i)=30d0*deg2rad
enddo
do i=2,nres-1
alph(i)=110d0*deg2rad
+ if (molnum(i).eq.2) alph(i)=30d0*deg2rad
enddo
do i=2,nres-1
omeg(i)=-120d0*deg2rad
+ if (molnum(i).eq.2) omeg(i)=60d0*deg2rad
if (itype(i,1).le.0) omeg(i)=-omeg(i)
enddo
call chainbuild
return
end subroutine molread
!-----------------------------------------------------------------------------
+ subroutine read_constr_homology
+ use control, only:init_int_table,homology_partition
+ use MD_data, only:iset
+! implicit none
+! include 'DIMENSIONS'
+!#ifdef MPI
+! include 'mpif.h'
+!#endif
+! include 'COMMON.SETUP'
+! include 'COMMON.CONTROL'
+! include 'COMMON.HOMOLOGY'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.MD'
+! include 'COMMON.QRESTR'
+! include 'COMMON.GEO'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.VAR'
+!
+
+! double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
+! & dist_cut
+! common /przechowalnia/ odl_temp(maxres,maxres,max_template),
+! & sigma_odl_temp(maxres,maxres,max_template)
+ character*2 kic2
+ character*24 model_ki_dist, model_ki_angle
+ character*500 controlcard
+ integer :: ki,i,ii,j,k,l
+ integer, dimension (:), allocatable :: ii_in_use
+ integer :: i_tmp,idomain_tmp,&
+ irec,ik,iistart,nres_temp
+! integer :: iset
+! external :: ilen
+ logical :: liiflag,lfirst
+ integer :: i01,i10
+!
+! FP - Nov. 2014 Temporary specifications for new vars
+!
+ real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
+ rescore3_tmp, dist_cut
+ real(kind=8), dimension (:,:),allocatable :: rescore
+ real(kind=8), dimension (:,:),allocatable :: rescore2
+ real(kind=8), dimension (:,:),allocatable :: rescore3
+ real(kind=8) :: distal
+ character*24 tpl_k_rescore
+ character*256 pdbfile
+
+! -----------------------------------------------------------------
+! Reading multiple PDB ref structures and calculation of retraints
+! not using pre-computed ones stored in files model_ki_{dist,angle}
+! FP (Nov., 2014)
+! -----------------------------------------------------------------
+!
+!
+! Alternative: reading from input
+ call card_concat(controlcard,.true.)
+ 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)
+ 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
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
+ write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
+ start_from_model=.false.
+ iranconf=(indpdb.le.0)
+ endif
+ if(start_from_model .and. (me.eq.king .or. .not. out1file))&
+ write(iout,*) 'START_FROM_MODELS is ON'
+! if(start_from_model .and. rest) then
+! if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+! write(iout,*) 'START_FROM_MODELS is OFF'
+! write(iout,*) 'remove restart keyword from input'
+! endif
+! endif
+ if (start_from_model) nmodel_start=constr_homology
+ if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
+ if (homol_nset.gt.1)then
+ call card_concat(controlcard,.true.)
+ read(controlcard,*) (waga_homology(i),i=1,homol_nset)
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+! write(iout,*) "iset homology_weight "
+ do i=1,homol_nset
+ write(iout,*) i,waga_homology(i)
+ enddo
+ endif
+ iset=mod(kolor,homol_nset)+1
+ else
+ iset=1
+ waga_homology(1)=1.0
+ endif
+
+!d write (iout,*) "nnt",nnt," nct",nct
+!d call flush(iout)
+
+ if (read_homol_frag) then
+ call read_klapaucjusz
+ else
+
+ lim_odl=0
+ lim_dih=0
+!
+! write(iout,*) 'nnt=',nnt,'nct=',nct
+!
+! do i = nnt,nct
+! do k=1,constr_homology
+! idomain(k,i)=0
+! enddo
+! enddo
+ idomain=0
+
+! ii=0
+! do i = nnt,nct-2
+! do j=i+2,nct
+! ii=ii+1
+! ii_in_use(ii)=0
+! enddo
+! enddo
+ ii_in_use=0
+ if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
+ if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
+
+ do k=1,constr_homology
+
+ read(inp,'(a)') pdbfile
+ pdbfiles_chomo(k)=pdbfile
+ if(me.eq.king .or. .not. out1file) &
+ 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
+! print *,'Begin reading pdb data'
+!
+! Files containing res sim or local scores (former containing sigmas)
+!
+
+ write(kic2,'(bz,i2.2)') k
+
+ tpl_k_rescore="template"//kic2//".sco"
+ write(iout,*) "tpl_k_rescore=",tpl_k_rescore
+ unres_pdb=.false.
+ nres_temp=nres
+ write(iout,*) "read2sigma",read2sigma
+
+ if (read2sigma) then
+ call readpdb_template(k)
+ else
+ call readpdb
+ endif
+ write(iout,*) "after readpdb"
+ if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
+ nres_chomo(k)=nres
+ nres=nres_temp
+ if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
+ if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
+ if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
+ if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
+ if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
+ if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
+ if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
+ if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
+ if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
+ if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
+ if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
+ if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
+ if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
+ if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
+! if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
+ if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
+ if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
+ if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
+ if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
+! if(.not.allocated(distance)) allocate(distance(constr_homology))
+! if(.not.allocated(distancek)) allocate(distancek(constr_homology))
+
+
+!
+! Distance restraints
+!
+! ... --> odl(k,ii)
+! Copy the coordinates from reference coordinates (?)
+ do i=1,2*nres_chomo(k)
+ do j=1,3
+ c(j,i)=cref(j,i,1)
+! write (iout,*) "c(",j,i,") =",c(j,i)
+ enddo
+ enddo
+!
+! From read_dist_constr (commented out 25/11/2014 <-> res sim)
+!
+! 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
+ if (.not. out1file .or. me.eq.king)&
+ 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
+
+! 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
+! 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)
+! 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.
+
+! write (iout,*) "k",k
+! write (iout,*) "i",i," j",j," constr_homology",
+! & 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
+ enddo
+ enddo
+ lim_odl=ii
+ endif
+! write (iout,*) "Distance restraints set"
+! call flush(iout)
+!
+! Theta, dihedral and SC retraints
+!
+ if (waga_angle.gt.0.0d0) then
+! open (ientin,file=tpl_k_sigma_dih,status='old')
+! do irec=1,maxres-3 ! loop for reading sigma_dih
+! read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
+! if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
+! sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
+! & sigma_dih(k,i+nnt-1)
+! enddo
+!1402 continue
+! 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?
+! read (ientin,*) sigma_dih(k,i) ! original variant
+! write (iout,*) "dih(",k,i,") =",dih(k,i)
+! write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+! & "rescore(",k,i-1,") =",rescore(k,i-1),
+! & "rescore(",k,i-2,") =",rescore(k,i-2),
+! & "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
+! if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
+! write (iout,*) "Raw sigmas for dihedral angle restraints"
+! write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
+! sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
+! rescore(k,i-2)*rescore(k,i-3) ! right expression ?
+! 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))
+! sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
+ enddo
+ lim_dih=nct-nnt-2
+ endif
+! write (iout,*) "Dihedral angle restraints set"
+! call flush(iout)
+
+ if (waga_theta.gt.0.0d0) then
+! open (ientin,file=tpl_k_sigma_theta,status='old')
+! do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
+! read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
+! sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
+! & sigma_theta(k,i+nnt-1)
+! enddo
+!1403 continue
+! close (ientin)
+
+ do i = nnt+2,nct ! right? without parallel.
+! do i = i=1,nres ! alternative for bounds acc to readpdb?
+! 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)
+! write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
+! write(iout,*) "rescore(",k,i,") =",rescore(k,i),
+! & "rescore(",k,i-1,") =",rescore(k,i-1),
+! & "rescore(",k,i-2,") =",rescore(k,i-2)
+! 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
+! 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))
+
+! sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
+! rescore(k,i-2) ! right expression ?
+! sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
+ enddo
+ endif
+! write (iout,*) "Angle restraints set"
+! call flush(iout)
+
+ if (waga_d.gt.0.0d0) then
+! open (ientin,file=tpl_k_sigma_d,status='old')
+! do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
+! read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
+! sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
+! & sigma_d(k,i+nnt-1)
+! enddo
+!1404 continue
+
+ do i = nnt,nct ! right? without parallel.
+! do i=2,nres-1 ! alternative for bounds acc to readpdb?
+! do i=loc_start,loc_end ! with FG parallel.
+ if (itype(i,1).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)
+! write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
+! write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
+! write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
+! 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))
+
+! sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
+! sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
+! read (ientin,*) sigma_d(k,i) ! 1st variant
+ enddo
+ endif
+ enddo
+! write (iout,*) "SC restraints set"
+! call flush(iout)
+!
+! remove distance restraints not used in any model from the list
+! shift data in all arrays
+!
+! write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
+ if (waga_dist.ne.0.0d0) then
+ ii=0
+ liiflag=.true.
+ lfirst=.true.
+ do i=nnt,nct-2
+ do j=i+2,nct
+ ii=ii+1
+! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+! & .and. distal.le.dist2_cut ) then
+! write (iout,*) "i",i," j",j," ii",ii
+! call flush(iout)
+ if (ii_in_use(ii).eq.0.and.liiflag.or. &
+ ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
+ liiflag=.false.
+ i10=ii
+ if (lfirst) then
+ lfirst=.false.
+ iistart=ii
+ else
+ if(i10.eq.lim_odl) i10=i10+1
+ do ki=0,i10-i01-1
+ ires_homo(iistart+ki)=ires_homo(ki+i01)
+ jres_homo(iistart+ki)=jres_homo(ki+i01)
+ ii_in_use(iistart+ki)=ii_in_use(ki+i01)
+ do k=1,constr_homology
+ odl(k,iistart+ki)=odl(k,ki+i01)
+ sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
+ l_homo(k,iistart+ki)=l_homo(k,ki+i01)
+ enddo
+ enddo
+ iistart=iistart+i10-i01
+ endif
+ endif
+ if (ii_in_use(ii).ne.0.and..not.liiflag) then
+ i01=ii
+ liiflag=.true.
+ endif
+ enddo
+ enddo
+ lim_odl=iistart-1
+ endif
+! write (iout,*) "Removing distances completed"
+! call flush(iout)
+ endif ! .not. klapaucjusz
+
+ if (constr_homology.gt.0) call homology_partition
+ write (iout,*) "After homology_partition"
+ call flush(iout)
+ if (constr_homology.gt.0) call init_int_table
+ write (iout,*) "After init_int_table"
+ call flush(iout)
+! endif ! .not. klapaucjusz
+! endif
+! if (constr_homology.gt.0) call homology_partition
+! write (iout,*) "After homology_partition"
+! call flush(iout)
+! if (constr_homology.gt.0) call init_int_table
+! write (iout,*) "After init_int_table"
+! call flush(iout)
+! write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
+! write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
+!
+! Print restraints
+!
+ if (.not.out_template_restr) return
+!d 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,'(3i7,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,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
+ (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,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
+ (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,'(i7,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
+ return
+ end subroutine read_constr_homology
+!-----------------------------------------------------------------------------
+ subroutine read_klapaucjusz
+ use energy_data
+ implicit none
+! include 'DIMENSIONS'
+!#ifdef MPI
+! include 'mpif.h'
+!#endif
+! include 'COMMON.SETUP'
+! include 'COMMON.CONTROL'
+! include 'COMMON.HOMOLOGY'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.MD'
+! include 'COMMON.GEO'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+ character(len=256) fragfile
+ integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,&
+ ii_in_use
+ integer, dimension(:,:), allocatable :: iresclust,inclust
+ integer :: nclust
+
+ character(len=2) :: kic2
+ character(len=24) :: model_ki_dist, model_ki_angle
+ character(len=500) :: controlcard
+ integer :: ki, i, j, jj,k, l, i_tmp,&
+ idomain_tmp,&
+ ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,&
+ i01,i10,nnt_chain,nct_chain
+ real(kind=8) :: distal
+ logical :: lprn = .true.
+ integer :: nres_temp
+! integer :: ilen
+! external :: ilen
+ logical :: liiflag,lfirst
+
+ real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut
+ real(kind=8), dimension (:,:), allocatable :: rescore
+ real(kind=8), dimension (:,:), allocatable :: rescore2
+ character(len=24) :: tpl_k_rescore
+ character(len=256) :: pdbfile
+
+!
+! For new homol impl
+!
+! include 'COMMON.VAR'
+!
+! write (iout,*) "READ_KLAPAUCJUSZ"
+! print *,"READ_KLAPAUCJUSZ"
+! call flush(iout)
+ call getenv("FRAGFILE",fragfile)
+ write (iout,*) "Opening", fragfile
+ call flush(iout)
+ open(ientin,file=fragfile,status="old",err=10)
+! write (iout,*) " opened"
+! call flush(iout)
+
+ sigma_theta=0.0
+ sigma_d=0.0
+ sigma_dih=0.0
+ l_homo = .false.
+
+ nres_temp=nres
+ itype_temp(:)=itype(:,1)
+ ii=0
+ lim_odl=0
+
+! write (iout,*) "Entering loop"
+! call flush(iout)
+
+ DO IGR = 1,NCHAIN_GROUP
+
+! write (iout,*) "igr",igr
+ call flush(iout)
+ read(ientin,*) constr_homology,nclust
+ if (start_from_model) then
+ nmodel_start=constr_homology
+ else
+ nmodel_start=0
+ endif
+
+ ii_old=lim_odl
+
+ ichain=iequiv(1,igr)
+ nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
+ nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
+! write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
+
+! Read pdb files
+ if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
+ if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
+ do k=1,constr_homology
+ read(ientin,'(a)') pdbfile
+ write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', &
+ pdbfile(:ilen(pdbfile))
+ pdbfiles_chomo(k)=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.
+! nres_temp=nres
+ call readpdb_template(k)
+ nres_chomo(k)=nres
+! nres=nres_temp
+ do i=1,nres
+ rescore(k,i)=0.2d0
+ rescore2(k,i)=1.0d0
+ enddo
+ enddo
+! 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
+!
+! Loop over clusters
+!
+ do l=1,nclust
+ do ll = 1,ninclust(l)
+
+ k = inclust(ll,l)
+! write (iout,*) "l",l," ll",ll," k",k
+ 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
+!
+! Distance restraints
+!
+! ... --> odl(k,ii)
+! Copy the coordinates from reference coordinates (?)
+! nres_temp=nres
+ nres=nres_chomo(k)
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=chomo(j,i,k)
+! write (iout,*) "c(",j,i,") =",c(j,i)
+ enddo
+ enddo
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+! nres=nres_temp
+ if (waga_dist.ne.0.0d0) then
+ ii=ii_old
+! do i = nnt,nct-2
+ do i = nnt_chain,nct_chain-2
+! do j=i+2,nct
+ do j=i+2,nct_chain
+
+ 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)
+! 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.
+
+! write (iout,*) "k",k
+! write (iout,*) "i",i," j",j," constr_homology",
+! & constr_homology
+ ires_homo(ii)=i+chain_border1(1,igr)-1
+ jres_homo(ii)=j+chain_border1(1,igr)-1
+ 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
+ enddo
+ enddo
+ lim_odl=ii
+ endif
+!
+! Theta, dihedral and SC retraints
+!
+ if (waga_angle.gt.0.0d0) then
+ do i = nnt_chain+3,nct_chain
+ iii=i+chain_border1(1,igr)-1
+ if (idomain(k,i).eq.0) then
+! sigma_dih(k,i)=0.0
+ cycle
+ endif
+ dih(k,iii)=phiref(i)
+ sigma_dih(k,iii)= &
+ (rescore(k,i)+rescore(k,i-1)+ &
+ rescore(k,i-2)+rescore(k,i-3))/4.0
+! write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
+! & " sigma_dihed",sigma_dih(k,i)
+ if (sigma_dih(k,iii).ne.0) &
+ sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
+ enddo
+! lim_dih=nct-nnt-2
+ endif
+
+ if (waga_theta.gt.0.0d0) then
+ do i = nnt_chain+2,nct_chain
+ iii=i+chain_border1(1,igr)-1
+ if (idomain(k,i).eq.0) then
+! sigma_theta(k,i)=0.0
+ cycle
+ endif
+ thetatpl(k,iii)=thetaref(i)
+ sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ &
+ rescore(k,i-2))/3.0
+ if (sigma_theta(k,iii).ne.0) &
+ sigma_theta(k,iii)=1.0d0/ &
+ (sigma_theta(k,iii)*sigma_theta(k,iii))
+ enddo
+ endif
+
+ if (waga_d.gt.0.0d0) then
+ do i = nnt_chain,nct_chain
+ iii=i+chain_border1(1,igr)-1
+ if (itype(i,1).eq.10) cycle
+ if (idomain(k,i).eq.0 ) then
+! sigma_d(k,i)=0.0
+ cycle
+ endif
+ xxtpl(k,iii)=xxref(i)
+ yytpl(k,iii)=yyref(i)
+ zztpl(k,iii)=zzref(i)
+ sigma_d(k,iii)=rescore(k,i)
+ if (sigma_d(k,iii).ne.0) &
+ sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
+! if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
+ enddo
+ endif
+ enddo ! l
+ enddo ! ll
+!
+! remove distance restraints not used in any model from the list
+! shift data in all arrays
+!
+! write (iout,*) "ii_old",ii_old
+ if (waga_dist.ne.0.0d0) then
+#ifdef DEBUG
+ write (iout,*) "Distance restraints from templates"
+ do iii=1,lim_odl
+ write(iout,'(4i5,100(2f8.2,1x,l1,4x))') &
+ iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), &
+ (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), &
+ ki=1,constr_homology)
+ enddo
+#endif
+ ii=ii_old
+ liiflag=.true.
+ lfirst=.true.
+ do i=nnt_chain,nct_chain-2
+ do j=i+2,nct_chain
+ ii=ii+1
+! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
+! & .and. distal.le.dist2_cut ) then
+! write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
+! call flush(iout)
+ if (ii_in_use(ii).eq.0.and.liiflag.or. &
+ ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
+ liiflag=.false.
+ i10=ii
+ if (lfirst) then
+ lfirst=.false.
+ iistart=ii
+ else
+ if(i10.eq.lim_odl) i10=i10+1
+ do ki=0,i10-i01-1
+ ires_homo(iistart+ki)=ires_homo(ki+i01)
+ jres_homo(iistart+ki)=jres_homo(ki+i01)
+ ii_in_use(iistart+ki)=ii_in_use(ki+i01)
+ do k=1,constr_homology
+ odl(k,iistart+ki)=odl(k,ki+i01)
+ sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
+ l_homo(k,iistart+ki)=l_homo(k,ki+i01)
+ enddo
+ enddo
+ iistart=iistart+i10-i01
+ endif
+ endif
+ if (ii_in_use(ii).ne.0.and..not.liiflag) then
+ i01=ii
+ liiflag=.true.
+ endif
+ enddo
+ enddo
+ lim_odl=iistart-1
+ endif
+
+ lll=lim_odl-ii_old
+
+ do i=2,nequiv(igr)
+
+ ichain=iequiv(i,igr)
+
+ do j=nnt_chain,nct_chain
+ jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
+ do k=1,constr_homology
+ dih(k,jj)=dih(k,j)
+ sigma_dih(k,jj)=sigma_dih(k,j)
+ thetatpl(k,jj)=thetatpl(k,j)
+ sigma_theta(k,jj)=sigma_theta(k,j)
+ xxtpl(k,jj)=xxtpl(k,j)
+ yytpl(k,jj)=yytpl(k,j)
+ zztpl(k,jj)=zztpl(k,j)
+ sigma_d(k,jj)=sigma_d(k,j)
+ enddo
+ enddo
+
+ jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
+! write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
+ do j=ii_old+1,lim_odl
+ ires_homo(j+lll)=ires_homo(j)+jj
+ jres_homo(j+lll)=jres_homo(j)+jj
+ do k=1,constr_homology
+ odl(k,j+lll)=odl(k,j)
+ sigma_odl(k,j+lll)=sigma_odl(k,j)
+ l_homo(k,j+lll)=l_homo(k,j)
+ enddo
+ enddo
+
+ ii_old=ii_old+lll
+ lim_odl=lim_odl+lll
+
+ enddo
+
+ ENDDO ! IGR
+
+ if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
+ nres=nres_temp
+ itype(:,1)=itype_temp(:)
+
+ return
+ 10 stop "Error in fragment file"
+ end subroutine read_klapaucjusz
+
+!-----------------------------------------------------------
+ subroutine seq2chains
+!c
+!c Split the total UNRES sequence, which has dummy residues separating
+!c the chains, into separate chains. The length of chain ichain is
+!c contained in chain_length(ichain), the first and last non-dummy
+!c residues are in chain_border(1,ichain) and chain_border(2,ichain),
+!c respectively. The lengths pertain to non-dummy residues only.
+!c
+! implicit none
+! include 'DIMENSIONS'
+ use energy_data, only:molnum,nchain,chain_length,ireschain
+ implicit none
+! integer ireschain(nres)
+ integer ii,ichain,i,j,mnum
+ logical new_chain
+ print *,"in seq2"
+ ichain=1
+ new_chain=.true.
+ if (.not.allocated(chain_length)) allocate(chain_length(50))
+ if (.not.allocated(chain_border)) allocate(chain_border(2,50))
+
+ chain_length(ichain)=0
+ ii=1
+ do while (ii.lt.nres)
+ write(iout,*) "in seq2chains",ii,new_chain
+ mnum=molnum(ii)
+ if (itype(ii,mnum).eq.ntyp1_molec(mnum)) then
+ if (.not.new_chain) then
+ new_chain=.true.
+ chain_border(2,ichain)=ii-1
+ ichain=ichain+1
+ chain_border(1,ichain)=ii+1
+ chain_length(ichain)=0
+ endif
+ else
+ if (new_chain) then
+ chain_border(1,ichain)=ii
+ new_chain=.false.
+ endif
+ chain_length(ichain)=chain_length(ichain)+1
+ endif
+ ii=ii+1
+ enddo
+ if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) then
+ ii=ii-1
+ else
+ chain_length(ichain)=chain_length(ichain)+1
+ endif
+ if (chain_length(ichain).gt.0) then
+ chain_border(2,ichain)=ii
+ nchain=ichain
+ else
+ nchain=ichain-1
+ endif
+ ireschain=0
+ do i=1,nchain
+ do j=chain_border(1,i),chain_border(2,i)
+ ireschain(j)=i
+ enddo
+ enddo
+ return
+ end subroutine
+!---------------------------------------------------------------------
+ subroutine chain_symmetry(npermchain,tabpermchain)
+!c
+!c Determine chain symmetry. nperm is the number of permutations and
+!c tabperchain contains the allowed permutations of the chains.
+!c
+! implicit none
+! include "DIMENSIONS"
+! include "COMMON.IOUNITS"
+ implicit none
+ integer itemp(50),&
+ npermchain,tabpermchain(50,5040),&
+ tabperm(50,5040),mapchain(50),&
+ iequiv(50,nres),iflag(nres)
+ integer i,j,k,l,ii,nchain_group,nequiv(50),iieq,&
+ nperm,npermc,ind,mnum
+ if (nchain.eq.1) then
+ npermchain=1
+ tabpermchain(1,1)=1
+!c print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
+ return
+ endif
+!c
+!c Look for equivalent chains
+#ifdef DEBUG
+ write(iout,*) "nchain",nchain
+ do i=1,nchain
+ write(iout,*) "chain",i," from",chain_border(1,i),&
+ " to",chain_border(2,i)
+ write(iout,*)&
+ "sequence ",(itype(j,molnum(j)),j=chain_border(1,i),chain_border(2,i))
+ enddo
+#endif
+ do i=1,nchain
+ iflag(i)=0
+ enddo
+ nchain_group=0
+ do i=1,nchain
+ if (iflag(i).gt.0) cycle
+ iflag(i)=1
+ nchain_group=nchain_group+1
+ iieq=1
+ iequiv(iieq,nchain_group)=i
+ do j=i+1,nchain
+ if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
+!c k=0
+!c do while(k.lt.chain_length(i) .and.
+!c & itype(chain_border(1,i)+k).eq.itype(chain_border(1,j)+k))
+ do k=0,chain_length(i)-1
+!c k=k+1
+ mnum=molnum(k)
+ if (itype(chain_border(1,i)+k,mnum).ne.&
+ itype(chain_border(1,j)+k,mnum)) exit
+ enddo
+ if (k.lt.chain_length(i)) cycle
+ iflag(j)=1
+ iieq=iieq+1
+ iequiv(iieq,nchain_group)=j
+ enddo
+ nequiv(nchain_group)=iieq
+ enddo
+ write(iout,*) "Number of equivalent chain groups:",nchain_group
+ write(iout,*) "Equivalent chain groups"
+ do i=1,nchain_group
+ write(iout,*) "group",i," #members",nequiv(i)," chains",&
+ (iequiv(j,i),j=1,nequiv(i))
+ enddo
+ ind=0
+ do i=1,nchain_group
+ do j=1,nequiv(i)
+ ind=ind+1
+ mapchain(ind)=iequiv(j,i)
+ enddo
+ enddo
+ write (iout,*) "mapchain"
+ do i=1,nchain
+ write (iout,*) i,mapchain(i)
+ enddo
+ ii=0
+ do i=1,nchain_group
+ call permut(nequiv(i),nperm,tabperm)
+ if (ii.eq.0) then
+ ii=nequiv(i)
+ npermchain=nperm
+ do j=1,nperm
+ do k=1,ii
+ tabpermchain(k,j)=iequiv(tabperm(k,j),i)
+ enddo
+ enddo
+ else
+ npermc=npermchain
+ npermchain=npermchain*nperm
+ ind=0
+ do k=1,nperm
+ do j=1,npermc
+ ind=ind+1
+ do l=1,ii
+ tabpermchain(l,ind)=tabpermchain(l,j)
+ enddo
+ do l=1,nequiv(i)
+ tabpermchain(ii+l,ind)=iequiv(tabperm(l,k),i)
+ enddo
+ enddo
+ enddo
+ ii=ii+nequiv(i)
+ endif
+ enddo
+ do i=1,npermchain
+ do j=1,nchain
+ itemp(mapchain(j))=tabpermchain(j,i)
+ enddo
+ do j=1,nchain
+ tabpermchain(j,i)=itemp(j)
+ enddo
+ enddo
+ write(iout,*) "Number of chain permutations",npermchain
+ write(iout,*) "Permutations"
+ do i=1,npermchain
+ write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
+ enddo
+ return
+ end
+!c---------------------------------------------------------------------
+ integer function tperm(i,iperm,tabpermchain)
+! implicit none
+! include 'DIMENSIONS'
+ integer i,iperm
+ integer tabpermchain(50,5040)
+ if (i.eq.0) then
+ tperm=0
+ else
+ tperm=tabpermchain(i,iperm)
+ endif
+ return
+ end function
+
!-----------------------------------------------------------------------------
end module io
projects / unres4.git / blobdiff