X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Freadrtns_CSA.F;h=afad609936acc69abfceeb6364c64c7bb4693145;hb=8bb1cdb2042196e72a125024b254c96a72d7c7ef;hp=0320484f901e1369ba0b67e78376c27e124fcc95;hpb=05c1789531ec3f44dfebf41b426cb1204329e7ff;p=unres.git diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index 0320484..afad609 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -149,6 +149,7 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0 call readi(controlcard,'SYM',symetr,1) + call readi(controlcard,'PERMUT',npermut,1) call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes @@ -167,6 +168,8 @@ c write (iout,'(a,f10.1)') 'Time limit (min):',timlim c endif call readi(controlcard,'NZ_START',nz_start,0) call readi(controlcard,'NZ_END',nz_end,0) + call readi(controlcard,'NRAN_START',nran_start,0) + write (iout,*) "nran_start",nran_start c call readi(controlcard,'IZ_SC',iz_sc,0) timlim=60.0D0*timlim safety = 60.0d0*safety @@ -187,13 +190,18 @@ c call readi(controlcard,'IZ_SC',iz_sc,0) pdbref=(index(controlcard,'PDBREF').gt.0) refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0) indpdb=index(controlcard,'PDBSTART') - extconf=(index(controlcard,'EXTCONF').gt.0) AFMlog=(index(controlcard,'AFM')) selfguide=(index(controlcard,'SELFGUIDE')) c print *,'AFMlog',AFMlog,selfguide,"KUPA" call readi(controlcard,'TUBEMOD',tubelog,0) c write (iout,*) TUBElog,"TUBEMODE" call readi(controlcard,'IPRINT',iprint,0) +c 6/22/2021 AL: alpha_GB: parameter to switch between the GB SC-SC +c interaction potential and the all-repulsive potential with singularity +c at zero site-site distance + call reada(controlcard,'ALPHA_GB',alpha_GB,1.0d-2) + alpha_GB1 = 1.0d0+1.0d0/alpha_GB + write (iout,*) "alpha_GB",alpha_GB," alpha_GB1",alpha_GB1 C SHIELD keyword sets if the shielding effect of side-chains is used C 0 denots no shielding is used all peptide are equally despite the C solvent accesible area @@ -296,6 +304,12 @@ cfmc modecalc=10 indphi=index(controlcard,'PHI') indback=index(controlcard,'BACK') iranconf=index(controlcard,'RAND_CONF') + start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0) + extconf=(index(controlcard,'EXTCONF').gt.0) + if (start_from_model) then + iranconf=0 + extconf=.false. + endif i2ndstr=index(controlcard,'USE_SEC_PRED') gradout=index(controlcard,'GRADOUT').gt.0 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0 @@ -476,6 +490,7 @@ C call readi(controlcard,"NSTEP",n_timestep,1000000) call readi(controlcard,"NTWE",ntwe,100) call readi(controlcard,"NTWX",ntwx,1000) + call readi(controlcard,"REST_FREQ",irest_freq,1000) call reada(controlcard,"DT",d_time,1.0d-1) call reada(controlcard,"DVMAX",dvmax,2.0d1) call reada(controlcard,"DAMAX",damax,1.0d1) @@ -551,6 +566,7 @@ c if performing umbrella sampling, fragments constrained are read from the frag & " A" write(iout,'(a60,i5)')"Frequency of updating interaction list", & imatupdate + write(iout,'(a60,i5)')"Restart writing frequency",irest_freq if (RESPA) then write (iout,'(2a,i4,a)') & "A-MTS algorithm used; initial time step for fast-varying", @@ -736,7 +752,7 @@ C integer ilen external ilen integer iperm,tperm - integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2 + integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp double precision sumv C C Read PDB structure if applicable @@ -831,7 +847,7 @@ c print '(20i4)',(itype(i),i=1,nres) do i=1,nres-1 write (iout,*) i,itype(i),itel(i) enddo - print *,'Call Read_Bridge.' +c print *,'Call Read_Bridge.' endif nnt=1 nct=nres @@ -844,7 +860,7 @@ cd print *,'NNT=',NNT,' NCT=',NCT chain_border1(1,i)=chain_border(1,i)-1 chain_border1(2,i)=chain_border(2,i)+1 enddo - chain_border1(1,nchain)=chain_border(1,nchain)-1 + 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 @@ -852,16 +868,19 @@ cd print *,'NNT=',NNT,' NCT=',NCT & chain_border(2,i),chain_border1(1,i),chain_border1(2,i) enddo call chain_symmetry(nchain,nres,itype,chain_border, - & chain_length,npermchain,tabpermchain) + & chain_length,npermchain,tabpermchain,nchain_group,nequiv, + & iequiv,mapchain) c do i=1,nres c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain), c & ii=1,npermchain) c enddo +#ifdef DEBUG write(iout,*) "residue permutations" do i=1,nres write(iout,*) i,(iperm(i,ii),ii=1,npermchain) enddo call flush(iout) +#endif if (itype(1).eq.ntyp1) nnt=2 if (itype(nres).eq.ntyp1) nct=nct-1 write (iout,*) "nnt",nnt," nct",nct @@ -870,9 +889,9 @@ c enddo if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and. & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then call init_dfa_vars - print*, 'init_dfa_vars finished!' +c print*, 'init_dfa_vars finished!' call read_dfa_info - print*, 'read_dfa_info finished!' +c print*, 'read_dfa_info finished!' endif #endif if (pdbref) then @@ -966,25 +985,34 @@ C & write (iout,*) 'FTORS',ftors call reada(weightcard,"WDIHC",wdihc,0.591D0) write (iout,*) "Weight of dihedral angle restraints",wdihc read(inp,'(9x,3f7.3)') +c & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct) & (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 + iconstr_dih=0 do i=nnt+3,nct if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then ndih_constr=ndih_constr+1 idih_constr(ndih_constr)=i + iconstr_dih(i)=ndih_constr 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 + if (sumv.gt.0.0d0) then + do j=1,3 + vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv + enddo + else + vpsipred(1,ndih_constr)=1.0d0 + vpsipred(2,ndih_constr)=0.0d0 + vpsipred(3,ndih_constr)=0.0d0 + endif phibound(1,ndih_constr)=phihel*deg2rad phibound(2,ndih_constr)=phibet*deg2rad sdihed(1,ndih_constr)=sigmahel*deg2rad @@ -1097,10 +1125,10 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) endif endif c print *, "A TU" - write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup +c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup call flush(iout) if (constr_dist.gt.0) call read_dist_constr - write (iout,*) "After read_dist_constr nhpb",nhpb +c write (iout,*) "After read_dist_constr nhpb",nhpb if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp call hpb_partition call NMRpeak_partition @@ -1120,8 +1148,11 @@ c print *, "A TU" endif write (iout,*) "calling read_saxs_consrtr",nsaxs if (nsaxs.gt.0) call read_saxs_constr - +c write (iout,*) "After read_saxs_constr" +c call flush(iout) if (constr_homology.gt.0) then +c write (iout,*) "Calling read_constr_homology" +c call flush(iout) call read_constr_homology if (indpdb.gt.0 .or. pdbref) then do i=1,2*nres @@ -1163,6 +1194,50 @@ c print *, "A TU" 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 +c 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 @@ -1172,18 +1247,24 @@ C endif & modecalc.ne.10) then C If input structure hasn't been supplied from the PDB file read or generate C initial geometry. - if (iranconf.eq.0 .and. .not. extconf) then + if (iranconf.eq.0 .and. .not. extconf .and. .not. + & start_from_model) then if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) & write (iout,'(a)') 'Initial geometry will be read in.' if (read_cart) then read(inp,'(8f10.5)',end=36,err=36) & ((c(l,k),l=1,3),k=1,nres), & ((c(l,k+nres),l=1,3),k=nnt,nct) - write (iout,*) "Exit READ_CART" + if (nnt.gt.1) c(:,nres+1)=c(:,1) + if (nct.lt.nres) c(:,2*nres)=c(:,nres) +c write (iout,*) "Exit READ_CART" c write (iout,'(8f10.5)') c & ((c(l,k),l=1,3),k=1,nres), c & ((c(l,k+nres),l=1,3),k=nnt,nct) call cartprint + do j=1,3 + dc(j,0)=c(j,1) + enddo do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) @@ -1243,7 +1324,7 @@ c return enddo call bond_regular call chainbuild_extconf - else + else if (.not. start_from_model) then if(me.eq.king.or..not.out1file) & write (iout,'(a)') 'Random-generated initial geometry.' call bond_regular @@ -1405,7 +1486,13 @@ C Read information about disulfide bridges. integer i,j C Read bridging residues. read (inp,*) ns,(iss(i),i=1,ns) - print *,'ns=',ns +c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine +c numbers + icys=0 + do i=1,ns + icys(iss(i))=i + enddo +c print *,'ns=',ns if(me.eq.king.or..not.out1file) & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns) C Check whether the specified bridging residues are cystines. @@ -1614,9 +1701,11 @@ C Generate CA distance constraints. include 'COMMON.CONTROL' include 'COMMON.DBASE' include 'COMMON.THREAD' + include 'COMMON.SPLITELE' include 'COMMON.TIME1' integer i,j,itype_pdb(maxres) common /pizda/ itype_pdb + double precision dd double precision dist character*2 iden cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct @@ -1627,11 +1716,14 @@ cd & ' nsup',nsup cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)), cd & ' seq_pdb', restyp(itype_pdb(i)) do j=i+2,nstart_sup+nsup-1 +c 5/24/2020 Adam: Cutoff included to reduce array size + dd = dist(i,j) + if (dd.gt.r_cut_int) cycle nhpb=nhpb+1 ihpb(nhpb)=i+nstart_seq-nstart_sup jhpb(nhpb)=j+nstart_seq-nstart_sup forcon(nhpb)=weidis - dhpb(nhpb)=dist(i,j) + dhpb(nhpb)=dd enddo enddo cd write (iout,'(a)') 'Distance constraints:' @@ -2265,6 +2357,7 @@ c1out open(iout,file=outname,status='unknown') #else if (me.eq.king .or. .not.out1file) & open(iout,file=outname,status='unknown') +#define DEBUG #ifdef DEBUG if (fg_rank.gt.0) then write (liczba,'(i3.3)') myrank/nfgtasks @@ -2273,6 +2366,7 @@ c1out open(iout,file=outname,status='unknown') & status='unknown') endif #endif +#undef DEBUG if(me.eq.king) then open(igeom,file=intname,status='unknown',access='append') open(ipdb,file=pdbname,status='unknown') @@ -2375,10 +2469,10 @@ c------------------------------------------------------------------------------ open(irest2,file=rest2name,status='unknown') read(irest2,*) totT,EK,potE,totE,t_bath totTafm=totT - do i=1,2*nres + do i=0,2*nres-1 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo - do i=1,2*nres + do i=0,2*nres-1 read(irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo if(usampl) then @@ -2468,7 +2562,7 @@ c print *, "wchodze" call readi(afmcard,"END",afmend,0) call reada(afmcard,"FORCE",forceAFMconst,0.0d0) call reada(afmcard,"VEL",velAFMconst,0.0d0) - print *,'FORCE=' ,forceAFMconst +c print *,'FORCE=' ,forceAFMconst CCCC NOW PROPERTIES FOR AFM distafminit=0.0d0 do i=1,3 @@ -2658,10 +2752,10 @@ c write (iout,*) "j",j," k",k endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", + & write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", + write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo @@ -2705,10 +2799,10 @@ c write (iout,*) "j",j," k",k endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", + & write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", + write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo @@ -2735,13 +2829,13 @@ c & fordepth_peak(nhpb_peak+1),npeak ipeak(2,npeak)=i #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ", + & write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak), & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak), & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak), & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak) #else - write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ", + write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak), & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak), & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak), @@ -2764,11 +2858,11 @@ c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) irestr_type(nhpb)=11 #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ", + & write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb) #else - write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ", + write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb) #endif @@ -2828,12 +2922,12 @@ c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1) endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", + & write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) #else - write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", + write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) @@ -2862,10 +2956,10 @@ C print *,"in else" endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", + & write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", + write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) #endif endif @@ -2908,13 +3002,13 @@ C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) dhpb(nhpb)=dist(i+nres*ii,j+nres*jj) #ifdef MPI if (.not.out1file .or. me.eq.king) then - write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", + write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) endif #else - write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", + write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) @@ -2978,11 +3072,13 @@ c & 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,ii_in_use(maxdim),i_tmp,idomain_tmp,irec, - & ik,iistart,iishift + integer ki,i,ii,j,k,l,ii_in_use(maxchain*maxdim),i_tmp, + & idomain_tmp, + & irec,ik,iistart,nres_temp integer ilen external ilen - logical liiflag + logical liiflag,lfirst + integer i01,i10 c c FP - Nov. 2014 Temporary specifications for new vars c @@ -2992,7 +3088,8 @@ c double precision, dimension (max_template,maxres) :: rescore2 double precision, dimension (max_template,maxres) :: rescore3 double precision distal - character*24 pdbfile,tpl_k_rescore + character*24 tpl_k_rescore + character*256 pdbfile c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints c not using pre-computed ones stored in files model_ki_{dist,angle} @@ -3008,6 +3105,7 @@ c Alternative: reading from input 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) + dist1cut=(index(controlcard,'DIST1CUT').gt.0) call readi(controlcard,"HOMOL_NSET",homol_nset,1) read2sigma=(index(controlcard,'READ2SIGMA').gt.0) start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0) @@ -3015,20 +3113,22 @@ c Alternative: reading from input 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 +c if(start_from_model .and. rest) then +c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then +c write(iout,*) 'START_FROM_MODELS is OFF' +c write(iout,*) 'remove restart keyword from input' +c endif +c endif + if (start_from_model) nmodel_start=constr_homology if (homol_nset.gt.1)then call card_concat(controlcard) 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 " +c write(iout,*) "iset homology_weight " do i=1,homol_nset write(iout,*) i,waga_homology(i) enddo @@ -3042,33 +3142,34 @@ c Alternative: reading from input cd write (iout,*) "nnt",nnt," nct",nct cd call flush(iout) + if (read_homol_frag) then + call read_klapaucjusz + else lim_odl=0 lim_dih=0 c c write(iout,*) 'nnt=',nnt,'nct=',nct c - do i = nnt,nct - do k=1,constr_homology - idomain(k,i)=0 - enddo - enddo - - ii=0 - do i = nnt,nct-2 - do j=i+2,nct - ii=ii+1 - ii_in_use(ii)=0 - enddo - enddo - - if (read_homol_frag) then - call read_klapaucjusz - else +c do i = nnt,nct +c do k=1,constr_homology +c idomain(k,i)=0 +c enddo +c enddo + idomain=0 +c ii=0 +c do i = nnt,nct-2 +c do j=i+2,nct +c ii=ii+1 +c ii_in_use(ii)=0 +c enddo +c enddo + ii_in_use=0 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)) @@ -3088,17 +3189,20 @@ c tpl_k_rescore="template"//kic2//".sco" unres_pdb=.false. + nres_temp=nres if (read2sigma) then call readpdb_template(k) else call readpdb endif + nres_chomo(k)=nres + nres=nres_temp c c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) - do i=1,2*nres + do i=1,2*nres_chomo(k) do j=1,3 c(j,i)=cref(j,i) c write (iout,*) "c(",j,i,") =",c(j,i) @@ -3144,7 +3248,19 @@ c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) z12=c(3,i)-c(3,j) distal=dsqrt(x12*x12+y12*y12+z12*z12) c write (iout,*) k,i,j,distal,dist2_cut - + if (dist1cut .and. k.gt.1) then + ii=ii+1 + if (l_homo(1,ii)) then + ii_in_use(ii)=1 + l_homo(k,ii)=.true. + ires_homo(ii)=i + jres_homo(ii)=j + odl(k,ii)=distal + sigma_odl(k,ii)=sigma_odl(1,ii) + else + l_homo(k,ii)=.false. + endif + else if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 & .and. distal.le.dist2_cut ) then @@ -3181,13 +3297,16 @@ c & constr_homology endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else - ii=ii+1 - l_homo(k,ii)=.false. +c ii=ii+1 +c l_homo(k,ii)=.false. endif + endif enddo enddo lim_odl=ii endif +c write (iout,*) "Distance restraints set" +c call flush(iout) c c Theta, dihedral and SC retraints c @@ -3228,6 +3347,8 @@ c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i) enddo lim_dih=nct-nnt-2 endif +c write (iout,*) "Dihedral angle restraints set" +c call flush(iout) if (waga_theta.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_theta,status='old') @@ -3263,6 +3384,8 @@ c rescore(k,i-2) ! right expression ? c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) enddo endif +c write (iout,*) "Angle restraints set" +c call flush(iout) if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') @@ -3298,51 +3421,66 @@ c read (ientin,*) sigma_d(k,i) ! 1st variant enddo endif enddo +c write (iout,*) "SC restraints set" +c call flush(iout) c c remove distance restraints not used in any model from the list c shift data in all arrays c +c 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 (ii_in_use(ii).eq.0.and.liiflag) then +c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 +c & .and. distal.le.dist2_cut ) then +c write (iout,*) "i",i," j",j," ii",ii +c 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. - iistart=ii + 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.or. - & .not.liiflag.and.ii.eq.lim_odl) then - if (ii.eq.lim_odl) then - iishift=ii-iistart+1 - else - iishift=ii-iistart - endif + if (ii_in_use(ii).ne.0.and..not.liiflag) then + i01=ii liiflag=.true. - do ki=iistart,lim_odl-iishift - ires_homo(ki)=ires_homo(ki+iishift) - jres_homo(ki)=jres_homo(ki+iishift) - ii_in_use(ki)=ii_in_use(ki+iishift) - do k=1,constr_homology - odl(k,ki)=odl(k,ki+iishift) - sigma_odl(k,ki)=sigma_odl(k,ki+iishift) - l_homo(k,ki)=l_homo(k,ki+iishift) - enddo - enddo - ii=ii-iishift - lim_odl=lim_odl-iishift endif enddo enddo + lim_odl=iistart-1 endif - +c write (iout,*) "Removing distances completed" +c call flush(iout) endif ! .not. klapaucjusz if (constr_homology.gt.0) call homology_partition +c write (iout,*) "After homology_partition" +c call flush(iout) if (constr_homology.gt.0) call init_int_table -c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end -c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end +c write (iout,*) "After init_int_table" +c call flush(iout) +c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end +c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end c c Print restraints c @@ -3351,26 +3489,26 @@ cd 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,'(3i5,100(2f8.2,1x,l1,4x))') + 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,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), + write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)), & (rad2deg*dih(ki,i), & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology) enddo write (iout,*) "Virtual-bond angle restraints from templates" do i=nnt+2,nct - write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), + write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)), & (rad2deg*thetatpl(ki,i), & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology) enddo write (iout,*) "SC restraints from templates" do i=nnt,nct - write(iout,'(i5,100(4f8.2,4x))') i, + 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 @@ -3536,37 +3674,78 @@ c---------------------------------------------------------------------- character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard - integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp, - & ik,ll,ii,kk,iistart,iishift,lim_xx + integer ki, i, j, jj,k, l, ii_in_use(maxdim_cont),i_tmp, + & idomain_tmp, + & ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr, + & i01,i10,nnt_chain,nct_chain + integer itype_temp(maxres) double precision distal logical lprn /.true./ + integer nres_temp integer ilen external ilen - logical liiflag + logical liiflag,lfirst c c double precision rescore_tmp,x12,y12,z12,rescore2_tmp double precision, dimension (max_template,maxres) :: rescore double precision, dimension (max_template,maxres) :: rescore2 - character*24 pdbfile,tpl_k_rescore + character*24 tpl_k_rescore + character*256 pdbfile c c For new homol impl c include 'COMMON.VAR' c +c write (iout,*) "READ_KLAPAUCJUSZ" +c print *,"READ_KLAPAUCJUSZ" +c call flush(iout) call getenv("FRAGFILE",fragfile) + write (iout,*) "Opening", fragfile + call flush(iout) open(ientin,file=fragfile,status="old",err=10) - read(ientin,*) constr_homology,nclust - l_homo = .false. +c write (iout,*) " opened" +c call flush(iout) + sigma_theta=0.0 sigma_d=0.0 sigma_dih=0.0 + l_homo = .false. + + nres_temp=nres + itype_temp=itype + ii=0 + lim_odl=0 + +c write (iout,*) "Entering loop" +c call flush(iout) + + DO IGR = 1,NCHAIN_GROUP + +c 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 +c write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain + c Read pdb files 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', @@ -3574,7 +3753,10 @@ c Read pdb files stop 34 continue unres_pdb=.false. +c nres_temp=nres call readpdb_template(k) + nres_chomo(k)=nres +c nres=nres_temp do i=1,nres rescore(k,i)=0.2d0 rescore2(k,i)=1.0d0 @@ -3593,6 +3775,7 @@ c do ll = 1,ninclust(l) k = inclust(ll,l) +c write (iout,*) "l",l," ll",ll," k",k do i=1,nres idomain(k,i)=0 enddo @@ -3608,6 +3791,8 @@ c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) +c nres_temp=nres + nres=nres_chomo(k) do i=1,2*nres do j=1,3 c(j,i)=chomo(j,i,k) @@ -3620,10 +3805,13 @@ c write (iout,*) "c(",j,i,") =",c(j,i) thetaref(i)=theta(i) phiref(i)=phi(i) enddo +c nres=nres_temp if (waga_dist.ne.0.0d0) then - ii=0 - do i = nnt,nct-2 - do j=i+2,nct + ii=ii_old +c do i = nnt,nct-2 + do i = nnt_chain,nct_chain-2 +c 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) @@ -3640,9 +3828,9 @@ c write (iout,*) k,i,j,distal,dist2_cut c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", -c & constr_homology - ires_homo(ii)=i - jres_homo(ii)=j +c & 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 @@ -3678,50 +3866,55 @@ c c Theta, dihedral and SC retraints c if (waga_angle.gt.0.0d0) then - do i = nnt+3,nct + do i = nnt_chain+3,nct_chain + iii=i+chain_border1(1,igr)-1 if (idomain(k,i).eq.0) then c sigma_dih(k,i)=0.0 cycle endif - dih(k,i)=phiref(i) - sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ + 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 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i), c & " sigma_dihed",sigma_dih(k,i) - if (sigma_dih(k,i).ne.0) - & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) + 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 +c lim_dih=nct-nnt-2 endif if (waga_theta.gt.0.0d0) then - do i = nnt+2,nct + do i = nnt_chain+2,nct_chain + iii=i+chain_border1(1,igr)-1 if (idomain(k,i).eq.0) then c sigma_theta(k,i)=0.0 cycle endif - thetatpl(k,i)=thetaref(i) - sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ + 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,i).ne.0) - & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + 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,nct + do i = nnt_chain,nct_chain + iii=i+chain_border1(1,igr)-1 if (itype(i).eq.10) cycle if (idomain(k,i).eq.0 ) then c sigma_d(k,i)=0.0 cycle endif - xxtpl(k,i)=xxref(i) - yytpl(k,i)=yyref(i) - zztpl(k,i)=zzref(i) - sigma_d(k,i)=rescore(k,i) - if (sigma_d(k,i).ne.0) - & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) - if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 + 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)) +c if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 enddo endif enddo ! l @@ -3730,41 +3923,101 @@ c c remove distance restraints not used in any model from the list c shift data in all arrays c +c write (iout,*) "ii_old",ii_old if (waga_dist.ne.0.0d0) then - ii=0 +#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. - do i=nnt,nct-2 - do j=i+2,nct + lfirst=.true. + do i=nnt_chain,nct_chain-2 + do j=i+2,nct_chain ii=ii+1 - if (ii_in_use(ii).eq.0.and.liiflag) then +c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 +c & .and. distal.le.dist2_cut ) then +c write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii) +c 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. - iistart=ii + 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.or. - & .not.liiflag.and.ii.eq.lim_odl) then - if (ii.eq.lim_odl) then - iishift=ii-iistart+1 - else - iishift=ii-iistart - endif + if (ii_in_use(ii).ne.0.and..not.liiflag) then + i01=ii liiflag=.true. - do ki=iistart,lim_odl-iishift - ires_homo(ki)=ires_homo(ki+iishift) - jres_homo(ki)=jres_homo(ki+iishift) - ii_in_use(ki)=ii_in_use(ki+iishift) - do k=1,constr_homology - odl(k,ki)=odl(k,ki+iishift) - sigma_odl(k,ki)=sigma_odl(k,ki+iishift) - l_homo(k,ki)=l_homo(k,ki+iishift) - enddo - enddo - ii=ii-iishift - lim_odl=lim_odl-iishift 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)) +c 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=itype_temp + return 10 stop "Error in fragment file" end