X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fmolread_zs.F;h=b39002c56791c98440a93a02ccebc6cc588e7c42;hb=34d3ad3987785642be58fb2f26557d3314215577;hp=885c57bfa9d730988c43f6d05d32ccb75ca30486;hpb=f690e8b70bab14132839afebf080d4a28363b226;p=unres.git diff --git a/source/wham/src-M/molread_zs.F b/source/wham/src-M/molread_zs.F index 885c57b..b39002c 100644 --- a/source/wham/src-M/molread_zs.F +++ b/source/wham/src-M/molread_zs.F @@ -22,6 +22,7 @@ C character*320 controlcard,ucase dimension itype_pdb(maxres) logical seq_comp + double precision secprob(3,maxdih_constr),phihel,phibet call card_concat(controlcard,.true.) call reada(controlcard,'SCAL14',scal14,0.4d0) call reada(controlcard,'SCALSCP',scalscp,1.0d0) @@ -52,54 +53,148 @@ C Convert sequence to numeric code write (iout,'(20i4)') (itype(i),i=1,nres) do i=1,nres-1 #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 itel(i)=2 endif enddo + write (iout,*) "ITEL" + do i=1,nres-1 + write (iout,*) i,itype(i),itel(i) + enddo call read_bridge + nnt=1 + nct=nres + if (itype(1).eq.ntyp1) nnt=2 + if (itype(nres).eq.ntyp1) nct=nct-1 + write(iout,*) 'NNT=',NNT,' NCT=',NCT if (with_dihed_constr) then read (inp,*) ndih_constr + write (iout,*) "ndih_constr",ndih_constr if (ndih_constr.gt.0) then - read (inp,*) ftors - write (iout,*) 'FTORS',ftors - read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) + raw_psipred=.false. +C read (inp,*) ftors +C write (iout,*) 'FTORS',ftors + 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.' + & 'There are',ndih_constr,' restraints on gamma angles.' do i=1,ndih_constr - write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i) + 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 + else if (ndih_constr.lt.0) then + raw_psipred=.true. + call card_concat(controlcard,.true.) + 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 - nnt=1 - nct=nres - if (itype(1).eq.21) nnt=2 - if (itype(nres).eq.21) nct=nct-1 - write(iout,*) 'NNT=',NNT,' NCT=',NCT + endif + 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 call setup_var call init_int_table if (ns.gt.0) then 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 @@ -111,7 +206,25 @@ C Convert sequence to numeric code & dhpb(i),ebr,forcon(i) enddo endif + endif write (iout,'(a)') + write (iout,*) "setting ss mask dyn_ss",dyn_ss + 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. + write (iout,*) "i",i," iss",iss(i), + & " mask",dyn_ss_mask(iss(i)) + enddo + endif return end c----------------------------------------------------------------------------- @@ -148,12 +261,14 @@ C Read bridging residues. write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,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?!!!' stop endif