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
 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         if (iti.eq.ntyp1) then
 957           ichain=ichain+1
 958           ires=0
 959           write (iunit,'(a)') 'TER'
 960         else
 961         ires=ires+1
 962         iatom=iatom+1
 963         ica(i)=iatom
 964         write (iunit,10) iatom,restyp(iti),chainid(ichain),&
 965            ires,(c(j,i),j=1,3),vtot(i)
 966         if (iti.ne.10) then
 967           iatom=iatom+1
 968           write (iunit,20) iatom,restyp(iti),chainid(ichain),&
 969             ires,(c(j,nres+i),j=1,3),&
 970             vtot(i+nres)
 971         endif
 972         endif
 973       enddo
 974       write (iunit,'(a)') 'TER'
 975       do i=nnt,nct-1
 976         if (itype(i).eq.ntyp1) cycle
 977         if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
 978           write (iunit,30) ica(i),ica(i+1)
 979         else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
 980           write (iunit,30) ica(i),ica(i+1),ica(i)+1
 981         else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
 982           write (iunit,30) ica(i),ica(i)+1
 983         endif
 984       enddo
 985       if (itype(nct).ne.10) then
 986         write (iunit,30) ica(nct),ica(nct)+1
 987       endif
 988       do i=1,nss
 989        if (dyn_ss) then
 990         write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
 991        else
 992         write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
 993        endif
 994       enddo
 995       write (iunit,'(a6)') 'ENDMDL'
 996   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
 997   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
 998   30  FORMAT ('CONECT',8I5)
 999       return
