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=101000!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       real, dimension(3) ::cbeg,cend
 439       integer i,j
 440       print *, "wchodze"
 441       call card_concat(afmcard,.true.)
 442       call readi(afmcard,"BEG",afmbeg,0)
 443       call readi(afmcard,"END",afmend,0)
 444       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
 445       call reada(afmcard,"VEL",velAFMconst,0.0d0)
 446       print *,'FORCE=' ,forceAFMconst
 447       cbeg=0.0d0
 448       cend=0.0d0
 449       if (afmbeg.eq.-1) then
 450         read(inp,*) nbegafmmat,(afmbegcentr(i),i=1,nbegafmmat)
 451         do i=1,nbegafmmat
 452          do j=1,3
 453           cbeg(j)=cbeg(j)+c(j,afmbegcentr(i))/nbegafmmat
 454          enddo
 455         enddo
 456       else
 457       do j=1,3
 458         cbeg(j)=c(j,afmend)
 459       enddo
 460       endif
 461       if (afmend.eq.-1) then
 462         read(inp,*) nendafmmat,(afmendcentr(i),i=1,nendafmmat)
 463         do i=1,nendafmmat
 464          do j=1,3
 465           cend(j)=cend(j)+c(j,afmendcentr(i))/nendafmmat
 466          enddo
 467         enddo
 468       else
 469         cend(j)=c(j,afmend)
 470       endif
 471 !------ NOW PROPERTIES FOR AFM
 472        distafminit=0.0d0
 473        do i=1,3
 474         distafminit=(cend(i)-cbeg(i))**2+distafminit
 475        enddo
 476         distafminit=dsqrt(distafminit)
 477         print *,'initdist',distafminit
 478       return
 479       end subroutine read_afminp
 480 !C------------------------------------------------------------
 481       subroutine read_afmnano
 482       use geometry_data
 483       use energy_data
 484       use control_data, only:out1file
 485       use MPI_data
 486       integer :: i
 487       real*8 :: cm
 488       read(inp,*,err=11) inanomove,(inanotab(i),i=1,inanomove),forcenanoconst
 489        cm=0.0d0
 490        do i=1,inanomove
 491          cm=cm+c(3,inanotab(i))
 492        enddo
 493        cm=cm/inanomove
 494         distnanoinit=cm-tubecenter(3)
 495       return
 496 11    write(iout,*)&
 497       "error in afmnano",&
 498        ", number of center, their index and forceconstance"
 499       stop
 500       return
 501       end subroutine read_afmnano
 502
 503 !-----------------------------------------------------------------------------
 504       subroutine read_dist_constr
 505       use MPI_data
 506 !     use control
 507       use geometry, only: dist
 508       use geometry_data
 509       use control_data
 510       use energy_data
 511 !      implicit real*8 (a-h,o-z)
 512 !      include 'DIMENSIONS'
 513 #ifdef MPI
 514       include 'mpif.h'
 515 #endif
 516 !      include 'COMMON.SETUP'
 517 !      include 'COMMON.CONTROL'
 518 !      include 'COMMON.CHAIN'
 519 !      include 'COMMON.IOUNITS'
 520 !      include 'COMMON.SBRIDGE'
 521       integer,dimension(2,100) :: ifrag_,ipair_
 522       real(kind=8),dimension(100) :: wfrag_,wpair_
 523       character(len=640) :: controlcard
 524
 525 !el local variables
 526       integer :: i,k,j,ii,jj,itemp,mnumkk,mnumjj,k1,j1
 527       integer :: nfrag_,npair_,ndist_
 528       real(kind=8) :: dist_cut,ddjk
 529
 530 !      write (iout,*) "Calling read_dist_constr"
 531 !      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 532 !      call flush(iout)
 533       if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
 534       if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
 535       if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
 536       if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim))
 537       if(.not.allocated(forcon)) allocate(forcon(maxdim))
 538       if(.not.allocated(fordepth)) allocate(fordepth(maxdim))
 539       if(.not.allocated(ibecarb)) allocate(ibecarb(maxdim))
 540       if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
 541       call gen_dist_constr2
 542       go to 1712
 543       endif
 544       call card_concat(controlcard,.true.)
 545       call readi(controlcard,"NFRAG",nfrag_,0)
 546       call readi(controlcard,"NPAIR",npair_,0)
 547       call readi(controlcard,"NDIST",ndist_,0)
 548       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 549       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 550       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 551       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 552       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 553       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 554 !      write (iout,*) "IFRAG"
 555       do i=1,nfrag_
 556         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 557       enddo
 558 !      write (iout,*) "IPAIR"
 559 !      do i=1,npair_
 560 !        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 561 !      enddo
 562 !      if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
 563 !      if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
 564 !      if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
 565 !      if(.not.allocated(forcon)) allocate(forcon(maxdim))
 566
 567       call flush(iout)
 568       do i=1,nfrag_
 569 !        if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 570 !        if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
 571 !          ifrag_(2,i)=nstart_sup+nsup-1
 572 !        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 573         call flush(iout)
 574         if (wfrag_(i).gt.0.0d0) then
 575         do j=ifrag_(1,i),ifrag_(2,i)-1
 576           do k=j+1,ifrag_(2,i)
 577           write (iout,*) "j",j," k",k,nres
 578           j1=j
 579           k1=k
 580           if (j.gt.nres) j1=j-nres
 581           if (k.gt.nres) k1=k-nres
 582           mnumkk=molnum(k1)
 583           mnumjj=molnum(j1)
 584
 585           if ((itype(k1,mnumkk).gt.ntyp_molec(mnumkk)).or.&
 586               (itype(j1,mnumjj).gt.ntyp_molec(mnumjj))) cycle
 587             write (iout,*) "j",j," k",k,itype(k1,mnumkk),itype(j1,mnumjj)
 588             ddjk=dist(j,k)
 589             if (constr_dist.eq.1) then
 590               nhpb=nhpb+1
 591               ihpb(nhpb)=j
 592               jhpb(nhpb)=k
 593               dhpb(nhpb)=ddjk
 594               forcon(nhpb)=wfrag_(i)
 595             else if (constr_dist.eq.2) then
 596               if (ddjk.le.dist_cut) then
 597                 nhpb=nhpb+1
 598                 ihpb(nhpb)=j
 599                 jhpb(nhpb)=k
 600                 dhpb(nhpb)=ddjk
 601                 forcon(nhpb)=wfrag_(i)
 602               endif
 603             else
 604               nhpb=nhpb+1
 605               ihpb(nhpb)=j
 606               jhpb(nhpb)=k
 607               dhpb(nhpb)=ddjk
 608               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 609             endif
 610 #ifdef MPI
 611             if (.not.out1file .or. me.eq.king) &
 612             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
 613              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 614 #else
 615             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
 616              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 617 #endif
 618           enddo
 619         enddo
 620         endif
 621       enddo
 622       do i=1,npair_
 623         if (wpair_(i).gt.0.0d0) then
 624         ii = ipair_(1,i)
 625         jj = ipair_(2,i)
 626         if (ii.gt.jj) then
 627           itemp=ii
 628           ii=jj
 629           jj=itemp
 630         endif
 631         do j=ifrag_(1,ii),ifrag_(2,ii)
 632           do k=ifrag_(1,jj),ifrag_(2,jj)
 633             nhpb=nhpb+1
 634             ihpb(nhpb)=j
 635             jhpb(nhpb)=k
 636             forcon(nhpb)=wpair_(i)
 637             dhpb(nhpb)=dist(j,k)
 638 #ifdef MPI
 639             if (.not.out1file .or. me.eq.king) &
 640             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 641              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 642 #else
 643             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 644              nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 645 #endif
 646           enddo
 647         enddo
 648         endif
 649       enddo
 650       do i=1,ndist_
 651 !        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
 652 !        if (forcon(nhpb+1).gt.0.0d0) then
 653 !          nhpb=nhpb+1
 654 !          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 655         if (constr_dist.eq.11) then
 656         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
 657           ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 658 !EL        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 659           fordepth(nhpb+1)=fordepth(nhpb+1)**(0.25d0)
 660           forcon(nhpb+1)=forcon(nhpb+1)**(0.25d0)
 661         else
 662 !C        print *,"in else"
 663         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
 664           ibecarb(i),forcon(nhpb+1)
 665         endif
 666         if (forcon(nhpb+1).gt.0.0d0) then
 667           nhpb=nhpb+1
 668           if (ibecarb(i).gt.0) then
 669             ihpb(i)=ihpb(i)+nres
 670             jhpb(i)=jhpb(i)+nres
 671           endif
 672           if (dhpb(nhpb).eq.0.0d0) &
 673             dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 674         endif
 675
 676 #ifdef MPI
 677           if (.not.out1file .or. me.eq.king) &
 678           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 679            nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 680 #else
 681           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
 682            nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 683 #endif
 684       enddo
 685  1712 continue
 686       call flush(iout)
 687       return
 688       end subroutine read_dist_constr
 689 !-----------------------------------------------------------------------------
 690       subroutine gen_dist_constr2
 691       use MPI_data
 692 !     use control
 693       use geometry, only: dist
 694       use geometry_data
 695       use control_data
 696       use energy_data
 697       integer :: i,j
 698       real(kind=8) :: distance
 699       if (constr_dist.eq.11) then
 700              do i=nstart_sup,nstart_sup+nsup-1
 701               do j=i+2,nstart_sup+nsup-1
 702                  distance=dist(i,j)
 703                  if (distance.le.15.0) then
 704                  nhpb=nhpb+1
 705                  ihpb(nhpb)=i+nstart_seq-nstart_sup
 706                  jhpb(nhpb)=j+nstart_seq-nstart_sup
 707                  forcon(nhpb)=sqrt(0.04*distance)
 708                  fordepth(nhpb)=sqrt(40.0/distance)
 709                  dhpb(nhpb)=distance-0.1d0
 710                  dhpb1(nhpb)=distance+0.1d0
 711
 712 #ifdef MPI
 713           if (.not.out1file .or. me.eq.king) &
 714           write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
 715           nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 716 #else
 717           write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
 718           nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 719 #endif
 720             endif
 721              enddo
 722            enddo
 723       else
 724       do i=nstart_sup,nstart_sup+nsup-1
 725         do j=i+2,nstart_sup+nsup-1
 726           nhpb=nhpb+1
 727           ihpb(nhpb)=i+nstart_seq-nstart_sup
 728           jhpb(nhpb)=j+nstart_seq-nstart_sup
 729           forcon(nhpb)=weidis
 730           dhpb(nhpb)=dist(i,j)
 731         enddo
 732       enddo
 733       endif
 734       return
 735       end subroutine gen_dist_constr2
 736
 737 !-----------------------------------------------------------------------------
 738 #ifdef WINIFL
 739       subroutine flush(iu)
 740       return
 741       end subroutine flush
 742 #endif
 743 #ifdef AIX
 744       subroutine flush(iu)
 745       call flush_(iu)
 746       return
 747       end subroutine flush
 748 #endif
 749 !-----------------------------------------------------------------------------
 750       subroutine copy_to_tmp(source)
 751
 752 !      include "DIMENSIONS"
 753 !      include "COMMON.IOUNITS"
 754       character*(*) :: source
 755       character(len=256) :: tmpfile
 756 !      integer ilen
 757 !el      external ilen
 758       logical :: ex
 759       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
 760       inquire(file=tmpfile,exist=ex)
 761       if (ex) then
 762         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
 763          " to temporary directory..."
 764         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
 765         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
 766       endif
 767       return
 768       end subroutine copy_to_tmp
 769 !-----------------------------------------------------------------------------
 770       subroutine move_from_tmp(source)
 771
 772 !      include "DIMENSIONS"
 773 !      include "COMMON.IOUNITS"
 774       character*(*) :: source
 775 !      integer ilen
 776 !el      external ilen
 777       write (*,*) "Moving ",source(:ilen(source)),&
 778        " from temporary directory to working directory"
 779       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
 780       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
 781       return
 782       end subroutine move_from_tmp
 783 !-----------------------------------------------------------------------------
 784 ! misc.f
 785 !-----------------------------------------------------------------------------
 786 ! $Date: 1994/10/12 17:24:21 $
 787 ! $Revision: 2.5 $
 788
 789       logical function find_arg(ipos,line,errflag)
 790
 791       integer, parameter :: maxlen = 80
 792       character(len=80) :: line
 793       character(len=1) :: empty=' ',equal='='
 794       logical :: errflag
 795       integer :: ipos
 796 ! This function returns .TRUE., if an argument follows keyword keywd; if so
 797 ! IPOS will point to the first non-blank character of the argument. Returns
 798 ! .FALSE., if no argument follows the keyword; in this case IPOS points
 799 ! to the first non-blank character of the next keyword.
 800
 801       do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
 802         ipos=ipos+1
 803       enddo
 804       errflag=.false.
 805       if (line(ipos:ipos).eq.equal) then
 806          find_arg=.true.
 807          ipos=ipos+1
 808          do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
 809            ipos=ipos+1
 810          enddo
 811          if (ipos.gt.maxlen) errflag=.true.
 812       else
 813          find_arg=.false.
 814       endif
 815
 816       return
 817       end function find_arg
 818 !-----------------------------------------------------------------------------
 819       logical function find_group(iunit,jout,key1)
 820
 821       character*(*) :: key1
 822       character(len=80) :: karta        !,ucase
 823       integer :: iunit,jout
 824       integer :: ll !,ilen
 825 !EL      external ilen
 826 !EL      logical lcom
 827       rewind (iunit)
 828       karta=' '
 829       ll=ilen(key1)
 830       do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
 831         read (iunit,'(a)',end=10) karta
 832       enddo
 833       write (jout,'(2a)') '> ',karta(1:78)
 834       find_group=.true.
 835       return
 836    10 find_group=.false.
 837       return
 838       end function find_group
 839 #endif
 840 !-----------------------------------------------------------------------------
 841       logical function iblnk(charc)
 842       character(len=1) :: charc
 843       integer :: n
 844       n = ichar(charc)
 845       iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
 846       return
 847       end function iblnk
 848 !-----------------------------------------------------------------------------
 849 #ifndef XDRFPDB
 850       integer function ilen(string)
 851       character*(*) ::  string
 852 !EL      logical :: iblnk
 853
 854       ilen = len(string)
 855 1     if ( ilen .gt. 0 ) then
 856          if ( iblnk( string(ilen:ilen) ) ) then
 857             ilen = ilen - 1
 858             goto 1
 859          endif
 860       endif
 861       return
 862       end function ilen
 863 !-----------------------------------------------------------------------------
 864       integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
 865       integer :: nkey,i,ikey,narg
 866       character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
 867 !      character(len=16) :: ucase
 868
 869       do i=1,narg
 870         if (ucase(keywd).eq.keywdset(i,ikey)) then
 871 ! Match found
 872           in_keywd_set=i
 873           return
 874         endif
 875       enddo
 876 ! No match to the allowed set of keywords if this point is reached.
 877       in_keywd_set=0
 878       return
 879       end function in_keywd_set
 880 !-----------------------------------------------------------------------------
 881       character function lcase(string)
 882       integer :: i, k, idiff
 883       character*(*) :: string
 884       character(len=1) :: c
 885       character(len=40) :: chtmp
 886 !
 887       i = len(lcase)
 888       k = len(string)
 889       if (i .lt. k) then
 890          k = i
 891          if (string(k+1:) .ne. ' ') then
 892             chtmp = string
 893          endif
 894       endif
 895       idiff = ichar('a') - ichar('A')
 896       lcase = string
 897       do 99 i = 1, k
 898          c = string(i:i)
 899          if (lge(c,'A') .and. lle(c,'Z')) then
 900             lcase(i:i) = char(ichar(c) + idiff)
 901          endif
 902    99 continue
 903       return
 904       end function lcase
 905 !-----------------------------------------------------------------------------
 906       logical function lcom(ipos,karta)
 907       character(len=80) :: karta
 908       character :: koment(2) = (/'!','#'/)
 909       integer :: ipos,i
 910
 911       lcom=.false.
 912       do i=1,2
 913         if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
 914       enddo
 915       return
 916       end function lcom
 917 !-----------------------------------------------------------------------------
 918       logical function lower_case(ch)
 919       character*(*) :: ch
 920       lower_case=(ch.ge.'a' .and. ch.le.'z')
 921       return
 922       end function lower_case
 923 !-----------------------------------------------------------------------------
 924       subroutine mykey(line,keywd,ipos,blankline,errflag)
 925
 926 ! This subroutine seeks a non-empty substring keywd in the string LINE.
 927 ! The substring begins with the first character different from blank and
 928 ! "=" encountered right to the pointer IPOS (inclusively) and terminates
 929 ! at the character left to the first blank or "=". When the subroutine is
 930 ! exited, the pointer IPOS is moved to the position of the terminator in LINE.
 931 ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
 932 ! only separators or the maximum length of the data line (80) has been reached.
 933 ! The logical variable ERRFLAG is set at .TRUE. if the string
 934 ! consists only from a "=".
 935       integer, parameter :: maxlen=80
 936       character(len=1) :: empty=' ',equal='=',comma=','
 937       character*(*) :: keywd
 938       character(len=80) :: line
 939       logical :: blankline,errflag      !EL,lcom
 940       integer :: ipos,istart,iend
 941       errflag=.false.
 942       do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
 943         ipos=ipos+1
 944       enddo
 945       if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
 946 ! At this point the rest of the input line turned out to contain only blanks
 947 ! or to be commented out.
 948         blankline=.true.
 949         return
 950       endif
 951       blankline=.false.
 952       istart=ipos
 953 ! Checks whether the current char is a separator.
 954       do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
 955        .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
 956         ipos=ipos+1
 957       enddo
 958       iend=ipos-1
 959 ! Error flag set to .true., if the length of the keyword was found less than 1.
 960       if (iend.lt.istart) then
 961         errflag=.true.
 962         return
 963       endif
 964       keywd=line(istart:iend)
 965       return
 966       end subroutine  mykey
 967 !-----------------------------------------------------------------------------
 968       subroutine numstr(inum,numm)
 969       character(len=10) :: huj='0123456789'
 970       character*(*) :: numm
 971       integer :: inumm,inum,inum1,inum2
 972       inumm=inum
 973       inum1=inumm/10
 974       inum2=inumm-10*inum1
 975       inumm=inum1
 976       numm(3:3)=huj(inum2+1:inum2+1)
 977       inum1=inumm/10
 978       inum2=inumm-10*inum1
 979       inumm=inum1
 980       numm(2:2)=huj(inum2+1:inum2+1)
 981       inum1=inumm/10
 982       inum2=inumm-10*inum1
 983       inumm=inum1
 984       numm(1:1)=huj(inum2+1:inum2+1)
 985       return
 986       end subroutine numstr
 987 #endif
 988 !-----------------------------------------------------------------------------
 989       function ucase(string)
 990       integer :: i, k, idiff
 991       character(*) :: string
 992       character(len=len(string)) :: ucase
 993       character(len=1) :: c
 994       character(len=40) :: chtmp
 995 !
 996       i = len(ucase)
 997       k = len(string)
 998       if (i .lt. k) then
 999          k = i
