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