X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio.F90;h=4750e579ee861416b994ecdb6493941cbf6c9971;hb=91b9f78d94b96277537615722323ebe03cc0a014;hp=2a6147e7984a62c36bd49c55a89d572e9cd7ab64;hpb=52d6371a0781ac6214af501666a2fbffbded3a40;p=unres4.git diff --git a/source/unres/io.F90 b/source/unres/io.F90 index 2a6147e..4750e57 100644 --- a/source/unres/io.F90 +++ b/source/unres/io.F90 @@ -531,14 +531,14 @@ 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" @@ -553,7 +553,7 @@ 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" @@ -561,21 +561,39 @@ 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 @@ -694,6 +712,7 @@ use MPI_data use compare, only:seq_comp,contact use control + implicit none ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI @@ -732,7 +751,7 @@ ! 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 @@ -741,7 +760,7 @@ 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 @@ -807,6 +826,8 @@ 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) @@ -852,7 +873,8 @@ 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,& @@ -976,11 +998,24 @@ 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)) @@ -1013,7 +1048,7 @@ enddo endif endif - + print *,"CZY TU DOCHODZE" if (indpdb.eq.0) then nres_molec(:)=0 allocate(sequence(maxres,5)) @@ -1047,16 +1082,20 @@ ! 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)(1: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 @@ -1099,23 +1138,37 @@ 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 - vbld(i)=7.0 - vbld_inv(i)=1.0/7.0 + 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 + 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) @@ -1186,7 +1239,9 @@ ! 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 @@ -1224,6 +1279,7 @@ 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)) @@ -1316,12 +1372,33 @@ 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 @@ -1346,7 +1423,7 @@ '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) & @@ -1362,7 +1439,7 @@ 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 @@ -1410,10 +1487,102 @@ 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 & @@ -1468,15 +1637,19 @@ 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 @@ -1611,5 +1784,1075 @@ 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