1000          if (string(k+1:) .ne. ' ') then
1001             chtmp = string
1002          endif
1003       endif
1004       idiff = ichar('a') - ichar('A')
1005       ucase = string
1006       do 99 i = 1, k
1007          c = string(i:i)
1008          if (lge(c,'a') .and. lle(c,'z')) then
1009             ucase(i:i) = char(ichar(c) - idiff)
1010          endif
1011    99 continue
1012       return
1013       end function ucase
1014 !-----------------------------------------------------------------------------
1015 ! geomout.F
1016 !-----------------------------------------------------------------------------
1017       subroutine pdbout(etot,tytul,iunit)
1018
1019       use geometry_data, only: c,nres,boxxsize,boxysize,boxzsize
1020       use energy_data
1021 !      use energy, only: to_box
1022 !     use control
1023       use compare_data
1024       use MD_data
1025 !      implicit real*8 (a-h,o-z)
1026 !      include 'DIMENSIONS'
1027 !      include 'COMMON.CHAIN'
1028 !      include 'COMMON.INTERACT'
1029 !      include 'COMMON.NAMES'
1030 !      include 'COMMON.IOUNITS'
1031 !      include 'COMMON.HEADER'
1032 !      include 'COMMON.SBRIDGE'
1033 !      include 'COMMON.DISTFIT'
1034 !      include 'COMMON.MD'
1035 !el      character(len=50) :: tytul
1036       character*(*) :: tytul
1037       character(len=1),dimension(58) :: chainid=  (/'A','B','C','D','E','F','G','H','I','J',&
1038          'K','L','M','O','Q','P','R','S','T','U','W','X','Y','Z','a','b','c','d','e','f',&
1039             'g','h','i','j','k','l','m','n','o','p','r','s','t','u','w','x','y','z',&
1040             '0','1','2','3','4','5','6','7','8','9'/)
1041 #ifdef XDRFPDB
1042       integer,dimension(maxres) :: ica  !(maxres)
1043 #else
1044       integer,dimension(nres) :: ica    !(maxres)
1045 #endif
1046       integer :: iti1
1047 !el  local variables
1048       integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
1049       real(kind=8) :: etot,xi,yi,zi
1050       integer :: nres2
1051       nres2=2*nres
1052 #ifdef XDRFPDB
1053       if(.not.allocated(molnum)) allocate(molnum(maxres))
1054 !      molnum(:)=mnumi(:)
1055       if(.not.allocated(vtot)) allocate(vtot(maxres*2)) !(maxres2)
1056       if(.not.allocated(istype)) allocate(istype(maxres))
1057 #else
1058       if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
1059 #endif
1060
1061       write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
1062 !model      write (iunit,'(a5,i6)') 'MODEL',1
1063 #ifndef XDRFPDB
1064       if (nhfrag.gt.0) then
1065        do j=1,nhfrag
1066         iti=itype(hfrag(1,j),1)
1067         itj=itype(hfrag(2,j),1)
1068         if (j.lt.10) then
1069            write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
1070                  'HELIX',j,'H',j,&
1071                  restyp(iti,1),hfrag(1,j)-1,&
1072                  restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1073         else
1074              write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
1075                  'HELIX',j,'H',j,&
1076                  restyp(iti,1),hfrag(1,j)-1,&
1077                  restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1078         endif
1079        enddo
1080       endif
1081
1082       if (nbfrag.gt.0) then
1083
1084        do j=1,nbfrag
1085
1086         iti=itype(bfrag(1,j),1)
1087         itj=itype(bfrag(2,j)-1,1)
1088
1089         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
1090                  'SHEET',1,'B',j,2,&
1091                  restyp(iti,1),bfrag(1,j)-1,&
1092                  restyp(itj,1),bfrag(2,j)-2,0
1093
1094         if (bfrag(3,j).gt.bfrag(4,j)) then
1095
1096          itk=itype(bfrag(3,j),1)
1097          itl=itype(bfrag(4,j)+1,1)
1098
1099          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)') &
1100                  'SHEET',2,'B',j,2,&
1101                  restyp(itl,1),bfrag(4,j),&
1102                  restyp(itk,1),bfrag(3,j)-1,-1,&
1103                  "N",restyp(itk,1),bfrag(3,j)-1,&
1104                  "O",restyp(iti,1),bfrag(1,j)-1
1105
1106         else
1107
1108          itk=itype(bfrag(3,j),1)
1109          itl=itype(bfrag(4,j)-1,1)
1110
1111
1112         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)') &
1113                  'SHEET',2,'B',j,2,&
1114                  restyp(itk,1),bfrag(3,j)-1,&
1115                  restyp(itl,1),bfrag(4,j)-2,1,&
1116                  "N",restyp(itk,1),bfrag(3,j)-1,&
1117                  "O",restyp(iti,1),bfrag(1,j)-1
1118
1119
1120
1121         endif
1122
1123        enddo
1124       endif
1125 #endif
1126       if (nss.gt.0) then
1127         do i=1,nss
1128          if (dyn_ss) then
1129           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1130                'SSBOND',i,'CYS',idssb(i)-nnt+1,&
1131                           'CYS',jdssb(i)-nnt+1
1132          else
1133           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1134                'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
1135                           'CYS',jhpb(i)-nnt+1-nres
1136          endif
1137         enddo
1138       endif
1139
1140       iatom=0
1141       ichain=1
1142       ires=0
1143       write(iout,*) "TUTUT"
1144       do i=nnt,nct
1145         write(iout,*), "coord",c(1,i),c(2,i),c(3,i)
1146         iti=itype(i,molnum(i))
1147         print *,i,molnum(i)
1148         if (molnum(i+1).eq.0) then
1149         iti1=ntyp1_molec(molnum(i))
1150         else
1151         iti1=itype(i+1,molnum(i+1))
1152         endif
1153         if ((iti.eq.ntyp1_molec(molnum(i))).and.(iti1.eq.ntyp1_molec(molnum(i)))) cycle
1154         if (i.lt.nnt) then
1155         if (iti.eq.ntyp1_molec(molnum(i)).and.(molnum(i+1).eq.5)) cycle
1156         endif
1157         if (iti.eq.ntyp1_molec(molnum(i))) then
1158           ichain=ichain+1
1159           if (ichain.gt.58) ichain=1
1160 !          ires=0
1161           write (iunit,'(a)') 'TER'
1162         else
1163         ires=ires+1
1164         iatom=iatom+1
1165         ica(i)=iatom
1166         if (molnum(i).eq.1) then
1167
1168         write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1169            ires,(c(j,i),j=1,3),vtot(i)
1170         elseif(molnum(i).eq.2) then
1171            if (istype(i).eq.0) istype(i)=1
1172         write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), &
1173           chainid(ichain),ires,(c(j,i),j=1,3),vtot(i)
1174         elseif (molnum(i).eq.4) then
1175         xi=c(1,i)
1176         yi=c(2,i)
1177         zi=c(3,i)
1178 !        xi=dmod(xi,boxxsize)
1179 !        if (xi.lt.0.0d0) xi=xi+boxxsize
1180 !        yi=dmod(yi,boxysize)
1181 !        if (yi.lt.0.0d0) yi=yi+boxysize
1182 !        zi=dmod(zi,boxzsize)
1183 !        if (zi.lt.0.0d0) zi=zi+boxzsize
1184         write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1185            ires,xi,yi,zi,vtot(i)
1186         else
1187         write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1188            ires,(c(j,i),j=1,3),vtot(i)
1189         endif
1190         if ((iti.ne.10).and.(molnum(i).ne.5)) then
1191           iatom=iatom+1
1192           if (molnum(i).eq.1) then
1193            write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1194             ires,(c(j,nres+i),j=1,3),&
1195             vtot(i+nres)
1196            else if (molnum(i).eq.2) then
1197            if (istype(i).eq.0) istype(i)=1
1198           write (iunit,50) iatom,sugartyp(istype(i)),restyp(iti,2), &
1199            chainid(ichain),ires,(c(j,nres+i),j=1,3),vtot(i+nres)
1200            endif
1201
1202         endif
1203         endif
1204       enddo
1205       write (iunit,'(a)') 'TER'
1206       do i=nnt,nct-1
1207         if (itype(i,1).eq.ntyp1) cycle
1208         if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
1209           write (iunit,30) ica(i),ica(i+1)
1210         else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1211           write (iunit,30) ica(i),ica(i+1),ica(i)+1
1212         else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1213           write (iunit,30) ica(i),ica(i)+1
1214         endif
1215       enddo
1216       if ((itype(nct,1).ne.10).and.(molnum(i).ne.5)) then
1217         write (iunit,30) ica(nct),ica(nct)+1
1218       endif
1219       do i=1,nss
1220        if (dyn_ss) then
1221         write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1222        else
1223         write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1224        endif
1225       enddo
1226       write (iunit,'(a6)') 'ENDMDL'
1227   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1228
1229   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1230   40  FORMAT ("ATOM",I7,"  C5' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1231   50  FORMAT ("ATOM",I7,"  C1' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1232
1233   30  FORMAT ('CONECT',8I5)
1234   60  FORMAT ('HETATM',I5,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1235
1236       return
1237       end subroutine pdbout
1238 #ifndef XDRFPDB
1239 !-----------------------------------------------------------------------------
1240       subroutine MOL2out(etot,tytul)
1241 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1242 ! format.
1243       use geometry_data, only: c
1244       use energy_data
1245 !      use control
1246 !      implicit real*8 (a-h,o-z)
1247 !      include 'DIMENSIONS'
1248 !      include 'COMMON.CHAIN'
1249 !      include 'COMMON.INTERACT'
1250 !      include 'COMMON.NAMES'
1251 !      include 'COMMON.IOUNITS'
1252 !      include 'COMMON.HEADER'
1253 !      include 'COMMON.SBRIDGE'
1254       character(len=32) :: tytul,fd
1255       character(len=3) :: zahl
1256       character(len=6) :: res_num,pom   !,ucase
1257
1258 !el  local variables
1259       integer :: i,j
1260       real(kind=8) :: etot
1261
1262 #ifdef AIX
1263       call fdate_(fd)
1264 #elif (defined CRAY)
1265       call date(fd)
1266 #else
1267       call fdate(fd)
1268 #endif
1269       write (imol2,'(a)') '#'
1270       write (imol2,'(a)') &
1271        '#         Creating user name:           unres'
1272       write (imol2,'(2a)') '#         Creation time:                ',&
1273        fd
1274       write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1275       write (imol2,'(a)') tytul
1276       write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1277       write (imol2,'(a)') 'SMALL'
1278       write (imol2,'(a)') 'USER_CHARGES'
1279       write (imol2,'(a)') '\@<TRIPOS>ATOM'
1280       do i=nnt,nct
1281         write (zahl,'(i3)') i
1282         pom=ucase(restyp(itype(i,1),1))
1283         res_num = pom(:3)//zahl(2:)
1284         write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1285       enddo
1286       write (imol2,'(a)') '\@<TRIPOS>BOND'
1287       do i=nnt,nct-1
1288         write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1289       enddo
1290       do i=1,nss
1291         write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1292       enddo
1293       write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1294       do i=nnt,nct
1295         write (zahl,'(i3)') i
1296         pom = ucase(restyp(itype(i,1),1))
1297         res_num = pom(:3)//zahl(2:)
1298         write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1299       enddo
1300   10  FORMAT (I7,' CA      ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1301   30  FORMAT (I7,1x,A,I14,' RESIDUE',I13,' ****  ****')
1302       return
1303       end subroutine MOL2out
1304 !-----------------------------------------------------------------------------
1305       subroutine intout
1306
1307       use geometry_data
1308       use energy_data, only: itype
1309 !      use control
1310 !      implicit real*8 (a-h,o-z)
1311 !      include 'DIMENSIONS'
1312 !      include 'COMMON.IOUNITS'
1313 !      include 'COMMON.CHAIN'
1314 !      include 'COMMON.VAR'
1315 !      include 'COMMON.LOCAL'
1316 !      include 'COMMON.INTERACT'
1317 !      include 'COMMON.NAMES'
1318 !      include 'COMMON.GEO'
1319 !      include 'COMMON.TORSION'
1320 !el  local variables
1321       integer :: i,iti
1322
1323       write (iout,'(/a)') 'Geometry of the virtual chain.'
1324       write (iout,'(7a)') '  Res  ','         d','     Theta',&
1325        '       Phi','       Dsc','     Alpha','      Omega'
1326       do i=1,nres
1327         iti=itype(i,1)
1328 !         print *,vbld(i),"vbld(i)",i
1329         write (iout,'(a3,i6,6f10.3)') restyp(iti,1),i,vbld(i),&
1330            rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1331            rad2deg*omeg(i)
1332       enddo
1333       return
1334       end subroutine intout
1335 !-----------------------------------------------------------------------------
1336 #ifdef CLUSTER
1337       subroutine briefout(it,ener,free)!,plik)
1338 #else
1339       subroutine briefout(it,ener)
1340 #endif
1341       use geometry_data
1342       use energy_data
1343 !      use control
1344 !      implicit real*8 (a-h,o-z)
1345 !      include 'DIMENSIONS'
1346 !      include 'COMMON.IOUNITS'
1347 !      include 'COMMON.CHAIN'
1348 !      include 'COMMON.VAR'
1349 !      include 'COMMON.LOCAL'
1350 !      include 'COMMON.INTERACT'
1351 !      include 'COMMON.NAMES'
1352 !      include 'COMMON.GEO'
1353 !      include 'COMMON.SBRIDGE'
1354 !     print '(a,i5)',intname,igeom
1355 !el  local variables
1356       integer :: i,it
1357       real(kind=8) :: ener,free
1358 !     character(len=80) :: plik
1359 !      integer :: iii
1360
1361 #if defined(AIX) || defined(PGI)
1362       open (igeom,file=intname,position='append')
1363 #else
1364       open (igeom,file=intname,access='append')
1365 #endif
1366 #ifdef WHAM_RUN
1367 !      iii=igeom
1368       igeom=iout
1369 #endif
1370       print *,nss
1371       IF (NSS.LE.9) THEN
1372 #ifdef CLUSTER
1373         WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1374       ELSE
1375         WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1376 #else
1377         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1378       ELSE
1379         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1380 #endif
1381         WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1382       ENDIF
1383 !     IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1384       WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1385       WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1386 !     if (nvar.gt.nphi+ntheta) then
1387         write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1388         write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1389 !     endif
1390       close(igeom)
1391   180 format (I5,F12.3,I2,9(1X,2I3))
1392   190 format (3X,11(1X,2I3))
1393   200 format (8F10.4)
1394       return
1395       end subroutine briefout
1396 !-----------------------------------------------------------------------------
1397 #ifdef WINIFL
1398       subroutine fdate(fd)
1399       character(len=32) :: fd
1400       write(fd,'(32x)')
1401       return
1402       end subroutine fdate
1403 #endif
1404 !-----------------------------------------------------------------------------
1405 #ifdef WHAM_RUN
1406       real(kind=8) function gyrate(jcon)
1407 #else
1408       real(kind=8) function gyrate()
1409 #endif
1410
1411       use geometry_data, only: c
1412 !      use geometry
1413       use energy_data
1414 !      implicit real*8 (a-h,o-z)
1415 !      include 'DIMENSIONS'
1416 !      include 'COMMON.INTERACT'
1417 !      include 'COMMON.CHAIN'
1418       real(kind=8) :: rg
1419       real(kind=8),dimension(3) :: cen
1420 !el  local variables
1421       integer :: i,j,jcon
1422
1423       do j=1,3
1424        cen(j)=0.0d0
1425       enddo
1426
1427       do i=nnt,nct
1428           do j=1,3
1429             cen(j)=cen(j)+c(j,i)
1430           enddo
1431       enddo
1432       do j=1,3
1433             cen(j)=cen(j)/dble(nct-nnt+1)
1434       enddo
1435       rg = 0.0d0
1436       do i = nnt, nct
1437         do j=1,3
1438          rg = rg + (c(j,i)-cen(j))**2
1439         enddo
1440       end do
1441 #ifdef WHAM_RUN
1442       gyrate = dsqrt(rg/dble(nct-nnt+1))
1443 #else
1444       gyrate = sqrt(rg/dble(nct-nnt+1))
1445 #endif
1446       return
1447       end function gyrate
1448 #ifdef WHAM_RUN
1449 !-----------------------------------------------------------------------------
1450 ! readrtns.F WHAM
1451       subroutine reads(rekord,lancuch,wartosc,default)
1452 !      implicit none
1453       character*(*) :: rekord,lancuch,wartosc,default
1454       character(len=80) :: aux
1455       integer :: lenlan,lenrec,iread,ireade
1456 !el      external ilen
1457 !el      logical iblnk
1458 !el      external iblnk
1459       lenlan=ilen(lancuch)
1460       lenrec=ilen(rekord)
1461       iread=index(rekord,lancuch(:lenlan)//"=")
1462 !      print *,"rekord",rekord," lancuch",lancuch
1463 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1464       if (iread.eq.0) then
1465         wartosc=default
1466         return
1467       endif
1468       iread=iread+lenlan+1
1469 !      print *,"iread",iread
1470 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1471       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1472         iread=iread+1
1473 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1474       enddo
1475 !      print *,"iread",iread
1476       if (iread.gt.lenrec) then
1477          wartosc=default
1478         return
1479       endif
1480       ireade=iread+1
1481 !      print *,"ireade",ireade
1482       do while (ireade.lt.lenrec .and. &
1483          .not.iblnk(rekord(ireade:ireade)))
1484         ireade=ireade+1
1485       enddo
1486       wartosc=rekord(iread:ireade)
1487       return
1488       end subroutine reads
1489 #endif
1490 !-----------------------------------------------------------------------------
1491 ! permut.F
1492 !-----------------------------------------------------------------------------
1493       subroutine permut(isym,nperm,tabperm)
1494 !c      integer maxperm,maxsym
1495 !c      parameter (maxperm=3628800)
1496 !c      parameter (maxsym=10)
1497 !      include "DIMENSIONS"
1498       integer n,a,tabperm,nperm,kkk,i,isym
1499 !      logical nextp
1500 !      external nextp
1501       dimension a(isym),tabperm(50,5040)
1502       n=isym
1503       nperm=1
1504       if (n.eq.1) then
1505         tabperm(1,1)=1
1506         return
1507       endif
1508       do i=2,n
1509         nperm=nperm*i
1510       enddo
1511       kkk=0
1512       do i=1,n
1513       a(i)=i
1514       enddo
1515    10 continue
1516 !c     print '(i3,2x,100i3)',kkk+1,(a(i),i=1,n)
1517       kkk=kkk+1
1518       do i=1,n
1519       tabperm(i,kkk)=a(i)
1520       enddo
1521       if(nextp(n,a)) go to 10
1522       return
1523       end subroutine
1524
1525
1526 !-----------------------------------------------------------------------------
1527       logical function nextp(n,a)
1528
1529       integer :: n,i,j,k,t
1530 !      logical :: nextp
1531       integer,dimension(n) :: a
1532       i=n-1
1533    10 if(a(i).lt.a(i+1)) go to 20
1534       i=i-1
1535       if(i.eq.0) go to 20
1536       go to 10
1537    20 j=i+1
1538       k=n
1539    30 t=a(j)
1540       a(j)=a(k)
1541       a(k)=t
1542       j=j+1
1543       k=k-1
1544       if(j.lt.k) go to 30
1545       j=i
1546       if(j.ne.0) go to 40
1547       nextp=.false.
1548       return
1549    40 j=j+1
1550       if(a(j).lt.a(i)) go to 40
1551       t=a(i)
1552       a(i)=a(j)
1553       a(j)=t
1554       nextp=.true.
1555       return
1556       end function nextp
1557 #endif
1558 !-----------------------------------------------------------------------------
1559 ! rescode.f
1560 !-----------------------------------------------------------------------------
1561       integer function rescode(iseq,nam,itype,molecule)
1562
1563 !      use io_base, only: ucase
1564 !      implicit real*8 (a-h,o-z)
1565 !      include 'DIMENSIONS'
1566 !      include 'COMMON.NAMES'
1567 !      include 'COMMON.IOUNITS'
1568       character(len=3) :: nam   !,ucase
1569       integer :: iseq,itype,i
1570       integer :: molecule
1571       print *,molecule,nam
1572       if (molecule.eq.1) then
1573       if (itype.eq.0) then
1574
1575       do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
1576         if (ucase(nam).eq.restyp(i,molecule)) then
1577           rescode=i
1578           return
1579         endif
1580       enddo
1581
1582       else
1583
1584       do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
1585         if (nam(1:1).eq.onelet(i)) then
1586           rescode=i
1587           return
1588         endif
1589       enddo
1590
1591       endif
1592       else if (molecule.eq.2) then
1593       do i=1,ntyp1_molec(molecule)
1594          print *,nam(1:1),restyp(i,molecule)(1:1)
1595          print *,nam(2:2),"read",i
1596         if (nam(2:2).eq.restyp(i,molecule)(1:1)) then
1597           rescode=i
1598           return
1599         endif
1600       enddo
1601       else if (molecule.eq.3) then
1602        write(iout,*) "SUGAR not yet implemented"
1603        stop
1604       else if (molecule.eq.4) then
1605        write(iout,*) "Explicit LIPID not yet implemented"
1606        stop
1607       else if (molecule.eq.5) then
1608       do i=-1,ntyp1_molec(molecule)
1609         print *,restyp(i,molecule)
1610         print *,i,restyp(i,molecule)(1:2),nam(1:2)
1611         if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then
1612           rescode=i
1613           return
1614         endif
1615       enddo
1616       else
1617        write(iout,*) "molecule not defined"
1618       endif
1619       write (iout,10) iseq,nam
1620       stop
1621    10 format ('**** Error - residue',i4,' has an unresolved name ',a3)
1622       end function rescode
1623       integer function sugarcode(sugar,ires)
1624       character sugar
1625       integer ires
1626       if (sugar.eq.'D') then
1627         sugarcode=1
1628       else if (sugar.eq.' ') then
1629         sugarcode=2
1630       else
1631         write (iout,*) 'UNKNOWN sugar type for residue',ires,' ',sugar
1632         stop
1633       endif
1634       return
1635       end function sugarcode
1636       integer function rescode_lip(res,ires)
1637       character(len=3):: res
1638       integer ires
1639        rescode_lip=0
1640        if  (res.eq.'NC3')  rescode_lip =4
1641        if  (res.eq.'NH3')  rescode_lip =2
1642        if  (res.eq.'PO4')  rescode_lip =3
1643        if  (res.eq.'GL1')  rescode_lip =12
1644        if  (res.eq.'GL2')  rescode_lip =12
1645        if  (res.eq.'C1A')  rescode_lip =18
1646        if  (res.eq.'C2A')  rescode_lip =18
1647        if  (res.eq.'C3A')  rescode_lip =18
1648        if  (res.eq.'C4A')  rescode_lip =18
1649        if  (res.eq.'C1B')  rescode_lip =18
1650        if  (res.eq.'C2B')  rescode_lip =18
1651        if  (res.eq.'C3B')  rescode_lip =18
1652        if  (res.eq.'C4B')  rescode_lip =18
1653        if  (res.eq.'D2A')  rescode_lip =16
1654        if  (res.eq.'D2B')  rescode_lip =16
1655
1656         if (rescode_lip.eq.0) write(iout,*) "UNKNOWN RESIDUE",ires,res
1657       return
1658       end function
1659 !-----------------------------------------------------------------------------
1660 !-----------------------------------------------------------------------------
1661       end module io_base