subroutine read_general_data(*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "DIMENSIONS.FREE" include "COMMON.TORSION" include "COMMON.INTERACT" include "COMMON.IOUNITS" include "COMMON.TIME1" include "COMMON.PROT" include "COMMON.PROTFILES" include "COMMON.CHAIN" include "COMMON.NAMES" include "COMMON.FFIELD" include "COMMON.ENEPS" include "COMMON.WEIGHTS" include "COMMON.FREE" include "COMMON.CONTROL" include "COMMON.ENERGIES" include "COMMON.SPLITELE" include "COMMON.SBRIDGE" include "COMMON.SHIELD" include "COMMON.SAXS" character*800 controlcard integer i,j,k,ii,n_ene_found integer ind,itype1,itype2,itypf,itypsc,itypp integer ilen external ilen character*16 ucase character*16 key external ucase double precision pi call card_concat(controlcard,.true.) call readi(controlcard,"N_ENE",n_ene,max_ene) if (n_ene.gt.max_ene) then write (iout,*) "Error: parameter out of range: N_ENE",n_ene, & max_ene return1 endif call readi(controlcard,"NPARMSET",nparmset,1) separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0 call readi(controlcard,"IPARMPRINT",iparmprint,1) write (iout,*) "PARMPRINT",iparmprint if (nparmset.gt.max_parm) then write (iout,*) "Error: parameter out of range: NPARMSET", & nparmset, Max_Parm return1 endif call readi(controlcard,"MAXIT",maxit,5000) call reada(controlcard,"FIMIN",fimin,1.0d-3) call readi(controlcard,"ENSEMBLES",ensembles,0) hamil_rep=index(controlcard,"HAMIL_REP").gt.0 write (iout,*) "Number of energy parameter sets",nparmset call multreadi(controlcard,"ISAMPL",isampl,nparmset,1) do i=1,nparmset if (isampl(i).eq.0) then write (iout,*) "ERROR: isampl is 0 for parmset",i call flush(iout) stop endif enddo write (iout,*) "MaxSlice",MaxSlice call readi(controlcard,"NSLICE",nslice,1) call flush(iout) if (nslice.gt.MaxSlice) then write (iout,*) "Error: parameter out of range: NSLICE",nslice, & MaxSlice return1 endif write (iout,*) "Frequency of storing conformations", & (isampl(i),i=1,nparmset) write (iout,*) "Maxit",maxit," Fimin",fimin call readi(controlcard,"NQ",nQ,1) if (nQ.gt.MaxQ) then write (iout,*) "Error: parameter out of range: NQ",nq, & maxq return1 endif indpdb=0 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) if (index(controlcard,"CLASSIFY").gt.0) indpdb=1 call reada(controlcard,"DELTA",delta,1.0d-2) call reada(controlcard,"TOLE",tole,1.0d-1) call readi(controlcard,"EINICHECK",einicheck,2) call reada(controlcard,"DELTRMS",deltrms,5.0d-2) call reada(controlcard,"DELTRGY",deltrgy,5.0d-2) call readi(controlcard,"RESCALE",rescale_mode,1) check_conf=index(controlcard,"NO_CHECK_CONF").eq.0 call readi(controlcard,'TORMODE',tor_mode,0) write(iout,*) "torsional and valence angle mode",tor_mode call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0) 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,25.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) write (iout,*) "Cutoff on interactions",r_cut write (iout,*) "lambda",rlamb call reada(controlcard,"LIPTHICK",lipthick,0.0d0) call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 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 call readi(controlcard,'SYM',symetr,1) write (iout,*) "DISTCHAINMAX",distchainmax write (iout,*) "delta",delta write (iout,*) "einicheck",einicheck write (iout,*) "rescale_mode",rescale_mode call flush(iout) bxfile=index(controlcard,"BXFILE").gt.0 cxfile=index(controlcard,"CXFILE").gt.0 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) & bxfile=.true. histfile=index(controlcard,"HISTFILE").gt.0 histout=index(controlcard,"HISTOUT").gt.0 entfile=index(controlcard,"ENTFILE").gt.0 zscfile=index(controlcard,"ZSCFILE").gt.0 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0 write (iout,*) "with_dihed_constr ",with_dihed_constr with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 write (iout,*) "with_theta_constr ",with_theta_constr call readi(controlcard,'SHIELD',shield_mode,0) write(iout,*) "shield_mode",shield_mode C endif call readi(controlcard,'TORMODE',tor_mode,0) write(iout,*) "torsional and valence angle mode",tor_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 C do i=1,ntyp C long_r_sidechain(i)=vbldsc0(1,i) C short_r_sidechain(i)=sigma0(i) C enddo buff_shield=1.0d0 endif call readi(controlcard,'CONSTR_DIST',constr_dist,0) call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) c if (constr_homology) tole=dmax1(tole,1.5d0) write (iout,*) "with_homology_constr ",with_dihed_constr, & " CONSTR_HOMOLOGY",constr_homology read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0 write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD write (iout,*) "out_template_restr",OUT_TEMPLATE_RESTR dyn_ss=index(controlcard,"DYN_SS").gt.0 adaptive = index(controlcard,"ADAPTIVE").gt.0 call readi(controlcard,'NSAXS',nsaxs,0) call readi(controlcard,'SAXS_MODE',saxs_mode,0) call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0) call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0) write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE", & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff return end c------------------------------------------------------------------------------ subroutine read_efree(*) C C Read molecular data C implicit none include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'DIMENSIONS.COMPAR' include 'DIMENSIONS.FREE' include 'COMMON.IOUNITS' include 'COMMON.TIME1' include 'COMMON.SBRIDGE' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.HEADER' include 'COMMON.GEO' include 'COMMON.FREE' character*320 controlcard,ucase integer iparm,ib,i,j,npars integer ilen external ilen if (hamil_rep) then npars=1 else npars=nParmSet endif do iparm=1,npars call card_concat(controlcard,.true.) call readi(controlcard,'NT',nT_h(iparm),1) write (iout,*) "iparm",iparm," nt",nT_h(iparm) call flush(iout) if (nT_h(iparm).gt.MaxT_h) then write (iout,*) "Error: parameter out of range: NT",nT_h(iparm), & MaxT_h return1 endif replica(iparm)=index(controlcard,"REPLICA").gt.0 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ", & replica(iparm)," umbrella ",umbrella(iparm), & " read_iset",read_iset(iparm) call flush(iout) do ib=1,nT_h(iparm) call card_concat(controlcard,.true.) call readi(controlcard,'NR',nR(ib,iparm),1) if (umbrella(iparm)) then nRR(ib,iparm)=1 else nRR(ib,iparm)=nR(ib,iparm) endif if (nR(ib,iparm).gt.MaxR) then write (iout,*) "Error: parameter out of range: NR", & nR(ib,iparm),MaxR return1 endif call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0) beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3) call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm), & 0.0d0) do i=1,nR(ib,iparm) call card_concat(controlcard,.true.) call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ, & 100.0d0) call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ, & 0.0d0) enddo enddo do ib=1,nT_h(iparm) write (iout,*) "ib",ib," beta_h", & 1.0d0/(0.001987*beta_h(ib,iparm)) write (iout,*) "nR",nR(ib,iparm) write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm)) do i=1,nR(ib,iparm) write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ), & "q0",(q0(j,i,ib,iparm),j=1,nQ) enddo call flush(iout) enddo enddo if (hamil_rep) then do iparm=2,nParmSet nT_h(iparm)=nT_h(1) do ib=1,nT_h(iparm) nR(ib,iparm)=nR(ib,1) if (umbrella(iparm)) then nRR(ib,iparm)=1 else nRR(ib,iparm)=nR(ib,1) endif beta_h(ib,iparm)=beta_h(ib,1) do i=1,nR(ib,iparm) f(i,ib,iparm)=f(i,ib,1) do j=1,nQ KH(j,i,ib,iparm)=KH(j,i,ib,1) Q0(j,i,ib,iparm)=Q0(j,i,ib,1) enddo enddo replica(iparm)=replica(1) umbrella(iparm)=umbrella(1) read_iset(iparm)=read_iset(1) enddo enddo endif return end c----------------------------------------------------------------------------- subroutine read_protein_data(*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "DIMENSIONS.FREE" #ifdef MPI include "mpif.h" integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE) include "COMMON.MPI" #endif include "COMMON.CHAIN" include "COMMON.IOUNITS" include "COMMON.PROT" include "COMMON.PROTFILES" include "COMMON.NAMES" include "COMMON.FREE" include "COMMON.OBCINKA" character*64 nazwa character*16000 controlcard integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars external ilen,iroof if (hamil_rep) then npars=1 else npars=nparmset endif do iparm=1,npars C Read names of files with conformation data. if (replica(iparm)) then nthr = 1 else nthr = nT_h(iparm) endif do ib=1,nthr do ii=1,nRR(ib,iparm) write (iout,*) "Parameter set",iparm," temperature",ib, & " window",ii call flush(iout) call card_concat(controlcard,.true.) write (iout,*) controlcard(:ilen(controlcard)) call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0) call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0) call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0) call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1) call readi(controlcard,"REC_END",rec_end(ii,ib,iparm), & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1) call reada(controlcard,"TIME_START", & time_start_collect(ii,ib,iparm),0.0d0) call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm), & 1.0d10) write (iout,*) "rec_start",rec_start(ii,ib,iparm), & " rec_end",rec_end(ii,ib,iparm) write (iout,*) "time_start",time_start_collect(ii,ib,iparm), & " time_end",time_end_collect(ii,ib,iparm) call flush(iout) if (replica(iparm)) then call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1) write (iout,*) "Number of trajectories",totraj(ii,iparm) call flush(iout) endif if (nfile_bin(ii,ib,iparm).lt.2 & .and. nfile_asc(ii,ib,iparm).eq.0 & .and. nfile_cx(ii,ib,iparm).eq.0) then write (iout,*) "Error - no action specified!" return1 endif if (nfile_bin(ii,ib,iparm).gt.0) then call card_concat(controlcard,.false.) call split_string(controlcard,protfiles(1,1,ii,ib,iparm), & maxfile_prot,nfile_bin(ii,ib,iparm)) #ifdef DEBUG write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm) write(iout,*) (protfiles(i,1,ii,ib,iparm), & i=1,nfile_bin(ii,ib,iparm)) #endif endif if (nfile_asc(ii,ib,iparm).gt.0) then call card_concat(controlcard,.false.) call split_string(controlcard,protfiles(1,2,ii,ib,iparm), & maxfile_prot,nfile_asc(ii,ib,iparm)) #ifdef DEBUG write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm) write(iout,*) (protfiles(i,2,ii,ib,iparm), & i=1,nfile_asc(ii,ib,iparm)) #endif else if (nfile_cx(ii,ib,iparm).gt.0) then call card_concat(controlcard,.false.) call split_string(controlcard,protfiles(1,2,ii,ib,iparm), & maxfile_prot,nfile_cx(ii,ib,iparm)) #ifdef DEBUG write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm) write(iout,*) (protfiles(i,2,ii,ib,iparm), & i=1,nfile_cx(ii,ib,iparm)) #endif endif call flush(iout) enddo enddo enddo return end c------------------------------------------------------------------------------- subroutine opentmp(islice,iunit,bprotfile_temp) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "DIMENSIONS.FREE" #ifdef MPI include "mpif.h" integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE) include "COMMON.MPI" #endif include "COMMON.IOUNITS" include "COMMON.PROTFILES" include "COMMON.PROT" include "COMMON.FREE" character*256 bprotfile_temp character*3 liczba,liczba2 character*2 liczba1 integer iunit,islice integer ilen,iroof external ilen,iroof logical lerr write (liczba1,'(bz,i2.2)') islice write (liczba,'(bz,i3.3)') me #ifdef MPI c write (iout,*) "separate_parset ",separate_parset, c & " myparm",myparm if (separate_parset) then write (liczba2,'(bz,i3.3)') myparm bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"// & prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1 open (iunit,file=bprotfile_temp,status="unknown", & form="unformatted",access="direct",recl=lenrec) else bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"// & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1 open (iunit,file=bprotfile_temp,status="unknown", & form="unformatted",access="direct",recl=lenrec) endif #else bprotfile_temp = scratchdir(:ilen(scratchdir))// & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1 open (iunit,file=bprotfile_temp,status="unknown", & form="unformatted",access="direct",recl=lenrec) #endif c write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp", c & bprotfile_temp c call flush(iout) return end c------------------------------------------------------------------------------- subroutine read_database(*) implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "DIMENSIONS.FREE" #ifdef MPI include "mpif.h" integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE) include "COMMON.MPI" #endif include "COMMON.CHAIN" include "COMMON.IOUNITS" include "COMMON.PROTFILES" include "COMMON.NAMES" include "COMMON.VAR" include "COMMON.GEO" include "COMMON.ENEPS" include "COMMON.PROT" include "COMMON.INTERACT" include "COMMON.FREE" include "COMMON.SBRIDGE" include "COMMON.OBCINKA" real*4 csingle(3,maxres2) character*256 nazwa,bprotfile_temp character*3 liczba character*2 liczba1 integer i,j,ii,jj(maxslice),k,kk(maxslice),l, & ll(maxslice),mm(maxslice),if integer nrec,nlines,iscor,iunit,islice double precision energ integer ilen,iroof external ilen,iroof double precision rmsdev,energia(0:max_ene),efree,eini,temp double precision prop(maxQ) integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff integer iparm,ib,iib,ir,nprop,nthr,npars double precision etot,time integer ixdrf,iret logical lerr,linit lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24 lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ lenrec=lenrec2+8 write (iout,*) "lenrec",lenrec," lenrec1",lenrec1, & " lenrec2",lenrec2 do i=1,nQ prop(i)=0.0d0 enddo do islice=1,nslice ll(islice)=0 mm(islice)=0 enddo write (iout,*) "nparmset",nparmset if (hamil_rep) then npars=1 else npars=nparmset endif do iparm=1,npars if (replica(iparm)) then nthr = 1 else nthr = nT_h(iparm) endif do ib=1,nthr do iR=1,nRR(ib,iparm) write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ do islice=1,nslice jj(islice)=0 kk(islice)=0 enddo IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN c Read conformations from binary DA files (one per batch) and write them to c a binary DA scratchfile. write (liczba,'(bz,i3.3)') me do if=1,nfile_bin(iR,ib,iparm) nazwa=protfiles(if,1,iR,ib,iparm) & (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx" open (ientin,file=nazwa,status="old",form="unformatted", & access="direct",recl=lenrec2,err=1111) ii=0 do islice=1,nslice call opentmp(islice,ientout,bprotfile_temp) call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice), & mm(islice),iR,ib,iparm) close(ientout) enddo close(ientin) enddo ENDIF ! NFILE_BIN>0 c IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN c Read conformations from multiple ASCII int files and write them to a binary c DA scratchfile. do if=1,nfile_asc(iR,ib,iparm) nazwa=protfiles(if,2,iR,ib,iparm) & (:ilen(protfiles(if,2,iR,ib,iparm)))//".x" open(unit=ientin,file=nazwa,status='old',err=1111) write(iout,*) "reading ",nazwa(:ilen(nazwa)) ii=0 call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm) enddo ! if ENDIF IF (NFILE_CX(iR,ib,iparm).gt.0) THEN c Read conformations from cx files and write them to a binary c DA scratchfile. do if=1,nfile_cx(iR,ib,iparm) nazwa=protfiles(if,2,iR,ib,iparm) & (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx" write(iout,*) "reading ",nazwa(:ilen(nazwa)) ii=0 print *,"Calling cxread" call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm, & *1111) close(ientout) write (iout,*) "exit cxread" call flush(iout) enddo ENDIF do islice=1,nslice stot(islice)=stot(islice)+jj(islice) enddo enddo enddo write (iout,*) "IPARM",iparm enddo if (nslice.eq.1) then #ifdef MPI write (liczba,'(bz,i3.3)') me bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"// & prefix(:ilen(prefix))//liczba//".xbin.tmp" #else bprotfile_temp = scratchdir(:ilen(scratchdir))// & "/"//prefix(:ilen(prefix))//".xbin.tmp" #endif write(iout,*) mm(1)," conformations read",ll(1), & " conformations written to ", & bprotfile_temp(:ilen(bprotfile_temp)) else do islice=1,nslice write (liczba1,'(bz,i2.2)') islice #ifdef MPI write (liczba,'(bz,i3.3)') me bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"// & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1 #else bprotfile_temp = scratchdir(:ilen(scratchdir))// & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1 #endif write(iout,*) mm(islice)," conformations read",ll(islice), & " conformations written to ", & bprotfile_temp(:ilen(bprotfile_temp)) enddo endif #ifdef MPI c Check if everyone has the same number of conformations c call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL, c & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR) maxslice_buff=maxslice call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER, & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR) lerr=.false. do i=0,nprocs-1 if (i.ne.me) then do islice=1,nslice if (stot(islice).ne.ntot_all(islice,i)) then write (iout,*) "Number of conformations at processor",i, & " differs from that at processor",me, & stot(islice),ntot_all(islice,i)," slice",islice lerr = .true. endif enddo endif enddo if (lerr) then write (iout,*) write (iout,*) "Numbers of conformations read by processors" write (iout,*) do i=0,nprocs-1 write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice) enddo write (iout,*) "Calculation terminated." call flush(iout) return1 endif do islice=1,nslice ntot(islice)=stot(islice) enddo return #endif 1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa)) call flush(iout) return1 end c------------------------------------------------------------------------------ subroutine card_concat(card,to_upper) implicit none include 'DIMENSIONS.ZSCOPT' include "COMMON.IOUNITS" character*(*) card character*80 karta,ucase logical to_upper integer ilen external ilen read (inp,'(a)') karta if (to_upper) karta=ucase(karta) card=' ' do while (karta(80:80).eq.'&') card=card(:ilen(card)+1)//karta(:79) read (inp,'(a)') karta if (to_upper) karta=ucase(karta) enddo card=card(:ilen(card)+1)//karta 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(:ilen(lancuch))//"=") if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*) wartosc return end c---------------------------------------------------------------------------- subroutine reada(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch character*80 aux double precision wartosc,default integer ilen,iread external ilen iread=index(rekord,lancuch(:ilen(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 multreada(rekord,lancuch,tablica,dim,default) implicit none integer dim,i double precision 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 reads(rekord,lancuch,wartosc,default) implicit none character*(*) rekord,lancuch,wartosc,default character*80 aux integer ilen,lenlan,lenrec,iread,ireade external ilen logical iblnk external iblnk lenlan=ilen(lancuch) lenrec=ilen(rekord) iread=index(rekord,lancuch(:lenlan)//"=") c print *,"rekord",rekord," lancuch",lancuch c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec if (iread.eq.0) then wartosc=default return endif iread=iread+lenlan+1 c print *,"iread",iread c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread)) do while (iread.le.lenrec .and. iblnk(rekord(iread:iread))) iread=iread+1 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread)) enddo c print *,"iread",iread if (iread.gt.lenrec) then wartosc=default return endif ireade=iread+1 c print *,"ireade",ireade do while (ireade.lt.lenrec .and. & .not.iblnk(rekord(ireade:ireade))) ireade=ireade+1 enddo wartosc=rekord(iread:ireade) return end c---------------------------------------------------------------------------- subroutine multreads(rekord,lancuch,tablica,dim,default) implicit none integer dim,i character*(*) rekord,lancuch,tablica(dim),default character*80 aux integer ilen,lenlan,lenrec,iread,ireade external ilen logical iblnk external iblnk do i=1,dim tablica(i)=default enddo lenlan=ilen(lancuch) lenrec=ilen(rekord) iread=index(rekord,lancuch(:lenlan)//"=") c print *,"rekord",rekord," lancuch",lancuch c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec if (iread.eq.0) return iread=iread+lenlan+1 do i=1,dim c print *,"iread",iread c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread)) do while (iread.le.lenrec .and. iblnk(rekord(iread:iread))) iread=iread+1 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread)) enddo c print *,"iread",iread if (iread.gt.lenrec) return ireade=iread+1 c print *,"ireade",ireade do while (ireade.lt.lenrec .and. & .not.iblnk(rekord(ireade:ireade))) ireade=ireade+1 enddo tablica(i)=rekord(iread:ireade) iread=ireade+1 enddo end c---------------------------------------------------------------------------- subroutine split_string(rekord,tablica,dim,nsub) implicit none integer dim,nsub,i,ii,ll,kk character*(*) tablica(dim) character*(*) rekord integer ilen external ilen do i=1,dim tablica(i)=" " enddo ii=1 ll = ilen(rekord) nsub=0 do i=1,dim C Find the start of term name kk = 0 do while (ii.le.ll .and. rekord(ii:ii).eq." ") ii = ii+1 enddo C Parse the name into TABLICA(i) until blank found do while (ii.le.ll .and. rekord(ii:ii).ne." ") kk = kk+1 tablica(i)(kk:kk)=rekord(ii:ii) ii = ii+1 enddo if (kk.gt.0) nsub=nsub+1 if (ii.gt.ll) return enddo return end c-------------------------------------------------------------------------------- integer function iroof(n,m) ii = n/m if (ii*m .lt. n) ii=ii+1 iroof = ii 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' include 'COMMON.INTERACT' integer ifrag_(2,100),ipair_(2,100) double precision wfrag_(100),wpair_(100) character*500 controlcard logical normalize,next integer restr_type double precision xlink(4,0:4) / c a b c sigma & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) restr_on_coord=.false. next=.true. DO WHILE (next) call card_concat(controlcard,.true.) next = index(controlcard,"NEXT").gt.0 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist) write (iout,*) "restr_type",restr_type call readi(controlcard,"NFRAG",nfrag_,0) 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 reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0) if (restr_type.eq.10) & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0) if (restr_type.eq.12) & call reada(controlcard,'SCAL_PEAK',scal_peak,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,*) "WBOLTZD",wboltzd write (iout,*) "SCAL_PEAK",scal_peak 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 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5) then nres0=nres read(inp,'(a)') pdbfile write (iout,*) & "Distance restraints will be constructed from structure ",pdbfile open(ipdbin,file=pdbfile,status='old',err=11) call readpdb(.true.) nres=nres0 close(ipdbin) endif 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) c call flush(iout) if (wfrag_(i).eq.0.0d0) cycle do j=ifrag_(1,i),ifrag_(2,i)-1 do k=j+1,ifrag_(2,i) c write (iout,*) "j",j," k",k ddjk=dist(j,k) if (restr_type.eq.1) then nhpb=nhpb+1 irestr_type(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 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(i) endif else if (restr_type.eq.3) then nhpb=nhpb+1 irestr_type(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.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) enddo enddo enddo do i=1,npair_ if (wpair_(i).eq.0.0d0) cycle 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) ddjk=dist(j,k) if (restr_type.eq.1) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wpair_(i) else if (constr_dist.eq.2) then if (ddjk.le.dist_cut) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wpair_(i) endif else if (restr_type.eq.3) then nhpb=nhpb+1 irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) endif write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) enddo enddo enddo c print *,ndist_ write (iout,*) "Distance restraints as read from input" do i=1,ndist_ if (restr_type.eq.12) then read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1), & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1), & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1), & fordepth_peak(nhpb_peak+1),npeak c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1), c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1), c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1), c & fordepth_peak(nhpb_peak+1),npeak if (forcon_peak(nhpb_peak+1).le.0.0d0.or. & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle nhpb_peak=nhpb_peak+1 irestr_type_peak(nhpb_peak)=12 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i ipeak(2,npeak)=i write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak), & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak), & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak), & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak) if (ibecarb_peak(nhpb_peak).eq.3) then jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres else if (ibecarb_peak(nhpb_peak).eq.2) then ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres else if (ibecarb_peak(nhpb_peak).eq.1) then ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres endif else if (restr_type.eq.11) then read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1), & dhpb1(nhpb+1),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 irestr_type(nhpb)=11 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb) c if (ibecarb(nhpb).gt.0) then c ihpb(nhpb)=ihpb(nhpb)+nres c jhpb(nhpb)=jhpb(nhpb)+nres c endif if (ibecarb(nhpb).eq.3) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.2) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.1) then ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+nres endif else if (restr_type.eq.10) then c Cross-lonk Markov-like potential call card_concat(controlcard,.true.) call readi(controlcard,"ILINK",ihpb(nhpb+1),0) call readi(controlcard,"JLINK",jhpb(nhpb+1),0) ibecarb(nhpb+1)=0 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle if (index(controlcard,"ZL").gt.0) then link_type=1 else if (index(controlcard,"ADH").gt.0) then link_type=2 else if (index(controlcard,"PDH").gt.0) then link_type=3 else if (index(controlcard,"DSS").gt.0) then link_type=4 else link_type=0 endif call reada(controlcard,"AXLINK",dhpb(nhpb+1), & xlink(1,link_type)) call reada(controlcard,"BXLINK",dhpb1(nhpb+1), & xlink(2,link_type)) call reada(controlcard,"CXLINK",fordepth(nhpb+1), & xlink(3,link_type)) call reada(controlcard,"SIGMA",forcon(nhpb+1), & xlink(4,link_type)) call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0) c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1), c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1) if (forcon(nhpb+1).le.0.0d0 .or. & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle nhpb=nhpb+1 irestr_type(nhpb)=10 c if (ibecarb(nhpb).gt.0) then c ihpb(nhpb)=ihpb(nhpb)+nres c jhpb(nhpb)=jhpb(nhpb)+nres c endif if (ibecarb(nhpb).eq.3) then jhpb(nhpb)=jhpb(nhpb)+nres else if (ibecarb(nhpb).eq.2) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.1) then ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+nres endif write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) else C print *,"in else" read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1), & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1) if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 if (dhpb1(nhpb).eq.0.0d0) then irestr_type(nhpb)=1 else irestr_type(nhpb)=2 endif c if (ibecarb(nhpb).gt.0) then c ihpb(nhpb)=ihpb(nhpb)+nres c jhpb(nhpb)=jhpb(nhpb)+nres c endif if (ibecarb(nhpb).eq.3) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.2) then ihpb(nhpb)=ihpb(nhpb)+nres else if (ibecarb(nhpb).eq.1) then ihpb(nhpb)=ihpb(nhpb)+nres jhpb(nhpb)=jhpb(nhpb)+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 (restr_type.eq.4) then write (iout,*) "The BFAC array" do i=nnt,nct write (iout,'(i5,f10.5)') i,bfac(i) enddo do i=nnt,nct if (itype(i).eq.ntyp1) cycle do j=nnt,i-1 if (itype(j).eq.ntyp1) cycle if (itype(i).eq.10) then iiend=0 else iiend=1 endif if (itype(j).eq.10) then jjend=0 else jjend=1 endif kk=0 do ii=0,iiend do jj=0,jjend nhpb=nhpb+1 irestr_type(nhpb)=1 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2) irestr_type(nhpb)=1 ibecarb(nhpb)=kk if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb) ihpb(nhpb)=i+nres*ii jhpb(nhpb)=j+nres*jj dhpb(nhpb)=dist(i+nres*ii,j+nres*jj) write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) kk=kk+1 enddo enddo enddo enddo endif if (restr_type.eq.5) then restr_on_coord=.true. do i=nnt,nct if (itype(i).eq.ntyp1) cycle bfac(i)=(scal_bfac/bfac(i))**2 enddo endif ENDDO ! next fordepthmax=0.0d0 if (normalize) then do i=nss+1,nhpb if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax) & fordepthmax=fordepth(i) enddo do i=nss+1,nhpb if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax enddo endif if (nhpb.gt.nss) then write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)') & "The following",nhpb-nss, & " distance restraints have been imposed:", & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V", & " score"," type" do i=nss+1,nhpb write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i), & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i), & irestr_type(i) enddo endif write (iout,*) call hpb_partition call flush(iout) return 11 write (iout,*)"read_dist_restr: error reading reference structure" stop end