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