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