subroutine molread(*) C C Read molecular data. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.TORCNSTR' include 'COMMON.CONTROL' include 'COMMON.SAXS' character*4 sequence(maxres) integer rescode,tperm double precision x(maxvar) 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 readi(controlcard,"NRES",nres,0) iscode=index(controlcard,"ONE_LETTER") if (nres.le.0) then write (iout,*) "Error: no residues in molecule" return1 endif if (nres.gt.maxres) then write (iout,*) "Error: too many residues",nres,maxres endif write(iout,*) 'nres=',nres C Read sequence of the protein if (iscode.gt.0) then read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres) else read (inp,'(20(1x,a3))') (sequence(i),i=1,nres) endif C Convert sequence to numeric code do i=1,nres itype(i)=rescode(i,sequence(i),iscode) enddo write (iout,*) "Numeric code:" write (iout,'(20i4)') (itype(i),i=1,nres) do i=1,nres-1 #ifdef PROCOR if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then #else if (itype(i).eq.ntyp1) then #endif itel(i)=0 #ifdef PROCOR else if (iabs(itype(i+1)).ne.20) then #else 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 call seq2chains(nres,itype,nchain,chain_length,chain_border, & ireschain) chain_border1(1,1)=1 chain_border1(2,1)=chain_border(2,1)+1 do i=2,nchain-1 chain_border1(1,i)=chain_border(1,i)-1 chain_border1(2,i)=chain_border(2,i)+1 enddo 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 write(iout,*)"chain",i,chain_length(i),chain_border(1,i), & chain_border(2,i),chain_border1(1,i),chain_border1(2,i) enddo call chain_symmetry(nchain,nres,itype,chain_border, & chain_length,npermchain,tabpermchain) write(iout,*) "ireschain permutations" do i=1,nres write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain), & ii=1,npermchain) enddo write(iout,*) "residue permutations" do i=1,nres write(iout,*) i,(iperm(i,ii),ii=1,npermchain) enddo 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 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,' restraints on gamma 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,.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 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 if (constr_homology.gt.0) then c write (iout,*) "About to call read_constr_homology" c call flush(iout) call read_constr_homology c write (iout,*) "Exit read_constr_homology" c call flush(iout) 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 #ifdef DEBUG call int_from_cart1(.false.) call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(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 #endif else homol_nset=0 endif call setup_var write (iout,*) "Calling init_int_table" 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 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 endif endif write (iout,'(a)') 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 write (iout,*) "calling read_saxs_consrtr",nsaxs if (nsaxs.gt.0) call read_saxs_constr return end c----------------------------------------------------------------------------- logical function seq_comp(itypea,itypeb,length) implicit none integer length,itypea(length),itypeb(length) integer i do i=1,length if (itypea(i).ne.itypeb(i)) then seq_comp=.false. return endif enddo seq_comp=.true. return end c----------------------------------------------------------------------------- subroutine read_bridge C Read information about disulfide bridges. implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' C Read bridging residues. read (inp,*) ns,(iss(i),i=1,ns) print *,'ns=',ns write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,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 Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then write (iout,'(2a,i3,a)') & '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(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' stop endif enddo 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 (nss.gt.0) then nhpb=nss C Check if the residues involved in bridges are in the specified list of C bridging residues. do i=1,nss do j=1,i-1 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then write (iout,'(a,i3,a)') 'Disulfide pair',i, & ' contains residues present in other pairs.' write (*,'(a,i3,a)') 'Disulfide pair',i, & ' contains residues present in other pairs.' stop endif enddo do j=1,ns if (ihpb(i).eq.iss(j)) goto 10 enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 10 continue do j=1,ns if (jhpb(i).eq.iss(j)) goto 20 enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 20 continue dhpb(i)=dbr forcon(i)=fbr enddo do i=1,nss ihpb(i)=ihpb(i)+nres jhpb(i)=jhpb(i)+nres enddo endif endif return end c------------------------------------------------------------------------------ subroutine read_angles(kanal,iscor,energ,iprot,*) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'COMMON.INTERACT' include 'COMMON.SBRIDGE' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' character*80 lineh read(kanal,'(a80)',end=10,err=10) lineh read(lineh(:5),*,err=8) ic read(lineh(6:),*,err=8) energ goto 9 8 ic=1 print *,'error, assuming e=1d10',lineh energ=1d10 nss=0 9 continue read(lineh(18:),*,end=10,err=10) nss IF (NSS.LT.9) THEN read (lineh(20:),*,end=10,err=10) & (IHPB(I),JHPB(I),I=1,NSS),iscor ELSE read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8) read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I), & I=9,NSS),iscor ENDIF c print *,"energy",energ," iscor",iscor read (kanal,*,err=10,end=10) (theta(i),i=3,nres) read (kanal,*,err=10,end=10) (phi(i),i=4,nres) read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1) read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1) do i=1,nres theta(i)=deg2rad*theta(i) phi(i)=deg2rad*phi(i) alph(i)=deg2rad*alph(i) omeg(i)=deg2rad*omeg(i) enddo return 10 return1 end c------------------------------------------------------------------------------- subroutine read_saxs_constr implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'DIMENSIONS.FREE' #ifdef MPI include 'mpif.h' #endif include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.SBRIDGE' include 'COMMON.SAXS' double precision cm(3) c read(inp,*) nsaxs write (iout,*) "Calling read_saxs nsaxs",nsaxs call flush(iout) if (saxs_mode.eq.0) then c SAXS distance distribution do i=1,nsaxs read(inp,*) distsaxs(i),Psaxs(i) enddo Cnorm = 0.0d0 do i=1,nsaxs Cnorm = Cnorm + Psaxs(i) enddo write (iout,*) "Cnorm",Cnorm do i=1,nsaxs Psaxs(i)=Psaxs(i)/Cnorm enddo write (iout,*) "Normalized distance distribution from SAXS" do i=1,nsaxs write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i) enddo Wsaxs0=0.0d0 do i=1,nsaxs Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i)) enddo write (iout,*) "Wsaxs0",Wsaxs0 else c SAXS "spheres". do i=1,nsaxs read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3) enddo do j=1,3 cm(j)=0.0d0 enddo do i=1,nsaxs do j=1,3 cm(j)=cm(j)+Csaxs(j,i) enddo enddo do j=1,3 cm(j)=cm(j)/nsaxs enddo do i=1,nsaxs do j=1,3 Csaxs(j,i)=Csaxs(j,i)-cm(j) enddo enddo write (iout,*) "SAXS sphere coordinates" do i=1,nsaxs write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3) enddo endif return end