module io_base !----------------------------------------------------------------------- use names use io_units implicit none !----------------------------------------------------------------------------- ! Max. number of AA residues integer,parameter :: maxres=101000!1200 ! Appr. max. number of interaction sites integer,parameter :: maxres2=2*maxres ! parameter (maxres6=6*maxres) ! parameter (mmaxres2=(maxres2*(maxres2+1)/2)) !----------------------------------------------------------------------------- ! Max. number of S-S bridges ! integer,parameter :: maxss=20 !----------------------------------------------------------------------------- ! Max. number of derivatives of virtual-bond and side-chain vectors in theta ! or phi. ! integer,parameter :: maxdim=(maxres-1)*(maxres-2)/2 !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains !----------------------------------------------------------------------------- ! readrtns_CSA.F !----------------------------------------------------------------------------- #ifndef XDRFPDB subroutine read_bridge ! Read information about disulfide bridges. use geometry_data, only: nres use energy_data use control_data, only:out1file use MPI_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif ! 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.DBASE' ! include 'COMMON.THREAD' ! include 'COMMON.TIME1' ! include 'COMMON.SETUP' !el local variables integer :: i,j,ierror print *,"ENTER READ" ! Read bridging residues. read (inp,*) ns write(iout,*) "ns",ns call flush(iout) if (ns.gt.0) then allocate(iss(ns)) read (inp,*) (iss(i),i=1,ns) endif ! print *,'ns=',ns if(me.eq.king.or..not.out1file) & write (iout,*) 'ns=',ns if (ns.gt.0) & write(iout,*) ' iss:',(iss(i),i=1,ns) ! Check whether the specified bridging residues are cystines. do i=1,ns write(iout,*) i,iss(i) if (itype(iss(i),1).ne.1) then if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') & 'Do you REALLY think that the residue ',& restyp(itype(iss(i),1),1),i,& ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') & 'Do you REALLY think that the residue ',& restyp(itype(iss(i),1),1),i,& ' can form a disulfide bridge?!!!' #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,ierror) stop #endif endif enddo if (dyn_ss) then if(.not.allocated(ihpb)) allocate(ihpb(ns)) if(.not.allocated(jhpb)) allocate(jhpb(ns)) endif ! Read preformed bridges. if (ns.gt.0) then read (inp,*) nss if (nss.gt.0) then if(.not.allocated(ihpb)) allocate(ihpb(nss)) if(.not.allocated(jhpb)) allocate(jhpb(nss)) read (inp,*) (ihpb(i),jhpb(i),i=1,nss) if(.not.allocated(dhpb)) allocate(dhpb(nss)) if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2 if(.not.allocated(dhpb1)) allocate(dhpb1(nss)) if(.not.allocated(ibecarb)) allocate(ibecarb(nss)) ! el Initialize the bridge array do i=1,nss dhpb(i)=0.0D0 enddo endif !-------------------- if(fg_rank.eq.0) & write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) if (nss.gt.0) then write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) nhpb=nss ! Check if the residues involved in bridges are in the specified list of ! 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 MPI call MPI_Finalize(MPI_COMM_WORLD,ierror) 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 dhpb(i)=d0cm forcon(i)=akcm enddo do i=1,nss ihpb(i)=ihpb(i)+nres jhpb(i)=jhpb(i)+nres enddo endif endif write(iout,*) "end read_bridge" return end subroutine read_bridge !----------------------------------------------------------------------------- subroutine read_x(kanal,*) use geometry_data use energy_data use geometry, only:int_from_cart1 ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.CONTROL' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! Read coordinates from input ! !el local variables integer :: l,k,j,i,kanal read(kanal,'(8f10.5)',end=10,err=10) & ((c(l,k),l=1,3),k=1,nres),& ((c(l,k+nres),l=1,3),k=nnt,nct) do j=1,3 c(j,nres+1)=c(j,1) c(j,2*nres)=c(j,nres) enddo call int_from_cart1(.false.) 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=nnt,nct if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then 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 endif enddo return 10 return 1 end subroutine read_x !----------------------------------------------------------------------------- subroutine read_threadbase use geometry_data use energy_data use compare_data ! implicit real*8 (a-h,o-z) ! 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.DBASE' ! include 'COMMON.THREAD' ! include 'COMMON.TIME1' !el local variables integer :: k,j,i ! Read pattern database for threading. read (icbase,*) nseq allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq) allocate(nres_base(3,nseq)) !(3,maxseq) allocate(str_nam(nseq)) !(maxseq) do i=1,nseq read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),& nres_base(2,i),nres_base(3,i) read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,& nres_base(1,i)) ! write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i), ! & nres_base(2,i),nres_base(3,i) ! write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1, ! & nres_base(1,i)) enddo close (icbase) if (weidis.eq.0.0D0) weidis=0.1D0 do i=nnt,nct do j=i+2,nct nhpb=nhpb+1 ihpb(nhpb)=i jhpb(nhpb)=j forcon(nhpb)=weidis enddo enddo read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl) write (iout,'(a,i5)') 'nexcl: ',nexcl write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl) return end subroutine read_threadbase !----------------------------------------------------------------------------- #ifdef WHAM_RUN !el subroutine read_angles(kanal,iscor,energ,iprot,*) subroutine read_angles(kanal,*) use geometry_data use energy_data ! 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(len=80) :: lineh integer :: iscor,iprot,ic real(kind=8) :: energ #else subroutine read_angles(kanal,*) use geometry_data ! use energy ! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.GEO' ! include 'COMMON.VAR' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.CONTROL' #endif ! Read angles from input ! !el local variables integer :: i,kanal #ifdef WHAM_RUN 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 ! print *,"energy",energ," iscor",iscor #endif 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 ! 9/7/01 avoid 180 deg valence angle if (theta(i).gt.179.99d0) theta(i)=179.99d0 ! theta(i)=deg2rad*theta(i) phi(i)=deg2rad*phi(i) alph(i)=deg2rad*alph(i) omeg(i)=deg2rad*omeg(i) enddo return 10 return 1 end subroutine read_angles !----------------------------------------------------------------------------- subroutine reada(rekord,lancuch,wartosc,default) ! implicit none character*(*) :: rekord,lancuch real(kind=8) :: wartosc,default integer :: iread !,ilen !el external ilen iread=index(rekord,lancuch) if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,err=10,end=10) wartosc return 10 wartosc=default return end subroutine reada !----------------------------------------------------------------------------- subroutine readi(rekord,lancuch,wartosc,default) ! implicit none character*(*) :: rekord,lancuch integer :: wartosc,default integer :: iread !,ilen !el external ilen iread=index(rekord,lancuch) if (iread.eq.0) then wartosc=default return endif iread=iread+ilen(lancuch)+1 read (rekord(iread:),*,err=10,end=10) wartosc return 10 wartosc=default return end subroutine readi !----------------------------------------------------------------------------- subroutine multreadi(rekord,lancuch,tablica,dim,default) ! implicit none integer :: dim,i integer :: tablica(dim),default character*(*) :: rekord,lancuch character(len=80) :: aux integer :: iread !,ilen !el 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 subroutine multreadi !----------------------------------------------------------------------------- subroutine multreada(rekord,lancuch,tablica,dim,default) ! implicit none integer :: dim,i real(kind=8) :: tablica(dim),default character*(*) :: rekord,lancuch character(len=80) :: aux integer :: iread !,ilen !el 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 subroutine multreada !----------------------------------------------------------------------------- subroutine card_concat(card,to_upper) ! dla UNRESA to_upper jest zawsze .true. ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' character(*) :: card character(len=80) :: karta !,ucase logical :: to_upper !el 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 subroutine card_concat !---------------------------------------------------------------------------- subroutine read_afminp use geometry_data use energy_data use control_data, only:out1file use MPI_data character*320 afmcard real, dimension(3) ::cbeg,cend integer i,j print *, "wchodze" call card_concat(afmcard,.true.) call readi(afmcard,"BEG",afmbeg,0) call readi(afmcard,"END",afmend,0) call reada(afmcard,"FORCE",forceAFMconst,0.0d0) call reada(afmcard,"VEL",velAFMconst,0.0d0) print *,'FORCE=' ,forceAFMconst cbeg=0.0d0 cend=0.0d0 if (afmbeg.eq.-1) then read(inp,*) nbegafmmat,(afmbegcentr(i),i=1,nbegafmmat) do i=1,nbegafmmat do j=1,3 cbeg(j)=cbeg(j)+c(j,afmbegcentr(i))/nbegafmmat enddo enddo else do j=1,3 cbeg(j)=c(j,afmend) enddo endif if (afmend.eq.-1) then read(inp,*) nendafmmat,(afmendcentr(i),i=1,nendafmmat) do i=1,nendafmmat do j=1,3 cend(j)=cend(j)+c(j,afmendcentr(i))/nendafmmat enddo enddo else cend(j)=c(j,afmend) endif !------ NOW PROPERTIES FOR AFM distafminit=0.0d0 do i=1,3 distafminit=(cend(i)-cbeg(i))**2+distafminit enddo distafminit=dsqrt(distafminit) print *,'initdist',distafminit return end subroutine read_afminp !C------------------------------------------------------------ subroutine read_afmnano use geometry_data use energy_data use control_data, only:out1file use MPI_data integer :: i real*8 :: cm read(inp,*,err=11) inanomove,(inanotab(i),i=1,inanomove),forcenanoconst cm=0.0d0 do i=1,inanomove cm=cm+c(3,inanotab(i)) enddo cm=cm/inanomove distnanoinit=cm-tubecenter(3) return 11 write(iout,*)& "error in afmnano",& ", number of center, their index and forceconstance" stop return end subroutine read_afmnano !----------------------------------------------------------------------------- subroutine read_dist_constr use MPI_data ! use control use geometry, only: dist use geometry_data use control_data use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif ! include 'COMMON.SETUP' ! include 'COMMON.CONTROL' ! include 'COMMON.CHAIN' ! include 'COMMON.IOUNITS' ! include 'COMMON.SBRIDGE' integer,dimension(2,100) :: ifrag_,ipair_ real(kind=8),dimension(100) :: wfrag_,wpair_ character(len=640) :: controlcard !el local variables integer :: i,k,j,ii,jj,itemp,mnumkk,mnumjj,k1,j1 integer :: nfrag_,npair_,ndist_ real(kind=8) :: dist_cut,ddjk ! write (iout,*) "Calling read_dist_constr" ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup ! call flush(iout) if(.not.allocated(ihpb)) allocate(ihpb(maxdim)) if(.not.allocated(jhpb)) allocate(jhpb(maxdim)) if(.not.allocated(dhpb)) allocate(dhpb(maxdim)) if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim)) if(.not.allocated(forcon)) allocate(forcon(maxdim)) if(.not.allocated(fordepth)) allocate(fordepth(maxdim)) if(.not.allocated(ibecarb)) allocate(ibecarb(maxdim)) if ((genconstr.gt.0).and.(constr_dist.eq.11)) then call gen_dist_constr2 go to 1712 endif call card_concat(controlcard,.true.) 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 ! if(.not.allocated(ihpb)) allocate(ihpb(maxdim)) ! if(.not.allocated(jhpb)) allocate(jhpb(maxdim)) ! if(.not.allocated(dhpb)) allocate(dhpb(maxdim)) ! if(.not.allocated(forcon)) allocate(forcon(maxdim)) 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 ! 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,nres j1=j k1=k if (j.gt.nres) j1=j-nres if (k.gt.nres) k1=k-nres mnumkk=molnum(k1) mnumjj=molnum(j1) if ((itype(k1,mnumkk).gt.ntyp_molec(mnumkk)).or.& (itype(j1,mnumjj).gt.ntyp_molec(mnumjj))) cycle write (iout,*) "j",j," k",k,itype(k1,mnumkk),itype(j1,mnumjj) 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 #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif 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) #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo enddo endif enddo do i=1,ndist_ ! read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) ! if (forcon(nhpb+1).gt.0.0d0) then ! nhpb=nhpb+1 ! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) 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) !EL fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) fordepth(nhpb+1)=fordepth(nhpb+1)**(0.25d0) forcon(nhpb+1)=forcon(nhpb+1)**(0.25d0) else !C print *,"in else" read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), & ibecarb(i),forcon(nhpb+1) endif 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 #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo 1712 continue call flush(iout) return end subroutine read_dist_constr !----------------------------------------------------------------------------- subroutine gen_dist_constr2 use MPI_data ! use control use geometry, only: dist use geometry_data use control_data use energy_data integer :: i,j real(kind=8) :: distance if (constr_dist.eq.11) then do i=nstart_sup,nstart_sup+nsup-1 do j=i+2,nstart_sup+nsup-1 distance=dist(i,j) if (distance.le.15.0) then nhpb=nhpb+1 ihpb(nhpb)=i+nstart_seq-nstart_sup jhpb(nhpb)=j+nstart_seq-nstart_sup forcon(nhpb)=sqrt(0.04*distance) fordepth(nhpb)=sqrt(40.0/distance) dhpb(nhpb)=distance-0.1d0 dhpb1(nhpb)=distance+0.1d0 #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif endif enddo enddo else do i=nstart_sup,nstart_sup+nsup-1 do j=i+2,nstart_sup+nsup-1 nhpb=nhpb+1 ihpb(nhpb)=i+nstart_seq-nstart_sup jhpb(nhpb)=j+nstart_seq-nstart_sup forcon(nhpb)=weidis dhpb(nhpb)=dist(i,j) enddo enddo endif return end subroutine gen_dist_constr2 !----------------------------------------------------------------------------- #ifdef WINIFL subroutine flush(iu) return end subroutine flush #endif #ifdef AIX subroutine flush(iu) call flush_(iu) return end subroutine flush #endif !----------------------------------------------------------------------------- subroutine copy_to_tmp(source) ! include "DIMENSIONS" ! include "COMMON.IOUNITS" character*(*) :: source character(len=256) :: tmpfile ! integer ilen !el external ilen logical :: ex tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source)) inquire(file=tmpfile,exist=ex) if (ex) then write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),& " to temporary directory..." write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir) endif return end subroutine copy_to_tmp !----------------------------------------------------------------------------- subroutine move_from_tmp(source) ! include "DIMENSIONS" ! include "COMMON.IOUNITS" character*(*) :: source ! integer ilen !el external ilen write (*,*) "Moving ",source(:ilen(source)),& " from temporary directory to working directory" write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir call system("/bin/mv "//source(:ilen(source))//" "//curdir) return end subroutine move_from_tmp !----------------------------------------------------------------------------- ! misc.f !----------------------------------------------------------------------------- ! $Date: 1994/10/12 17:24:21 $ ! $Revision: 2.5 $ logical function find_arg(ipos,line,errflag) integer, parameter :: maxlen = 80 character(len=80) :: line character(len=1) :: empty=' ',equal='=' logical :: errflag integer :: ipos ! This function returns .TRUE., if an argument follows keyword keywd; if so ! IPOS will point to the first non-blank character of the argument. Returns ! .FALSE., if no argument follows the keyword; in this case IPOS points ! to the first non-blank character of the next keyword. do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen) ipos=ipos+1 enddo errflag=.false. if (line(ipos:ipos).eq.equal) then find_arg=.true. ipos=ipos+1 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen) ipos=ipos+1 enddo if (ipos.gt.maxlen) errflag=.true. else find_arg=.false. endif return end function find_arg !----------------------------------------------------------------------------- logical function find_group(iunit,jout,key1) character*(*) :: key1 character(len=80) :: karta !,ucase integer :: iunit,jout integer :: ll !,ilen !EL external ilen !EL logical lcom rewind (iunit) karta=' ' ll=ilen(key1) do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta)) read (iunit,'(a)',end=10) karta enddo write (jout,'(2a)') '> ',karta(1:78) find_group=.true. return 10 find_group=.false. return end function find_group #endif !----------------------------------------------------------------------------- logical function iblnk(charc) character(len=1) :: charc integer :: n n = ichar(charc) iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ') return end function iblnk !----------------------------------------------------------------------------- #ifndef XDRFPDB integer function ilen(string) character*(*) :: string !EL logical :: iblnk ilen = len(string) 1 if ( ilen .gt. 0 ) then if ( iblnk( string(ilen:ilen) ) ) then ilen = ilen - 1 goto 1 endif endif return end function ilen !----------------------------------------------------------------------------- integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset) integer :: nkey,i,ikey,narg character(len=16) :: keywd,keywdset(1:nkey,0:nkey) ! character(len=16) :: ucase do i=1,narg if (ucase(keywd).eq.keywdset(i,ikey)) then ! Match found in_keywd_set=i return endif enddo ! No match to the allowed set of keywords if this point is reached. in_keywd_set=0 return end function in_keywd_set !----------------------------------------------------------------------------- character function lcase(string) integer :: i, k, idiff character*(*) :: string character(len=1) :: c character(len=40) :: chtmp ! i = len(lcase) k = len(string) if (i .lt. k) then k = i if (string(k+1:) .ne. ' ') then chtmp = string endif endif idiff = ichar('a') - ichar('A') lcase = string do 99 i = 1, k c = string(i:i) if (lge(c,'A') .and. lle(c,'Z')) then lcase(i:i) = char(ichar(c) + idiff) endif 99 continue return end function lcase !----------------------------------------------------------------------------- logical function lcom(ipos,karta) character(len=80) :: karta character :: koment(2) = (/'!','#'/) integer :: ipos,i lcom=.false. do i=1,2 if (karta(ipos:ipos).eq.koment(i)) lcom=.true. enddo return end function lcom !----------------------------------------------------------------------------- logical function lower_case(ch) character*(*) :: ch lower_case=(ch.ge.'a' .and. ch.le.'z') return end function lower_case !----------------------------------------------------------------------------- subroutine mykey(line,keywd,ipos,blankline,errflag) ! This subroutine seeks a non-empty substring keywd in the string LINE. ! The substring begins with the first character different from blank and ! "=" encountered right to the pointer IPOS (inclusively) and terminates ! at the character left to the first blank or "=". When the subroutine is ! exited, the pointer IPOS is moved to the position of the terminator in LINE. ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains ! only separators or the maximum length of the data line (80) has been reached. ! The logical variable ERRFLAG is set at .TRUE. if the string ! consists only from a "=". integer, parameter :: maxlen=80 character(len=1) :: empty=' ',equal='=',comma=',' character*(*) :: keywd character(len=80) :: line logical :: blankline,errflag !EL,lcom integer :: ipos,istart,iend errflag=.false. do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen)) ipos=ipos+1 enddo if (ipos.gt.maxlen .or. lcom(ipos,line) ) then ! At this point the rest of the input line turned out to contain only blanks ! or to be commented out. blankline=.true. return endif blankline=.false. istart=ipos ! Checks whether the current char is a separator. do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal & .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen) ipos=ipos+1 enddo iend=ipos-1 ! Error flag set to .true., if the length of the keyword was found less than 1. if (iend.lt.istart) then errflag=.true. return endif keywd=line(istart:iend) return end subroutine mykey !----------------------------------------------------------------------------- subroutine numstr(inum,numm) character(len=10) :: huj='0123456789' character*(*) :: numm integer :: inumm,inum,inum1,inum2 inumm=inum inum1=inumm/10 inum2=inumm-10*inum1 inumm=inum1 numm(3:3)=huj(inum2+1:inum2+1) inum1=inumm/10 inum2=inumm-10*inum1 inumm=inum1 numm(2:2)=huj(inum2+1:inum2+1) inum1=inumm/10 inum2=inumm-10*inum1 inumm=inum1 numm(1:1)=huj(inum2+1:inum2+1) return end subroutine numstr #endif !----------------------------------------------------------------------------- function ucase(string) integer :: i, k, idiff character(*) :: string character(len=len(string)) :: ucase character(len=1) :: c character(len=40) :: chtmp ! i = len(ucase) k = len(string) if (i .lt. k) then k = i if (string(k+1:) .ne. ' ') then chtmp = string endif endif idiff = ichar('a') - ichar('A') ucase = string do 99 i = 1, k c = string(i:i) if (lge(c,'a') .and. lle(c,'z')) then ucase(i:i) = char(ichar(c) - idiff) endif 99 continue return end function ucase !----------------------------------------------------------------------------- ! geomout.F !----------------------------------------------------------------------------- subroutine pdbout(etot,tytul,iunit) use geometry_data, only: c,nres,boxxsize,boxysize,boxzsize use energy_data ! use control use compare_data use MD_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' ! include 'COMMON.IOUNITS' ! include 'COMMON.HEADER' ! include 'COMMON.SBRIDGE' ! include 'COMMON.DISTFIT' ! include 'COMMON.MD' !el character(len=50) :: tytul character*(*) :: tytul character(len=1),dimension(58) :: chainid= (/'A','B','C','D','E','F','G','H','I','J',& 'K','L','M','O','Q','P','R','S','T','U','W','X','Y','Z','a','b','c','d','e','f',& 'g','h','i','j','k','l','m','n','o','p','r','s','t','u','w','x','y','z',& '0','1','2','3','4','5','6','7','8','9'/) #ifdef XDRFPDB integer,dimension(maxres) :: ica !(maxres) #else integer,dimension(nres) :: ica !(maxres) #endif integer :: iti1 !el local variables integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit real(kind=8) :: etot,xi,yi,zi integer :: nres2 nres2=2*nres #ifdef XDRFPDB if(.not.allocated(molnum)) allocate(molnum(maxres)) ! molnum(:)=mnumi(:) if(.not.allocated(vtot)) allocate(vtot(maxres*2)) !(maxres2) if(.not.allocated(istype)) allocate(istype(maxres)) #else if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2) #endif write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot !model write (iunit,'(a5,i6)') 'MODEL',1 #ifndef XDRFPDB if (nhfrag.gt.0) then do j=1,nhfrag iti=itype(hfrag(1,j),1) itj=itype(hfrag(2,j),1) if (j.lt.10) then write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') & 'HELIX',j,'H',j,& restyp(iti,1),hfrag(1,j)-1,& restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j) else write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') & 'HELIX',j,'H',j,& restyp(iti,1),hfrag(1,j)-1,& restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j) endif enddo endif if (nbfrag.gt.0) then do j=1,nbfrag iti=itype(bfrag(1,j),1) itj=itype(bfrag(2,j)-1,1) write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') & 'SHEET',1,'B',j,2,& restyp(iti,1),bfrag(1,j)-1,& restyp(itj,1),bfrag(2,j)-2,0 if (bfrag(3,j).gt.bfrag(4,j)) then itk=itype(bfrag(3,j),1) itl=itype(bfrag(4,j)+1,1) write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') & 'SHEET',2,'B',j,2,& restyp(itl,1),bfrag(4,j),& restyp(itk,1),bfrag(3,j)-1,-1,& "N",restyp(itk,1),bfrag(3,j)-1,& "O",restyp(iti,1),bfrag(1,j)-1 else itk=itype(bfrag(3,j),1) itl=itype(bfrag(4,j)-1,1) write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') & 'SHEET',2,'B',j,2,& restyp(itk,1),bfrag(3,j)-1,& restyp(itl,1),bfrag(4,j)-2,1,& "N",restyp(itk,1),bfrag(3,j)-1,& "O",restyp(iti,1),bfrag(1,j)-1 endif enddo endif #endif if (nss.gt.0) then do i=1,nss if (dyn_ss) then write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') & 'SSBOND',i,'CYS',idssb(i)-nnt+1,& 'CYS',jdssb(i)-nnt+1 else write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') & 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,& 'CYS',jhpb(i)-nnt+1-nres endif enddo endif iatom=0 ichain=1 ires=0 write(iout,*) "TUTUT" do i=nnt,nct ! write(iout,*), "coord",c(1,i),c(2,i),c(3,i) iti=itype(i,molnum(i)) print *,i,molnum(i) if (molnum(i+1).eq.0) then iti1=ntyp1_molec(molnum(i)) else iti1=itype(i+1,molnum(i+1)) endif if ((iti.eq.ntyp1_molec(molnum(i))).and.(iti1.eq.ntyp1_molec(molnum(i)))) cycle if (i.lt.nnt) then if (iti.eq.ntyp1_molec(molnum(i)).and.(molnum(i+1).eq.5)) cycle endif if (iti.eq.ntyp1_molec(molnum(i))) then ichain=ichain+1 if (ichain.gt.58) ichain=1 ! ires=0 write (iunit,'(a)') 'TER' else ires=ires+1 iatom=iatom+1 ica(i)=iatom if (molnum(i).eq.1) then write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,(c(j,i),j=1,3),vtot(i) elseif(molnum(i).eq.2) then if (istype(i).eq.0) istype(i)=1 write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), & chainid(ichain),ires,(c(j,i),j=1,3),vtot(i) elseif (molnum(i).eq.4) then xi=c(1,i) yi=c(2,i) zi=c(3,i) ! xi=dmod(xi,boxxsize) ! if (xi.lt.0.0d0) xi=xi+boxxsize ! yi=dmod(yi,boxysize) ! if (yi.lt.0.0d0) yi=yi+boxysize ! zi=dmod(zi,boxzsize) ! if (zi.lt.0.0d0) zi=zi+boxzsize write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,xi,yi,zi,vtot(i) elseif (molnum(i).eq.5) then xi=c(1,i) yi=c(2,i) zi=c(3,i) xi=dmod(xi,boxxsize) if (xi.lt.0.0d0) xi=xi+boxxsize yi=dmod(yi,boxysize) if (yi.lt.0.0d0) yi=yi+boxysize zi=dmod(zi,boxzsize) if (zi.lt.0.0d0) zi=zi+boxzsize write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,xi,yi,zi,vtot(i) else write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),& ires,(c(j,i),j=1,3),vtot(i) endif if ((iti.ne.10).and.(molnum(i).ne.5)) then iatom=iatom+1 if (molnum(i).eq.1) then write (iunit,20) iatom,restyp(iti,1),chainid(ichain),& ires,(c(j,nres+i),j=1,3),& vtot(i+nres) else if (molnum(i).eq.2) then if (istype(i).eq.0) istype(i)=1 write (iunit,50) iatom,sugartyp(istype(i)),restyp(iti,2), & chainid(ichain),ires,(c(j,nres+i),j=1,3),vtot(i+nres) endif endif endif enddo write (iunit,'(a)') 'TER' do i=nnt,nct-1 if (molnum(i).eq.5) cycle if (itype(i,1).eq.ntyp1) cycle if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then write (iunit,30) ica(i),ica(i+1) else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then write (iunit,30) ica(i),ica(i+1),ica(i)+1 else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then write (iunit,30) ica(i),ica(i)+1 endif enddo if ((itype(nct,1).ne.10).and.(molnum(i).ne.5)) then write (iunit,30) ica(nct),ica(nct)+1 endif do i=1,nss if (dyn_ss) then write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1 else write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1 endif enddo write (iunit,'(a6)') 'ENDMDL' 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3) 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3) 40 FORMAT ("ATOM",I7," C5' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3) 50 FORMAT ("ATOM",I7," C1' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3) 30 FORMAT ('CONECT',8I5) 60 FORMAT ('HETATM',I5,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3) return end subroutine pdbout #ifndef XDRFPDB !----------------------------------------------------------------------------- subroutine MOL2out(etot,tytul) ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2 ! format. use geometry_data, only: c use energy_data ! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' ! include 'COMMON.IOUNITS' ! include 'COMMON.HEADER' ! include 'COMMON.SBRIDGE' character(len=32) :: tytul,fd character(len=3) :: zahl character(len=6) :: res_num,pom !,ucase !el local variables integer :: i,j real(kind=8) :: etot #ifdef AIX call fdate_(fd) #elif (defined CRAY) call date(fd) #else call fdate(fd) #endif write (imol2,'(a)') '#' write (imol2,'(a)') & '# Creating user name: unres' write (imol2,'(2a)') '# Creation time: ',& fd write (imol2,'(/a)') '\@MOLECULE' write (imol2,'(a)') tytul write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0 write (imol2,'(a)') 'SMALL' write (imol2,'(a)') 'USER_CHARGES' write (imol2,'(a)') '\@ATOM' do i=nnt,nct write (zahl,'(i3)') i pom=ucase(restyp(itype(i,1),1)) res_num = pom(:3)//zahl(2:) write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0 enddo write (imol2,'(a)') '\@BOND' do i=nnt,nct-1 write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1 enddo do i=1,nss write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1 enddo write (imol2,'(a)') '\@SUBSTRUCTURE' do i=nnt,nct write (zahl,'(i3)') i pom = ucase(restyp(itype(i,1),1)) res_num = pom(:3)//zahl(2:) write (imol2,30) i-nnt+1,res_num,i-nnt+1,0 enddo 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****') 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****') return end subroutine MOL2out !----------------------------------------------------------------------------- subroutine intout use geometry_data use energy_data, only: itype ! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.VAR' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' ! include 'COMMON.GEO' ! include 'COMMON.TORSION' !el local variables integer :: i,iti write (iout,'(/a)') 'Geometry of the virtual chain.' write (iout,'(7a)') ' Res ',' d',' Theta',& ' Phi',' Dsc',' Alpha',' Omega' do i=1,nres iti=itype(i,1) ! print *,vbld(i),"vbld(i)",i write (iout,'(a3,i6,6f10.3)') restyp(iti,1),i,vbld(i),& rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),& rad2deg*omeg(i) enddo return end subroutine intout !----------------------------------------------------------------------------- #ifdef CLUSTER subroutine briefout(it,ener,free)!,plik) #else subroutine briefout(it,ener) #endif use geometry_data use energy_data ! use control ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.IOUNITS' ! include 'COMMON.CHAIN' ! include 'COMMON.VAR' ! include 'COMMON.LOCAL' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' ! include 'COMMON.GEO' ! include 'COMMON.SBRIDGE' ! print '(a,i5)',intname,igeom !el local variables integer :: i,it real(kind=8) :: ener,free ! character(len=80) :: plik ! integer :: iii #if defined(AIX) || defined(PGI) open (igeom,file=intname,position='append') #else open (igeom,file=intname,access='append') #endif #ifdef WHAM_RUN ! iii=igeom igeom=iout #endif print *,nss IF (NSS.LE.9) THEN #ifdef CLUSTER WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS) ELSE WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8) #else WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS) ELSE WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9) #endif WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS) ENDIF ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES) WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES) ! if (nvar.gt.nphi+ntheta) then write (igeom,200) (rad2deg*alph(i),i=2,nres-1) write (igeom,200) (rad2deg*omeg(i),i=2,nres-1) ! endif close(igeom) 180 format (I5,F12.3,I2,9(1X,2I3)) 190 format (3X,11(1X,2I3)) 200 format (8F10.4) return end subroutine briefout !----------------------------------------------------------------------------- #ifdef WINIFL subroutine fdate(fd) character(len=32) :: fd write(fd,'(32x)') return end subroutine fdate #endif !----------------------------------------------------------------------------- #ifdef WHAM_RUN real(kind=8) function gyrate(jcon) #else real(kind=8) function gyrate() #endif use geometry_data, only: c ! use geometry use energy_data ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.INTERACT' ! include 'COMMON.CHAIN' real(kind=8) :: rg real(kind=8),dimension(3) :: cen !el local variables integer :: i,j,jcon do j=1,3 cen(j)=0.0d0 enddo do i=nnt,nct do j=1,3 cen(j)=cen(j)+c(j,i) enddo enddo do j=1,3 cen(j)=cen(j)/dble(nct-nnt+1) enddo rg = 0.0d0 do i = nnt, nct do j=1,3 rg = rg + (c(j,i)-cen(j))**2 enddo end do #ifdef WHAM_RUN gyrate = dsqrt(rg/dble(nct-nnt+1)) #else gyrate = sqrt(rg/dble(nct-nnt+1)) #endif return end function gyrate #ifdef WHAM_RUN !----------------------------------------------------------------------------- ! readrtns.F WHAM subroutine reads(rekord,lancuch,wartosc,default) ! implicit none character*(*) :: rekord,lancuch,wartosc,default character(len=80) :: aux integer :: lenlan,lenrec,iread,ireade !el external ilen !el logical iblnk !el external iblnk lenlan=ilen(lancuch) lenrec=ilen(rekord) iread=index(rekord,lancuch(:lenlan)//"=") ! print *,"rekord",rekord," lancuch",lancuch ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec if (iread.eq.0) then wartosc=default return endif iread=iread+lenlan+1 ! print *,"iread",iread ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread)) do while (iread.le.lenrec .and. iblnk(rekord(iread:iread))) iread=iread+1 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread)) enddo ! print *,"iread",iread if (iread.gt.lenrec) then wartosc=default return endif ireade=iread+1 ! print *,"ireade",ireade do while (ireade.lt.lenrec .and. & .not.iblnk(rekord(ireade:ireade))) ireade=ireade+1 enddo wartosc=rekord(iread:ireade) return end subroutine reads #endif !----------------------------------------------------------------------------- ! permut.F !----------------------------------------------------------------------------- subroutine permut(isym,nperm,tabperm) !c integer maxperm,maxsym !c parameter (maxperm=3628800) !c parameter (maxsym=10) ! include "DIMENSIONS" integer n,a,tabperm,nperm,kkk,i,isym ! logical nextp ! external nextp dimension a(isym),tabperm(50,5040) n=isym nperm=1 if (n.eq.1) then tabperm(1,1)=1 return endif do i=2,n nperm=nperm*i enddo kkk=0 do i=1,n a(i)=i enddo 10 continue !c print '(i3,2x,100i3)',kkk+1,(a(i),i=1,n) kkk=kkk+1 do i=1,n tabperm(i,kkk)=a(i) enddo if(nextp(n,a)) go to 10 return end subroutine !----------------------------------------------------------------------------- logical function nextp(n,a) integer :: n,i,j,k,t ! logical :: nextp integer,dimension(n) :: a i=n-1 10 if(a(i).lt.a(i+1)) go to 20 i=i-1 if(i.eq.0) go to 20 go to 10 20 j=i+1 k=n 30 t=a(j) a(j)=a(k) a(k)=t j=j+1 k=k-1 if(j.lt.k) go to 30 j=i if(j.ne.0) go to 40 nextp=.false. return 40 j=j+1 if(a(j).lt.a(i)) go to 40 t=a(i) a(i)=a(j) a(j)=t nextp=.true. return end function nextp #endif !----------------------------------------------------------------------------- ! rescode.f !----------------------------------------------------------------------------- integer function rescode(iseq,nam,itype,molecule) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.NAMES' ! include 'COMMON.IOUNITS' character(len=3) :: nam !,ucase integer :: iseq,itype,i integer :: molecule print *,molecule,nam if (molecule.eq.1) then if (itype.eq.0) then do i=-ntyp1_molec(molecule),ntyp1_molec(molecule) if (ucase(nam).eq.restyp(i,molecule)) then rescode=i return endif enddo else do i=-ntyp1_molec(molecule),ntyp1_molec(molecule) if (nam(1:1).eq.onelet(i)) then rescode=i return endif enddo endif else if (molecule.eq.2) then do i=1,ntyp1_molec(molecule) print *,nam(1:1),restyp(i,molecule)(1:1) print *,nam(2:2),"read",i if (nam(2:2).eq.restyp(i,molecule)(1:1)) then rescode=i return endif enddo else if (molecule.eq.3) then write(iout,*) "SUGAR not yet implemented" stop else if (molecule.eq.4) then write(iout,*) "Explicit LIPID not yet implemented" stop else if (molecule.eq.5) then do i=-1,ntyp1_molec(molecule) print *,restyp(i,molecule) print *,i,restyp(i,molecule)(1:2),nam(1:2) if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then rescode=i return endif enddo else write(iout,*) "molecule not defined" endif write (iout,10) iseq,nam stop 10 format ('**** Error - residue',i4,' has an unresolved name ',a3) end function rescode integer function sugarcode(sugar,ires) character sugar integer ires if (sugar.eq.'D') then sugarcode=1 else if (sugar.eq.' ') then sugarcode=2 else write (iout,*) 'UNKNOWN sugar type for residue',ires,' ',sugar stop endif return end function sugarcode integer function rescode_lip(res,ires) character(len=3):: res integer ires rescode_lip=0 if (res.eq.'NC3') rescode_lip =4 if (res.eq.'NH3') rescode_lip =2 if (res.eq.'PO4') rescode_lip =3 if (res.eq.'GL1') rescode_lip =12 if (res.eq.'GL2') rescode_lip =12 if (res.eq.'C1A') rescode_lip =18 if (res.eq.'C2A') rescode_lip =18 if (res.eq.'C3A') rescode_lip =18 if (res.eq.'C4A') rescode_lip =18 if (res.eq.'C1B') rescode_lip =18 if (res.eq.'C2B') rescode_lip =18 if (res.eq.'C3B') rescode_lip =18 if (res.eq.'C4B') rescode_lip =18 if (res.eq.'D2A') rescode_lip =16 if (res.eq.'D2B') rescode_lip =16 if (rescode_lip.eq.0) write(iout,*) "UNKNOWN RESIDUE",ires,res return end function !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module io_base