X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Freadrtns.F;h=4e0ee5642574068b3668b5a6c1f8bfb508e5b7cc;hb=34d3ad3987785642be58fb2f26557d3314215577;hp=a723920328aa803a509793eed0543cf63d196edb;hpb=f690e8b70bab14132839afebf080d4a28363b226;p=unres.git diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index a723920..4e0ee56 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -15,30 +15,92 @@ 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) + energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) + call readi(controlcard,'TORMODE',tor_mode,0) call readi(controlcard,'NRES',nres,0) 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 + 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 call readi(controlcard,'PDBOUT',outpdb,0) call readi(controlcard,'MOL2OUT',outmol2,0) refstr=(index(controlcard,'REFSTR').gt.0) - write (iout,*) "REFSTR",refstr pdbref=(index(controlcard,'PDBREF').gt.0) + refstr = refstr .or. pdbref + write (iout,*) "REFSTR",refstr," PDBREF",pdbref 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) - call readi(controlcard,'NCUT',ncut,1) + call readi(controlcard,'NCUT',ncut,0) + if (ncut.gt.0) then + call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) + nclust=0 + else + call readi(controlcard,'NCLUST',nclust,5) + endif call readi(controlcard,'SYM',symetr,1) write (iout,*) 'sym', symetr call readi(controlcard,'NSTART',nstart,0) @@ -49,10 +111,9 @@ C lgrp=(index(controlcard,'LGRP').gt.0) caonly=(index(controlcard,'CA_ONLY').gt.0) print_dist=(index(controlcard,'PRINT_DIST').gt.0) - call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) call readi(controlcard,'IOPT',iopt,2) lside = index(controlcard,"SIDE").gt.0 - efree = index(controlcard,"EFREE").gt.0 + lefree = index(controlcard,"EFREE").gt.0 call readi(controlcard,'NTEMP',nT,1) write (iout,*) "nT",nT call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0) @@ -87,16 +148,20 @@ C include 'COMMON.CONTROL' include 'COMMON.CONTACTS' include 'COMMON.TIME1' + include 'COMMON.TORCNSTR' + include 'COMMON.SHIELD' #ifdef MPL include 'COMMON.INFO' #endif character*4 sequence(maxres) - character*800 weightcard + character*800 weightcard,controlcard integer rescode double precision x(maxvar) + double precision phihel,phibet,sigmahel,sigmabet,sumv, + & secprob(3,maxres) integer itype_pdb(maxres) logical seq_comp - integer i,j + integer i,j,kkk,i1,i2,it1,it2 C C Body C @@ -115,6 +180,7 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'WTURN4',wturn4,1.0D0) call reada(weightcard,'WTURN6',wturn6,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) + call reada(weightcard,'WSCCOR',wsccor,1.0D0) call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) call reada(weightcard,'WTORD',wtor_d,1.0D0) @@ -125,6 +191,54 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) if (index(weightcard,'SOFT').gt.0) ipot=6 + call reada(weightcard,"D0CM",d0cm,3.78d0) + call reada(weightcard,"AKCM",akcm,15.1d0) + call reada(weightcard,"AKTH",akth,11.0d0) + call reada(weightcard,"AKCT",akct,12.0d0) + call reada(weightcard,"V1SS",v1ss,-1.08d0) + 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,'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) + call reada(weightcard,"DTRISS",dtriss,1.001D0) + write (iout,*) "ATRISS=", atriss + write (iout,*) "BTRISS=", btriss + write (iout,*) "CTRISS=", ctriss + write (iout,*) "DTRISS=", dtriss + 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 + write (iout,'(/a)') "Disulfide bridge parameters:" + write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr + write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm + write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct + write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, + & ' v3ss:',v3ss + C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) if (wcorr4.gt.0.0d0) wcorr=wcorr4 @@ -146,9 +260,10 @@ C 12/1/95 Added weight for the multi-body term WCORR weights(16)=wvdwpp weights(17)=wbond weights(18)=scal14 + weights(19)=wsccor write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3, - & wturn4,wturn6 + & wturn4,wturn6,wsccor 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ @@ -166,7 +281,9 @@ C 12/1/95 Added weight for the multi-body term WCORR & 'WCORR6= ',f10.6,' (multi-body 6th order)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ - & 'WTURN6= ',f10.6,' (turns, 6th order)') + & 'WTURN6= ',f10.6,' (turns, 6th order)'/ + & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)') + if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', & 'between contact pairs of peptide groups' @@ -190,7 +307,7 @@ C 12/1/95 Added weight for the multi-body term WCORR enddo call flush(iout) - print *,'indpdb=',indpdb,' pdbref=',pdbref +c print *,'indpdb=',indpdb,' pdbref=',pdbref C Read sequence if not taken from the pdb file. if (iscode.gt.0) then @@ -202,20 +319,20 @@ C Convert sequence to numeric code do i=1,nres itype(i)=rescode(i,sequence(i),iscode) enddo - print *,nres - print '(20i4)',(itype(i),i=1,nres) +c print *,nres +c print '(20i4)',(itype(i),i=1,nres) do i=1,nres #ifdef PROCOR - if (itype(i).eq.21 .or. itype(i+1).eq.21) then + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then #else - if (itype(i).eq.21) then + if (itype(i).eq.ntyp1) then #endif itel(i)=0 #ifdef PROCOR - else if (itype(i+1).ne.20) then + else if (iabs(itype(i+1)).ne.20) then #else - else if (itype(i).ne.20) then + else if (iabs(itype(i)).ne.20) then #endif itel(i)=1 else @@ -227,17 +344,125 @@ C Convert sequence to numeric code write (iout,*) i,itype(i),itel(i) enddo - print *,'Call Read_Bridge.' +c print *,'Call Read_Bridge.' call read_bridge +C this fragment reads diheadral constrains nnt=1 nct=nres - print *,'NNT=',NNT,' NCT=',NCT - if (itype(1).eq.21) nnt=2 - if (itype(nres).eq.21) nct=nct-1 +c print *,'NNT=',NNT,' NCT=',NCT + if (itype(1).eq.ntyp1) nnt=2 + if (itype(nres).eq.ntyp1) nct=nct-1 if (nstart.lt.nnt) nstart=nnt if (nend.gt.nct .or. nend.eq.0) nend=nct write (iout,*) "nstart",nstart," nend",nend nres0=nres + if (with_dihed_constr) then + + read (inp,*) ndih_constr + if (ndih_constr.gt.0) then + raw_psipred=.false. +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 + else if (ndih_constr.lt.0) then + raw_psipred=.true. + call card_concat(controlcard) + call reada(controlcard,"PHIHEL",phihel,50.0D0) + call reada(controlcard,"PHIBET",phibet,180.0D0) + call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0) + call reada(controlcard,"SIGMABET",sigmabet,40.0d0) + call reada(controlcard,"WDIHC",wdihc,0.591d0) + write (iout,*) "Weight of the dihedral restraint term",wdihc + read(inp,'(9x,3f7.3)') + & (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 + 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 + 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 + phibound(1,ndih_constr)=phihel*deg2rad + phibound(2,ndih_constr)=phibet*deg2rad + sdihed(1,ndih_constr)=sigmahel*deg2rad + sdihed(2,ndih_constr)=sigmabet*deg2rad + endif + enddo + write (iout,*) + & 'There are',ndih_constr, + & ' bimodal restraints on gamma angles.' + do i=1,ndih_constr + write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i, + & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2, + & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1, + & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg, + & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg, + & (vpsipred(j,i),j=1,3) + 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 + c if (pdbref) then c read(inp,'(a)') pdbfile c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile @@ -299,13 +524,53 @@ c endif nstart_sup=nnt nstart_seq=nnt nsup=nct-nnt+1 + kkk=1 do i=1,2*nres do j=1,3 - cref(j,i)=c(j,i) + cref(j,i,kkk)=c(j,i) enddo enddo endif - call contact(.true.,ncont_ref,icont_ref) +c call contact(.true.,ncont_ref,icont_ref) + endif + 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" + 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) + 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 +c Read distance restraints + if (constr_dist.gt.0) then + call read_dist_constr + call hpb_partition endif return end @@ -346,15 +611,17 @@ 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 print *,'ns=',ns C Check whether the specified bridging residues are cystines. do i=1,ns - if (itype(iss(i)).ne.1) then + if (iabs(itype(iss(i))).ne.1) then 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 MPL call mp_stopall(error_msg) @@ -395,8 +662,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 @@ -477,6 +744,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' @@ -556,8 +842,10 @@ C Get parameter filenames and open the parameter files. open (irotam,file=rotname,status='old') call getenv('TORPAR',torname) open (itorp,file=torname,status='old') +#ifndef NEWCORR call getenv('TORDPAR',tordname) open (itordp,file=tordname,status='old') +#endif call getenv('FOURIER',fouriername) open (ifourier,file=fouriername,status='old') call getenv('ELEPAR',elename) @@ -568,6 +856,8 @@ 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') #ifndef OLDSCP C C 8/9/01 In the newest version SCp interaction constants are read from a file @@ -578,3 +868,163 @@ C #endif return end +c-------------------------------------------------------------------------- + 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 normalize + 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) + normalize = index(controlcard,"NORMALIZE").gt.0 + 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 + write (iout,*) "Distance restraints as read from input" + 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) +c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) + if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle + nhpb=nhpb+1 + write (iout,'(a,4i5,2f8.2,2f10.5)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), + & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb) + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + else +C print *,"in else" + 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 + 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 + write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) + endif +C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) +C if (forcon(nhpb+1).gt.0.0d0) then +C nhpb=nhpb+1 +C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + enddo + if (constr_dist.eq.11 .and. normalize) then + fordepthmax=fordepth(1) + do i=2,nhpb + if (fordepth(i).gt.fordepthmax) fordepthmax=fordepth(i) + enddo + do i=1,nhpb + fordepth(i)=fordepth(i)/fordepthmax + enddo + write (iout,'(/a/4a5,2a8,2a10)') + & "Normalized Lorenzian-like distance restraints", + & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V" + do i=1,nhpb + write (iout,'(4i5,2f8.2,2f10.5)')i,ihpb(i),jhpb(i),ibecarb(i), + & dhpb(i),dhpb1(i),forcon(i),fordepth(i) + enddo + endif + write (iout,*) + call hpb_partition + call flush(iout) + return + end