subroutine read_control C C Read molecular data C implicit none include 'DIMENSIONS' include 'sizesclu.dat' include 'COMMON.IOUNITS' include 'COMMON.TIME1' include 'COMMON.SBRIDGE' include 'COMMON.CONTROL' include 'COMMON.CLUSTER' include 'COMMON.CHAIN' include 'COMMON.HEADER' include 'COMMON.FFIELD' include 'COMMON.FREE' character*320 controlcard,ucase #ifdef MPL include 'COMMON.INFO' #endif integer i read (INP,'(a80)') titel call card_concat(controlcard) call readi(controlcard,'NRES',nres,0) call readi(controlcard,'RESCALE',rescale_mode,2) 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) dyn_ss=(index(controlcard,'DYN_SS').gt.0) iscode=index(controlcard,'ONE_LETTER') tree=(index(controlcard,'MAKE_TREE').gt.0) 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,0) call readi(controlcard,'NCLUST',nclust,5) call readi(controlcard,'NSTART',nstart,0) call readi(controlcard,'NEND',nend,0) call reada(controlcard,'ECUT',ecut,10.0d0) call reada(controlcard,'PROB',prob_limit,0.99d0) write (iout,*) "Probability limit",prob_limit lgrp=(index(controlcard,'LGRP').gt.0) caonly=(index(controlcard,'CA_ONLY').gt.0) print_dist=(index(controlcard,'PRINT_DIST').gt.0) if (ncut.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 call readi(controlcard,'NTEMP',nT,1) write (iout,*) "nT",nT call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0) write (iout,*) "nT",nT write (iout,*) 'beta_h',(beta_h(i),i=1,nT) do i=1,nT beta_h(i)=1.0d0/(1.987D-3*beta_h(i)) enddo write (iout,*) 'beta_h',(beta_h(i),i=1,nT) lprint_cart=index(controlcard,"PRINT_CART") .gt.0 lprint_int=index(controlcard,"PRINT_INT") .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 call flush(iout) if (min_var) iopt=1 return end c-------------------------------------------------------------------------- subroutine molread C C Read molecular data. C implicit none include 'DIMENSIONS' 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.HEADER' include 'COMMON.CONTROL' include 'COMMON.CONTACTS' include 'COMMON.TIME1' include 'COMMON.TORCNSTR' #ifdef MPL include 'COMMON.INFO' #endif character*4 sequence(maxres) character*800 weightcard integer rescode double precision x(maxvar) integer itype_pdb(maxres) logical seq_comp integer i,j C C Body C C Read weights of the subsequent energy terms. call card_concat(weightcard) call reada(weightcard,'WSC',wsc,1.0d0) call reada(weightcard,'WLONG',wsc,wsc) call reada(weightcard,'WSCP',wscp,1.0d0) call reada(weightcard,'WELEC',welec,1.0D0) call reada(weightcard,'WVDWPP',wvdwpp,welec) call reada(weightcard,'WEL_LOC',wel_loc,1.0D0) call reada(weightcard,'WCORR4',wcorr4,0.0D0) call reada(weightcard,'WCORR5',wcorr5,0.0D0) call reada(weightcard,'WCORR6',wcorr6,0.0D0) call reada(weightcard,'WTURN3',wturn3,1.0D0) call reada(weightcard,'WTURN4',wturn4,1.0D0) call reada(weightcard,'WTURN6',wturn6,1.0D0) call reada(weightcard,'WSTRAIN',wstrain,1.0D0) call reada(weightcard,'WBOND',wbond,1.0D0) call reada(weightcard,'WTOR',wtor,1.0D0) 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,'SCAL14',scal14,0.4D0) call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'WSCCOR',wsccor,1.0D0) 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) if (index(weightcard,'SOFT').gt.0) ipot=6 C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) 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,'(a,f10.2)') 'S-S depth: ',ss_depth 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 write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1) if (wcorr4.gt.0.0d0) wcorr=wcorr4 weights(1)=wsc weights(2)=wscp weights(3)=welec weights(4)=wcorr weights(5)=wcorr5 weights(6)=wcorr6 weights(7)=wel_loc weights(8)=wturn3 weights(9)=wturn4 weights(10)=wturn6 weights(11)=wang weights(12)=wscloc weights(13)=wtor weights(14)=wtor_d weights(15)=wstrain weights(16)=wvdwpp weights(17)=wbond weights(18)=scal14 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3, & wturn4,wturn6,wsccor 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ & 'WELEC= ',f10.6,' (p-p electr)'/ & 'WVDWPP= ',f10.6,' (p-p VDW)'/ & 'WBOND= ',f10.6,' (stretching)'/ & 'WANG= ',f10.6,' (bending)'/ & 'WSCLOC= ',f10.6,' (SC local)'/ & 'WTOR= ',f10.6,' (torsional)'/ & 'WTORD= ',f10.6,' (double torsional)'/ & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ & 'WCORR4= ',f10.6,' (multi-body 4th order)'/ & 'WCORR5= ',f10.6,' (multi-body 5th order)'/ & '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)'/ & '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' write (iout,'(2(a,f5.3/))') & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr, & 'Range of quenching the correlation terms:',2*delt_corr else if (wcorr.gt.0.0d0) then write (iout,'(/2a/)') 'Hydrogen-bonding correlation ', & 'between contact pairs of peptide groups' endif write (iout,'(a,f8.3)') & 'Scaling factor of 1,4 SC-p interactions:',scal14 write (iout,'(a,f8.3)') & 'General scaling factor of SC-p interactions:',scalscp r0_corr=cutoff_corr-delt_corr do i=1,20 aad(i,1)=scalscp*aad(i,1) aad(i,2)=scalscp*aad(i,2) bad(i,1)=scalscp*bad(i,1) bad(i,2)=scalscp*bad(i,2) enddo call flush(iout) print *,'indpdb=',indpdb,' pdbref=',pdbref C Read sequence if not taken from the pdb file. 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 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 print *,nres 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 #else if (itype(i).eq.21) then #endif itel(i)=0 #ifdef PROCOR else if (itype(i+1).ne.20) then #else else if (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 print *,'Call Read_Bridge.' call read_bridge if (with_dihed_constr) then read (inp,*) 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) write (iout,*) & 'There are',ndih_constr,' constraints on phi angles.' do i=1,ndih_constr write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i) enddo do i=1,ndih_constr phi0(i)=deg2rad*phi0(i) drange(i)=deg2rad*drange(i) enddo endif endif nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT if (itype(1).eq.21) nnt=2 if (itype(nres).eq.21) 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 c if (pdbref) then c read(inp,'(a)') pdbfile c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile c open(ipdbin,file=pdbfile,status='old',err=33) c goto 34 c 33 write (iout,'(a)') 'Error opening PDB file.' c stop c 34 continue c print *,'Begin reading pdb data' c call readpdb c print *,'Finished reading pdb data' c write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup c do i=1,nres c itype_pdb(i)=itype(i) c enddo c close (ipdbin) c write (iout,'(a,i3)') 'nsup=',nsup c nstart_seq=nnt c if (nsup.le.(nct-nnt+1)) then c do i=0,nct-nnt+1-nsup c if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then c nstart_seq=nnt+i c goto 111 c endif c enddo c write (iout,'(a)') c & 'Error - sequences to be superposed do not match.' c stop c else c do i=0,nsup-(nct-nnt+1) c if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) c & then c nstart_sup=nstart_sup+i c nsup=nct-nnt+1 c goto 111 c endif c enddo c write (iout,'(a)') c & 'Error - sequences to be superposed do not match.' c endif c 111 continue c write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup, c & ' nstart_seq=',nstart_seq c endif call init_int_table call setup_var write (iout,*) "molread: REFSTR",refstr if (refstr) then if (.not.pdbref) then call read_angles(inp,*38) goto 39 38 write (iout,'(a)') 'Error reading reference structure.' #ifdef MPL call mp_stopall(Error_Msg) #else stop 'Error reading reference structure' #endif 39 call chainbuild nstart_sup=nnt nstart_seq=nnt nsup=nct-nnt+1 do i=1,2*nres do j=1,3 cref(j,i)=c(j,i) enddo enddo endif call contact(print_contact_map,ncont_ref,icont_ref) endif c Read distance restraints if (constr_dist.gt.0) then call read_dist_constr call hpb_partition endif 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 none include 'DIMENSIONS' 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.HEADER' include 'COMMON.CONTROL' include 'COMMON.TIME1' #ifdef MPL include 'COMMON.INFO' #endif integer i,j C Read bridging residues. read (inp,*) ns,(iss(i),i=1,ns) print *,'ns=',ns 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(iss(i)),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') & 'Do you REALLY think that the residue ',restyp(iss(i)),i, & ' can form a disulfide bridge?!!!' #ifdef MPL call mp_stopall(error_msg) #else stop #endif endif enddo C Read preformed bridges. if (ns.gt.0) then read (inp,*) nss,(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.' #ifdef MPL call mp_stopall(error_msg) #else stop #endif 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 c dhpb(i)=dbr c forcon(i)=fbr enddo do i=1,nss ihpb(i)=ihpb(i)+nres jhpb(i)=jhpb(i)+nres enddo 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. c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU" enddo endif print *, "Leaving brigde read" return end c---------------------------------------------------------------------------- subroutine read_angles(kanal,*) implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' integer i,kanal 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 reada(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch double precision wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch) if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*) wartosc return end c---------------------------------------------------------------------------- subroutine multreada(rekord,lancuch,tablica,dim,default) implicit none integer dim,i double precision tablica(dim),default character*(*) rekord,lancuch integer ilen,iread external ilen do i=1,dim tablica(i)=default enddo iread=index(rekord,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 readi(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch integer wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch) if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 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' include 'COMMON.IOUNITS' character*(*) card character*80 karta,ucase external ilen read (inp,'(a)') karta karta=ucase(karta) card=' ' do while (karta(80:80).eq.'&') card=card(:ilen(card)+1)//karta(:79) read (inp,'(a)') karta karta=ucase(karta) enddo card=card(:ilen(card)+1)//karta return end c---------------------------------------------------------------------------- subroutine openunits implicit none include 'DIMENSIONS' #ifdef MPI include "mpif.h" character*3 liczba include "COMMON.MPI" #endif include 'COMMON.IOUNITS' include 'COMMON.CONTROL' integer lenpre,lenpot,ilen external ilen character*16 cformat,cprint character*16 ucase integer lenint,lenout call getenv('INPUT',prefix) call getenv('OUTPUT',prefout) call getenv('INTIN',prefintin) call getenv('COORD',cformat) call getenv('PRINTCOOR',cprint) call getenv('SCRATCHDIR',scratchdir) from_bx=.true. from_cx=.false. if (index(ucase(cformat),'CX').gt.0) then from_cx=.true. from_bx=.false. endif from_cart=.true. lenpre=ilen(prefix) lenout=ilen(prefout) lenint=ilen(prefintin) C Get the names and open the input files open (inp,file=prefix(:ilen(prefix))//'.inp',status='old') #ifdef MPI write (liczba,'(bz,i3.3)') me outname=prefout(:lenout)//'_clust.out_'//liczba #else outname=prefout(:lenout)//'_clust.out' #endif if (from_bx) then intinname=prefintin(:lenint)//'.bx' else if (from_cx) then intinname=prefintin(:lenint)//'.cx' else intinname=prefintin(:lenint)//'.int' endif rmsname=prefintin(:lenint)//'.rms' open (jplot,file=prefout(:ilen(prefout))//'.tex', & status='unknown') open (jrms,file=rmsname,status='unknown') open(iout,file=outname,status='unknown') C Get parameter filenames and open the parameter files. call getenv('BONDPAR',bondname) open (ibond,file=bondname,status='old') call getenv('THETPAR',thetname) open (ithep,file=thetname,status='old') call getenv('ROTPAR',rotname) open (irotam,file=rotname,status='old') call getenv('TORPAR',torname) open (itorp,file=torname,status='old') call getenv('TORDPAR',tordname) open (itordp,file=tordname,status='old') call getenv('FOURIER',fouriername) open (ifourier,file=fouriername,status='old') call getenv('ELEPAR',elename) open (ielep,file=elename,status='old') call getenv('SIDEPAR',sidename) open (isidep,file=sidename,status='old') call getenv('SIDEP',sidepname) open (isidep1,file=sidepname,status="old") call getenv('SCCORPAR',sccorname) open (isccor,file=sccorname,status="old") #ifndef OLDSCP C C 8/9/01 In the newest version SCp interaction constants are read from a file C Use -DOLDSCP to use hard-coded constants instead. C call getenv('SCPPAR',scpname) open (iscpp,file=scpname,status='old') #endif return end c------------------------------------------------------------------------------- subroutine read_dist_constr implicit real*8 (a-h,o-z) include 'DIMENSIONS' 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 c write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup call card_concat(controlcard) c call flush(iout) c write (iout,'(a)') controlcard 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) 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 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 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_ 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 enddo 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) return end