!-----------------------------------------------------------------------------
use clust_data
use io_units
+ use geometry, only:int_from_cart,sc_loc_geom
+ use io_config,only:readpdb,readpdb_template
! use names
use io_base !, only: ilen
use geometry_data, only: nres,c,nres_molec,boxxsize,boxysize,boxzsize
real(kind=8) :: temper,curr_dist,emin,qpart,boltz,&
ave_dim,amax_dim,emin1
-
+ if (allocated(tempfac)) deallocate(tempfac)
allocate(tempfac(2,nres))
do i=1,64
RETURN
END SUBROUTINE WRTCLUST
!------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
subroutine ave_coord(igr)
use control_data, only:lside
allocate(nss_all(maxstr_proc)) !(maxstr_proc)
allocate(ihpb_all(maxss,maxstr_proc),jhpb_all(maxss,maxstr_proc))!(maxss,maxstr_proc)
allocate(iscore(maxconf)) !(maxconf)
-
+ nss_all=0
+ ihpb_all=0
+ jhpb_all=0
#ifdef CHUJ
ICON=1
! call flush(iout)
if (iret.eq.0) goto 101
call xdrfint(ixdrf, nss, iret)
-! write (iout,*) "iret",iret
-! write (iout,*) "nss",nss
+ write (iout,*) "iret",iret
+! write (iout,*) "nss",nss,i,"TUTU"
call flush(iout)
if (iret.eq.0) goto 101
do k=1,nss
if (iret.eq.0) goto 101
call xdrfint(ixdrf, jhpb(k), iret)
if (iret.eq.0) goto 101
+! write(iout,*), ihpb(k),jhpb(k),"TUTU"
enddo
call xdrffloat(ixdrf,reini,iret)
if (iret.eq.0) goto 101
+! write(iout,*) "TUTU", reini
call xdrffloat(ixdrf,refree,iret)
+! write(iout,*) "TUTU", refree
if (iret.eq.0) goto 101
call xdrffloat(ixdrf,rmsdev,iret)
if (iret.eq.0) goto 101
call int_from_cart1(.false.)
do j=nnt+1,nct-1
mnum=molnum(j)
- write (iout,*) "Check atom",j
- if (mnum.eq.5) cycle
+ write (iout,*) "Check atom",j,itype(j,mnum),vbld(j)
+ if (mnum.ne.1) cycle
if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
- if (itype(j+1,molnum(j+1)).eq.ntyp1_molec(molnum(j+1))) cycle
+ if (itype(j-1,molnum(j-1)).eq.ntyp1_molec(molnum(j-1))) cycle
if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then
if (j.gt.2) then
if (itel(j).ne.0 .and. itel(j-1).ne.0) then
- write (iout,*) "Conformation",jjj,jj+1
+ write (iout,*) "Conformation",jjj,jj+1,itype(j-1,molnum(j-1)),itype(j,mnum)
write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
chalen
write (iout,*) "The Cartesian geometry is:"
do j=nnt,nct
mnum=molnum(j)
itj=itype(j,mnum)
- if (mnum.eq.5) cycle
+ if (mnum.ne.1) cycle
if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) &
.and. (vbld(nres+j)-dsc(iabs(itj))) &
.gt.2.0d0) then
!------------------------------------------------------------------------------
subroutine daread_ccoords(istart_conf,iend_conf)
+ use energy_data, only: dyn_ss
+
! implicit none
! include "DIMENSIONS"
! include "sizesclu.dat"
write (iout,*) "Reading binary file, record",iii," ii",ii
call flush(iout)
#endif
+ if (dyn_ss) then
+ read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
+ ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
+ entfac(ii),rmstb(ii)
+ else
read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
- ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
- nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
- entfac(ii),rmstb(ii)
+ ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
+ nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
+ entfac(ii),rmstb(ii)
+ endif
+
#ifdef DEBUG
write (iout,*) ii,iii,ij,entfac(ii)
write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
! implicit none
! include "DIMENSIONS"
! include "sizesclu.dat"
+ use energy_data, only: dyn_ss
+
#ifdef MPI
use MPI_data
include "mpif.h"
write (iout,*) "Writing binary file, record",iii," ii",ii
call flush(iout)
#endif
+ if (dyn_ss) then
write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
- ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
- nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
- entfac(ii),rmstb(ii)
+ ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
+ entfac(ii),rmstb(ii)
+ else
+ write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
+ ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
+ nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
+ entfac(ii),rmstb(ii)
+ endif
+
#ifdef DEBUG
write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
!
! Read molecular data
!
- use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
+ use energy_data, only: rescale_mode,distchainmax,ipot,constr_homology,&
+ with_dihed_constr,read_homol_frag,out_template_coord,&
+ out_template_restr
use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
iscode,symetr,punch_dist,print_dist,nstart,nend,&
caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,&
- r_cut_ele,nclust,tor_mode,scelemode
+ r_cut_ele,nclust,tor_mode,scelemode,r_cut_mart,r_cut_ang,&
+ rlamb_mart,constr_dist
+ use geometry_data, only:bordliptop,bordlipbot,&
+ bufliptop,buflipbot,lipthick,lipbufthick
! implicit none
! include 'DIMENSIONS'
! include 'sizesclu.dat'
call card_concat(controlcard,.true.)
call readi(controlcard,'NRES',nres_molec(1),0)
- call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
- call readi(controlcard,"NRES_CAT",nres_molec(5),0)
+ call readi(controlcard,'NRES_NUCL',nres_molec(2),0)
+ call readi(controlcard,'NRES_CAT',nres_molec(5),0)
nres=0
do i=1,5
nres=nres_molec(i)+nres
enddo
+ print *,"TU",nres_molec(:)
! call alloc_clust_arrays
allocate(rcutoff(max_cut+1)) !(max_cut+1)
call reada(controlcard,'BOXX',boxxsize,100.0d0)
call reada(controlcard,'BOXY',boxysize,100.0d0)
call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+ call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+ call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+ if (lipthick.gt.0.0d0) then
+ bordliptop=(boxzsize+lipthick)/2.0
+ bordlipbot=bordliptop-lipthick
+ if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
+ write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+ buflipbot=bordlipbot+lipbufthick
+ bufliptop=bordliptop-lipbufthick
+ if ((lipbufthick*2.0d0).gt.lipthick) &
+ write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+ endif !lipthick.gt.0
+ write(iout,*) "bordliptop=",bordliptop
+ write(iout,*) "bordlipbot=",bordlipbot
+ write(iout,*) "bufliptop=",bufliptop
+ write(iout,*) "buflipbot=",buflipbot
+
call readi(controlcard,'NCLUST',nclust,5)
! ions=index(controlcard,"IONS").gt.0
call readi(controlcard,"SCELEMODE",scelemode,0)
call readi(controlcard,'TORMODE',tor_mode,0)
!C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
write(iout,*) "torsional and valence angle mode",tor_mode
- call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+ call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
write(iout,*) "R_CUT_ELE=",r_cut_ele
+ call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
+ call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
+ call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
call readi(controlcard,'NEND',nend,0)
call reada(controlcard,'ECUT',ecut,10.0d0)
call reada(controlcard,'PROB',prob_limit,0.99d0)
write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
lprint_cart=index(controlcard,"PRINT_CART") .gt.0
lprint_int=index(controlcard,"PRINT_INT") .gt.0
+ call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+ call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
+!c if (constr_homology) tole=dmax1(tole,1.5d0)
+ write (iout,*) "with_homology_constr ",with_dihed_constr,&
+ " CONSTR_HOMOLOGY",constr_homology
+ read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
+ out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
+ out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
+ write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD
if (min_var) iopt=1
return
end subroutine read_control
!-----------------------------------------------------------------------------
+ subroutine read_constr_homology
+ use MPI_data
+ use energy_data
+ use geometry_data
+ use control_data
+ 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 molread
!
! Read molecular data.
!
- use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
+ use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc,&
+ deg2rad,phibound,crefjlee,cref,rad2deg
use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
! wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
! wturn3,wturn4,wturn6,wvdwpp,weights
use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
- indpdb
+ indpdb,constr_dist,raw_psipred, with_theta_constr
use geometry, only: chainbuild,alloc_geo_arrays
use energy, only: alloc_ener_arrays
use io_base, only: rescode
! include 'COMMON.INFO'
!#endif
character(len=4) :: sequence(nres) !(maxres)
- character(len=800) :: weightcard
+ character(len=800) :: weightcard,controlcard
! integer rescode
real(kind=8) :: x(6*nres) !(maxvar)
+ real(kind=8) :: ssscale
+ real(kind=8) :: phihel,phibet,sigmahel,sigmabet,sumv,&
+ secprob(3,maxres)
integer :: itype_pdb(nres) !(maxres)
! logical seq_comp
integer :: i,j,kkk,mnum
call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
call reada(weightcard,'TEMP0',temp0,300.0d0) !!! el
if (index(weightcard,'SOFT').gt.0) ipot=6
+ call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
+ call reada(weightcard,'WELPP',welpp,0.0d0)
+ call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
+ call reada(weightcard,'WELPSB',welpsb,0.0D0)
+ call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
+ call reada(weightcard,'WELSB',welsb,0.0D0)
+ call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
+ 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_NUC',wcorr3_nucl,0.0D0)
! 12/1/95 Added weight for the multi-body term WCORR
call reada(weightcard,'WCORRH',wcorr,1.0D0)
+ call reada(weightcard,'WSCBASE',wscbase,0.0D0)
+ call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
+ call reada(weightcard,'WSCPHO',wscpho,0.0d0)
+ call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
+ call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
+ call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0)
+ call reada(weightcard,"D0CM",d0cm,3.78d0)
+ call reada(weightcard,"AKCM",akcm,15.1d0)
+ call reada(weightcard,"AKTH",akth,11.0d0)
+ call reada(weightcard,"AKCT",akct,12.0d0)
+ call reada(weightcard,"V1SS",v1ss,-1.08d0)
+ call reada(weightcard,"V2SS",v2ss,7.61d0)
+ call reada(weightcard,"V3SS",v3ss,13.7d0)
+ call reada(weightcard,"EBR",ebr,-5.50D0)
+ call reada(weightcard,"ATRISS",atriss,0.301D0)
+ call reada(weightcard,"BTRISS",btriss,0.021D0)
+ call reada(weightcard,"CTRISS",ctriss,1.001D0)
+ call reada(weightcard,"DTRISS",dtriss,1.001D0)
+ call reada(weightcard,"SSSCALE",ssscale,1.0D0)
+ dyn_ss=(index(weightcard,'DYN_SS').gt.0)
+ call reada(weightcard,"LIPSCALE",lipscale,1.0D0)
+ call reada(weightcard,"WTL",wliptran,1.0D0)
+
+ call reada(weightcard,"HT",Ht,0.0D0)
+ if (dyn_ss) then
+ ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
+ Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
+ akcm=akcm/wsc*ssscale
+ akth=akth/wsc*ssscale
+ akct=akct/wsc*ssscale
+ v1ss=v1ss/wsc*ssscale
+ v2ss=v2ss/wsc*ssscale
+ v3ss=v3ss/wsc*ssscale
+ else
+ ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+ endif
+
if (wcorr4.gt.0.0d0) wcorr=wcorr4
weights(1)=wsc
weights(2)=wscp
weights(18)=scal14
!el weights(19)=wsccor !!!!!!!!!!!!!!!!
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)=wscpho
+ weights(48)=wpeppho
+ weights(49)=wpeppho
+ weights(50)=wcatnucl
+ weights(56)=wcat_tran
+ weights(22)=wliptran
+
write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
wturn4,wturn6,wsccor
else
read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
endif
+ print *,nres_molec(:),nres
! Convert sequence to numeric code
do i=1,nres_molec(1)
mnum=1
do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
mnum=2
molnum(i)=2
+ write (iout,*),i,sequence(i)
itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
+
enddo
do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
mnum=5
! write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
! & ' nstart_seq=',nstart_seq
! endif
-write(iout,*)"przed ini_int_tab"
+ if (with_dihed_constr) then
+
+ read (inp,*) ndih_constr
+ if (ndih_constr.gt.0) then
+ raw_psipred=.false.
+!C read (inp,*) ftors
+!C write (iout,*) 'FTORS',ftors
+!C ftors is the force constant for torsional quartic constrains
+ read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),&
+ i=1,ndih_constr)
+ write (iout,*)&
+ 'There are',ndih_constr,' constraints on phi angles.'
+ do i=1,ndih_constr
+ write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),&
+ ftors(i)
+ enddo
+ do i=1,ndih_constr
+ phi0(i)=deg2rad*phi0(i)
+ drange(i)=deg2rad*drange(i)
+ enddo
+ else if (ndih_constr.lt.0) then
+ raw_psipred=.true.
+ call card_concat(controlcard,.true.)
+ call reada(controlcard,"PHIHEL",phihel,50.0D0)
+ call reada(controlcard,"PHIBET",phibet,180.0D0)
+ call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
+ call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
+ call reada(controlcard,"WDIHC",wdihc,0.591d0)
+ write (iout,*) "Weight of the dihedral restraint term",wdihc
+ read(inp,'(9x,3f7.3)')&
+ (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
+ write (iout,*) "The secprob array"
+ do i=nnt,nct
+ write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
+ enddo
+ ndih_constr=0
+ do i=nnt+3,nct
+ if (itype(i-3,molnum(i-3)).ne.ntyp1 .and. itype(i-2,molnum(i-2)).ne.ntyp1&
+ .and. itype(i-1,molnum(i-1)).ne.ntyp1 .and. itype(i,molnum(i)).ne.ntyp1) then
+ ndih_constr=ndih_constr+1
+ idih_constr(ndih_constr)=i
+ sumv=0.0d0
+ do j=1,3
+ vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
+ sumv=sumv+vpsipred(j,ndih_constr)
+ enddo
+ do j=1,3
+ vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
+ enddo
+ phibound(1,ndih_constr)=phihel*deg2rad
+ phibound(2,ndih_constr)=phibet*deg2rad
+ sdihed(1,ndih_constr)=sigmahel*deg2rad
+ sdihed(2,ndih_constr)=sigmabet*deg2rad
+ endif
+ enddo
+ write (iout,*)&
+ 'There are',ndih_constr,&
+ ' bimodal restraints on gamma angles.'
+ do i=1,ndih_constr
+ write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,&
+ restyp(itype(idih_constr(i)-2,molnum(idih_constr(i)-2)),&
+ molnum(idih_constr(i)-2)),idih_constr(i)-2,&
+ restyp(itype(idih_constr(i)-1,molnum(idih_constr(i)-1)),&
+ molnum(idih_constr(i)-2)),idih_constr(i)-1,&
+ phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,&
+ phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,&
+ (vpsipred(j,i),j=1,3)
+ enddo
+
+ endif ! endif ndif_constr.gt.0
+ endif ! with_dihed_constr
+ if (with_theta_constr) then
+!C with_theta_constr is keyword allowing for occurance of theta constrains
+ read (inp,*) ntheta_constr
+!C ntheta_constr is the number of theta constrains
+ if (ntheta_constr.gt.0) then
+!C read (inp,*) ftors
+ read (inp,*) (itheta_constr(i),theta_constr0(i),&
+ theta_drange(i),for_thet_constr(i),&
+ i=1,ntheta_constr)
+!C the above code reads from 1 to ntheta_constr
+!C itheta_constr(i) residue i for which is theta_constr
+!C theta_constr0 the global minimum value
+!C theta_drange is range for which there is no energy penalty
+!C for_thet_constr is the force constant for quartic energy penalty
+!C E=k*x**4
+ write (iout,*)&
+ 'There are',ntheta_constr,' constraints on phi angles.'
+ do i=1,ntheta_constr
+ write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),&
+ theta_drange(i),&
+ for_thet_constr(i)
+ enddo
+!C endif
+ do i=1,ntheta_constr
+ theta_constr0(i)=deg2rad*theta_constr0(i)
+ theta_drange(i)=deg2rad*theta_drange(i)
+ enddo
+!C if(me.eq.king.or..not.out1file)
+!C & write (iout,*) 'FTORS',ftors
+!C do i=1,ntheta_constr
+!C ii = itheta_constr(i)
+!C thetabound(1,ii) = phi0(i)-drange(i)
+!C thetabound(2,ii) = phi0(i)+drange(i)
+!C enddo
+ endif ! ntheta_constr.gt.0
+ endif! with_theta_constr
+ if (constr_homology.gt.0) then
+!c write (iout,*) "About to call read_constr_homology"
+!c call flush(iout)
+ call read_constr_homology
+!c write (iout,*) "Exit read_constr_homology"
+!c call flush(iout)
+ if (indpdb.gt.0 .or. pdbref) then
+ kkk=1
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=crefjlee(j,i)
+ cref(j,i,kkk)=crefjlee(j,i)
+ enddo
+ enddo
+ endif
+#ifdef DEBUG
+ write (iout,*) "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),j=1,3),&
+ (cref(j,i+nres),j=1,3)
+ enddo
+#endif
+#ifdef DEBUG
+ call int_from_cart1(.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(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
+#endif
+ else
+ homol_nset=0
+ endif
+
+
+ write(iout,*)"przed ini_int_tab"
call init_int_table
-write(iout,*)"po ini_int_tab"
-write(iout,*)"przed setup var"
+ write(iout,*)"po ini_int_tab"
+ write(iout,*)"przed setup var"
call setup_var
-write(iout,*)"po setup var"
+ write(iout,*)"po setup var"
+ if (ns.gt.0.and.dyn_ss) then
+ do i=nss+1,nhpb
+ ihpb(i-nss)=ihpb(i)
+ jhpb(i-nss)=jhpb(i)
+ forcon(i-nss)=forcon(i)
+ dhpb(i-nss)=dhpb(i)
+ enddo
+ nhpb=nhpb-nss
+ nss=0
+ link_start=1
+ link_end=nhpb
+! call hpb_partition
+ do i=1,ns
+ dyn_ss_mask(iss(i))=.true.
+ enddo
+ endif
+
write (iout,*) "molread: REFSTR",refstr
if (refstr) then
if (.not.pdbref) then
open (isidep1,file=sidepname,status="old")
call getenv('SCCORPAR',sccorname)
open (isccor,file=sccorname,status="old")
+ call getenv('BONDPAR_NUCL',bondname_nucl)
+ open (ibond_nucl,file=bondname_nucl,status='old')
call getenv('THETPAR_NUCL',thetname_nucl)
open (ithep_nucl,file=thetname_nucl,status='old')
call getenv('ROTPAR_NUCL',rotname_nucl)
open (itube,file=tubename,status='old')
call getenv('IONPAR',ionname)
open (iion,file=ionname,status='old')
+ call getenv('IONPAR_NUCL',ionnuclname)
+ open (iionnucl,file=ionnuclname,status='old')
+ call getenv('IONPAR_TRAN',iontranname)
+ open (iiontran,file=iontranname,status='old')
+
#ifndef OLDSCP
!
!
call getenv('SCPPAR',scpname)
open (iscpp,file=scpname,status='old')
+ call getenv('SCPPAR_NUCL',scpname_nucl)
+ open (iscpp_nucl,file=scpname_nucl,status='old')
+
#endif
return
end subroutine openunits
subroutine pdboutC(etot,rmsd,tytul)
use energy_data, only: ihpb,jhpb,itype,molnum
+ use energy, only:boxshift,to_box
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.CONTROL'
character(len=50) :: tytul
character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
'G','H','I','J'/)
- integer :: ica(nres)
+ integer :: ica(nres),k,iti1,ki
real(kind=8) :: etot,rmsd
integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
- real(kind=8) :: boxxxx(3)
+ real(kind=8) :: boxxxx(3),cbeg(3),Rdist(3)
write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
' ENERGY ',etot,' RMS ',rmsd
iatom=0
ichain=1
ires=0
- boxxxshift(1)=int(c(1,nnt)/boxxsize)
+! boxxxshift(1)=int(c(1,nnt)/boxxsize)
! if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
- boxxxshift(2)=int(c(2,nnt)/boxzsize)
+! boxxxshift(2)=int(c(2,nnt)/boxzsize)
! if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
- boxxxshift(3)=int(c(3,nnt)/boxzsize)
+! boxxxshift(3)=int(c(3,nnt)/boxzsize)
! if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
-
+! do i=1,3
+! if (boxxxshift(i).gt.2) boxxxshift(i)=2
+! if (boxxxshift(i).lt.-2) boxxxshift(i)=-2
+!
+! enddo
+ do j=1,3
+ cbeg(j)=c(j,nnt)
+ enddo
+ call to_box(cbeg(1),cbeg(2),cbeg(3))
boxxxx(1)=boxxsize
boxxxx(2)=boxysize
boxxxx(3)=boxzsize
do i=nnt,nct
mnum=molnum(i)
iti=itype(i,mnum)
+ if (i.ne.nct) then
+ iti1=itype(i+1,mnum)
+ mnum1=molnum(i+1)
+ if ((iti.eq.ntyp1_molec(mnum)).and.(iti1.eq.ntyp1_molec(mnum1))) cycle
+ endif
if (iti.eq.ntyp1_molec(mnum)) then
ichain=ichain+1
ires=0
write (ipdb,'(a)') 'TER'
else
ires=ires+1
+! if (molnum(i).ge.1) then
+ if (i.le.3) then
+ ki=2
+ else
+ if (itype(i,mnum).ne.ntyp1_molec(mnum)) then
+ ki=ki+1
+ endif
+ endif
+! endif
+! print *,"tu2",i,k
+ if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
+ if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
+ do j=1,3
+ cbeg(j)=c(j,ki)
+ enddo
iatom=iatom+1
ica(i)=iatom
- if (mnum.eq.5) then
- do j=1,3
- if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then
- c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j)
- else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then
- c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j)
- else
- c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j)
- endif
- enddo
- endif
+ call to_box(c(1,i),c(2,i),c(3,i))
+
+ DO k=1,3
+ Rdist(k) = boxshift(c(k,i) - cbeg(k),boxxxx(k))
+ c(k,i)=cbeg(k)+Rdist(k)
+ ENDDO
+ if ((itype(i,molnum(i)).ne.10).and.(molnum(i).ne.5)) then
+ call to_box(c(1,i+nres),c(2,i+nres),c(3,i+nres))
+ DO k=1,3
+ Rdist(k) = boxshift(c(k,i+nres) - cbeg(k),boxxxx(k))
+ c(k,i+nres)=cbeg(k)+Rdist(k)
+ ENDDO
+ endif
write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
if ((iti.ne.10).and.(mnum.ne.5)) then
return
end subroutine cartout
!------------------------------------------------------------------------------
-! subroutine alloc_clust_arrays(n_conf)
-
-! integer :: n_conf
-!COMMON.CLUSTER
-! common /clu/
-! allocate(diss(maxdist)) !(maxdist)
-!el allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
-! allocatable :: enetb !(max_ene,maxstr_proc)
-!el allocate(entfac(maxconf)) !(maxconf)
-! allocatable :: totfree_gr !(maxgr)
-!el allocate(rcutoff(max_cut+1)) !(max_cut+1)
-! common /clu1/
-! allocatable :: licz,iass !(maxgr)
-! allocatable :: nconf !(maxgr,maxingr)
-! allocatable :: iass_tot !(maxgr,max_cut)
-! allocatable :: list_conf !(maxconf)
-! common /alles/
-!el allocatable :: allcart !(3,maxres2,maxstr_proc)
-!el allocate(rmstb(maxconf)) !(maxconf)
-!el allocate(mult(nres)) !(maxres)
-!el allocatable :: nss_all !(maxstr_proc)
-!el allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
-! allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
-!COMMON.TEMPFAC
-! common /factemp/
-! allocatable :: tempfac !(2,maxres)
-!COMMON.FREE
-! common /free/
-!el allocate(beta_h(maxT)) !(maxT)
-
-! end subroutine alloc_clust_arrays
-!-----------------------------------------------------------------------------
+ subroutine read_klapaucjusz
+ use energy_data
+ use control_data, only:unres_pdb
+ use geometry_data, only:theta,thetaref,phi,phiref,&
+ xxref,yyref,zzref
+ 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
+
!-----------------------------------------------------------------------------
end module io_clust
projects / unres4.git / blobdiff