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 control
1022       use compare_data
1023       use MD_data
1024 !      implicit real*8 (a-h,o-z)
1025 !      include 'DIMENSIONS'
1026 !      include 'COMMON.CHAIN'
1027 !      include 'COMMON.INTERACT'
1028 !      include 'COMMON.NAMES'
1029 !      include 'COMMON.IOUNITS'
1030 !      include 'COMMON.HEADER'
1031 !      include 'COMMON.SBRIDGE'
1032 !      include 'COMMON.DISTFIT'
1033 !      include 'COMMON.MD'
1034 !el      character(len=50) :: tytul
1035       character*(*) :: tytul
1036       character(len=1),dimension(58) :: chainid=  (/'A','B','C','D','E','F','G','H','I','J',&
1037          'K','L','M','O','Q','P','R','S','T','U','W','X','Y','Z','a','b','c','d','e','f',&
1038             'g','h','i','j','k','l','m','n','o','p','r','s','t','u','w','x','y','z',&
1039             '0','1','2','3','4','5','6','7','8','9'/)
1040 #ifdef XDRFPDB
1041       integer,dimension(maxres) :: ica  !(maxres)
1042 #else
1043       integer,dimension(nres) :: ica    !(maxres)
1044 #endif
1045       integer :: iti1
1046 !el  local variables
1047       integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
1048       real(kind=8) :: etot,xi,yi,zi
1049       integer :: nres2
1050       nres2=2*nres
1051 #ifdef XDRFPDB
1052       if(.not.allocated(molnum)) allocate(molnum(maxres))
1053 !      molnum(:)=mnumi(:)
1054       if(.not.allocated(vtot)) allocate(vtot(maxres*2)) !(maxres2)
1055       if(.not.allocated(istype)) allocate(istype(maxres))
1056 #else
1057       if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
1058 #endif
1059
1060       write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
1061 !model      write (iunit,'(a5,i6)') 'MODEL',1
1062 #ifndef XDRFPDB
1063       if (nhfrag.gt.0) then
1064        do j=1,nhfrag
1065         iti=itype(hfrag(1,j),1)
1066         itj=itype(hfrag(2,j),1)
1067         if (j.lt.10) then
1068            write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
1069                  'HELIX',j,'H',j,&
1070                  restyp(iti,1),hfrag(1,j)-1,&
1071                  restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1072         else
1073              write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
1074                  'HELIX',j,'H',j,&
1075                  restyp(iti,1),hfrag(1,j)-1,&
1076                  restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1077         endif
1078        enddo
1079       endif
1080
1081       if (nbfrag.gt.0) then
1082
1083        do j=1,nbfrag
1084
1085         iti=itype(bfrag(1,j),1)
1086         itj=itype(bfrag(2,j)-1,1)
1087
1088         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
1089                  'SHEET',1,'B',j,2,&
1090                  restyp(iti,1),bfrag(1,j)-1,&
1091                  restyp(itj,1),bfrag(2,j)-2,0
1092
1093         if (bfrag(3,j).gt.bfrag(4,j)) then
1094
1095          itk=itype(bfrag(3,j),1)
1096          itl=itype(bfrag(4,j)+1,1)
1097
1098          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)') &
1099                  'SHEET',2,'B',j,2,&
1100                  restyp(itl,1),bfrag(4,j),&
1101                  restyp(itk,1),bfrag(3,j)-1,-1,&
1102                  "N",restyp(itk,1),bfrag(3,j)-1,&
1103                  "O",restyp(iti,1),bfrag(1,j)-1
1104
1105         else
1106
1107          itk=itype(bfrag(3,j),1)
1108          itl=itype(bfrag(4,j)-1,1)
1109
1110
1111         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)') &
1112                  'SHEET',2,'B',j,2,&
1113                  restyp(itk,1),bfrag(3,j)-1,&
1114                  restyp(itl,1),bfrag(4,j)-2,1,&
1115                  "N",restyp(itk,1),bfrag(3,j)-1,&
1116                  "O",restyp(iti,1),bfrag(1,j)-1
1117
1118
1119
1120         endif
1121
1122        enddo
1123       endif
1124 #endif
1125       if (nss.gt.0) then
1126         do i=1,nss
1127          if (dyn_ss) then
1128           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1129                'SSBOND',i,'CYS',idssb(i)-nnt+1,&
1130                           'CYS',jdssb(i)-nnt+1
1131          else
1132           write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1133                'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
1134                           'CYS',jhpb(i)-nnt+1-nres
1135          endif
1136         enddo
1137       endif
1138
1139       iatom=0
1140       ichain=1
1141       ires=0
1142       write(iout,*) "TUTUT"
1143       do i=nnt,nct
1144 !        write(iout,*), "coord",c(1,i),c(2,i),c(3,i)
1145         iti=itype(i,molnum(i))
1146         print *,i,molnum(i)
1147         if (molnum(i+1).eq.0) then
1148         iti1=ntyp1_molec(molnum(i))
1149         else
1150         iti1=itype(i+1,molnum(i+1))
1151         endif
1152         if ((iti.eq.ntyp1_molec(molnum(i))).and.(iti1.eq.ntyp1_molec(molnum(i)))) cycle
1153         if (i.lt.nnt) then
1154         if (iti.eq.ntyp1_molec(molnum(i)).and.(molnum(i+1).eq.5)) cycle
1155         endif
1156         if (iti.eq.ntyp1_molec(molnum(i))) then
1157           ichain=ichain+1
1158           if (ichain.gt.58) ichain=1
1159 !          ires=0
1160           write (iunit,'(a)') 'TER'
1161         else
1162         ires=ires+1
1163         iatom=iatom+1
1164         ica(i)=iatom
1165         if (molnum(i).eq.1) then
1166
1167         write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1168            ires,(c(j,i),j=1,3),vtot(i)
1169         elseif(molnum(i).eq.2) then
1170            if (istype(i).eq.0) istype(i)=1
1171         write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), &
1172           chainid(ichain),ires,(c(j,i),j=1,3),vtot(i)
1173         elseif (molnum(i).eq.4) then
1174         xi=c(1,i)
1175         yi=c(2,i)
1176         zi=c(3,i)
1177 !        xi=dmod(xi,boxxsize)
1178 !        if (xi.lt.0.0d0) xi=xi+boxxsize
1179 !        yi=dmod(yi,boxysize)
1180 !        if (yi.lt.0.0d0) yi=yi+boxysize
1181 !        zi=dmod(zi,boxzsize)
1182 !        if (zi.lt.0.0d0) zi=zi+boxzsize
1183         write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1184            ires,xi,yi,zi,vtot(i)
1185         elseif (molnum(i).eq.5) then
1186         xi=c(1,i)
1187         yi=c(2,i)
1188         zi=c(3,i)
1189         xi=dmod(xi,boxxsize)
1190         if (xi.lt.0.0d0) xi=xi+boxxsize
1191         yi=dmod(yi,boxysize)
1192         if (yi.lt.0.0d0) yi=yi+boxysize
1193         zi=dmod(zi,boxzsize)
1194         if (zi.lt.0.0d0) zi=zi+boxzsize
1195         write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1196            ires,xi,yi,zi,vtot(i)
1197         else
1198
1199         write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1200            ires,(c(j,i),j=1,3),vtot(i)
1201         endif
1202         if ((iti.ne.10).and.(molnum(i).ne.5)) then
1203           iatom=iatom+1
1204           if (molnum(i).eq.1) then
1205            write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1206             ires,(c(j,nres+i),j=1,3),&
1207             vtot(i+nres)
1208            else if (molnum(i).eq.2) then
1209            if (istype(i).eq.0) istype(i)=1
1210           write (iunit,50) iatom,sugartyp(istype(i)),restyp(iti,2), &
1211            chainid(ichain),ires,(c(j,nres+i),j=1,3),vtot(i+nres)
1212            endif
1213
1214         endif
1215         endif
1216       enddo
1217       write (iunit,'(a)') 'TER'
1218       do i=nnt,nct-1
1219         if (molnum(i).eq.5) cycle
1220         if (itype(i,1).eq.ntyp1) cycle
1221         if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
1222           write (iunit,30) ica(i),ica(i+1)
1223         else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1224           write (iunit,30) ica(i),ica(i+1),ica(i)+1
1225         else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1226           write (iunit,30) ica(i),ica(i)+1
1227         endif
1228       enddo
1229       if ((itype(nct,1).ne.10).and.(molnum(i).ne.5)) then
1230         write (iunit,30) ica(nct),ica(nct)+1
1231       endif
1232       do i=1,nss
1233        if (dyn_ss) then
1234         write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1235        else
1236         write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1237        endif
1238       enddo
1239       write (iunit,'(a6)') 'ENDMDL'
1240   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1241
1242   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1243   40  FORMAT ("ATOM",I7,"  C5' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1244   50  FORMAT ("ATOM",I7,"  C1' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1245
1246   30  FORMAT ('CONECT',8I5)
1247   60  FORMAT ('HETATM',I5,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1248
1249       return
1250       end subroutine pdbout
1251 #ifndef XDRFPDB
1252 !-----------------------------------------------------------------------------
1253       subroutine MOL2out(etot,tytul)
1254 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1255 ! format.
1256       use geometry_data, only: c
1257       use energy_data
1258 !      use control
1259 !      implicit real*8 (a-h,o-z)
1260 !      include 'DIMENSIONS'
1261 !      include 'COMMON.CHAIN'
1262 !      include 'COMMON.INTERACT'
1263 !      include 'COMMON.NAMES'
1264 !      include 'COMMON.IOUNITS'
1265 !      include 'COMMON.HEADER'
1266 !      include 'COMMON.SBRIDGE'
1267       character(len=32) :: tytul,fd
1268       character(len=3) :: zahl
1269       character(len=6) :: res_num,pom   !,ucase
1270
1271 !el  local variables
1272       integer :: i,j
1273       real(kind=8) :: etot
1274
1275 #ifdef AIX
1276       call fdate_(fd)
1277 #elif (defined CRAY)
1278       call date(fd)
1279 #else
1280       call fdate(fd)
1281 #endif
1282       write (imol2,'(a)') '#'
1283       write (imol2,'(a)') &
1284        '#         Creating user name:           unres'
1285       write (imol2,'(2a)') '#         Creation time:                ',&
1286        fd
1287       write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1288       write (imol2,'(a)') tytul
1289       write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1290       write (imol2,'(a)') 'SMALL'
1291       write (imol2,'(a)') 'USER_CHARGES'
1292       write (imol2,'(a)') '\@<TRIPOS>ATOM'
1293       do i=nnt,nct
1294         write (zahl,'(i3)') i
1295         pom=ucase(restyp(itype(i,1),1))
1296         res_num = pom(:3)//zahl(2:)
1297         write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1298       enddo
1299       write (imol2,'(a)') '\@<TRIPOS>BOND'
1300       do i=nnt,nct-1
1301         write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1302       enddo
1303       do i=1,nss
1304         write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1305       enddo
1306       write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1307       do i=nnt,nct
1308         write (zahl,'(i3)') i
1309         pom = ucase(restyp(itype(i,1),1))
1310         res_num = pom(:3)//zahl(2:)
1311         write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1312       enddo
1313   10  FORMAT (I7,' CA      ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1314   30  FORMAT (I7,1x,A,I14,' RESIDUE',I13,' ****  ****')
1315       return
1316       end subroutine MOL2out
1317 !-----------------------------------------------------------------------------
1318       subroutine intout
1319
1320       use geometry_data
1321       use energy_data, only: itype
1322 !      use control
1323 !      implicit real*8 (a-h,o-z)
1324 !      include 'DIMENSIONS'
1325 !      include 'COMMON.IOUNITS'
1326 !      include 'COMMON.CHAIN'
1327 !      include 'COMMON.VAR'
1328 !      include 'COMMON.LOCAL'
1329 !      include 'COMMON.INTERACT'
1330 !      include 'COMMON.NAMES'
1331 !      include 'COMMON.GEO'
1332 !      include 'COMMON.TORSION'
1333 !el  local variables
1334       integer :: i,iti
1335
1336       write (iout,'(/a)') 'Geometry of the virtual chain.'
1337       write (iout,'(7a)') '  Res  ','         d','     Theta',&
1338        '       Phi','       Dsc','     Alpha','      Omega'
1339       do i=1,nres
1340         iti=itype(i,1)
1341 !         print *,vbld(i),"vbld(i)",i
1342         write (iout,'(a3,i6,6f10.3)') restyp(iti,1),i,vbld(i),&
1343            rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1344            rad2deg*omeg(i)
1345       enddo
1346       return
1347       end subroutine intout
1348 !-----------------------------------------------------------------------------
1349 #ifdef CLUSTER
1350       subroutine briefout(it,ener,free)!,plik)
1351 #else
1352       subroutine briefout(it,ener)
1353 #endif
1354       use geometry_data
1355       use energy_data
1356 !      use control
1357 !      implicit real*8 (a-h,o-z)
1358 !      include 'DIMENSIONS'
1359 !      include 'COMMON.IOUNITS'
1360 !      include 'COMMON.CHAIN'
1361 !      include 'COMMON.VAR'
1362 !      include 'COMMON.LOCAL'
1363 !      include 'COMMON.INTERACT'
1364 !      include 'COMMON.NAMES'
1365 !      include 'COMMON.GEO'
1366 !      include 'COMMON.SBRIDGE'
1367 !     print '(a,i5)',intname,igeom
1368 !el  local variables
1369       integer :: i,it
1370       real(kind=8) :: ener,free
1371 !     character(len=80) :: plik
1372 !      integer :: iii
1373
1374 #if defined(AIX) || defined(PGI)
1375       open (igeom,file=intname,position='append')
1376 #else
1377       open (igeom,file=intname,access='append')
1378 #endif
1379 #ifdef WHAM_RUN
1380 !      iii=igeom
1381       igeom=iout
1382 #endif
1383       print *,nss
1384       IF (NSS.LE.9) THEN
1385 #ifdef CLUSTER
1386         WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1387       ELSE
1388         WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1389 #else
1390         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1391       ELSE
1392         WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1393 #endif
1394         WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1395       ENDIF
1396 !     IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1397       WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1398       WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1399 !     if (nvar.gt.nphi+ntheta) then
1400         write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1401         write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1402 !     endif
1403       close(igeom)
1404   180 format (I5,F12.3,I2,9(1X,2I3))
1405   190 format (3X,11(1X,2I3))
1406   200 format (8F10.4)
1407       return
1408       end subroutine briefout
1409 !-----------------------------------------------------------------------------
1410 #ifdef WINIFL
1411       subroutine fdate(fd)
1412       character(len=32) :: fd
1413       write(fd,'(32x)')
1414       return
1415       end subroutine fdate
1416 #endif
1417 !-----------------------------------------------------------------------------
1418 #ifdef WHAM_RUN
1419       real(kind=8) function gyrate(jcon)
1420 #else
1421       real(kind=8) function gyrate()
1422 #endif
1423
1424       use geometry_data, only: c
1425 !      use geometry
1426       use energy_data
1427 !      implicit real*8 (a-h,o-z)
1428 !      include 'DIMENSIONS'
1429 !      include 'COMMON.INTERACT'
1430 !      include 'COMMON.CHAIN'
1431       real(kind=8) :: rg
1432       real(kind=8),dimension(3) :: cen
1433 !el  local variables
1434       integer :: i,j,jcon
1435
1436       do j=1,3
1437        cen(j)=0.0d0
1438       enddo
1439
1440       do i=nnt,nct
1441           do j=1,3
1442             cen(j)=cen(j)+c(j,i)
1443           enddo
1444       enddo
1445       do j=1,3
1446             cen(j)=cen(j)/dble(nct-nnt+1)
1447       enddo
1448       rg = 0.0d0
1449       do i = nnt, nct
1450         do j=1,3
1451          rg = rg + (c(j,i)-cen(j))**2
1452         enddo
1453       end do
1454 #ifdef WHAM_RUN
1455       gyrate = dsqrt(rg/dble(nct-nnt+1))
1456 #else
1457       gyrate = sqrt(rg/dble(nct-nnt+1))
1458 #endif
1459       return
1460       end function gyrate
1461 #ifdef WHAM_RUN
1462 !-----------------------------------------------------------------------------
1463 ! readrtns.F WHAM
1464       subroutine reads(rekord,lancuch,wartosc,default)
1465 !      implicit none
1466       character*(*) :: rekord,lancuch,wartosc,default
1467       character(len=80) :: aux
1468       integer :: lenlan,lenrec,iread,ireade
1469 !el      external ilen
1470 !el      logical iblnk
1471 !el      external iblnk
1472       lenlan=ilen(lancuch)
1473       lenrec=ilen(rekord)
1474       iread=index(rekord,lancuch(:lenlan)//"=")
1475 !      print *,"rekord",rekord," lancuch",lancuch
1476 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1477       if (iread.eq.0) then
1478         wartosc=default
1479         return
1480       endif
1481       iread=iread+lenlan+1
1482 !      print *,"iread",iread
1483 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1484       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1485         iread=iread+1
1486 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1487       enddo
1488 !      print *,"iread",iread
1489       if (iread.gt.lenrec) then
1490          wartosc=default
1491         return
1492       endif
1493       ireade=iread+1
1494 !      print *,"ireade",ireade
1495       do while (ireade.lt.lenrec .and. &
1496          .not.iblnk(rekord(ireade:ireade)))
1497         ireade=ireade+1
1498       enddo
1499       wartosc=rekord(iread:ireade)
1500       return
1501       end subroutine reads
1502 #endif
1503 !-----------------------------------------------------------------------------
1504 ! permut.F
1505 !-----------------------------------------------------------------------------
1506       subroutine permut(isym,nperm,tabperm)
1507 !c      integer maxperm,maxsym
1508 !c      parameter (maxperm=3628800)
1509 !c      parameter (maxsym=10)
1510 !      include "DIMENSIONS"
1511       integer n,a,tabperm,nperm,kkk,i,isym
1512 !      logical nextp
1513 !      external nextp
1514       dimension a(isym),tabperm(50,5040)
1515       n=isym
1516       nperm=1
1517       if (n.eq.1) then
1518         tabperm(1,1)=1
1519         return
1520       endif
1521       do i=2,n
1522         nperm=nperm*i
1523       enddo
1524       kkk=0
1525       do i=1,n
1526       a(i)=i
1527       enddo
1528    10 continue
1529 !c     print '(i3,2x,100i3)',kkk+1,(a(i),i=1,n)
1530       kkk=kkk+1
1531       do i=1,n
1532       tabperm(i,kkk)=a(i)
1533       enddo
1534       if(nextp(n,a)) go to 10
1535       return
1536       end subroutine
1537
1538
1539 !-----------------------------------------------------------------------------
1540       logical function nextp(n,a)
1541
1542       integer :: n,i,j,k,t
1543 !      logical :: nextp
1544       integer,dimension(n) :: a
1545       i=n-1
1546    10 if(a(i).lt.a(i+1)) go to 20
1547       i=i-1
1548       if(i.eq.0) go to 20
1549       go to 10
1550    20 j=i+1
1551       k=n
1552    30 t=a(j)
1553       a(j)=a(k)
1554       a(k)=t
1555       j=j+1
1556       k=k-1
1557       if(j.lt.k) go to 30
1558       j=i
1559       if(j.ne.0) go to 40
1560       nextp=.false.
1561       return
1562    40 j=j+1
1563       if(a(j).lt.a(i)) go to 40
1564       t=a(i)
1565       a(i)=a(j)
1566       a(j)=t
1567       nextp=.true.
1568       return
1569       end function nextp
1570 #endif
1571 !-----------------------------------------------------------------------------
1572 ! rescode.f
1573 !-----------------------------------------------------------------------------
1574       integer function rescode(iseq,nam,itype,molecule)
1575
1576 !      implicit real*8 (a-h,o-z)
1577 !      include 'DIMENSIONS'
1578 !      include 'COMMON.NAMES'
1579 !      include 'COMMON.IOUNITS'
1580       character(len=3) :: nam   !,ucase
1581       integer :: iseq,itype,i
1582       integer :: molecule
1583       print *,molecule,nam
1584       if (molecule.eq.1) then
1585       if (itype.eq.0) then
1586
1587       do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
1588         if (ucase(nam).eq.restyp(i,molecule)) then
1589           rescode=i
1590           return
1591         endif
1592       enddo
1593
1594       else
1595
1596       do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
1597         if (nam(1:1).eq.onelet(i)) then
1598           rescode=i
1599           return
1600         endif
1601       enddo
1602
1603       endif
1604       else if (molecule.eq.2) then
1605       do i=1,ntyp1_molec(molecule)
1606          print *,nam(1:1),restyp(i,molecule)(1:1)
1607          print *,nam(2:2),"read",i
1608         if (nam(2:2).eq.restyp(i,molecule)(1:1)) then
1609           rescode=i
1610           return
1611         endif
1612       enddo
1613       else if (molecule.eq.3) then
1614        write(iout,*) "SUGAR not yet implemented"
1615        stop
1616       else if (molecule.eq.4) then
1617        write(iout,*) "Explicit LIPID not yet implemented"
1618        stop
1619       else if (molecule.eq.5) then
1620       do i=-1,ntyp1_molec(molecule)
1621         print *,restyp(i,molecule)
1622         print *,i,restyp(i,molecule)(1:2),nam(1:2)
1623         if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then
1624           rescode=i
1625           return
1626         endif
1627       enddo
1628       else
1629        write(iout,*) "molecule not defined"
1630       endif
1631       write (iout,10) iseq,nam
1632       stop
1633    10 format ('**** Error - residue',i4,' has an unresolved name ',a3)
1634       end function rescode
1635       integer function sugarcode(sugar,ires)
1636       character sugar
1637       integer ires
1638       if (sugar.eq.'D') then
1639         sugarcode=1
1640       else if (sugar.eq.' ') then
1641         sugarcode=2
1642       else
1643         write (iout,*) 'UNKNOWN sugar type for residue',ires,' ',sugar
1644         stop
1645       endif
1646       return
1647       end function sugarcode
1648       integer function rescode_lip(res,ires)
1649       character(len=3):: res
1650       integer ires
1651        rescode_lip=0
1652        if  (res.eq.'NC3')  rescode_lip =4
1653        if  (res.eq.'NH3')  rescode_lip =2
1654        if  (res.eq.'PO4')  rescode_lip =3
1655        if  (res.eq.'GL1')  rescode_lip =12
1656        if  (res.eq.'GL2')  rescode_lip =12
1657        if  (res.eq.'C1A')  rescode_lip =18
1658        if  (res.eq.'C2A')  rescode_lip =18
1659        if  (res.eq.'C3A')  rescode_lip =18
1660        if  (res.eq.'C4A')  rescode_lip =18
1661        if  (res.eq.'C1B')  rescode_lip =18
1662        if  (res.eq.'C2B')  rescode_lip =18
1663        if  (res.eq.'C3B')  rescode_lip =18
1664        if  (res.eq.'C4B')  rescode_lip =18
1665        if  (res.eq.'D2A')  rescode_lip =16
1666        if  (res.eq.'D2B')  rescode_lip =16
1667
1668         if (rescode_lip.eq.0) write(iout,*) "UNKNOWN RESIDUE",ires,res
1669       return
1670       end function
1671 !-----------------------------------------------------------------------------
1672 !-----------------------------------------------------------------------------
1673       end module io_base