source/unres/io_base.f90

   1       module io_base
   2 !-----------------------------------------------------------------------
   3       use names
   4       use io_units
   5       implicit none
   6 !-----------------------------------------------------------------------------
   7 ! Max. number of AA residues
   8       integer,parameter :: maxres=6000!1200
   9 ! Appr. max. number of interaction sites
  10       integer,parameter :: maxres2=2*maxres
  11 !      parameter (maxres6=6*maxres)
  12 !      parameter (mmaxres2=(maxres2*(maxres2+1)/2))
  13 !-----------------------------------------------------------------------------
  14 ! Max. number of S-S bridges
  15 !      integer,parameter :: maxss=20
  16 !-----------------------------------------------------------------------------
  17 ! Max. number of derivatives of virtual-bond and side-chain vectors in theta
  18 ! or phi.
  19 !      integer,parameter :: maxdim=(maxres-1)*(maxres-2)/2
  20 !-----------------------------------------------------------------------------
  21 !
  22 !
  23 !-----------------------------------------------------------------------------
  24       contains
  25 !-----------------------------------------------------------------------------
  26 ! readrtns_CSA.F
  27 !-----------------------------------------------------------------------------
  28       subroutine read_bridge
  29 ! Read information about disulfide bridges.
  30       use geometry_data, only: nres
  31       use energy_data
  32       use control_data, only:out1file
  33       use MPI_data
  34 !      implicit real*8 (a-h,o-z)
  35 !      include 'DIMENSIONS'
  36 #ifdef MPI
  37       include 'mpif.h'
  38 #endif
  39 !      include 'COMMON.IOUNITS'
  40 !      include 'COMMON.GEO'
  41 !      include 'COMMON.VAR'
  42 !      include 'COMMON.INTERACT'
  43 !      include 'COMMON.LOCAL'
  44 !      include 'COMMON.NAMES'
  45 !      include 'COMMON.CHAIN'
  46 !      include 'COMMON.FFIELD'
  47 !      include 'COMMON.SBRIDGE'
  48 !      include 'COMMON.HEADER'
  49 !      include 'COMMON.CONTROL'
  50 !      include 'COMMON.DBASE'
  51 !      include 'COMMON.THREAD'
  52 !      include 'COMMON.TIME1'
  53 !      include 'COMMON.SETUP'
  54 !el local variables
  55       integer :: i,j,ierror
  56
  57 ! Read bridging residues.
  58       read (inp,*) ns
  59       if (ns.gt.0) then
  60         allocate(iss(ns))
  61         read (inp,*) (iss(i),i=1,ns)
  62       endif
  63
  64 !      print *,'ns=',ns
  65       if(me.eq.king.or..not.out1file) &
  66         write (iout,*) 'ns=',ns
  67       if (ns.gt.0) &
  68         write(iout,*) ' iss:',(iss(i),i=1,ns)
  69 ! Check whether the specified bridging residues are cystines.
  70       do i=1,ns
  71         if (itype(iss(i)).ne.1) then
  72           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
  73          'Do you REALLY think that the residue ',&
  74           restyp(itype(iss(i))),i,&
  75          ' can form a disulfide bridge?!!!'
  76           write (*,'(2a,i3,a)') &
  77          'Do you REALLY think that the residue ',&
  78           restyp(itype(iss(i))),i,&
  79          ' can form a disulfide bridge?!!!'
  80 #ifdef MPI
  81          call MPI_Finalize(MPI_COMM_WORLD,ierror)
  82          stop
  83 #endif
  84         endif
  85       enddo
  86 ! Read preformed bridges.
  87       if (ns.gt.0) then
  88         read (inp,*) nss
  89       if (nss.gt.0) then
  90         if(.not.allocated(ihpb)) allocate(ihpb(nss))
  91         if(.not.allocated(jhpb)) allocate(jhpb(nss))
  92         read (inp,*) (ihpb(i),jhpb(i),i=1,nss)
  93
  94         if(.not.allocated(dhpb)) allocate(dhpb(nss))
  95         if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2
  96         if(.not.allocated(dhpb1)) allocate(dhpb1(nss))
  97         if(.not.allocated(ibecarb)) allocate(ibecarb(nss))
  98 ! el Initialize the bridge array
  99         do i=1,nss
 100           dhpb(i)=0.0D0
 101         enddo
 102       endif
 103 !--------------------
 104       if(fg_rank.eq.0) &
 105        write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 106       if (nss.gt.0) then
 107        write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 108         nhpb=nss
 109 ! Check if the residues involved in bridges are in the specified list of
 110 ! bridging residues.
 111         do i=1,nss
 112           do j=1,i-1
 113             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) &
 114             .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 115               write (iout,'(a,i3,a)') 'Disulfide pair',i,&
 116             ' contains residues present in other pairs.'
 117               write (*,'(a,i3,a)') 'Disulfide pair',i,&
 118             ' contains residues present in other pairs.'
 119 #ifdef MPI
 120               call MPI_Finalize(MPI_COMM_WORLD,ierror)
 121               stop
 122 #endif
 123             endif
 124           enddo
 125           do j=1,ns
 126             if (ihpb(i).eq.iss(j)) goto 10
 127           enddo
 128           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 129    10     continue
 130           do j=1,ns
 131             if (jhpb(i).eq.iss(j)) goto 20
 132           enddo
 133           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 134    20     continue
 135           dhpb(i)=d0cm
 136           forcon(i)=akcm
 137         enddo
 138         do i=1,nss
 139           ihpb(i)=ihpb(i)+nres
 140           jhpb(i)=jhpb(i)+nres
 141         enddo
 142       endif
 143       endif
 144 !      write(iout,*) "end read_bridge"
 145       return
 146       end subroutine read_bridge
 147 !-----------------------------------------------------------------------------
 148       subroutine read_x(kanal,*)
 149
 150      use geometry_data
 151      use energy_data
 152      use geometry, only:int_from_cart1
 153 !      implicit real*8 (a-h,o-z)
 154 !      include 'DIMENSIONS'
 155 !      include 'COMMON.GEO'
 156 !      include 'COMMON.VAR'
 157 !      include 'COMMON.CHAIN'
 158 !      include 'COMMON.IOUNITS'
 159 !      include 'COMMON.CONTROL'
 160 !      include 'COMMON.LOCAL'
 161 !      include 'COMMON.INTERACT'
 162 ! Read coordinates from input
 163 !
 164 !el local variables
 165       integer :: l,k,j,i,kanal
 166
 167       read(kanal,'(8f10.5)',end=10,err=10) &
 168         ((c(l,k),l=1,3),k=1,nres),&
 169         ((c(l,k+nres),l=1,3),k=nnt,nct)
 170       do j=1,3
 171         c(j,nres+1)=c(j,1)
 172         c(j,2*nres)=c(j,nres)
 173       enddo
 174       call int_from_cart1(.false.)
 175       do i=1,nres-1
 176         do j=1,3
 177           dc(j,i)=c(j,i+1)-c(j,i)
 178           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 179         enddo
 180       enddo
 181       do i=nnt,nct
 182         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 183           do j=1,3
 184             dc(j,i+nres)=c(j,i+nres)-c(j,i)
 185             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 186           enddo
 187         endif
 188       enddo
 189
 190       return
 191    10 return 1
 192       end subroutine read_x
 193 !-----------------------------------------------------------------------------
 194       subroutine read_threadbase
 195
 196       use geometry_data
 197       use energy_data
 198       use compare_data
 199 !      implicit real*8 (a-h,o-z)
 200 !      include 'DIMENSIONS'
 201 !      include 'COMMON.IOUNITS'
 202 !      include 'COMMON.GEO'
 203 !      include 'COMMON.VAR'
 204 !      include 'COMMON.INTERACT'
 205 !      include 'COMMON.LOCAL'
 206 !      include 'COMMON.NAMES'
 207 !      include 'COMMON.CHAIN'
 208 !      include 'COMMON.FFIELD'
 209 !      include 'COMMON.SBRIDGE'
 210 !      include 'COMMON.HEADER'
 211 !      include 'COMMON.CONTROL'
 212 !      include 'COMMON.DBASE'
 213 !      include 'COMMON.THREAD'
 214 !      include 'COMMON.TIME1'
 215
 216 !el local variables
 217       integer :: k,j,i
 218
 219 ! Read pattern database for threading.
 220       read (icbase,*) nseq
 221       allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq)
 222       allocate(nres_base(3,nseq)) !(3,maxseq)
 223       allocate(str_nam(nseq)) !(maxseq)
 224       do i=1,nseq
 225         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),&
 226          nres_base(2,i),nres_base(3,i)
 227         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,&
 228          nres_base(1,i))
 229 !       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
 230 !    &   nres_base(2,i),nres_base(3,i)
 231 !       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
 232 !    &   nres_base(1,i))
 233       enddo
 234       close (icbase)
 235       if (weidis.eq.0.0D0) weidis=0.1D0
 236       do i=nnt,nct
 237         do j=i+2,nct
 238           nhpb=nhpb+1
 239           ihpb(nhpb)=i
 240           jhpb(nhpb)=j
 241           forcon(nhpb)=weidis
 242         enddo
 243       enddo
 244       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
 245       write (iout,'(a,i5)') 'nexcl: ',nexcl
 246       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
 247       return
 248       end subroutine read_threadbase
 249 !-----------------------------------------------------------------------------
 250 #ifdef WHAM_RUN
 251 !el      subroutine read_angles(kanal,iscor,energ,iprot,*)
 252       subroutine read_angles(kanal,*)
 253
 254       use geometry_data
 255       use energy_data
 256 !      implicit real*8 (a-h,o-z)
 257 !      include 'DIMENSIONS'
 258 !      include 'DIMENSIONS.ZSCOPT'
 259 !      include 'COMMON.INTERACT'
 260 !      include 'COMMON.SBRIDGE'
 261 !      include 'COMMON.GEO'
 262 !      include 'COMMON.VAR'
 263 !      include 'COMMON.CHAIN'
 264 !      include 'COMMON.IOUNITS'
 265       character(len=80) :: lineh
 266       integer :: iscor,iprot,ic
 267       real(kind=8) :: energ
 268 #else
 269       subroutine read_angles(kanal,*)
 270
 271       use geometry_data
 272    !  use energy
 273    !  use control
 274 !      implicit real*8 (a-h,o-z)
 275 !      include 'DIMENSIONS'
 276 !      include 'COMMON.GEO'
 277 !      include 'COMMON.VAR'
 278 !      include 'COMMON.CHAIN'
 279 !      include 'COMMON.IOUNITS'
 280 !      include 'COMMON.CONTROL'
 281 #endif
 282 ! Read angles from input
 283 !
 284 !el local variables
 285       integer :: i,kanal
 286 #ifdef WHAM_RUN
 287       read(kanal,'(a80)',end=10,err=10) lineh
 288       read(lineh(:5),*,err=8) ic
 289       read(lineh(6:),*,err=8) energ
 290       goto 9
 291     8 ic=1
 292       print *,'error, assuming e=1d10',lineh
 293       energ=1d10
 294       nss=0
 295     9 continue
 296       read(lineh(18:),*,end=10,err=10) nss
 297       IF (NSS.LT.9) THEN
 298         read (lineh(20:),*,end=10,err=10) &
 299         (IHPB(I),JHPB(I),I=1,NSS),iscor
 300       ELSE
 301         read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
 302         read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),&
 303           I=9,NSS),iscor
 304       ENDIF
 305 !      print *,"energy",energ," iscor",iscor
 306 #endif
 307        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 308        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 309        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 310        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 311
 312        do i=1,nres
 313 ! 9/7/01 avoid 180 deg valence angle
 314         if (theta(i).gt.179.99d0) theta(i)=179.99d0
 315 !
 316         theta(i)=deg2rad*theta(i)
 317         phi(i)=deg2rad*phi(i)
 318         alph(i)=deg2rad*alph(i)
 319         omeg(i)=deg2rad*omeg(i)
 320        enddo
 321       return
 322    10 return 1
 323       end subroutine read_angles
 324 !-----------------------------------------------------------------------------
 325       subroutine reada(rekord,lancuch,wartosc,default)
 326
 327 !      implicit none
 328       character*(*) :: rekord,lancuch
 329       real(kind=8) :: wartosc,default
 330       integer :: iread !,ilen
 331 !el      external ilen
 332       iread=index(rekord,lancuch)
 333       if (iread.eq.0) then
 334         wartosc=default
 335         return
 336       endif
 337       iread=iread+ilen(lancuch)+1
 338       read (rekord(iread:),*,err=10,end=10) wartosc
 339       return
 340   10  wartosc=default
 341       return
 342       end subroutine reada
 343 !-----------------------------------------------------------------------------
 344       subroutine readi(rekord,lancuch,wartosc,default)
 345
 346 !      implicit none
 347       character*(*) :: rekord,lancuch
 348       integer :: wartosc,default
 349       integer :: iread !,ilen
 350 !el      external ilen
 351       iread=index(rekord,lancuch)
 352       if (iread.eq.0) then
 353         wartosc=default
 354         return
 355       endif
 356       iread=iread+ilen(lancuch)+1
 357       read (rekord(iread:),*,err=10,end=10) wartosc
 358       return
 359   10  wartosc=default
 360       return
 361       end subroutine readi
 362 !-----------------------------------------------------------------------------
 363       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 364
 365 !      implicit none
 366       integer :: dim,i
 367       integer :: tablica(dim),default
 368       character*(*) :: rekord,lancuch
 369       character(len=80) :: aux
 370       integer :: iread !,ilen
 371 !el      external ilen
 372       do i=1,dim
 373         tablica(i)=default
 374       enddo
 375       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 376       if (iread.eq.0) return
 377       iread=iread+ilen(lancuch)+1
 378       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 379    10 return
 380       end subroutine multreadi
 381 !-----------------------------------------------------------------------------
 382       subroutine multreada(rekord,lancuch,tablica,dim,default)
 383
 384 !      implicit none
 385       integer :: dim,i
 386       real(kind=8) :: tablica(dim),default
 387       character*(*) :: rekord,lancuch
 388       character(len=80) :: aux
 389       integer :: iread !,ilen
 390 !el      external ilen
 391       do i=1,dim
 392         tablica(i)=default
 393       enddo
 394       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 395       if (iread.eq.0) return
 396       iread=iread+ilen(lancuch)+1
 397       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 398    10 return
 399       end subroutine multreada
 400 !-----------------------------------------------------------------------------
 401       subroutine card_concat(card,to_upper)
 402
 403 ! dla UNRESA to_upper jest zawsze .true.
 404 !      implicit real*8 (a-h,o-z)
 405 !      include 'DIMENSIONS'
 406 !      include 'COMMON.IOUNITS'
 407       character(*) :: card
 408       character(len=80) :: karta        !,ucase
 409       logical :: to_upper
 410 !el      external ilen
 411       read (inp,'(a)') karta
 412       if (to_upper) karta=ucase(karta)
 413       card=' '
 414       do while (karta(80:80).eq.'&')
 415         card=card(:ilen(card)+1)//karta(:79)
 416         read (inp,'(a)') karta
 417         if (to_upper) karta=ucase(karta)
 418       enddo
 419       card=card(:ilen(card)+1)//karta
 420       return
 421       end subroutine card_concat
 422 !-----------------------------------------------------------------------------
 423       subroutine read_dist_constr
 424       use MPI_data
 425    !  use control
 426       use geometry, only: dist
 427       use geometry_data
 428       use control_data
 429       use energy_data
 430 !      implicit real*8 (a-h,o-z)
 431 !      include 'DIMENSIONS'
 432 #ifdef MPI
 433       include 'mpif.h'
 434 #endif
 435 !      include 'COMMON.SETUP'
 436 !      include 'COMMON.CONTROL'
 437 !      include 'COMMON.CHAIN'
 438 !      include 'COMMON.IOUNITS'
 439 !      include 'COMMON.SBRIDGE'
 440       integer,dimension(2,100) :: ifrag_,ipair_
 441       real(kind=8),dimension(100) :: wfrag_,wpair_
 442       character(len=640) :: controlcard
 443
 444 !el local variables
 445       integer :: i,k,j,ddjk,ii,jj,itemp
 446       integer :: nfrag_,npair_,ndist_
 447       real(kind=8) :: dist_cut
 448
 449 !      write (iout,*) "Calling read_dist_constr"
 450 !      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 451 !      call flush(iout)
 452       call card_concat(controlcard,.true.)
 453       call readi(controlcard,"NFRAG",nfrag_,0)
 454       call readi(controlcard,"NPAIR",npair_,0)
 455       call readi(controlcard,"NDIST",ndist_,0)
 456       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 457       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 458       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 459       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 460       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 461 !      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 462 !      write (iout,*) "IFRAG"
 463 !      do i=1,nfrag_
 464 !        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 465 !      enddo
 466 !      write (iout,*) "IPAIR"
 467 !      do i=1,npair_
 468 !        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 469 !      enddo
 470       if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
 471       if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
 472       if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
 473       if(.not.allocated(forcon)) allocate(forcon(maxdim))
 474
 475       call flush(iout)
 476       do i=1,nfrag_
 477         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 478         if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
 479           ifrag_(2,i)=nstart_sup+nsup-1
 480 !        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 481         call flush(iout)
 482         if (wfrag_(i).gt.0.0d0) then
 483         do j=ifrag_(1,i),ifrag_(2,i)-1
 484           do k=j+1,ifrag_(2,i)
 485 !            write (iout,*) "j",j," k",k
 486             ddjk=dist(j,k)
 487             if (constr_dist.eq.1) then
 488               nhpb=nhpb+1
 489               ihpb(nhpb)=j
 490               jhpb(nhpb)=k
 491               dhpb(nhpb)=ddjk
 492               forcon(nhpb)=wfrag_(i)
 493             else if (constr_dist.eq.2) then
 494               if (ddjk.le.dist_cut) then
 495                 nhpb=nhpb+1
 496                 ihpb(nhpb)=j
 497                 jhpb(nhpb)=k
 498                 dhpb(nhpb)=ddjk
 499                 forcon(nhpb)=wfrag_(i)
 500               endif
 501             else
 502               nhpb=nhpb+1
 503               ihpb(nhpb)=j
 504               jhpb(nhpb)=k
 505               dhpb(nhpb)=ddjk
 506               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 507             endif
 508 #ifdef MPI
 509             if (.not.out1file .or. me.eq.king) &
 510             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
 511              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 512 #else
 513             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
 514              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 515 #endif
 516           enddo
 517         enddo
 518         endif
 519       enddo
 520       do i=1,npair_
 521         if (wpair_(i).gt.0.0d0) then
 522         ii = ipair_(1,i)
 523         jj = ipair_(2,i)
 524         if (ii.gt.jj) then
 525           itemp=ii
 526           ii=jj
 527           jj=itemp
 528         endif
 529         do j=ifrag_(1,ii),ifrag_(2,ii)
 530           do k=ifrag_(1,jj),ifrag_(2,jj)
 531             nhpb=nhpb+1
 532             ihpb(nhpb)=j
 533             jhpb(nhpb)=k
 534             forcon(nhpb)=wpair_(i)
 535             dhpb(nhpb)=dist(j,k)
 536 #ifdef MPI
 537             if (.not.out1file .or. me.eq.king) &
 538             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 539              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 540 #else
 541             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 542              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 543 #endif
 544           enddo
 545         enddo
 546         endif
 547       enddo
 548       do i=1,ndist_
 549         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
 550         if (forcon(nhpb+1).gt.0.0d0) then
 551           nhpb=nhpb+1
 552           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 553 #ifdef MPI
 554           if (.not.out1file .or. me.eq.king) &
 555           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 556            nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 557 #else
 558           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 559            nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 560 #endif
 561         endif
 562       enddo
 563       call flush(iout)
 564       return
 565       end subroutine read_dist_constr
 566 !-----------------------------------------------------------------------------
 567 #ifdef WINIFL
 568       subroutine flush(iu)
 569       return
 570       end subroutine flush
 571 #endif
 572 #ifdef AIX
 573       subroutine flush(iu)
 574       call flush_(iu)
 575       return
 576       end subroutine flush
 577 #endif
 578 !-----------------------------------------------------------------------------
 579       subroutine copy_to_tmp(source)
 580
 581 !      include "DIMENSIONS"
 582 !      include "COMMON.IOUNITS"
 583       character*(*) :: source
 584       character(len=256) :: tmpfile
 585 !      integer ilen
 586 !el      external ilen
 587       logical :: ex
 588       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
 589       inquire(file=tmpfile,exist=ex)
 590       if (ex) then
 591         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
 592          " to temporary directory..."
 593         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
 594         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
 595       endif
 596       return
 597       end subroutine copy_to_tmp
 598 !-----------------------------------------------------------------------------
 599       subroutine move_from_tmp(source)
 600
 601 !      include "DIMENSIONS"
 602 !      include "COMMON.IOUNITS"
 603       character*(*) :: source
 604 !      integer ilen
 605 !el      external ilen
 606       write (*,*) "Moving ",source(:ilen(source)),&
 607        " from temporary directory to working directory"
 608       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
 609       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
 610       return
 611       end subroutine move_from_tmp
 612 !-----------------------------------------------------------------------------
 613 ! misc.f
 614 !-----------------------------------------------------------------------------
 615 ! $Date: 1994/10/12 17:24:21 $
 616 ! $Revision: 2.5 $
 617
 618       logical function find_arg(ipos,line,errflag)
 619
 620       integer, parameter :: maxlen = 80
 621       character(len=80) :: line
 622       character(len=1) :: empty=' ',equal='='
 623       logical :: errflag
 624       integer :: ipos
 625 ! This function returns .TRUE., if an argument follows keyword keywd; if so
 626 ! IPOS will point to the first non-blank character of the argument. Returns
 627 ! .FALSE., if no argument follows the keyword; in this case IPOS points
 628 ! to the first non-blank character of the next keyword.
 629
 630       do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
 631         ipos=ipos+1
 632       enddo
 633       errflag=.false.
 634       if (line(ipos:ipos).eq.equal) then
 635          find_arg=.true.
 636          ipos=ipos+1
 637          do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
 638            ipos=ipos+1
 639          enddo
 640          if (ipos.gt.maxlen) errflag=.true.
 641       else
 642          find_arg=.false.
 643       endif
 644
 645       return
 646       end function find_arg
 647 !-----------------------------------------------------------------------------
 648       logical function find_group(iunit,jout,key1)
 649
 650       character*(*) :: key1
 651       character(len=80) :: karta        !,ucase
 652       integer :: iunit,jout
 653       integer :: ll !,ilen
 654 !EL      external ilen
 655 !EL      logical lcom
 656       rewind (iunit)
 657       karta=' '
 658       ll=ilen(key1)
 659       do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
 660         read (iunit,'(a)',end=10) karta
 661       enddo
 662       write (jout,'(2a)') '> ',karta(1:78)
 663       find_group=.true.
 664       return
 665    10 find_group=.false.
 666       return
 667       end function find_group
 668 !-----------------------------------------------------------------------------
 669       logical function iblnk(charc)
 670       character(len=1) :: charc
 671       integer :: n
 672       n = ichar(charc)
 673       iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
 674       return
 675       end function iblnk
 676 !-----------------------------------------------------------------------------
 677       integer function ilen(string)
 678       character*(*) ::  string
 679 !EL      logical :: iblnk
 680
 681       ilen = len(string)
 682 1     if ( ilen .gt. 0 ) then
 683          if ( iblnk( string(ilen:ilen) ) ) then
 684             ilen = ilen - 1
 685             goto 1
 686          endif
 687       endif
 688       return
 689       end function ilen
 690 !-----------------------------------------------------------------------------
 691       integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
 692       integer :: nkey,i,ikey,narg
 693       character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
 694 !      character(len=16) :: ucase
 695
 696       do i=1,narg
 697         if (ucase(keywd).eq.keywdset(i,ikey)) then
 698 ! Match found
 699           in_keywd_set=i
 700           return
 701         endif
 702       enddo
 703 ! No match to the allowed set of keywords if this point is reached.
 704       in_keywd_set=0
 705       return
 706       end function in_keywd_set
 707 !-----------------------------------------------------------------------------
 708       character function lcase(string)
 709       integer :: i, k, idiff
 710       character*(*) :: string
 711       character(len=1) :: c
 712       character(len=40) :: chtmp
 713 !
 714       i = len(lcase)
 715       k = len(string)
 716       if (i .lt. k) then
 717          k = i
 718          if (string(k+1:) .ne. ' ') then
 719             chtmp = string
 720          endif
 721       endif
 722       idiff = ichar('a') - ichar('A')
 723       lcase = string
 724       do 99 i = 1, k
 725          c = string(i:i)
 726          if (lge(c,'A') .and. lle(c,'Z')) then
 727             lcase(i:i) = char(ichar(c) + idiff)
 728          endif
 729    99 continue
 730       return
 731       end function lcase
 732 !-----------------------------------------------------------------------------
 733       logical function lcom(ipos,karta)
 734       character(len=80) :: karta
 735       character :: koment(2) = (/'!','#'/)
 736       integer :: ipos,i
 737
 738       lcom=.false.
 739       do i=1,2
 740         if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
 741       enddo
 742       return
 743       end function lcom
 744 !-----------------------------------------------------------------------------
 745       logical function lower_case(ch)
 746       character*(*) :: ch
 747       lower_case=(ch.ge.'a' .and. ch.le.'z')
 748       return
 749       end function lower_case
 750 !-----------------------------------------------------------------------------
 751       subroutine mykey(line,keywd,ipos,blankline,errflag)
 752
 753 ! This subroutine seeks a non-empty substring keywd in the string LINE.
 754 ! The substring begins with the first character different from blank and
 755 ! "=" encountered right to the pointer IPOS (inclusively) and terminates
 756 ! at the character left to the first blank or "=". When the subroutine is
 757 ! exited, the pointer IPOS is moved to the position of the terminator in LINE.
 758 ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
 759 ! only separators or the maximum length of the data line (80) has been reached.
 760 ! The logical variable ERRFLAG is set at .TRUE. if the string
 761 ! consists only from a "=".
 762       integer, parameter :: maxlen=80
 763       character(len=1) :: empty=' ',equal='=',comma=','
 764       character*(*) :: keywd
 765       character(len=80) :: line
 766       logical :: blankline,errflag      !EL,lcom
 767       integer :: ipos,istart,iend
 768       errflag=.false.
 769       do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
 770         ipos=ipos+1
 771       enddo
 772       if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
 773 ! At this point the rest of the input line turned out to contain only blanks
 774 ! or to be commented out.
 775         blankline=.true.
 776         return
 777       endif
 778       blankline=.false.
 779       istart=ipos
 780 ! Checks whether the current char is a separator.
 781       do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
 782        .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
 783         ipos=ipos+1
 784       enddo
 785       iend=ipos-1
 786 ! Error flag set to .true., if the length of the keyword was found less than 1.
 787       if (iend.lt.istart) then
 788         errflag=.true.
 789         return
 790       endif
 791       keywd=line(istart:iend)
 792       return
 793       end subroutine  mykey
 794 !-----------------------------------------------------------------------------
 795       subroutine numstr(inum,numm)
 796       character(len=10) :: huj='0123456789'
 797       character*(*) :: numm
 798       integer :: inumm,inum,inum1,inum2
 799       inumm=inum
 800       inum1=inumm/10
 801       inum2=inumm-10*inum1
 802       inumm=inum1
 803       numm(3:3)=huj(inum2+1:inum2+1)
 804       inum1=inumm/10
 805       inum2=inumm-10*inum1
 806       inumm=inum1
 807       numm(2:2)=huj(inum2+1:inum2+1)
 808       inum1=inumm/10
 809       inum2=inumm-10*inum1
 810       inumm=inum1
 811       numm(1:1)=huj(inum2+1:inum2+1)
 812       return
 813       end subroutine numstr
 814 !-----------------------------------------------------------------------------
 815       function ucase(string)
 816       integer :: i, k, idiff
 817       character(*) :: string
 818       character(len=len(string)) :: ucase
 819       character(len=1) :: c
 820       character(len=40) :: chtmp
 821 !
 822       i = len(ucase)
 823       k = len(string)
 824       if (i .lt. k) then
 825          k = i
 826          if (string(k+1:) .ne. ' ') then
 827             chtmp = string
 828          endif
 829       endif
 830       idiff = ichar('a') - ichar('A')
 831       ucase = string
 832       do 99 i = 1, k
 833          c = string(i:i)
 834          if (lge(c,'a') .and. lle(c,'z')) then
 835             ucase(i:i) = char(ichar(c) - idiff)
 836          endif
 837    99 continue
 838       return
 839       end function ucase
 840 !-----------------------------------------------------------------------------
 841 ! geomout.F
 842 !-----------------------------------------------------------------------------
 843       subroutine pdbout(etot,tytul,iunit)
 844
 845       use geometry_data, only: c,nres
 846       use energy_data
 847    !  use control
 848       use compare_data
 849       use MD_data
 850 !      implicit real*8 (a-h,o-z)
 851 !      include 'DIMENSIONS'
 852 !      include 'COMMON.CHAIN'
 853 !      include 'COMMON.INTERACT'
 854 !      include 'COMMON.NAMES'
 855 !      include 'COMMON.IOUNITS'
 856 !      include 'COMMON.HEADER'
 857 !      include 'COMMON.SBRIDGE'
 858 !      include 'COMMON.DISTFIT'
 859 !      include 'COMMON.MD'
 860 !el      character(len=50) :: tytul
 861       character*(*) :: tytul
 862       character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
 863       integer,dimension(nres) :: ica    !(maxres)
 864        integer iti1
 865 !el  local variables
 866       integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
 867       real(kind=8) :: etot
 868       integer :: nres2
 869       nres2=2*nres
 870
 871       if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
 872
 873       write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
 874 !model      write (iunit,'(a5,i6)') 'MODEL',1
 875       if (nhfrag.gt.0) then
 876        do j=1,nhfrag
 877         iti=itype(hfrag(1,j))
 878         itj=itype(hfrag(2,j))
 879         if (j.lt.10) then
 880            write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
 881                  'HELIX',j,'H',j,&
 882                  restyp(iti),hfrag(1,j)-1,&
 883                  restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
 884         else
 885              write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
 886                  'HELIX',j,'H',j,&
 887                  restyp(iti),hfrag(1,j)-1,&
 888                  restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
 889         endif
 890        enddo
 891       endif
 892
 893       if (nbfrag.gt.0) then
 894
 895        do j=1,nbfrag
 896
 897         iti=itype(bfrag(1,j))
 898         itj=itype(bfrag(2,j)-1)
 899
 900         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
 901                  'SHEET',1,'B',j,2,&
 902                  restyp(iti),bfrag(1,j)-1,&
 903                  restyp(itj),bfrag(2,j)-2,0
 904
 905         if (bfrag(3,j).gt.bfrag(4,j)) then
 906
 907          itk=itype(bfrag(3,j))
 908          itl=itype(bfrag(4,j)+1)
 909
 910          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)') &
 911                  'SHEET',2,'B',j,2,&
 912                  restyp(itl),bfrag(4,j),&
 913                  restyp(itk),bfrag(3,j)-1,-1,&
 914                  "N",restyp(itk),bfrag(3,j)-1,&
 915                  "O",restyp(iti),bfrag(1,j)-1
 916
 917         else
 918
 919          itk=itype(bfrag(3,j))
 920          itl=itype(bfrag(4,j)-1)
 921
 922
 923         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)') &
 924                  'SHEET',2,'B',j,2,&
 925                  restyp(itk),bfrag(3,j)-1,&
 926                  restyp(itl),bfrag(4,j)-2,1,&
 927                  "N",restyp(itk),bfrag(3,j)-1,&
 928                  "O",restyp(iti),bfrag(1,j)-1
 929
 930
 931
 932         endif
 933
 934        enddo
 935       endif
 936
 937       if (nss.gt.0) then
 938         do i=1,nss
 939          if (dyn_ss) then
 940           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
 941                'SSBOND',i,'CYS',idssb(i)-nnt+1,&
 942                           'CYS',jdssb(i)-nnt+1
 943          else
 944           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
 945                'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
 946                           'CYS',jhpb(i)-nnt+1-nres
 947          endif
 948         enddo
 949       endif
 950
 951       iatom=0
 952       ichain=1
 953       ires=0
 954       do i=nnt,nct
 955         iti=itype(i)
 956         iti1=itype(i+1)
 957         if ((iti.eq.ntyp1).and.(iti1.eq.ntyp1)) cycle
 958         if (iti.eq.ntyp1) then
 959           ichain=ichain+1
 960 !          ires=0
 961           write (iunit,'(a)') 'TER'
 962         else
 963         ires=ires+1
 964         iatom=iatom+1
 965         ica(i)=iatom
 966         write (iunit,10) iatom,restyp(iti),chainid(ichain),&
 967            ires,(c(j,i),j=1,3),vtot(i)
 968         if (iti.ne.10) then
 969           iatom=iatom+1
 970           write (iunit,20) iatom,restyp(iti),chainid(ichain),&
 971             ires,(c(j,nres+i),j=1,3),&
 972             vtot(i+nres)
 973         endif
 974         endif
 975       enddo
 976       write (iunit,'(a)') 'TER'
 977       do i=nnt,nct-1
 978         if (itype(i).eq.ntyp1) cycle
 979         if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
 980           write (iunit,30) ica(i),ica(i+1)
 981         else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
 982           write (iunit,30) ica(i),ica(i+1),ica(i)+1
 983         else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
 984           write (iunit,30) ica(i),ica(i)+1
 985         endif
 986       enddo
 987       if (itype(nct).ne.10) then
 988         write (iunit,30) ica(nct),ica(nct)+1
 989       endif
 990       do i=1,nss
 991        if (dyn_ss) then
 992         write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
 993        else
 994         write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
 995        endif
 996       enddo
 997       write (iunit,'(a6)') 'ENDMDL'
 998   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
 999   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1000   30  FORMAT ('CONECT',8I5)
