X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=6b59d69374adb84fd4bc34735694dcc46fa943dd;hb=96d29f8250d30733e8a1f624f918388ef4f4050f;hp=988983ec337ae70ee0a972c5d7d76288375ce56f;hpb=4866688f3356059b9a933ec2be6da24e6784fa9f;p=unres.git diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index 988983e..6b59d69 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -8,6 +8,7 @@ include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' logical file_exist C Read force-field parameters except weights call parmread @@ -44,7 +45,8 @@ C Print restraint information &write (iout,'(a,i5,a)') "The following",nhpb-nss, & " distance constraints have been imposed" do i=nss+1,nhpb - write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i) + write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i), + & ibecarb(i),dhpb(i),dhpb1(i),forcon(i) enddo #ifdef MPI endif @@ -95,6 +97,7 @@ c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file C Set up the time limit (caution! The time must be input in minutes!) read_cart=index(controlcard,'READ_CART').gt.0 call readi(controlcard,'CONSTR_DIST',constr_dist,0) + call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes @@ -128,6 +131,7 @@ C Set up the time limit (caution! The time must be input in minutes!) sideadd=(index(controlcard,'SIDEADD').gt.0) energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) outpdb=(index(controlcard,'PDBOUT').gt.0) + outx=(index(controlcard,'XOUT').gt.0) outmol2=(index(controlcard,'MOL2OUT').gt.0) pdbref=(index(controlcard,'PDBREF').gt.0) refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0) @@ -137,7 +141,7 @@ C Set up the time limit (caution! The time must be input in minutes!) call readi(controlcard,'MAXGEN',maxgen,10000) call readi(controlcard,'MAXOVERLAP',maxoverlap,1000) call readi(controlcard,"KDIAG",kdiag,0) - call readi(controlcard,"RESCALE_MODE",rescale_mode,1) + call readi(controlcard,"RESCALE_MODE",rescale_mode,2) if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) & write (iout,*) "RESCALE_MODE",rescale_mode split_ene=index(controlcard,'SPLIT_ENE').gt.0 @@ -371,7 +375,6 @@ C endif call reada(controlcard,"Q_NP",Q_np,0.1d0) usampl = index(controlcard,"USAMPL").gt.0 - mdpdb = index(controlcard,"MDPDB").gt.0 call reada(controlcard,"T_BATH",t_bath,300.0d0) call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1) @@ -385,6 +388,11 @@ C large = index(controlcard,"LARGE").gt.0 print_compon = index(controlcard,"PRINT_COMPON").gt.0 rattle = index(controlcard,"RATTLE").gt.0 + preminim = index(controlcard,"PREMINIM").gt.0 + if (preminim) then + dccart=(index(controlcard,'CART').gt.0) + call read_minim + endif c if performing umbrella sampling, fragments constrained are read from the fragment file nset=0 if(usampl) then @@ -425,6 +433,10 @@ c if performing umbrella sampling, fragments constrained are read from the frag write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx if (rattle) write (iout,'(a60)') & "Rattle algorithm used to constrain the virtual bonds" + if (preminim .or. iranconf.gt.0) then + write (iout,'(a60)') + & "Initial structure will be energy-minimized" + endif endif reset_fricmat=1000 if (lang.gt.0) then @@ -626,6 +638,7 @@ C common /pizda/ itype_pdb logical seq_comp,fail double precision energia(0:n_ene) + integer ilen external ilen C @@ -748,6 +761,10 @@ C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WTORD',wtor_d,1.0D0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) + call reada(weightcard,'WDFAD',wdfa_dist,0.0d0) + call reada(weightcard,'WDFAT',wdfa_tor,0.0d0) + call reada(weightcard,'WDFAN',wdfa_nei,0.0d0) + call reada(weightcard,'WDFAB',wdfa_beta,0.0d0) call reada(weightcard,'SCAL14',scal14,0.4D0) call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) @@ -777,11 +794,16 @@ C 12/1/95 Added weight for the multi-body term WCORR weights(18)=scal14 weights(21)=wsccor endif + weights(25)=wdfa_dist + weights(26)=wdfa_tor + weights(27)=wdfa_nei + weights(28)=wdfa_beta 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, - & wturn4,wturn6 + & wturn4,wturn6, + & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ @@ -800,7 +822,11 @@ C 12/1/95 Added weight for the multi-body term WCORR & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ - & 'WTURN6= ',f10.6,' (turns, 6th order)') + & 'WTURN6= ',f10.6,' (turns, 6th order)'/ + & 'WDFA_D= ',f10.6,' (DFA, distance)' / + & 'WDFA_T= ',f10.6,' (DFA, torsional)' / + & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' / + & 'WDFA_B= ',f10.6,' (DFA, beta formation)') if(me.eq.king.or..not.out1file)then if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', @@ -828,7 +854,8 @@ C 12/1/95 Added weight for the multi-body term WCORR if(me.eq.king.or..not.out1file) & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, - & wturn4,wturn6 + & wturn4,wturn6, + & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta 22 format (/'Energy-term weights (scaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ @@ -847,7 +874,11 @@ C 12/1/95 Added weight for the multi-body term WCORR & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ - & 'WTURN6= ',f10.6,' (turns, 6th order)') + & 'WTURN6= ',f10.6,' (turns, 6th order)'/ + & 'WDFA_D= ',f10.6,' (DFA, distance)' / + & 'WDFA_T= ',f10.6,' (DFA, torsional)' / + & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' / + & 'WDFA_B= ',f10.6,' (DFA, beta formation)') if(me.eq.king.or..not.out1file) & write (iout,*) "Reference temperature for weights calculation:", & temp0 @@ -859,12 +890,36 @@ C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) + dyn_ss=(index(weightcard,'DYN_SS').gt.0) + do i=1,maxres + dyn_ss_mask(i)=.false. + enddo + do i=1,maxres-1 + do j=i+1,maxres + dyn_ssbond_ij(i,j)=1.0d300 + enddo + enddo + call reada(weightcard,"HT",Ht,0.0D0) + if (dyn_ss) then + ss_depth=ebr/wsc-0.25*eps(1,1) + Ht=Ht/wsc-0.25*eps(1,1) + akcm=akcm*wstrain/wsc + akth=akth*wstrain/wsc + akct=akct*wstrain/wsc + v1ss=v1ss*wstrain/wsc + v2ss=v2ss*wstrain/wsc + v3ss=v3ss*wstrain/wsc + else + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + endif + if(me.eq.king.or..not.out1file) then write (iout,*) "Parameters of the SS-bond potential:" write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth, & " AKCT",akct write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss - write (iout,*) "EBR",ebr + write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth + write (iout,*)" HT",Ht print *,'indpdb=',indpdb,' pdbref=',pdbref endif if (indpdb.gt.0 .or. pdbref) then @@ -879,6 +934,17 @@ C 12/1/95 Added weight for the multi-body term WCORR 34 continue c print *,'Begin reading pdb data' call readpdb + 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 c print *,'Finished reading pdb data' if(me.eq.king.or..not.out1file) & write (iout,'(a,i3,a,i3)')'nsup=',nsup, @@ -892,6 +958,9 @@ c print *,'Finished reading pdb data' call contact(.false.,ncont_ref,icont_ref,co) if (sideadd) then +C Following 2 lines for diagnostics; comment out if not needed + write (iout,*) "Before sideadd" + call intout if(me.eq.king.or..not.out1file) & write(iout,*)'Adding sidechains' maxsi=1000 @@ -908,7 +977,12 @@ c print *,'Finished reading pdb data' & i,' after ',nsi,' trials' endif enddo +C 10/03/12 Adam: Recalculate coordinates with new side chain positions + call chainbuild endif +C Following 2 lines for diagnostics; comment out if not needed +c write (iout,*) "After sideadd" +c call intout endif if (indpdb.eq.0) then C Read sequence if not taken from the pdb file. @@ -923,6 +997,20 @@ C Convert sequence to numeric code do i=1,nres itype(i)=rescode(i,sequence(i),iscode) enddo + if (itype(2).eq.10.and.itype(1).eq.ntyp1) then + write (iout,*) + & "Glycine is the first full residue, initial dummy deleted" + do i=1,nres + itype(i)=itype(i+1) + enddo + nres=nres-1 + endif + if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then + write (iout,*) + & "Glycine is the last full residue, terminal dummy deleted" + nres=nres-1 + endif + C Assign initial virtual bond lengths do i=2,nres vbld(i)=vbl @@ -1006,6 +1094,23 @@ C 8/13/98 Set limits to generating the dihedral angles cd print *,'NNT=',NNT,' NCT=',NCT if (itype(1).eq.21) nnt=2 if (itype(nres).eq.21) nct=nct-1 + +C Bartek:READ init_vars +C Initialize variables! +C Juyong:READ read_info +C READ fragment information!! +C both routines should be in dfa.F file!! + +#ifdef DFA + 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!' + call read_dfa_info + print*, 'read_dfa_info finished!' + endif +#endif +C if (pdbref) then if(me.eq.king.or..not.out1file) & write (iout,'(a,i3)') 'nsup=',nsup @@ -1075,11 +1180,6 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) enddo call contact(.true.,ncont_ref,icont_ref,co) endif -c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup - call flush(iout) - if (constr_dist.gt.0) call read_dist_constr -c write (iout,*) "After read_dist_constr nhpb",nhpb - call hpb_partition if(me.eq.king.or..not.out1file) & write (iout,*) 'Contact order:',co if (pdbref) then @@ -1096,6 +1196,60 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb enddo endif endif +c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup + if (constr_dist.gt.0) then + call read_dist_constr + endif + + + if (constr_homology.gt.0) then + 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)=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 + 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 + endif + + + if (nhpb.gt.0) call hpb_partition +c write (iout,*) "After read_dist_constr nhpb",nhpb +c call flush(iout) if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and. & modecalc.ne.10) then @@ -1105,6 +1259,8 @@ C initial geometry. 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,*) time,potE,uconst,t_bath, + & nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn) 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) @@ -1123,7 +1279,7 @@ C initial geometry. enddo endif enddo - return +c return else call read_angles(inp,*36) endif @@ -1207,6 +1363,8 @@ C Generate distance constraints, if the PDB structure is to be regularized. if (nthread.gt.0) then call read_threadbase endif + write (iout,*) "READRTNS: Calling setup_var" + call flush(iout) call setup_var if (me.eq.king .or. .not. out1file) & call intout @@ -1214,18 +1372,35 @@ C Generate distance constraints, if the PDB structure is to be regularized. write (iout,'(/a,i3,a)') & 'The chain contains',ns,' disulfide-bridging cysteines.' write (iout,'(20i4)') (iss(i),i=1,ns) + if (dyn_ss) then + write(iout,*)"Running with dynamic disulfide-bond formation" + else write (iout,'(/a/)') 'Pre-formed links are:' do i=1,nss i1=ihpb(i)-nres i2=jhpb(i)-nres it1=itype(i1) it2=itype(i2) - if (me.eq.king.or..not.out1file) - & write (iout,'(2a,i3,3a,i3,a,3f10.3)') + write (iout,'(2a,i3,3a,i3,a,3f10.3)') & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i), & ebr,forcon(i) enddo write (iout,'(a)') + endif + endif + 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 + call hpb_partition + do i=1,ns + dyn_ss_mask(iss(i))=.true. + enddo endif if (i2ndstr.gt.0) call secstrp2dihc c call geom_to_var(nvar,x) @@ -1289,10 +1464,12 @@ C Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') - & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & 'Do you REALLY think that the residue ', + & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') - & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & 'Do you REALLY think that the residue ', + & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,ierror) @@ -1303,7 +1480,8 @@ C Check whether the specified bridging residues are cystines. C Read preformed bridges. if (ns.gt.0) then read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss) - write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) + if(fg_rank.eq.0) + & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) if (nss.gt.0) then nhpb=nss C Check if the residues involved in bridges are in the specified list of @@ -1810,6 +1988,9 @@ c---------------------------------------------------------------------------- call readi(minimcard,'MINFUN',minfun,maxmin) call reada(minimcard,'TOLF',tolf,1.0D-2) call reada(minimcard,'RTOLF',rtolf,1.0D-4) + print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1) + print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1) + print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1) write (iout,'(/80(1h*)/20x,a/80(1h*))') & 'Options in energy minimization:' write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') @@ -2080,38 +2261,38 @@ C Get parameter filenames and open the parameter files. open (isidep,file=sidename,status='old') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old', - & readonly) + &action='read') open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown') C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown') C Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) - open (ibond,file=bondname,status='old',readonly) + open (ibond,file=bondname,status='old',action='read') call getenv_loc('THETPAR',thetname) - open (ithep,file=thetname,status='old',readonly) + open (ithep,file=thetname,status='old',action='read') #ifndef CRYST_THETA call getenv_loc('THETPARPDB',thetname_pdb) print *,"thetname_pdb ",thetname_pdb - open (ithep_pdb,file=thetname_pdb,status='old',readonly) + open (ithep_pdb,file=thetname_pdb,status='old',action='read') print *,ithep_pdb," opened" #endif call getenv_loc('ROTPAR',rotname) - open (irotam,file=rotname,status='old',readonly) + open (irotam,file=rotname,status='old',action='read') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) - open (irotam_pdb,file=rotname_pdb,status='old',readonly) + open (irotam_pdb,file=rotname_pdb,status='old',action='read') #endif call getenv_loc('TORPAR',torname) - open (itorp,file=torname,status='old',readonly) + open (itorp,file=torname,status='old',action='read') call getenv_loc('TORDPAR',tordname) - open (itordp,file=tordname,status='old',readonly) + open (itordp,file=tordname,status='old',action='read') call getenv_loc('SCCORPAR',sccorname) - open (isccor,file=sccorname,status='old',readonly) + open (isccor,file=sccorname,status='old',action='read') call getenv_loc('FOURIER',fouriername) - open (ifourier,file=fouriername,status='old',readonly) + open (ifourier,file=fouriername,status='old',action='read') call getenv_loc('ELEPAR',elename) - open (ielep,file=elename,status='old',readonly) + open (ielep,file=elename,status='old',action='read') call getenv_loc('SIDEPAR',sidename) - open (isidep,file=sidename,status='old',readonly) + open (isidep,file=sidename,status='old',action='read') #endif #ifndef OLDSCP C @@ -2126,7 +2307,7 @@ C #elif (defined G77) open (iscpp,file=scpname,status='old') #else - open (iscpp,file=scpname,status='old',readonly) + open (iscpp,file=scpname,status='old',action='read') #endif #endif call getenv_loc('PATTERN',patname) @@ -2137,7 +2318,7 @@ C #elif (defined G77) open (icbase,file=patname,status='old') #else - open (icbase,file=patname,status='old',readonly) + open (icbase,file=patname,status='old',action='read') #endif #ifdef MPI C Open output file only for CG processes @@ -2162,6 +2343,8 @@ c print *,"Processor",myrank," fg_rank",fg_rank & //'.pdb' mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// & liczba(:ilen(liczba))//'.mol2' + cartname=prefix(:lenpre)//'_'//pot(:lenpot)// + & liczba(:ilen(liczba))//'.x' statname=prefix(:lenpre)//'_'//pot(:lenpot)// & liczba(:ilen(liczba))//'.stat' if (lentmp.gt.0) @@ -2180,6 +2363,7 @@ c print *,"Processor",myrank," fg_rank",fg_rank intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int' pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb' mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2' + cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x' statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat' if (lentmp.gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) @@ -2313,6 +2497,7 @@ c------------------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.MD' + include 'COMMON.CONTROL' open(irest2,file=rest2name,status='unknown') read(irest2,*) totT,EK,potE,totE,t_bath do i=1,2*nres @@ -2321,7 +2506,7 @@ c------------------------------------------------------------------------------- do i=1,2*nres read(irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo - if(usampl) then + if(usampl.or.homol_nset.gt.1) then read (irest2,*) iset endif close(irest2) @@ -2339,33 +2524,36 @@ c------------------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.MD' include 'COMMON.CONTROL' + integer iset1 read(inp,*) nset,nfrag,npair,nfrag_back if(me.eq.king.or..not.out1file) & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair, & " nfrag_back",nfrag_back - do iset=1,nset - read(inp,*) mset(iset) + do iset1=1,nset + read(inp,*) mset(iset1) do i=1,nfrag - read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset), - & qinfrag(i,iset) + read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1), + & qinfrag(i,iset1) if(me.eq.king.or..not.out1file) - & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset), - & ifrag(2,i,iset), qinfrag(i,iset) + & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1), + & ifrag(2,i,iset1), qinfrag(i,iset1) enddo do i=1,npair - read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset), - & qinpair(i,iset) + read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1), + & qinpair(i,iset1) if(me.eq.king.or..not.out1file) - & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset), - & ipair(2,i,iset), qinpair(i,iset) + & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1), + & ipair(2,i,iset1), qinpair(i,iset1) enddo do i=1,nfrag_back - read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset), - & wfrag_back(3,i,iset), - & ifrag_back(1,i,iset),ifrag_back(2,i,iset) + read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1), + & wfrag_back(3,i,iset1), + & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1) if(me.eq.king.or..not.out1file) - & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset), - & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset) + & write(iout,*) "A",i,wfrag_back(1,i,iset1), + & wfrag_back(2,i,iset1), + & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1), + & ifrag_back(2,i,iset1) enddo enddo return @@ -2393,6 +2581,7 @@ c call flush(iout) call readi(controlcard,"NPAIR",npair_,0) call readi(controlcard,"NDIST",ndist_,0) call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) + call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0) call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0) call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0) call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0) @@ -2406,6 +2595,18 @@ c write (iout,*) "IPAIR" c do i=1,npair_ c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) c enddo + if (.not.refstr .and. nfrag.gt.0) then + write (iout,*) + & "ERROR: no reference structure to compute distance restraints" + write (iout,*) + & "Restraints must be specified explicitly (NDIST=number)" + stop + endif + if (nfrag.lt.2 .and. npair.gt.0) then + write (iout,*) "ERROR: Less than 2 fragments specified", + & " but distance restraints between pairs requested" + stop + endif call flush(iout) do i=1,nfrag_ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup @@ -2416,7 +2617,7 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) if (wfrag_(i).gt.0.0d0) then do j=ifrag_(1,i),ifrag_(2,i)-1 do k=j+1,ifrag_(2,i) - write (iout,*) "j",j," k",k +c write (iout,*) "j",j," k",k ddjk=dist(j,k) if (constr_dist.eq.1) then nhpb=nhpb+1 @@ -2480,24 +2681,437 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) endif enddo do i=1,ndist_ - read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1) if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 - dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + if (dhpb(nhpb).eq.0.0d0) + & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + endif + enddo #ifdef MPI - if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", - & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) -#else - write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", - & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + if (.not.out1file .or. me.eq.king) then #endif - endif + do i=1,nhpb + write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ", + & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i) enddo call flush(iout) +#ifdef MPI + endif +#endif return end c------------------------------------------------------------------------------- + + subroutine read_constr_homology + + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.MD' + include 'COMMON.GEO' + include 'COMMON.INTERACT' +c +c For new homol impl +c + include 'COMMON.VAR' +c + +c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d, +c & dist_cut +c common /przechowalnia/ odl_temp(maxres,maxres,max_template), +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, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp + logical lprn /.true./ + integer ilen + external ilen +c +c FP - Nov. 2014 Temporary specifications for new vars +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 +c ----------------------------------------------------------------- +c Reading multiple PDB ref structures and calculation of retraints +c not using pre-computed ones stored in files model_ki_{dist,angle} +c FP (Nov., 2014) +c ----------------------------------------------------------------- +c +c +c Alternative: reading from input + call card_concat(controlcard) + 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 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 + write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA' + start_from_model=.false. + endif + if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON' + if(start_from_model .and. rest) then + write(iout,*) 'START_FROM_MODELS is OFF' + write(iout,*) 'remove restart keyword from input' + endif + 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 " + 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 + +cd write (iout,*) "nnt",nnt," nct",nct +cd call flush(iout) + + + lim_odl=0 + lim_dih=0 +c +c New +c + lim_theta=0 + lim_xx=0 +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 + + do k=1,constr_homology + + read(inp,'(a)') pdbfile +c Next stament causes error upon compilation (?) +c if(me.eq.king.or. .not. out1file) +c write (iout,'(2a)') 'PDB data will be read from file ', +c & pdbfile(:ilen(pdbfile)) + 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 +c print *,'Begin reading pdb data' +c +c Files containing res sim or local scores (former containing sigmas) +c + + write(kic2,'(bz,i2.2)') k + + tpl_k_rescore="template"//kic2//".sco" + + unres_pdb=.false. + if (read2sigma) then + call readpdb_template(k) + else + call readpdb + endif +c +c Distance restraints +c +c ... --> odl(k,ii) +C Copy the coordinates from reference coordinates (?) + do i=1,2*nres + do j=1,3 + c(j,i)=cref(j,i) +c write (iout,*) "c(",j,i,") =",c(j,i) + enddo + enddo +c +c From read_dist_constr (commented out 25/11/2014 <-> res sim) +c +c 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,maxdim ! loop for reading res sim + if (read2sigma) then + read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_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 + else + idomain(k,irec)=1 + read (ientin,*,end=1401) rescore_tmp + +c 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 +c 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) + + + 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. + +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 + 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 +c +c Theta, dihedral and SC retraints +c + if (waga_angle.gt.0.0d0) then +c open (ientin,file=tpl_k_sigma_dih,status='old') +c do irec=1,maxres-3 ! loop for reading sigma_dih +c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for? +c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right? +c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity +c & sigma_dih(k,i+nnt-1) +c enddo +c1402 continue +c 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? +c read (ientin,*) sigma_dih(k,i) ! original variant +c write (iout,*) "dih(",k,i,") =",dih(k,i) +c write(iout,*) "rescore(",k,i,") =",rescore(k,i), +c & "rescore(",k,i-1,") =",rescore(k,i-1), +c & "rescore(",k,i-2,") =",rescore(k,i-2), +c & "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 +c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0 +c write (iout,*) "Raw sigmas for dihedral angle restraints" +c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i) +c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* +c rescore(k,i-2)*rescore(k,i-3) ! right expression ? +c Instead of res sim other local measure of b/b str reliability possible + sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) +c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i) + enddo + lim_dih=nct-nnt-2 + endif + + if (waga_theta.gt.0.0d0) then +c open (ientin,file=tpl_k_sigma_theta,status='old') +c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds? +c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for? +c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity +c & sigma_theta(k,i+nnt-1) +c enddo +c1403 continue +c close (ientin) + + do i = nnt+2,nct ! right? without parallel. +c do i = i=1,nres ! alternative for bounds acc to readpdb? +c 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) +c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i) +c write(iout,*) "rescore(",k,i,") =",rescore(k,i), +c & "rescore(",k,i-1,") =",rescore(k,i-1), +c & "rescore(",k,i-2,") =",rescore(k,i-2) +c 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 +c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 + sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + +c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* +c rescore(k,i-2) ! right expression ? +c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) + enddo + endif + lim_theta=nct-nnt-1 + + if (waga_d.gt.0.0d0) then +c open (ientin,file=tpl_k_sigma_d,status='old') +c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds? +c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for? +c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity +c & sigma_d(k,i+nnt-1) +c enddo +c1404 continue + + do i = nnt,nct ! right? without parallel. +c do i=2,nres-1 ! alternative for bounds acc to readpdb? +c do i=loc_start,loc_end ! with FG parallel. + if (itype(i).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) +c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i) +c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i) +c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i) +c write(iout,*) "rescore(",k,i,") =",rescore(k,i) + sigma_d(k,i)=rescore(k,i) ! right expression ? + sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) + +c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ? +c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i) +c read (ientin,*) sigma_d(k,i) ! 1st variant + if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right? + enddo + lim_xx=nct-nnt+1 + endif + enddo +c +c remove distance restraints not used in any model from the list +c shift data in all arrays +c + if (waga_dist.ne.0.0d0) then + ii=0 + do i=nnt,nct-2 + do j=i+2,nct + ii=ii+1 + if (ii_in_use(ii).eq.0) then + do ki=ii,lim_odl-1 + ires_homo(ki)=ires_homo(ki+1) + jres_homo(ki)=jres_homo(ki+1) + ii_in_use(ki)=ii_in_use(ki+1) + do k=1,constr_homology + odl(k,ki)=odl(k,ki+1) + sigma_odl(k,ki)=sigma_odl(k,ki+1) + l_homo(k,ki)=l_homo(k,ki+1) + enddo + enddo + ii=ii-1 + lim_odl=lim_odl-1 + endif + enddo + enddo + endif + if (constr_homology.gt.0) call homology_partition + if (constr_homology.gt.0) call init_int_table +cd write (iout,*) "homology_partition: lim_theta= ",lim_theta, +cd & "lim_xx=",lim_xx +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 + if (.not.lprn) return +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,10(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,lim_dih + write (iout,'(i5,10(2f8.2,4x))') 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,lim_theta + write (iout,'(i5,10(2f8.2,4x))') 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,lim_xx + write(iout,'(i5,10(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 +c ----------------------------------------------------------------- + return + end +c---------------------------------------------------------------------- + #ifdef WINIFL subroutine flush(iu) return @@ -2509,6 +3123,7 @@ c------------------------------------------------------------------------------- return end #endif + c------------------------------------------------------------------------------ subroutine copy_to_tmp(source) include "DIMENSIONS"