1000       end subroutine pdbout
1001 !-----------------------------------------------------------------------------
1002       subroutine MOL2out(etot,tytul)
1003 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1004 ! format.
1005       use geometry_data, only: c
1006       use energy_data
1007    !  use control
1008 !      implicit real*8 (a-h,o-z)
1009 !      include 'DIMENSIONS'
1010 !      include 'COMMON.CHAIN'
1011 !      include 'COMMON.INTERACT'
1012 !      include 'COMMON.NAMES'
1013 !      include 'COMMON.IOUNITS'
1014 !      include 'COMMON.HEADER'
1015 !      include 'COMMON.SBRIDGE'
1016       character(len=32) :: tytul,fd
1017       character(len=3) :: zahl
1018       character(len=6) :: res_num,pom   !,ucase
1019
1020 !el  local variables
1021       integer :: i,j
1022       real(kind=8) :: etot
1023
1024 #ifdef AIX
1025       call fdate_(fd)
1026 #elif (defined CRAY)
1027       call date(fd)
1028 #else
1029       call fdate(fd)
1030 #endif
1031       write (imol2,'(a)') '#'
1032       write (imol2,'(a)') &
1033        '#         Creating user name:           unres'
1034       write (imol2,'(2a)') '#         Creation time:                ',&
1035        fd
1036       write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1037       write (imol2,'(a)') tytul
1038       write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1039       write (imol2,'(a)') 'SMALL'
1040       write (imol2,'(a)') 'USER_CHARGES'
1041       write (imol2,'(a)') '\@<TRIPOS>ATOM'
1042       do i=nnt,nct
1043         write (zahl,'(i3)') i
1044         pom=ucase(restyp(itype(i)))
1045         res_num = pom(:3)//zahl(2:)
1046         write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1047       enddo
1048       write (imol2,'(a)') '\@<TRIPOS>BOND'
1049       do i=nnt,nct-1
1050         write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1051       enddo
1052       do i=1,nss
1053         write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1054       enddo
1055       write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1056       do i=nnt,nct
1057         write (zahl,'(i3)') i
1058         pom = ucase(restyp(itype(i)))
1059         res_num = pom(:3)//zahl(2:)
1060         write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1061       enddo
1062   10  FORMAT (I7,' CA      ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1063   30  FORMAT (I7,1x,A,I14,' RESIDUE',I13,' ****  ****')
1064       return
1065       end subroutine MOL2out
1066 !-----------------------------------------------------------------------------
1067       subroutine intout
1068
1069       use geometry_data
1070       use energy_data, only: itype
1071    !  use control
1072 !      implicit real*8 (a-h,o-z)
1073 !      include 'DIMENSIONS'
1074 !      include 'COMMON.IOUNITS'
1075 !      include 'COMMON.CHAIN'
1076 !      include 'COMMON.VAR'
1077 !      include 'COMMON.LOCAL'
1078 !      include 'COMMON.INTERACT'
1079 !      include 'COMMON.NAMES'
1080 !      include 'COMMON.GEO'
1081 !      include 'COMMON.TORSION'
1082 !el  local variables
1083       integer :: i,iti
1084
1085       write (iout,'(/a)') 'Geometry of the virtual chain.'
1086       write (iout,'(7a)') '  Res  ','         d','     Theta',&
1087        '       Phi','       Dsc','     Alpha','      Omega'
1088       do i=1,nres
1089         iti=itype(i)
1090         write (iout,'(a3,i4,6f10.3)') restyp(iti),i,vbld(i),&
1091            rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1092            rad2deg*omeg(i)
1093       enddo
1094       return
1095       end subroutine intout
1096 !-----------------------------------------------------------------------------
1097 #ifdef CLUSTER
1098       subroutine briefout(it,ener,free)!,plik)
1099 #else
1100       subroutine briefout(it,ener)
1101 #endif
1102       use geometry_data
1103       use energy_data
1104    !  use control
1105 !      implicit real*8 (a-h,o-z)
1106 !      include 'DIMENSIONS'
1107 !      include 'COMMON.IOUNITS'
1108 !      include 'COMMON.CHAIN'
1109 !      include 'COMMON.VAR'
1110 !      include 'COMMON.LOCAL'
1111 !      include 'COMMON.INTERACT'
1112 !      include 'COMMON.NAMES'
1113 !      include 'COMMON.GEO'
1114 !      include 'COMMON.SBRIDGE'
1115 !     print '(a,i5)',intname,igeom
1116 !el  local variables
1117       integer :: i,it
1118       real(kind=8) :: ener,free
1119 !     character(len=80) :: plik
1120 !      integer :: iii
1121
1122 #if defined(AIX) || defined(PGI)
1123       open (igeom,file=intname,position='append')
1124 #else
1125       open (igeom,file=intname,access='append')
1126 #endif
1127 #ifdef WHAM_RUN
1128 !      iii=igeom
1129       igeom=iout
1130 #endif
1131       IF (NSS.LE.9) THEN
1132 #ifdef CLUSTER
1133         WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1134       ELSE
1135         WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1136 #else
1137         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1138       ELSE
1139         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1140 #endif
1141         WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1142       ENDIF
1143 !     IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1144       WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1145       WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1146 !     if (nvar.gt.nphi+ntheta) then
1147         write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1148         write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1149 !     endif
1150       close(igeom)
1151   180 format (I5,F12.3,I2,9(1X,2I3))
1152   190 format (3X,11(1X,2I3))
1153   200 format (8F10.4)
1154       return
1155       end subroutine briefout
1156 !-----------------------------------------------------------------------------
1157 #ifdef WINIFL
1158       subroutine fdate(fd)
1159       character(len=32) :: fd
1160       write(fd,'(32x)')
1161       return
1162       end subroutine fdate
1163 #endif
1164 !-----------------------------------------------------------------------------
1165 #ifdef WHAM_RUN
1166       real(kind=8) function gyrate(jcon)
1167 #else
1168       real(kind=8) function gyrate()
1169 #endif
1170
1171       use geometry_data, only: c
1172    !  use geometry
1173       use energy_data
1174 !      implicit real*8 (a-h,o-z)
1175 !      include 'DIMENSIONS'
1176 !      include 'COMMON.INTERACT'
1177 !      include 'COMMON.CHAIN'
1178       real(kind=8) :: rg
1179       real(kind=8),dimension(3) :: cen
1180 !el  local variables
1181       integer :: i,j,jcon
1182
1183       do j=1,3
1184        cen(j)=0.0d0
1185       enddo
1186
1187       do i=nnt,nct
1188           do j=1,3
1189             cen(j)=cen(j)+c(j,i)
1190           enddo
1191       enddo
1192       do j=1,3
1193             cen(j)=cen(j)/dble(nct-nnt+1)
1194       enddo
1195       rg = 0.0d0
1196       do i = nnt, nct
1197         do j=1,3
1198          rg = rg + (c(j,i)-cen(j))**2
1199         enddo
1200       end do
1201 #ifdef WHAM_RUN
1202       gyrate = dsqrt(rg/dble(nct-nnt+1))
1203 #else
1204       gyrate = sqrt(rg/dble(nct-nnt+1))
1205 #endif
1206       return
1207       end function gyrate
1208 #ifdef WHAM_RUN
1209 !-----------------------------------------------------------------------------
1210 ! readrtns.F WHAM
1211       subroutine reads(rekord,lancuch,wartosc,default)
1212 !      implicit none
1213       character*(*) :: rekord,lancuch,wartosc,default
1214       character(len=80) :: aux
1215       integer :: lenlan,lenrec,iread,ireade
1216 !el      external ilen
1217 !el      logical iblnk
1218 !el      external iblnk
1219       lenlan=ilen(lancuch)
1220       lenrec=ilen(rekord)
1221       iread=index(rekord,lancuch(:lenlan)//"=")
1222 !      print *,"rekord",rekord," lancuch",lancuch
1223 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1224       if (iread.eq.0) then
1225         wartosc=default
1226         return
1227       endif
1228       iread=iread+lenlan+1
1229 !      print *,"iread",iread
1230 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1231       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1232         iread=iread+1
1233 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1234       enddo
1235 !      print *,"iread",iread
1236       if (iread.gt.lenrec) then
1237          wartosc=default
1238         return
1239       endif
1240       ireade=iread+1
1241 !      print *,"ireade",ireade
1242       do while (ireade.lt.lenrec .and. &
1243          .not.iblnk(rekord(ireade:ireade)))
1244         ireade=ireade+1
1245       enddo
1246       wartosc=rekord(iread:ireade)
1247       return
1248       end subroutine reads
1249 #endif
1250 !-----------------------------------------------------------------------------
1251 ! permut.F
1252 !-----------------------------------------------------------------------------
1253       subroutine permut(isym)
1254
1255       use geometry_data, only: tabperm
1256 !      implicit real*8 (a-h,o-z)
1257 !      include 'DIMENSIONS'
1258 !      include 'COMMON.LOCAL'
1259 !      include 'COMMON.VAR'
1260 !      include 'COMMON.CHAIN'
1261 !      include 'COMMON.INTERACT'
1262 !      include 'COMMON.IOUNITS'
1263 !      include 'COMMON.GEO'
1264 !      include 'COMMON.NAMES'
1265 !      include 'COMMON.CONTROL'
1266
1267       integer :: n,isym
1268 !      logical nextp
1269 !el      external nextp
1270       integer,dimension(isym) :: a
1271 !      parameter(n=symetr)
1272 !el local variables
1273       integer :: kkk,i
1274
1275       n=isym
1276       if (n.eq.1) then
1277         tabperm(1,1)=1
1278         return
1279       endif
1280       kkk=0
1281       do i=1,n
1282       a(i)=i
1283       enddo
1284    10 print *,(a(i),i=1,n)
1285       kkk=kkk+1
1286       do i=1,n
1287       tabperm(kkk,i)=a(i)
1288 !      write (iout,*) "tututu", kkk
1289       enddo
1290       if(nextp(n,a)) go to 10
1291       return
1292       end subroutine permut
1293 !-----------------------------------------------------------------------------
1294       logical function nextp(n,a)
1295
1296       integer :: n,i,j,k,t
1297 !      logical :: nextp
1298       integer,dimension(n) :: a
1299       i=n-1
1300    10 if(a(i).lt.a(i+1)) go to 20
1301       i=i-1
1302       if(i.eq.0) go to 20
1303       go to 10
1304    20 j=i+1
1305       k=n
1306    30 t=a(j)
1307       a(j)=a(k)
1308       a(k)=t
1309       j=j+1
1310       k=k-1
1311       if(j.lt.k) go to 30
1312       j=i
1313       if(j.ne.0) go to 40
1314       nextp=.false.
1315       return
1316    40 j=j+1
1317       if(a(j).lt.a(i)) go to 40
1318       t=a(i)
1319       a(i)=a(j)
1320       a(j)=t
1321       nextp=.true.
1322       return
1323       end function nextp
1324 !-----------------------------------------------------------------------------
1325 !-----------------------------------------------------------------------------
1326       end module io_base