X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Freadrtns.F;h=21a2b7658ec0f659c85e02198f22913424684012;hb=e12c8db731e266c6e38b2f5883f3ae7c15bbfb98;hp=bb1e950541325cad99c29be2250c94cd5d7ad5dc;hpb=f7ed201184254bb7be2333857d07a78183752d45;p=unres.git diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index bb1e950..21a2b76 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -15,12 +15,14 @@ C include 'COMMON.FFIELD' include 'COMMON.FREE' include 'COMMON.INTERACT' + include "COMMON.SPLITELE" + include 'COMMON.SHIELD' character*320 controlcard,ucase #ifdef MPL include 'COMMON.INFO' #endif - integer i - + integer i,i1,i2,it1,it2 + double precision pi read (INP,'(a80)') titel call card_concat(controlcard) @@ -28,6 +30,63 @@ C call readi(controlcard,'RESCALE',rescale_mode,2) call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0) write (iout,*) "DISTCHAINMAX",distchainmax +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) +c Cutoff range for interactions + call reada(controlcard,"R_CUT",r_cut,15.0d0) + call reada(controlcard,"LAMBDA",rlamb,0.3d0) + call reada(controlcard,"LIPTHICK",lipthick,0.0d0) + call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) + if (lipthick.gt.0.0d0) then + bordliptop=(boxzsize+lipthick)/2.0 + bordlipbot=bordliptop-lipthick +C endif + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) + & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" + buflipbot=bordlipbot+lipbufthick + bufliptop=bordliptop-lipbufthick + if ((lipbufthick*2.0d0).gt.lipthick) + &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" + endif + write(iout,*) "bordliptop=",bordliptop + write(iout,*) "bordlipbot=",bordlipbot + write(iout,*) "bufliptop=",bufliptop + write(iout,*) "buflipbot=",buflipbot +C Shielding mode + call readi(controlcard,'SHIELD',shield_mode,0) + write (iout,*) "SHIELD MODE",shield_mode + call readi(controlcard,'TUBEMOD',tubelog,0) + write (iout,*) TUBElog,"TUBEMODE" + if (shield_mode.gt.0) then + pi=3.141592d0 +C VSolvSphere the volume of solving sphere +C print *,pi,"pi" +C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +C there will be no distinction between proline peptide group and normal peptide +C group in case of shielding parameters + VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 + write (iout,*) VSolvSphere,VSolvSphere_div +C long axis of side chain + do i=1,ntyp + long_r_sidechain(i)=vbldsc0(1,i) + short_r_sidechain(i)=sigma0(i) + enddo + buff_shield=1.0d0 + endif + if (TUBElog.gt.0) then + call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) + call reada(controlcard,"YTUBE",tubecenter(2),0.0d0) + call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0) + call reada(controlcard,"RTUBE",tubeR0,0.0d0) + call reada(controlcard,"TUBETOP",bordtubetop,boxzsize) + call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0) + call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0) + buftubebot=bordtubebot+tubebufthick + buftubetop=bordtubetop-tubebufthick + endif call readi(controlcard,'PDBOUT',outpdb,0) call readi(controlcard,'MOL2OUT',outmol2,0) refstr=(index(controlcard,'REFSTR').gt.0) @@ -35,6 +94,13 @@ C pdbref=(index(controlcard,'PDBREF').gt.0) iscode=index(controlcard,'ONE_LETTER') tree=(index(controlcard,'MAKE_TREE').gt.0) + with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0 + call readi(controlcard,'CONSTR_DIST',constr_dist,0) + write (iout,*) "with_dihed_constr ",with_dihed_constr, + & " CONSTR_DIST",constr_dist + with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 + write (iout,*) "with_theta_constr ",with_theta_constr + call flush(iout) min_var=(index(controlcard,'MINVAR').gt.0) plot_tree=(index(controlcard,'PLOT_TREE').gt.0) punch_dist=(index(controlcard,'PUNCH_DIST').gt.0) @@ -87,6 +153,8 @@ C include 'COMMON.CONTROL' include 'COMMON.CONTACTS' include 'COMMON.TIME1' + include 'COMMON.TORCNSTR' + include 'COMMON.SHIELD' #ifdef MPL include 'COMMON.INFO' #endif @@ -96,13 +164,15 @@ C double precision x(maxvar) integer itype_pdb(maxres) logical seq_comp - integer i,j,kkk + integer i,j,kkk,i1,i2,it1,it2 C C Body C C Read weights of the subsequent energy terms. call card_concat(weightcard) - call reada(weightcard,'WSC',wsc,1.0d0) + write(iout,*) weightcard +C call reada(weightcard,'WSC',wsc,1.0d0) + write(iout,*) wsc call reada(weightcard,'WLONG',wsc,wsc) call reada(weightcard,'WSCP',wscp,1.0d0) call reada(weightcard,'WELEC',welec,1.0D0) @@ -134,6 +204,12 @@ C Read weights of the subsequent energy terms. call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) + call reada(weightcard,'WSHIELD',wshield,1.0d0) + write(iout,*) 'WSHIELD',wshield + call reada(weightcard,'WTUBE',wtube,0.0d0) + write(iout,*) 'WTUBE',wtube + call reada(weightcard,"LIPSCALE",lipscale,1.3D0) + call reada(weightcard,'WLT',wliptran,0.0D0) call reada(weightcard,"ATRISS",atriss,0.301D0) call reada(weightcard,"BTRISS",btriss,0.021D0) call reada(weightcard,"CTRISS",ctriss,1.001D0) @@ -142,6 +218,10 @@ C Read weights of the subsequent energy terms. write (iout,*) "BTRISS=", btriss write (iout,*) "CTRISS=", ctriss write (iout,*) "DTRISS=", dtriss + if (shield_mode.gt.0) then + lipscale=0.0d0 + write (iout,*) "WARNING:liscale not used in shield mode" + endif dyn_ss=(index(weightcard,'DYN_SS').gt.0) do i=1,maxres dyn_ss_mask(i)=.false. @@ -278,6 +358,66 @@ C Convert sequence to numeric code print *,'Call Read_Bridge.' call read_bridge +C this fragment reads diheadral constrains + if (with_dihed_constr) then + + read (inp,*) ndih_constr + if (ndih_constr.gt.0) then +C read (inp,*) ftors +C write (iout,*) 'FTORS',ftors +C ftors is the force constant for torsional quartic constrains + read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), + & i=1,ndih_constr) + write (iout,*) + & 'There are',ndih_constr,' constraints on phi angles.' + do i=1,ndih_constr + write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), + & ftors(i) + enddo + do i=1,ndih_constr + phi0(i)=deg2rad*phi0(i) + drange(i)=deg2rad*drange(i) + enddo + endif ! endif ndif_constr.gt.0 + endif ! with_dihed_constr + if (with_theta_constr) then +C with_theta_constr is keyword allowing for occurance of theta constrains + read (inp,*) ntheta_constr +C ntheta_constr is the number of theta constrains + if (ntheta_constr.gt.0) then +C read (inp,*) ftors + read (inp,*) (itheta_constr(i),theta_constr0(i), + & theta_drange(i),for_thet_constr(i), + & i=1,ntheta_constr) +C the above code reads from 1 to ntheta_constr +C itheta_constr(i) residue i for which is theta_constr +C theta_constr0 the global minimum value +C theta_drange is range for which there is no energy penalty +C for_thet_constr is the force constant for quartic energy penalty +C E=k*x**4 +C if(me.eq.king.or..not.out1file)then + write (iout,*) + & 'There are',ntheta_constr,' constraints on phi angles.' + do i=1,ntheta_constr + write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), + & theta_drange(i), + & for_thet_constr(i) + enddo +C endif + do i=1,ntheta_constr + theta_constr0(i)=deg2rad*theta_constr0(i) + theta_drange(i)=deg2rad*theta_drange(i) + enddo +C if(me.eq.king.or..not.out1file) +C & write (iout,*) 'FTORS',ftors +C do i=1,ntheta_constr +C ii = itheta_constr(i) +C thetabound(1,ii) = phi0(i)-drange(i) +C thetabound(2,ii) = phi0(i)+drange(i) +C enddo + endif ! ntheta_constr.gt.0 + endif! with_theta_constr + nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT @@ -357,9 +497,9 @@ c endif endif call contact(.true.,ncont_ref,icont_ref) endif - if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then - write (iout,'(/a,i3,a)') - & 'The chain contains',ns,' disulfide-bridging cysteines.' + if (ns.gt.0) then +C write (iout,'(/a,i3,a)') +C & '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" @@ -391,6 +531,11 @@ c endif dyn_ss_mask(iss(i))=.true. enddo endif +c Read distance restraints + if (constr_dist.gt.0) then + call read_dist_constr + call hpb_partition + endif return end c----------------------------------------------------------------------------- @@ -481,8 +626,8 @@ C bridging residues. enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 20 continue - dhpb(i)=dbr - forcon(i)=fbr +C dhpb(i)=dbr +C forcon(i)=fbr enddo do i=1,nss ihpb(i)=ihpb(i)+nres @@ -563,6 +708,25 @@ c---------------------------------------------------------------------------- read (rekord(iread:),*) wartosc return end +C---------------------------------------------------------------------- + subroutine multreadi(rekord,lancuch,tablica,dim,default) + implicit none + integer dim,i + integer tablica(dim),default + character*(*) rekord,lancuch + character*80 aux + integer ilen,iread + external ilen + do i=1,dim + tablica(i)=default + enddo + iread=index(rekord,lancuch(:ilen(lancuch))//"=") + if (iread.eq.0) return + iread=iread+ilen(lancuch)+1 + read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) + 10 return + end + c---------------------------------------------------------------------------- subroutine card_concat(card) include 'DIMENSIONS' @@ -654,6 +818,11 @@ C Get parameter filenames and open the parameter files. open (isidep1,file=sidepname,status="old") call getenv('SCCORPAR',sccorname) open (isccor,file=sccorname,status="old") + call getenv('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old') + call getenv('TUBEPAR',tubename) + open (itube,file=tubename,status='old') + #ifndef OLDSCP C C 8/9/01 In the newest version SCp interaction constants are read from a file @@ -664,3 +833,131 @@ C #endif return end + subroutine read_dist_constr + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.SBRIDGE' + integer ifrag_(2,100),ipair_(2,100) + double precision wfrag_(100),wpair_(100) + character*500 controlcard + logical lprn /.true./ + write (iout,*) "Calling read_dist_constr" +C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup +C call flush(iout) + write(iout,*) "TU sie wywalam?" + call card_concat(controlcard) + write (iout,*) controlcard + call flush(iout) + call readi(controlcard,"NFRAG",nfrag_,0) + call readi(controlcard,"NPAIR",npair_,0) + call readi(controlcard,"NDIST",ndist_,0) + call reada(controlcard,'DIST_CUT',dist_cut,5.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) + call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0) + write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ + write (iout,*) "IFRAG" + do i=1,nfrag_ + write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) + enddo + write (iout,*) "IPAIR" + do i=1,npair_ + write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) + enddo + call flush(iout) + do i=1,nfrag_ + if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup + if (ifrag_(2,i).gt.nstart_sup+nsup-1) + & ifrag_(2,i)=nstart_sup+nsup-1 +c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) + call flush(iout) + 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 + ddjk=dist(j,k) + if (constr_dist.eq.1) then + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + else if (constr_dist.eq.2) then + if (ddjk.le.dist_cut) then + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + endif + else + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) + endif + if (lprn) + & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + enddo + enddo + endif + enddo + do i=1,npair_ + if (wpair_(i).gt.0.0d0) then + ii = ipair_(1,i) + jj = ipair_(2,i) + if (ii.gt.jj) then + itemp=ii + ii=jj + jj=itemp + endif + do j=ifrag_(1,ii),ifrag_(2,ii) + do k=ifrag_(1,jj),ifrag_(2,jj) + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + forcon(nhpb)=wpair_(i) + dhpb(nhpb)=dist(j,k) + write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + enddo + enddo + endif + enddo + do i=1,ndist_ + if (constr_dist.eq.11) then + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) + fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) +C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", +C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + else + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) + endif + if (forcon(nhpb+1).gt.0.0d0) then + nhpb=nhpb+1 + 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)) +C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + endif +C endif + enddo + call hpb_partition + call flush(iout) + return + end