X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Freadrtns_CSA.F;h=d76b29e0bf84b2232c41dc8c0975b47ccd4ce87a;hb=57038e4bdff4cc9534106b25bfbd4b9a844d47fd;hp=3c6fb51bde7ed0ec3ab7a35230ec3b43f8c2a68e;hpb=48ae9e01d2dd6571fa2cca6c704dc04f86e5fd7b;p=unres.git diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index 3c6fb51..d76b29e 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -171,13 +171,14 @@ c call readi(controlcard,'IZ_SC',iz_sc,0) timlim=60.0D0*timlim safety = 60.0d0*safety modecalc=0 + call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100) call reada(controlcard,"T_BATH",t_bath,300.0d0) minim=(index(controlcard,'MINIMIZE').gt.0) dccart=(index(controlcard,'CART').gt.0) overlapsc=(index(controlcard,'OVERLAP').gt.0) overlapsc=.not.overlapsc - searchsc=(index(controlcard,'NOSEARCHSC').gt.0) - searchsc=.not.searchsc + searchsc=(index(controlcard,'SEARCHSC').gt.0 .and. + & index(controlcard,'NOSEARCHSC').eq.0) sideadd=(index(controlcard,'SIDEADD').gt.0) energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) mremd_dec=(index(controlcard,'MREMD_DEC').gt.0) @@ -186,7 +187,6 @@ 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" @@ -295,6 +295,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 @@ -304,9 +310,16 @@ C Reading the dimensions of box in x,y,z coordinates call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) + write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize c Cutoff range for interactions - call reada(controlcard,"R_CUT",r_cut,15.0d0) + call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0) + call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) + write (iout,*) "Cutoff on interactions",r_cut_int + write (iout,*) + & "Cutoff in switching short and long range interactions in RESPA", + & r_cut_respa + write (iout,*) "lambda in switch function",rlamb call reada(controlcard,"LIPTHICK",lipthick,0.0d0) call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) if (lipthick.gt.0.0d0) then @@ -539,12 +552,16 @@ c if performing umbrella sampling, fragments constrained are read from the frag & "Initial time step of numerical integration:",d_time, & " natural units" write (iout,'(60x,f10.5,a)') d_time*48.9," fs" + write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int, + & " A" + write(iout,'(a60,i5)')"Frequency of updating interaction list", + & imatupdate if (RESPA) then write (iout,'(2a,i4,a)') & "A-MTS algorithm used; initial time step for fast-varying", & " short-range forces split into",ntime_split," steps." write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff", - & r_cut," lambda",rlamb + & r_cut_respa," lambda",rlamb endif write (iout,'(2a,f10.5)') & "Maximum acceleration threshold to reduce the time step", @@ -724,7 +741,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 @@ -819,7 +836,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 @@ -832,7 +849,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 @@ -858,9 +875,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 @@ -1085,10 +1102,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 @@ -1151,6 +1168,49 @@ 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 + call readpdb + 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,2i5,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 @@ -1160,14 +1220,15 @@ 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" +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) @@ -1393,7 +1454,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. @@ -1568,12 +1635,11 @@ C Set up variable list. nphi=nres-3 nvar=ntheta+nphi nside=0 - write (iout,*) "SETUP_VAR ialph" do i=2,nres-1 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then - nside=nside+1 + nside=nside+1 ialph(i,1)=nvar+nside - ialph(nside,2)=i + ialph(nside,2)=i endif enddo if (indphi.gt.0) then @@ -1583,7 +1649,6 @@ C Set up variable list. else nvar=nvar+2*nside endif - write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1) return end c---------------------------------------------------------------------------- @@ -1604,9 +1669,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 @@ -1617,11 +1684,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:' @@ -2365,10 +2435,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 @@ -2458,7 +2528,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 @@ -2969,10 +3039,11 @@ c & sigma_odl_temp(maxres,maxres,max_template) 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 + & 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 @@ -3005,15 +3076,17 @@ 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) @@ -3078,17 +3151,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) @@ -3178,6 +3254,8 @@ c & constr_homology enddo lim_odl=ii endif +c write (iout,*) "Distance restraints set" +c call flush(iout) c c Theta, dihedral and SC retraints c @@ -3218,6 +3296,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') @@ -3253,6 +3333,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') @@ -3288,51 +3370,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 @@ -3530,6 +3627,7 @@ c---------------------------------------------------------------------- & ik,ll,ii,kk,iistart,iishift,lim_xx double precision distal logical lprn /.true./ + integer nres_temp integer ilen external ilen logical liiflag @@ -3564,7 +3662,10 @@ c Read pdb files 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 @@ -3598,6 +3699,8 @@ c Distance restraints c c ... --> odl(k,ii) C 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) @@ -3605,11 +3708,12 @@ c write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo call int_from_cart(.true.,.false.) - call sc_loc_geom(.true.) + 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=0 do i = nnt,nct-2