1001       return
1002       end subroutine pdbout
1003 !-----------------------------------------------------------------------------
1004       subroutine MOL2out(etot,tytul)
1005 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1006 ! format.
1007       use geometry_data, only: c
1008       use energy_data
1009    !  use control
1010 !      implicit real*8 (a-h,o-z)
1011 !      include 'DIMENSIONS'
1012 !      include 'COMMON.CHAIN'
1013 !      include 'COMMON.INTERACT'
1014 !      include 'COMMON.NAMES'
1015 !      include 'COMMON.IOUNITS'
1016 !      include 'COMMON.HEADER'
1017 !      include 'COMMON.SBRIDGE'
1018       character(len=32) :: tytul,fd
1019       character(len=3) :: zahl
1020       character(len=6) :: res_num,pom   !,ucase
1021
1022 !el  local variables
1023       integer :: i,j
1024       real(kind=8) :: etot
1025
1026 #ifdef AIX
1027       call fdate_(fd)
1028 #elif (defined CRAY)
1029       call date(fd)
1030 #else
1031       call fdate(fd)
1032 #endif
1033       write (imol2,'(a)') '#'
1034       write (imol2,'(a)') &
1035        '#         Creating user name:           unres'
1036       write (imol2,'(2a)') '#         Creation time:                ',&
1037        fd
1038       write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1039       write (imol2,'(a)') tytul
1040       write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1041       write (imol2,'(a)') 'SMALL'
1042       write (imol2,'(a)') 'USER_CHARGES'
1043       write (imol2,'(a)') '\@<TRIPOS>ATOM'
1044       do i=nnt,nct
1045         write (zahl,'(i3)') i
1046         pom=ucase(restyp(itype(i)))
1047         res_num = pom(:3)//zahl(2:)
1048         write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1049       enddo
1050       write (imol2,'(a)') '\@<TRIPOS>BOND'
1051       do i=nnt,nct-1
1052         write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1053       enddo
1054       do i=1,nss
1055         write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1056       enddo
1057       write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1058       do i=nnt,nct
1059         write (zahl,'(i3)') i
1060         pom = ucase(restyp(itype(i)))
1061         res_num = pom(:3)//zahl(2:)
1062         write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1063       enddo
1064   10  FORMAT (I7,' CA      ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1065   30  FORMAT (I7,1x,A,I14,' RESIDUE',I13,' ****  ****')
1066       return
1067       end subroutine MOL2out
1068 !-----------------------------------------------------------------------------
1069       subroutine intout
1070
1071       use geometry_data
1072       use energy_data, only: itype
1073    !  use control
1074 !      implicit real*8 (a-h,o-z)
1075 !      include 'DIMENSIONS'
1076 !      include 'COMMON.IOUNITS'
1077 !      include 'COMMON.CHAIN'
1078 !      include 'COMMON.VAR'
1079 !      include 'COMMON.LOCAL'
1080 !      include 'COMMON.INTERACT'
1081 !      include 'COMMON.NAMES'
1082 !      include 'COMMON.GEO'
1083 !      include 'COMMON.TORSION'
1084 !el  local variables
1085       integer :: i,iti
1086
1087       write (iout,'(/a)') 'Geometry of the virtual chain.'
1088       write (iout,'(7a)') '  Res  ','         d','     Theta',&
1089        '       Phi','       Dsc','     Alpha','      Omega'
1090       do i=1,nres
1091         iti=itype(i)
1092         write (iout,'(a3,i4,6f10.3)') restyp(iti),i,vbld(i),&
1093            rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1094            rad2deg*omeg(i)
1095       enddo
1096       return
1097       end subroutine intout
1098 !-----------------------------------------------------------------------------
1099 #ifdef CLUSTER
1100       subroutine briefout(it,ener,free)!,plik)
1101 #else
1102       subroutine briefout(it,ener)
1103 #endif
1104       use geometry_data
1105       use energy_data
1106    !  use control
1107 !      implicit real*8 (a-h,o-z)
1108 !      include 'DIMENSIONS'
1109 !      include 'COMMON.IOUNITS'
1110 !      include 'COMMON.CHAIN'
1111 !      include 'COMMON.VAR'
1112 !      include 'COMMON.LOCAL'
1113 !      include 'COMMON.INTERACT'
1114 !      include 'COMMON.NAMES'
1115 !      include 'COMMON.GEO'
1116 !      include 'COMMON.SBRIDGE'
1117 !     print '(a,i5)',intname,igeom
1118 !el  local variables
1119       integer :: i,it
1120       real(kind=8) :: ener,free
1121 !     character(len=80) :: plik
1122 !      integer :: iii
1123
1124 #if defined(AIX) || defined(PGI)
1125       open (igeom,file=intname,position='append')
1126 #else
1127       open (igeom,file=intname,access='append')
1128 #endif
1129 #ifdef WHAM_RUN
1130 !      iii=igeom
1131       igeom=iout
1132 #endif
1133       IF (NSS.LE.9) THEN
1134 #ifdef CLUSTER
1135         WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1136       ELSE
1137         WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1138 #else
1139         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1140       ELSE
1141         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1142 #endif
1143         WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1144       ENDIF
1145 !     IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1146       WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1147       WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1148 !     if (nvar.gt.nphi+ntheta) then
1149         write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1150         write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1151 !     endif
1152       close(igeom)
1153   180 format (I5,F12.3,I2,9(1X,2I3))
1154   190 format (3X,11(1X,2I3))
1155   200 format (8F10.4)
1156       return
1157       end subroutine briefout
1158 !-----------------------------------------------------------------------------
1159 #ifdef WINIFL
1160       subroutine fdate(fd)
1161       character(len=32) :: fd
1162       write(fd,'(32x)')
1163       return
1164       end subroutine fdate
1165 #endif
1166 !-----------------------------------------------------------------------------
1167 #ifdef WHAM_RUN
1168       real(kind=8) function gyrate(jcon)
1169 #else
1170       real(kind=8) function gyrate()
1171 #endif
1172
1173       use geometry_data, only: c
1174    !  use geometry
1175       use energy_data
1176 !      implicit real*8 (a-h,o-z)
1177 !      include 'DIMENSIONS'
1178 !      include 'COMMON.INTERACT'
1179 !      include 'COMMON.CHAIN'
1180       real(kind=8) :: rg
1181       real(kind=8),dimension(3) :: cen
1182 !el  local variables
1183       integer :: i,j,jcon
1184
1185       do j=1,3
1186        cen(j)=0.0d0
1187       enddo
1188
1189       do i=nnt,nct
1190           do j=1,3
1191             cen(j)=cen(j)+c(j,i)
1192           enddo
1193       enddo
1194       do j=1,3
1195             cen(j)=cen(j)/dble(nct-nnt+1)
1196       enddo
1197       rg = 0.0d0
1198       do i = nnt, nct
1199         do j=1,3
1200          rg = rg + (c(j,i)-cen(j))**2
1201         enddo
1202       end do
1203 #ifdef WHAM_RUN
1204       gyrate = dsqrt(rg/dble(nct-nnt+1))
1205 #else
1206       gyrate = sqrt(rg/dble(nct-nnt+1))
1207 #endif
1208       return
1209       end function gyrate
1210 #ifdef WHAM_RUN
1211 !-----------------------------------------------------------------------------
1212 ! readrtns.F WHAM
1213       subroutine reads(rekord,lancuch,wartosc,default)
1214 !      implicit none
1215       character*(*) :: rekord,lancuch,wartosc,default
1216       character(len=80) :: aux
1217       integer :: lenlan,lenrec,iread,ireade
1218 !el      external ilen
1219 !el      logical iblnk
1220 !el      external iblnk
1221       lenlan=ilen(lancuch)
1222       lenrec=ilen(rekord)
1223       iread=index(rekord,lancuch(:lenlan)//"=")
1224 !      print *,"rekord",rekord," lancuch",lancuch
1225 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1226       if (iread.eq.0) then
1227         wartosc=default
1228         return
1229       endif
1230       iread=iread+lenlan+1
1231 !      print *,"iread",iread
1232 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1233       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1234         iread=iread+1
1235 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1236       enddo
1237 !      print *,"iread",iread
1238       if (iread.gt.lenrec) then
1239          wartosc=default
1240         return
1241       endif
1242       ireade=iread+1
1243 !      print *,"ireade",ireade
1244       do while (ireade.lt.lenrec .and. &
1245          .not.iblnk(rekord(ireade:ireade)))
1246         ireade=ireade+1
1247       enddo
1248       wartosc=rekord(iread:ireade)
1249       return
1250       end subroutine reads
1251 #endif
1252 !-----------------------------------------------------------------------------
1253 ! permut.F
1254 !-----------------------------------------------------------------------------
1255       subroutine permut(isym)
1256
1257       use geometry_data, only: tabperm
1258 !      implicit real*8 (a-h,o-z)
1259 !      include 'DIMENSIONS'
1260 !      include 'COMMON.LOCAL'
1261 !      include 'COMMON.VAR'
1262 !      include 'COMMON.CHAIN'
1263 !      include 'COMMON.INTERACT'
1264 !      include 'COMMON.IOUNITS'
1265 !      include 'COMMON.GEO'
1266 !      include 'COMMON.NAMES'
1267 !      include 'COMMON.CONTROL'
1268
1269       integer :: n,isym
1270 !      logical nextp
1271 !el      external nextp
1272       integer,dimension(isym) :: a
1273 !      parameter(n=symetr)
1274 !el local variables
1275       integer :: kkk,i
1276
1277       n=isym
1278       if (n.eq.1) then
1279         tabperm(1,1)=1
1280         return
1281       endif
1282       kkk=0
1283       do i=1,n
1284       a(i)=i
1285       enddo
1286    10 print *,(a(i),i=1,n)
1287       kkk=kkk+1
1288       do i=1,n
1289       tabperm(kkk,i)=a(i)
1290 !      write (iout,*) "tututu", kkk
1291       enddo
1292       if(nextp(n,a)) go to 10
1293       return
1294       end subroutine permut
1295 !-----------------------------------------------------------------------------
1296       logical function nextp(n,a)
1297
1298       integer :: n,i,j,k,t
1299 !      logical :: nextp
1300       integer,dimension(n) :: a
1301       i=n-1
1302    10 if(a(i).lt.a(i+1)) go to 20
1303       i=i-1
1304       if(i.eq.0) go to 20
1305       go to 10
1306    20 j=i+1
1307       k=n
1308    30 t=a(j)
1309       a(j)=a(k)
1310       a(k)=t
1311       j=j+1
1312       k=k-1
1313       if(j.lt.k) go to 30
1314       j=i
1315       if(j.ne.0) go to 40
1316       nextp=.false.
1317       return
1318    40 j=j+1
1319       if(a(j).lt.a(i)) go to 40
1320       t=a(i)
1321       a(i)=a(j)
1322       a(j)=t
1323       nextp=.true.
1324       return
1325       end function nextp
1326 !-----------------------------------------------------------------------------
1327 !-----------------------------------------------------------------------------
1328       end module io_base