source/unres/io.f90

   1       module io
   2 !-----------------------------------------------------------------------
   3       use io_units
   4       use names
   5       use io_base
   6       use io_config
   7       implicit none
   8 !-----------------------------------------------------------------------------
   9 !
  10 !
  11 !-----------------------------------------------------------------------------
  12       contains
  13 !-----------------------------------------------------------------------------
  14 ! bank.F    io_csa
  15 !-----------------------------------------------------------------------------
  16       subroutine write_csa_pdb(var,ene,nft,ik,iw_pdb)
  17
  18       use csa_data
  19       use geometry_data, only:nres,nvar
  20       use geometry, only:var_to_geom,chainbuild
  21       use compare, only:secondary2
  22 !      implicit real*8 (a-h,o-z)
  23 !      include 'DIMENSIONS'
  24 !      include 'COMMON.CSA'
  25 !      include 'COMMON.BANK'
  26 !      include 'COMMON.VAR'
  27 !      include 'COMMON.IOUNITS'
  28 !      include 'COMMON.MINIM'
  29 !      include 'COMMON.SETUP'
  30 !      include 'COMMON.GEO'
  31 !      include 'COMMON.CHAIN'
  32 !      include 'COMMON.LOCAL'
  33 !      include 'COMMON.INTERACT'
  34 !      include 'COMMON.NAMES'
  35 !      include 'COMMON.SBRIDGE'
  36       integer :: lenpre,lenpot  !,ilen
  37 !el      external ilen
  38       real(kind=8),dimension(nvar) :: var       !(maxvar)       (maxvar=6*maxres)
  39       character(len=50) :: titelloc
  40       character(len=3) :: zahl
  41       real(kind=8),dimension(mxch*(mxch+1)/2+1) :: ene
  42 !el local variables
  43       integer :: nft,ik,iw_pdb
  44
  45       nmin_csa=nmin_csa+1
  46       if(ene(1).lt.eglob_csa) then
  47         eglob_csa=ene(1)
  48         nglob_csa=nglob_csa+1
  49         call numstr(nglob_csa,zahl)
  50
  51         call var_to_geom(nvar,var)
  52         call chainbuild
  53         call secondary2(.false.)
  54
  55         lenpre=ilen(prefix)
  56         open(icsa_pdb,file=prefix(:lenpre)//'@'//zahl//'.pdb')
  57
  58         if (iw_pdb.eq.1) then
  59           write(titelloc,'(a2,i3,a3,i9,a3,i6)') &
  60           'GM',nglob_csa,' e ',nft,' m ',nmin_csa
  61         else
  62           write(titelloc,'(a2,i3,a3,i9,a3,i6,a5,f5.2,a5,f5.1)') &
  63          'GM',nglob_csa,' e ',nft,' m ',nmin_csa,' rms ',&
  64                rmsn(ik),' %NC ',pncn(ik)*100
  65         endif
  66         call pdbout(eglob_csa,titelloc,icsa_pdb)
  67         close(icsa_pdb)
  68       endif
  69
  70       return
  71       end subroutine write_csa_pdb
  72 !-----------------------------------------------------------------------------
  73 ! csa.f         io_csa
  74 !-----------------------------------------------------------------------------
  75       subroutine from_pdb(n,idum)
  76 ! This subroutine stores the UNRES int variables generated from
  77 ! subroutine readpdb into the 1st conformation of in dihang_in.
  78 ! Subsequent n-1 conformations of dihang_in have identical values
  79 ! of theta and phi as the 1st conformation but random values for
  80 ! alph and omeg.
  81 ! The array cref (also generated from subroutine readpdb) is stored
  82 ! to crefjlee to be used for rmsd calculation in CSA, if necessary.
  83
  84       use csa_data
  85       use geometry_data
  86       use random, only: ran1
  87 !      implicit real*8 (a-h,o-z)
  88 !      include 'DIMENSIONS'
  89 !      include 'COMMON.IOUNITS'
  90 !      include 'COMMON.CHAIN'
  91 !      include 'COMMON.VAR'
  92 !      include 'COMMON.BANK'
  93 !      include 'COMMON.GEO'
  94 !el local variables
  95       integer :: n,idum,m,i,j,k,kk,kkk
  96       real(kind=8) :: e
  97
  98       m=1
  99       do j=2,nres-1
 100         dihang_in(1,j,1,m)=theta(j+1)
 101         dihang_in(2,j,1,m)=phi(j+2)
 102         dihang_in(3,j,1,m)=alph(j)
 103         dihang_in(4,j,1,m)=omeg(j)
 104       enddo
 105       dihang_in(2,nres-1,1,k)=0.0d0
 106
 107       do m=2,n
 108        do k=2,nres-1
 109         dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
 110         dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
 111         if(dabs(dihang_in(3,k,1,1)).gt.1.d-6) then
 112          dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
 113          dihang_in(3,k,1,m)=dihang_in(3,k,1,m)*deg2rad
 114         endif
 115         if(dabs(dihang_in(4,k,1,1)).gt.1.d-6) then
 116          dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
 117          dihang_in(4,k,1,m)=dihang_in(4,k,1,m)*deg2rad
 118         endif
 119        enddo
 120       enddo
 121
 122 ! Store cref to crefjlee (they are in COMMON.CHAIN).
 123       do k=1,2*nres
 124        do kk=1,3
 125         kkk=1
 126         crefjlee(kk,k)=cref(kk,k,kkk)
 127        enddo
 128       enddo
 129
 130       open(icsa_native_int,file=csa_native_int,status="old")
 131       do m=1,n
 132          write(icsa_native_int,*) m,e
 133          write(icsa_native_int,200) &
 134           (dihang_in(1,k,1,m)*rad2deg,k=2,nres-1)
 135          write(icsa_native_int,200) &
 136           (dihang_in(2,k,1,m)*rad2deg,k=2,nres-2)
 137          write(icsa_native_int,200) &
 138           (dihang_in(3,k,1,m)*rad2deg,k=2,nres-1)
 139          write(icsa_native_int,200) &
 140           (dihang_in(4,k,1,m)*rad2deg,k=2,nres-1)
 141       enddo
 142
 143       do k=1,nres
 144        write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
 145       enddo
 146       close(icsa_native_int)
 147
 148   200 format (8f10.4)
 149
 150       return
 151       end subroutine from_pdb
 152 !-----------------------------------------------------------------------------
 153       subroutine from_int(n,mm,idum)
 154
 155       use csa_data
 156       use geometry_data
 157       use energy_data
 158       use geometry, only:chainbuild,gen_side
 159       use energy, only:etotal
 160       use compare
 161 !      implicit real*8 (a-h,o-z)
 162 !      include 'DIMENSIONS'
 163 !      include 'COMMON.IOUNITS'
 164 !      include 'COMMON.CHAIN'
 165 !      include 'COMMON.VAR'
 166 !      include 'COMMON.INTERACT'
 167 !      include 'COMMON.BANK'
 168 !      include 'COMMON.GEO'
 169 !      include 'COMMON.CONTACTS'
 170 !      integer ilen
 171 !el      external ilen
 172       logical :: fail
 173       real(kind=8),dimension(0:n_ene) :: energia
 174 !el local variables
 175       integer :: n,mm,idum,i,ii,j,m,k,kk,maxcount_fail,icount_fail,maxsi
 176       real(kind=8) :: co
 177
 178       open(icsa_native_int,file=csa_native_int,status="old")
 179       read (icsa_native_int,*)
 180       call read_angles(icsa_native_int,*10)
 181       goto 11
 182    10 write (iout,'(2a)') "CHUJ NASTAPIL - error in ",&
 183         csa_native_int(:ilen(csa_native_int))
 184    11 continue
 185       call intout
 186       do j=2,nres-1
 187         dihang_in(1,j,1,1)=theta(j+1)
 188         dihang_in(2,j,1,1)=phi(j+2)
 189         dihang_in(3,j,1,1)=alph(j)
 190         dihang_in(4,j,1,1)=omeg(j)
 191       enddo
 192       dihang_in(2,nres-1,1,1)=0.0d0
 193
 194 !         read(icsa_native_int,*) ind,e
 195 !         read(icsa_native_int,200) (dihang_in(1,k,1,1),k=2,nres-1)
 196 !         read(icsa_native_int,200) (dihang_in(2,k,1,1),k=2,nres-2)
 197 !         read(icsa_native_int,200) (dihang_in(3,k,1,1),k=2,nres-1)
 198 !         read(icsa_native_int,200) (dihang_in(4,k,1,1),k=2,nres-1)
 199 !         dihang_in(2,nres-1,1,1)=0.d0
 200
 201          maxsi=100
 202          maxcount_fail=100
 203
 204          do m=mm+2,n
 205 !          do k=2,nres-1
 206 !           dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
 207 !           dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
 208 !           if(abs(dihang_in(3,k,1,1)).gt.1.d-3) then
 209 !            dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
 210 !           endif
 211 !           if(abs(dihang_in(4,k,1,1)).gt.1.d-3) then
 212 !            dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
 213 !           endif
 214 !          enddo
 215 !           call intout
 216            fail=.true.
 217
 218            icount_fail=0
 219
 220            DO WHILE (FAIL .AND. ICOUNT_FAIL .LE. MAXCOUNT_FAIL)
 221
 222            do i=nnt,nct
 223              if (itype(i,1).ne.10) then
 224 !d             print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
 225                fail=.true.
 226                ii=0
 227                do while (fail .and. ii .le. maxsi)
 228                  call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail)
 229                  ii = ii+1
 230                enddo
 231              endif
 232            enddo
 233            call chainbuild
 234            call etotal(energia)
 235            fail = (energia(0).ge.1.0d20)
 236            icount_fail=icount_fail+1
 237
 238            ENDDO
 239
 240            if (icount_fail.gt.maxcount_fail) then
 241              write (iout,*) &
 242              'Failed to generate non-overlaping near-native conf.',&
 243              m
 244            endif
 245
 246            do j=2,nres-1
 247              dihang_in(1,j,1,m)=theta(j+1)
 248              dihang_in(2,j,1,m)=phi(j+2)
 249              dihang_in(3,j,1,m)=alph(j)
 250              dihang_in(4,j,1,m)=omeg(j)
 251            enddo
 252            dihang_in(2,nres-1,1,m)=0.0d0
 253          enddo
 254
 255 !      do m=1,n
 256 !        write(icsa_native_int,*) m,e
 257 !         write(icsa_native_int,200) (dihang_in(1,k,1,m),k=2,nres-1)
 258 !         write(icsa_native_int,200) (dihang_in(2,k,1,m),k=2,nres-2)
 259 !         write(icsa_native_int,200) (dihang_in(3,k,1,m),k=2,nres-1)
 260 !         write(icsa_native_int,200) (dihang_in(4,k,1,m),k=2,nres-1)
 261 !      enddo
 262 !     close(icsa_native_int)
 263
 264 !      do m=mm+2,n
 265 !       do i=1,4
 266 !        do j=2,nres-1
 267 !         dihang_in(i,j,1,m)=dihang_in(i,j,1,m)*deg2rad
 268 !        enddo
 269 !       enddo
 270 !      enddo
 271
 272       call dihang_to_c(dihang_in(1,1,1,1))
 273
 274 ! Store c to cref (they are in COMMON.CHAIN).
 275       do k=1,2*nres
 276        do kk=1,3
 277         crefjlee(kk,k)=c(kk,k)
 278        enddo
 279       enddo
 280
 281       call contact(.true.,ncont_ref,icont_ref,co)
 282
 283 !      do k=1,nres
 284 !       write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
 285 !      enddo
 286       close(icsa_native_int)
 287
 288   200 format (8f10.4)
 289
 290       return
 291       end subroutine from_int
 292 !-----------------------------------------------------------------------------
 293       subroutine dihang_to_c(aarray)
 294
 295       use geometry_data
 296       use csa_data
 297       use geometry, only:chainbuild
 298 !      implicit real*8 (a-h,o-z)
 299 !      include 'DIMENSIONS'
 300 !      include 'COMMON.CSA'
 301 !      include 'COMMON.BANK'
 302 !      include 'COMMON.CHAIN'
 303 !      include 'COMMON.GEO'
 304 !      include 'COMMON.VAR'
 305       integer :: i
 306       real(kind=8),dimension(mxang,nres,mxch) :: aarray !(mxang,maxres,mxch)
 307
 308 !     do i=4,nres
 309 !      phi(i)=dihang_in(1,i-2,1,1)
 310 !     enddo
 311       do i=2,nres-1
 312        theta(i+1)=aarray(1,i,1)
 313        phi(i+2)=aarray(2,i,1)
 314        alph(i)=aarray(3,i,1)
 315        omeg(i)=aarray(4,i,1)
 316       enddo
 317
 318       call chainbuild
 319
 320       return
 321       end subroutine dihang_to_c
 322 !-----------------------------------------------------------------------------
 323 ! geomout.F
 324 !-----------------------------------------------------------------------------
 325 #ifdef NOXDR
 326       subroutine cartout(time)
 327 #else
 328       subroutine cartoutx(time)
 329 #endif
 330       use geometry_data, only: c,nres
 331       use energy_data
 332       use MD_data, only: potE,t_bath
 333 !      implicit real*8 (a-h,o-z)
 334 !      include 'DIMENSIONS'
 335 !      include 'COMMON.CHAIN'
 336 !      include 'COMMON.INTERACT'
 337 !      include 'COMMON.NAMES'
 338 !      include 'COMMON.IOUNITS'
 339 !      include 'COMMON.HEADER'
 340 !      include 'COMMON.SBRIDGE'
 341 !      include 'COMMON.DISTFIT'
 342 !      include 'COMMON.MD'
 343       real(kind=8) :: time
 344 !el  local variables
 345       integer :: j,k,i
 346
 347 #if defined(AIX) || defined(PGI)
 348       open(icart,file=cartname,position="append")
 349 #else
 350       open(icart,file=cartname,access="append")
 351 #endif
 352       write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 353       if (dyn_ss) then
 354        write (icart,'(i4,$)') &
 355          nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)
 356       else
 357        write (icart,'(i4,$)') &
 358          nss,(ihpb(j),jhpb(j),j=1,nss)
 359        endif
 360        write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,&
 361        (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),&
 362        (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 363       write (icart,'(8f10.5)') &
 364        ((c(k,j),k=1,3),j=1,nres),&
 365        ((c(k,j+nres),k=1,3),j=nnt,nct)
 366       close(icart)
 367       return
 368
 369 #ifdef NOXDR
 370       end subroutine cartout
 371 #else
 372       end subroutine cartoutx
 373 #endif
 374 !-----------------------------------------------------------------------------
 375 #ifndef NOXDR
 376       subroutine cartout(time)
 377 !      implicit real*8 (a-h,o-z)
 378 !      include 'DIMENSIONS'
 379       use geometry_data, only: c,nres
 380       use energy_data
 381       use MD_data, only: potE,t_bath
 382 #ifdef MPI
 383       use MPI_data
 384       include 'mpif.h'
 385 !      include 'COMMON.SETUP'
 386 #else
 387       integer,parameter :: me=0
 388 #endif
 389 !      include 'COMMON.CHAIN'
 390 !      include 'COMMON.INTERACT'
 391 !      include 'COMMON.NAMES'
 392 !      include 'COMMON.IOUNITS'
 393 !      include 'COMMON.HEADER'
 394 !      include 'COMMON.SBRIDGE'
 395 !      include 'COMMON.DISTFIT'
 396 !      include 'COMMON.MD'
 397       real(kind=8) :: time
 398       integer :: iret,itmp
 399       real(kind=4) :: prec
 400       real(kind=4),dimension(3,2*nres+2) :: xcoord      !(3,maxres2+2)  (maxres2=2*maxres
 401 !el  local variables
 402       integer :: j,i,ixdrf
 403
 404 #ifdef AIX
 405       call xdrfopen_(ixdrf,cartname, "a", iret)
 406       call xdrffloat_(ixdrf, real(time), iret)
 407       call xdrffloat_(ixdrf, real(potE), iret)
 408       call xdrffloat_(ixdrf, real(uconst), iret)
 409       call xdrffloat_(ixdrf, real(uconst_back), iret)
 410       call xdrffloat_(ixdrf, real(t_bath), iret)
 411       call xdrfint_(ixdrf, nss, iret)
 412       do j=1,nss
 413        if (dyn_ss) then
 414         call xdrfint_(ixdrf, idssb(j)+nres, iret)
 415         call xdrfint_(ixdrf, jdssb(j)+nres, iret)
 416        else
 417         call xdrfint_(ixdrf, ihpb(j), iret)
 418         call xdrfint_(ixdrf, jhpb(j), iret)
 419        endif
 420       enddo
 421       call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
 422       do i=1,nfrag
 423         call xdrffloat_(ixdrf, real(qfrag(i)), iret)
 424       enddo
 425       do i=1,npair
 426         call xdrffloat_(ixdrf, real(qpair(i)), iret)
 427       enddo
 428       do i=1,nfrag_back
 429         call xdrffloat_(ixdrf, real(utheta(i)), iret)
 430         call xdrffloat_(ixdrf, real(ugamma(i)), iret)
 431         call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
 432       enddo
 433 #else
 434       call xdrfopen(ixdrf,cartname, "a", iret)
 435       call xdrffloat(ixdrf, real(time), iret)
 436       call xdrffloat(ixdrf, real(potE), iret)
 437       call xdrffloat(ixdrf, real(uconst), iret)
 438       call xdrffloat(ixdrf, real(uconst_back), iret)
 439       call xdrffloat(ixdrf, real(t_bath), iret)
 440       call xdrfint(ixdrf, nss, iret)
 441       do j=1,nss
 442        if (dyn_ss) then
 443         call xdrfint(ixdrf, idssb(j)+nres, iret)
 444         call xdrfint(ixdrf, jdssb(j)+nres, iret)
 445        else
 446         call xdrfint(ixdrf, ihpb(j), iret)
 447         call xdrfint(ixdrf, jhpb(j), iret)
 448        endif
 449       enddo
 450       call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
 451       do i=1,nfrag
 452         call xdrffloat(ixdrf, real(qfrag(i)), iret)
 453       enddo
 454       do i=1,npair
 455         call xdrffloat(ixdrf, real(qpair(i)), iret)
 456       enddo
 457       do i=1,nfrag_back
 458         call xdrffloat(ixdrf, real(utheta(i)), iret)
 459         call xdrffloat(ixdrf, real(ugamma(i)), iret)
 460         call xdrffloat(ixdrf, real(uscdiff(i)), iret)
 461       enddo
 462 #endif
 463       prec=10000.0
 464       do i=1,nres
 465        do j=1,3
 466         xcoord(j,i)=c(j,i)
 467        enddo
 468       enddo
 469       do i=nnt,nct
 470        do j=1,3
 471         xcoord(j,nres+i-nnt+1)=c(j,i+nres)
 472        enddo
 473       enddo
 474
 475       itmp=nres+nct-nnt+1
 476 #ifdef AIX
 477       call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
 478       call xdrfclose_(ixdrf, iret)
 479 #else
 480       call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
 481       call xdrfclose(ixdrf, iret)
 482 #endif
 483       return
 484       end subroutine cartout
 485 #endif
 486 !-----------------------------------------------------------------------------
 487       subroutine statout(itime)
 488
 489       use energy_data
 490       use control_data
 491       use MD_data
 492       use MPI_data
 493       use compare, only:rms_nac_nnc
 494 !      implicit real*8 (a-h,o-z)
 495 !      include 'DIMENSIONS'
 496 !      include 'COMMON.CONTROL'
 497 !      include 'COMMON.CHAIN'
 498 !      include 'COMMON.INTERACT'
 499 !      include 'COMMON.NAMES'
 500 !      include 'COMMON.IOUNITS'
 501 !      include 'COMMON.HEADER'
 502 !      include 'COMMON.SBRIDGE'
 503 !      include 'COMMON.DISTFIT'
 504 !      include 'COMMON.MD'
 505 !      include 'COMMON.REMD'
 506 !      include 'COMMON.SETUP'
 507       integer :: itime
 508       real(kind=8),dimension(0:n_ene) :: energia
 509 !      double precision gyrate
 510 !el      external gyrate
 511 !el      common /gucio/ cm
 512       character(len=256) :: line1,line2
 513       character(len=4) :: format1,format2
 514       character(len=30) :: format
 515 !el  local variables
 516       integer :: i,ii1,ii2,j
 517       real(kind=8) :: rms,frac,frac_nn,co,distance
 518
 519 #ifdef AIX
 520       if(itime.eq.0) then
 521        open(istat,file=statname,position="append")
 522       endif
 523 #else
 524 #ifdef PGI
 525       open(istat,file=statname,position="append")
 526 #else
 527       open(istat,file=statname,access="append")
 528 #endif
 529 #endif
 530        if (AFMlog.gt.0) then
 531        if (refstr) then
 532          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 533           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')&
 534                itime,totT,EK,potE,totE,&
 535                rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
 536                potEcomp(23),me
 537           format1="a133"
 538          else
 539 !C          print *,'A CHUJ',potEcomp(23)
 540           write (line1,'(i10,f15.2,7f12.3,i5,$)') &
 541                 itime,totT,EK,potE,totE,&
 542                 kinetic_T,t_bath,gyrate(),&
 543                 potEcomp(23),me
 544           format1="a114"
 545         endif
 546        else if (selfguide.gt.0) then
 547        distance=0.0
 548        do j=1,3
 549        distance=distance+(c(j,afmend)-c(j,afmbeg))**2
 550        enddo
 551        distance=dsqrt(distance)
 552        if (refstr) then
 553          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 554           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, &
 555          f9.3,i5,$)') &
 556                itime,totT,EK,potE,totE,&
 557                rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
 558                distance,potEcomp(23),me
 559           format1="a133"
 560 !C          print *,"CHUJOWO"
 561          else
 562 !C          print *,'A CHUJ',potEcomp(23)
 563           write (line1,'(i10,f15.2,8f12.3,i5,$)')&
 564                 itime,totT,EK,potE,totE, &
 565                 kinetic_T,t_bath,gyrate(),&
 566                 distance,potEcomp(23),me
 567           format1="a114"
 568         endif
 569        else
 570        if (refstr) then
 571          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 572           write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') &
 573                 itime,totT,EK,potE,totE,&
 574                 rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
 575           format1="a133"
 576         else
 577           write (line1,'(i10,f15.2,7f12.3,i5,$)') &
 578                  itime,totT,EK,potE,totE,&
 579                  amax,kinetic_T,t_bath,gyrate(),me
 580           format1="a114"
 581         endif
 582         ENDIF
 583         if(usampl.and.totT.gt.eq_time) then
 584            write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,&
 585             (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),&
 586             (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 587            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair &
 588                    +21*nfrag_back
 589         else
 590            format2="a001"
 591            line2=' '
 592         endif
 593         if (print_compon) then
 594           if(itime.eq.0) then
 595            write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
 596                                                            ",20a12)"
 597            write (istat,format) "#","",&
 598             (ename(print_order(i)),i=1,nprint_ene)
 599           endif
 600           write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
 601                                                            ",20f12.3)"
 602           write (istat,format) line1,line2,&
 603             (potEcomp(print_order(i)),i=1,nprint_ene)
 604         else
 605           write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
 606           write (istat,format) line1,line2
 607         endif
 608 #if defined(AIX)
 609         call flush(istat)
 610 #else
 611         close(istat)
 612 #endif
 613       return
 614       end subroutine  statout
 615 !-----------------------------------------------------------------------------
 616 ! readrtns_CSA.F
 617 !-----------------------------------------------------------------------------
 618       subroutine readrtns
 619
 620       use control_data
 621       use energy_data
 622       use MPI_data
 623       use muca_md, only:read_muca
 624 !      implicit real*8 (a-h,o-z)
 625 !      include 'DIMENSIONS'
 626 #ifdef MPI
 627       include 'mpif.h'
 628 #endif
 629 !      include 'COMMON.SETUP'
 630 !      include 'COMMON.CONTROL'
 631 !      include 'COMMON.SBRIDGE'
 632 !      include 'COMMON.IOUNITS'
 633       logical :: file_exist
 634       integer :: i
 635 ! Read force-field parameters except weights
 636       call parmread
 637 ! Read job setup parameters
 638       call read_control
 639 ! Read control parameters for energy minimzation if required
 640       if (minim) call read_minim
 641 ! Read MCM control parameters if required
 642       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
 643 ! Read MD control parameters if reqjuired
 644       if (modecalc.eq.12) call read_MDpar
 645 ! Read MREMD control parameters if required
 646       if (modecalc.eq.14) then
 647          call read_MDpar
 648          call read_REMDpar
 649       endif
 650 ! Read MUCA control parameters if required
 651       if (lmuca) call read_muca
 652 ! Read CSA control parameters if required (from fort.40 if exists
 653 ! otherwise from general input file)
 654       if (modecalc.eq.8) then
 655        inquire (file="fort.40",exist=file_exist)
 656        if (.not.file_exist) call csaread
 657       endif
 658 !fmc      if (modecalc.eq.10) call mcmfread
 659 ! Read molecule information, molecule geometry, energy-term weights, and
 660 ! restraints if requested
 661       call molread
 662 ! Print restraint information
 663 #ifdef MPI
 664       if (.not. out1file .or. me.eq.king) then
 665 #endif
 666       if (nhpb.gt.nss) &
 667       write (iout,'(a,i5,a)') "The following",nhpb-nss,&
 668        " distance constraints have been imposed"
 669       do i=nss+1,nhpb
 670         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
 671       enddo
 672 #ifdef MPI
 673       endif
 674 #endif
 675 !      print *,"Processor",myrank," leaves READRTNS"
 676 !      write(iout,*) "end readrtns"
 677       return
 678       end subroutine readrtns
 679 !-----------------------------------------------------------------------------
 680       subroutine molread
 681 !
 682 ! Read molecular data.
 683 !
 684 !      use control, only: ilen
 685       use control_data
 686       use geometry_data
 687       use energy_data
 688       use energy
 689       use compare_data
 690       use MD_data, only: t_bath
 691       use MPI_data
 692       use compare, only:seq_comp,contact
 693       use control
 694 !      implicit real*8 (a-h,o-z)
 695 !      include 'DIMENSIONS'
 696 #ifdef MPI
 697       include 'mpif.h'
 698       integer :: error_msg,ierror,ierr,ierrcode
 699 #endif
 700 !      include 'COMMON.IOUNITS'
 701 !      include 'COMMON.GEO'
 702 !      include 'COMMON.VAR'
 703 !      include 'COMMON.INTERACT'
 704 !      include 'COMMON.LOCAL'
 705 !      include 'COMMON.NAMES'
 706 !      include 'COMMON.CHAIN'
 707 !      include 'COMMON.FFIELD'
 708 !      include 'COMMON.SBRIDGE'
 709 !      include 'COMMON.HEADER'
 710 !      include 'COMMON.CONTROL'
 711 !      include 'COMMON.DBASE'
 712 !      include 'COMMON.THREAD'
 713 !      include 'COMMON.CONTACTS'
 714 !      include 'COMMON.TORCNSTR'
 715 !      include 'COMMON.TIME1'
 716 !      include 'COMMON.BOUNDS'
 717 !      include 'COMMON.MD'
 718 !      include 'COMMON.SETUP'
 719       character(len=4),dimension(:,:),allocatable :: sequence   !(maxres,maxmolecules)
 720 !      integer :: rescode
 721 !      double precision x(maxvar)
 722       character(len=256) :: pdbfile
 723       character(len=320) :: weightcard
 724       character(len=80) :: weightcard_t!,ucase
 725 !      integer,dimension(:),allocatable :: itype_pdb    !(maxres)
 726 !      common /pizda/ itype_pdb
 727       logical :: fail   !seq_comp,
 728       real(kind=8) :: energia(0:n_ene)
 729 !      integer ilen
 730 !el      external ilen
 731 !el local varables
 732       integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2
 733
 734       real(kind=8),dimension(3,maxres2+2) :: c_alloc
 735       real(kind=8),dimension(3,0:maxres2) :: dc_alloc
 736       integer,dimension(maxres) :: itype_alloc
 737
 738       integer :: iti,nsi,maxsi,itrial,itmp
 739       real(kind=8) :: wlong,scalscp,co,ssscale
 740       allocate(weights(n_ene))
 741 !-----------------------------
 742       allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
 743       allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
 744       allocate(itype(maxres,5)) !(maxres)
 745       allocate(istype(maxres))
 746 !
 747 ! Zero out tables.
 748 !
 749       c(:,:)=0.0D0
 750       dc(:,:)=0.0D0
 751       itype(:,:)=0
 752 !-----------------------------
 753 !
 754 ! Body
 755 !
 756 ! Read weights of the subsequent energy terms.
 757       call card_concat(weightcard,.true.)
 758       call reada(weightcard,'WLONG',wlong,1.0D0)
 759       call reada(weightcard,'WSC',wsc,wlong)
 760       call reada(weightcard,'WSCP',wscp,wlong)
 761       call reada(weightcard,'WELEC',welec,1.0D0)
 762       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 763       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 764       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 765       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 766       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 767       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 768       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 769       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 770       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 771       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 772       call reada(weightcard,'WBOND',wbond,1.0D0)
 773       call reada(weightcard,'WTOR',wtor,1.0D0)
 774       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 775       call reada(weightcard,'WSHIELD',wshield,0.05D0)
 776       call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
 777       call reada(weightcard,'WLT',wliptran,1.0D0)
 778       call reada(weightcard,'WTUBE',wtube,1.0d0)
 779       call reada(weightcard,'WANG',wang,1.0D0)
 780       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 781       call reada(weightcard,'SCAL14',scal14,0.4D0)
 782       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 783       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 784       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 785       call reada(weightcard,'TEMP0',temp0,300.0d0)
 786       if (index(weightcard,'SOFT').gt.0) ipot=6
 787 ! 12/1/95 Added weight for the multi-body term WCORR
 788       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 789       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 790       weights(1)=wsc
 791       weights(2)=wscp
 792       weights(3)=welec
 793       weights(4)=wcorr
 794       weights(5)=wcorr5
 795       weights(6)=wcorr6
 796       weights(7)=wel_loc
 797       weights(8)=wturn3
 798       weights(9)=wturn4
 799       weights(10)=wturn6
 800       weights(11)=wang
 801       weights(12)=wscloc
 802       weights(13)=wtor
 803       weights(14)=wtor_d
 804       weights(15)=wstrain
 805       weights(16)=wvdwpp
 806       weights(17)=wbond
 807       weights(18)=scal14
 808       weights(21)=wsccor
 809       if(me.eq.king.or..not.out1file) &
 810        write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 811         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 812         wturn4,wturn6
 813    10 format (/'Energy-term weights (unscaled):'// &
 814        'WSCC=   ',f10.6,' (SC-SC)'/ &
 815        'WSCP=   ',f10.6,' (SC-p)'/ &
 816        'WELEC=  ',f10.6,' (p-p electr)'/ &
 817        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 818        'WBOND=  ',f10.6,' (stretching)'/ &
 819        'WANG=   ',f10.6,' (bending)'/ &
 820        'WSCLOC= ',f10.6,' (SC local)'/ &
 821        'WTOR=   ',f10.6,' (torsional)'/ &
 822        'WTORD=  ',f10.6,' (double torsional)'/ &
 823        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 824        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 825        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 826        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 827        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 828        'WSCCOR= ',f10.6,' (back-scloc correlation)'/ &
 829        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 830        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 831        'WTURN6= ',f10.6,' (turns, 6th order)')
 832       if(me.eq.king.or..not.out1file)then
 833        if (wcorr4.gt.0.0d0) then
 834         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
 835          'between contact pairs of peptide groups'
 836         write (iout,'(2(a,f5.3/))') &
 837         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
 838         'Range of quenching the correlation terms:',2*delt_corr
 839        else if (wcorr.gt.0.0d0) then
 840         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
 841          'between contact pairs of peptide groups'
 842        endif
 843        write (iout,'(a,f8.3)') &
 844         'Scaling factor of 1,4 SC-p interactions:',scal14
 845        write (iout,'(a,f8.3)') &
 846         'General scaling factor of SC-p interactions:',scalscp
 847       endif
 848       r0_corr=cutoff_corr-delt_corr
 849       do i=1,ntyp
 850         aad(i,1)=scalscp*aad(i,1)
 851         aad(i,2)=scalscp*aad(i,2)
 852         bad(i,1)=scalscp*bad(i,1)
 853         bad(i,2)=scalscp*bad(i,2)
 854       enddo
 855       call rescale_weights(t_bath)
 856       if(me.eq.king.or..not.out1file) &
 857        write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 858         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 859         wturn4,wturn6
 860    22 format (/'Energy-term weights (scaled):'// &
 861        'WSCC=   ',f10.6,' (SC-SC)'/ &
 862        'WSCP=   ',f10.6,' (SC-p)'/ &
 863        'WELEC=  ',f10.6,' (p-p electr)'/ &
 864        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 865        'WBOND=  ',f10.6,' (stretching)'/ &
 866        'WANG=   ',f10.6,' (bending)'/ &
 867        'WSCLOC= ',f10.6,' (SC local)'/ &
 868        'WTOR=   ',f10.6,' (torsional)'/ &
 869        'WTORD=  ',f10.6,' (double torsional)'/ &
 870        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 871        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 872        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 873        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 874        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 875        'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ &
 876        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 877        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 878        'WTURN6= ',f10.6,' (turns, 6th order)')
 879       if(me.eq.king.or..not.out1file) &
 880        write (iout,*) "Reference temperature for weights calculation:",&
 881         temp0
 882       call reada(weightcard,"D0CM",d0cm,3.78d0)
 883       call reada(weightcard,"AKCM",akcm,15.1d0)
 884       call reada(weightcard,"AKTH",akth,11.0d0)
 885       call reada(weightcard,"AKCT",akct,12.0d0)
 886       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 887       call reada(weightcard,"V2SS",v2ss,7.61d0)
 888       call reada(weightcard,"V3SS",v3ss,13.7d0)
 889       call reada(weightcard,"EBR",ebr,-5.50D0)
 890       call reada(weightcard,"ATRISS",atriss,0.301D0)
 891       call reada(weightcard,"BTRISS",btriss,0.021D0)
 892       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 893       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 894       call reada(weightcard,"SSSCALE",ssscale,1.0D0)
 895       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 896
 897       call reada(weightcard,"HT",Ht,0.0D0)
 898       if (dyn_ss) then
 899        ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
 900         Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
 901         akcm=akcm/wsc*ssscale
 902         akth=akth/wsc*ssscale
 903         akct=akct/wsc*ssscale
 904         v1ss=v1ss/wsc*ssscale
 905         v2ss=v2ss/wsc*ssscale
 906         v3ss=v3ss/wsc*ssscale
 907       else
 908         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 909       endif
 910
 911       if(me.eq.king.or..not.out1file) then
 912        write (iout,*) "Parameters of the SS-bond potential:"
 913        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,&
 914        " AKCT",akct
 915        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 916        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 917        write (iout,*)" HT",Ht
 918        print *,'indpdb=',indpdb,' pdbref=',pdbref
 919       endif
 920       if (indpdb.gt.0 .or. pdbref) then
 921         read(inp,'(a)') pdbfile
 922         if(me.eq.king.or..not.out1file) &
 923          write (iout,'(2a)') 'PDB data will be read from file ',&
 924          pdbfile(:ilen(pdbfile))
 925         open(ipdbin,file=pdbfile,status='old',err=33)
 926         goto 34
 927   33    write (iout,'(a)') 'Error opening PDB file.'
 928         stop
 929   34    continue
 930 !        print *,'Begin reading pdb data'
 931         call readpdb
 932 !        print *,'Finished reading pdb data'
 933         if(me.eq.king.or..not.out1file) &
 934          write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
 935          ' nstart_sup=',nstart_sup !,"ergwergewrgae"
 936 !el        if(.not.allocated(itype_pdb))
 937         allocate(itype_pdb(nres))
 938         do i=1,nres
 939           itype_pdb(i)=itype(i,1)
 940         enddo
 941         close (ipdbin)
 942         nnt=nstart_sup
 943         nct=nstart_sup+nsup-1
 944 !el        if(.not.allocated(icont_ref))
 945         allocate(icont_ref(2,12*nres)) ! maxcont=12*maxres
 946         call contact(.false.,ncont_ref,icont_ref,co)
 947
 948         if (sideadd) then
 949          if(me.eq.king.or..not.out1file) &
 950           write(iout,*)'Adding sidechains'
 951          maxsi=1000
 952          do i=2,nres-1
 953           iti=itype(i,1)
 954           if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then
 955             nsi=0
 956             fail=.true.
 957             do while (fail.and.nsi.le.maxsi)
 958               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 959               nsi=nsi+1
 960             enddo
 961             if(fail) write(iout,*)'Adding sidechain failed for res ',&
 962                     i,' after ',nsi,' trials'
 963           endif
 964          enddo
 965         endif
 966       endif
 967
 968       if (indpdb.eq.0) then
 969       nres_molec(:)=0
 970         allocate(sequence(maxres,5))
 971
 972      if (protein) then
 973 ! Read sequence if not taken from the pdb file.
 974         molec=1
 975         read (inp,*) nres_molec(molec)
 976 !        print *,'nres=',nres
 977         if (iscode.gt.0) then
 978           read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec))
 979         else
 980           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
 981         endif
 982 ! Convert sequence to numeric code
 983
 984         do i=1,nres_molec(molec)
 985           itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
 986         enddo
 987        endif
 988        if (nucleic) then
 989 ! Read sequence if not taken from the pdb file.
 990         molec=2
 991         read (inp,*) nres_molec(molec)
 992 !        print *,'nres=',nres
 993 !        allocate(sequence(maxres,5))
 994 !        if (iscode.gt.0) then
 995           read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
 996 ! Convert sequence to numeric code
 997
 998         do i=1,nres_molec(molec)
 999           istype(i)=sugarcode(sequence(i,molec)(1:1),i)
1000           itype(i,1)=rescode(i,sequence(i,molec)(2:2),iscode,molec)
1001         enddo
1002        endif
1003
1004        if (ions) then
1005 ! Read sequence if not taken from the pdb file.
1006         molec=5
1007         read (inp,*) nres_molec(molec)
1008 !        print *,'nres=',nres
1009           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1010 ! Convert sequence to numeric code
1011         print *,nres_molec(molec)
1012         do i=1,nres_molec(molec)
1013           itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
1014         enddo
1015        endif
1016        nres=0
1017        do i=1,5
1018         nres=nres+nres_molec(i)
1019         print *,nres_molec(i)
1020        enddo
1021
1022 ! Assign initial virtual bond lengths
1023 !elwrite(iout,*) "test_alloc"
1024         if(.not.allocated(vbld)) allocate(vbld(2*nres))
1025 !elwrite(iout,*) "test_alloc"
1026         if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
1027 !elwrite(iout,*) "test_alloc"
1028         do i=2,nres
1029           vbld(i)=vbl
1030           vbld_inv(i)=vblinv
1031         enddo
1032         do i=2,nres-1
1033           vbld(i+nres)=dsc(iabs(itype(i,1)))
1034           vbld_inv(i+nres)=dsc_inv(iabs(itype(i,1)))
1035 !          write (iout,*) "i",i," itype",itype(i,1),
1036 !     &      " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
1037         enddo
1038       endif
1039 !      print *,nres
1040 !      print '(20i4)',(itype(i,1),i=1,nres)
1041 !----------------------------
1042 !el reallocate tables
1043 !      do i=1,maxres2
1044 !        do j=1,3
1045 !          c_alloc(j,i)=c(j,i)
1046 !          dc_alloc(j,i)=dc(j,i)
1047 !        enddo
1048 !      enddo
1049 !      do i=1,maxres
1050 !elwrite(iout,*) "itype",i,itype(i,1)
1051 !        itype_alloc(i)=itype(i,1)
1052 !      enddo
1053
1054 !      deallocate(c)
1055 !      deallocate(dc)
1056 !      deallocate(itype)
1057 !      allocate(c(3,2*nres+4))
1058 !      allocate(dc(3,0:2*nres+2))
1059 !      allocate(itype(nres+2))
1060       allocate(itel(nres+2))
1061       itel(:)=0
1062
1063 !      do i=1,2*nres+2
1064 !        do j=1,3
1065 !          c(j,i)=c_alloc(j,i)
1066 !          dc(j,i)=dc_alloc(j,i)
1067 !        enddo
1068 !      enddo
1069 !      do i=1,nres+2
1070 !        itype(i,1)=itype_alloc(i)
1071 !        itel(i)=0
1072 !      enddo
1073 !--------------------------
1074       do i=1,nres
1075 #ifdef PROCOR
1076         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then
1077 #else
1078         if (itype(i,1).eq.ntyp1) then
1079 #endif
1080           itel(i)=0
1081 #ifdef PROCOR
1082         else if (iabs(itype(i+1,1)).ne.20) then
1083 #else
1084         else if (iabs(itype(i,1)).ne.20) then
1085 #endif
1086           itel(i)=1
1087         else
1088           itel(i)=2
1089         endif
1090       enddo
1091       if(me.eq.king.or..not.out1file)then
1092        write (iout,*) "ITEL"
1093        do i=1,nres-1
1094          write (iout,*) i,itype(i,1),itel(i)
1095        enddo
1096        print *,'Call Read_Bridge.'
1097       endif
1098       call read_bridge
1099 !--------------------------------
1100 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1101       call alloc_geo_arrays
1102       call alloc_ener_arrays
1103 !--------------------------------
1104 ! 8/13/98 Set limits to generating the dihedral angles
1105       do i=1,nres
1106         phibound(1,i)=-pi
1107         phibound(2,i)=pi
1108       enddo
1109       read (inp,*) ndih_constr
1110       if (ndih_constr.gt.0) then
1111         allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1112         allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1113
1114         read (inp,*) ftors
1115         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1116         if(me.eq.king.or..not.out1file)then
1117          write (iout,*) &
1118          'There are',ndih_constr,' constraints on phi angles.'
1119          do i=1,ndih_constr
1120           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1121          enddo
1122         endif
1123         do i=1,ndih_constr
1124           phi0(i)=deg2rad*phi0(i)
1125           drange(i)=deg2rad*drange(i)
1126         enddo
1127         if(me.eq.king.or..not.out1file) &
1128          write (iout,*) 'FTORS',ftors
1129         do i=1,ndih_constr
1130           ii = idih_constr(i)
1131           phibound(1,ii) = phi0(i)-drange(i)
1132           phibound(2,ii) = phi0(i)+drange(i)
1133         enddo
1134       endif
1135       if (with_theta_constr) then
1136 !C with_theta_constr is keyword allowing for occurance of theta constrains
1137       read (inp,*) ntheta_constr
1138 !C ntheta_constr is the number of theta constrains
1139       if (ntheta_constr.gt.0) then
1140 !C        read (inp,*) ftors
1141         allocate(itheta_constr(ntheta_constr))
1142         allocate(theta_constr0(ntheta_constr))
1143         allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr))
1144         read (inp,*) (itheta_constr(i),theta_constr0(i), &
1145        theta_drange(i),for_thet_constr(i), &
1146        i=1,ntheta_constr)
1147 !C the above code reads from 1 to ntheta_constr
1148 !C itheta_constr(i) residue i for which is theta_constr
1149 !C theta_constr0 the global minimum value
1150 !C theta_drange is range for which there is no energy penalty
1151 !C for_thet_constr is the force constant for quartic energy penalty
1152 !C E=k*x**4
1153         if(me.eq.king.or..not.out1file)then
1154          write (iout,*) &
1155         'There are',ntheta_constr,' constraints on phi angles.'
1156          do i=1,ntheta_constr
1157           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), &
1158          theta_drange(i), &
1159          for_thet_constr(i)
1160          enddo
1161         endif
1162         do i=1,ntheta_constr
1163           theta_constr0(i)=deg2rad*theta_constr0(i)
1164           theta_drange(i)=deg2rad*theta_drange(i)
1165         enddo
1166 !C        if(me.eq.king.or..not.out1file)
1167 !C     &   write (iout,*) 'FTORS',ftors
1168 !C        do i=1,ntheta_constr
1169 !C          ii = itheta_constr(i)
1170 !C          thetabound(1,ii) = phi0(i)-drange(i)
1171 !C          thetabound(2,ii) = phi0(i)+drange(i)
1172 !C        enddo
1173       endif ! ntheta_constr.gt.0
1174       endif! with_theta_constr
1175
1176       nnt=1
1177 #ifdef MPI
1178       if (me.eq.king) then
1179 #endif
1180        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1181        do i=1,nres
1182          write (iout,'(a3,i5,2f10.1)') &
1183          restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1184        enddo
1185 #ifdef MP
1186       endif
1187 #endif
1188       nct=nres
1189 !d      print *,'NNT=',NNT,' NCT=',NCT
1190       if (itype(1,1).eq.ntyp1) nnt=2
1191       if (itype(nres,1).eq.ntyp1) nct=nct-1
1192       if (pdbref) then
1193         if(me.eq.king.or..not.out1file) &
1194          write (iout,'(a,i3)') 'nsup=',nsup
1195         nstart_seq=nnt
1196         if (nsup.le.(nct-nnt+1)) then
1197           do i=0,nct-nnt+1-nsup
1198             if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
1199               nstart_seq=nnt+i
1200               goto 111
1201             endif
1202           enddo
1203           write (iout,'(a)') &
1204                   'Error - sequences to be superposed do not match.'
1205           stop
1206         else
1207           do i=0,nsup-(nct-nnt+1)
1208             if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1209             then
1210               nstart_sup=nstart_sup+i
1211               nsup=nct-nnt+1
1212               goto 111
1213             endif
1214           enddo
1215           write (iout,'(a)') &
1216                   'Error - sequences to be superposed do not match.'
1217         endif
1218   111   continue
1219         if (nsup.eq.0) nsup=nct-nnt
1220         if (nstart_sup.eq.0) nstart_sup=nnt
1221         if (nstart_seq.eq.0) nstart_seq=nnt
1222         if(me.eq.king.or..not.out1file) &
1223          write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1224                        ' nstart_seq=',nstart_seq !,"242343453254"
1225       endif
1226 !--- Zscore rms -------
1227       if (nz_start.eq.0) nz_start=nnt
1228       if (nz_end.eq.0 .and. nsup.gt.0) then
1229         nz_end=nnt+nsup-1
1230       else if (nz_end.eq.0) then
1231         nz_end=nct
1232       endif
1233       if(me.eq.king.or..not.out1file)then
1234        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1235        write (iout,*) 'IZ_SC=',iz_sc
1236       endif
1237 !----------------------
1238       call init_int_table
1239       if (refstr) then
1240         if (.not.pdbref) then
1241           call read_angles(inp,*38)
1242           goto 39
1243    38     write (iout,'(a)') 'Error reading reference structure.'
1244 #ifdef MPI
1245           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1246           stop 'Error reading reference structure'
1247 #endif
1248    39     call chainbuild
1249           call setup_var
1250 !zscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1251           nstart_sup=nnt
1252           nstart_seq=nnt
1253           nsup=nct-nnt+1
1254           kkk=1
1255           do i=1,2*nres
1256             do j=1,3
1257               cref(j,i,kkk)=c(j,i)
1258             enddo
1259           enddo
1260           call contact(.true.,ncont_ref,icont_ref,co)
1261         endif
1262 !        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1263         call flush(iout)
1264         if (constr_dist.gt.0) call read_dist_constr
1265         write (iout,*) "After read_dist_constr nhpb",nhpb
1266         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1267         call hpb_partition
1268         if(me.eq.king.or..not.out1file) &
1269          write (iout,*) 'Contact order:',co
1270         if (pdbref) then
1271         if(me.eq.king.or..not.out1file) &
1272          write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1273         do i=1,ncont_ref
1274           do j=1,2
1275             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1276           enddo
1277           if(me.eq.king.or..not.out1file) &
1278            write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1279            icont_ref(1,i),' ',&
1280            restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1281         enddo
1282         endif
1283       endif
1284       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1285           .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1286           modecalc.ne.10) then
1287 ! If input structure hasn't been supplied from the PDB file read or generate
1288 ! initial geometry.
1289         if (iranconf.eq.0 .and. .not. extconf) then
1290           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1291            write (iout,'(a)') 'Initial geometry will be read in.'
1292           if (read_cart) then
1293             read(inp,'(8f10.5)',end=36,err=36) &
1294              ((c(l,k),l=1,3),k=1,nres),&
1295              ((c(l,k+nres),l=1,3),k=nnt,nct)
1296             write (iout,*) "Exit READ_CART"
1297             write (iout,'(8f10.5)') &
1298              ((c(l,k),l=1,3),k=1,nres),&
1299              ((c(l,k+nres),l=1,3),k=nnt,nct)
1300             call int_from_cart1(.true.)
1301             write (iout,*) "Finish INT_TO_CART"
1302             do i=1,nres-1
1303               do j=1,3
1304                 dc(j,i)=c(j,i+1)-c(j,i)
1305                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1306               enddo
1307             enddo
1308             do i=nnt,nct
1309               if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1310                 do j=1,3
1311                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1312                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1313                 enddo
1314               endif
1315             enddo
1316             return
1317           else
1318             call read_angles(inp,*36)
1319           endif
1320           goto 37
1321    36     write (iout,'(a)') 'Error reading angle file.'
1322 #ifdef MPI
1323           call mpi_finalize( MPI_COMM_WORLD,IERR )
1324 #endif
1325           stop 'Error reading angle file.'
1326    37     continue
1327         else if (extconf) then
1328          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1329           write (iout,'(a)') 'Extended chain initial geometry.'
1330          do i=3,nres
1331           theta(i)=90d0*deg2rad
1332          enddo
1333          do i=4,nres
1334           phi(i)=180d0*deg2rad
1335          enddo
1336          do i=2,nres-1
1337           alph(i)=110d0*deg2rad
1338          enddo
1339 !elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
1340          do i=2,nres-1
1341           omeg(i)=-120d0*deg2rad
1342           if (itype(i,1).le.0) omeg(i)=-omeg(i)
1343          enddo
1344         else
1345           if(me.eq.king.or..not.out1file) &
1346            write (iout,'(a)') 'Random-generated initial geometry.'
1347
1348
1349 #ifdef MPI
1350           if (me.eq.king  .or. fg_rank.eq.0 .and. &
1351                  ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1352 #endif
1353             do itrial=1,100
1354               itmp=1
1355               call gen_rand_conf(itmp,*30)
1356               goto 40
1357    30         write (iout,*) 'Failed to generate random conformation',&
1358                 ', itrial=',itrial
1359               write (*,*) 'Processor:',me,&
1360                 ' Failed to generate random conformation',&
1361                 ' itrial=',itrial
1362               call intout
1363
1364 #ifdef AIX
1365               call flush_(iout)
1366 #else
1367               call flush(iout)
1368 #endif
1369             enddo
1370             write (iout,'(a,i3,a)') 'Processor:',me,&
1371               ' error in generating random conformation.'
1372             write (*,'(a,i3,a)') 'Processor:',me,&
1373               ' error in generating random conformation.'
1374             call flush(iout)
1375 #ifdef MPI
1376             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1377    40       continue
1378           endif
1379 #else
1380           do itrial=1,100
1381             itmp=1
1382             call gen_rand_conf(itmp,*335)
1383             goto 40
1384   335       write (iout,*) 'Failed to generate random conformation',&
1385               ', itrial=',itrial
1386             write (*,*) 'Failed to generate random conformation',&
1387               ', itrial=',itrial
1388           enddo
1389           write (iout,'(a,i3,a)') 'Processor:',me,&
1390             ' error in generating random conformation.'
1391           write (*,'(a,i3,a)') 'Processor:',me,&
1392             ' error in generating random conformation.'
1393           stop
1394    40     continue
1395 #endif
1396         endif
1397       elseif (modecalc.eq.4) then
1398         read (inp,'(a)') intinname
1399         open (intin,file=intinname,status='old',err=333)
1400         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1401         write (iout,'(a)') 'intinname',intinname
1402         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1403         goto 334
1404   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1405 #ifdef MPI
1406         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1407 #endif
1408         stop 'Error opening angle file.'
1409   334   continue
1410
1411       endif
1412 ! Generate distance constraints, if the PDB structure is to be regularized.
1413       if (nthread.gt.0) then
1414         call read_threadbase
1415       endif
1416       call setup_var
1417 !elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
1418       if (me.eq.king .or. .not. out1file) &
1419        call intout
1420 !elwrite (iout,*)"alph(i)*rad2deg",(alph(i)*rad2deg, i=1,nres)
1421       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1422         write (iout,'(/a,i3,a)') &
1423         'The chain contains',ns,' disulfide-bridging cysteines.'
1424         write (iout,'(20i4)') (iss(i),i=1,ns)
1425        if (dyn_ss) then
1426           write(iout,*)"Running with dynamic disulfide-bond formation"
1427        else
1428         write (iout,'(/a/)') 'Pre-formed links are:'
1429         do i=1,nss
1430           i1=ihpb(i)-nres
1431           i2=jhpb(i)-nres
1432           it1=itype(i1,1)
1433           it2=itype(i2,1)
1434           if (me.eq.king.or..not.out1file) &
1435           write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1436           restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1437           ebr,forcon(i)
1438         enddo
1439         write (iout,'(a)')
1440        endif
1441       endif
1442       if (ns.gt.0.and.dyn_ss) then
1443           do i=nss+1,nhpb
1444             ihpb(i-nss)=ihpb(i)
1445             jhpb(i-nss)=jhpb(i)
1446             forcon(i-nss)=forcon(i)
1447             dhpb(i-nss)=dhpb(i)
1448           enddo
1449           nhpb=nhpb-nss
1450           nss=0
1451           call hpb_partition
1452           do i=1,ns
1453             dyn_ss_mask(iss(i))=.true.
1454           enddo
1455       endif
1456       if (i2ndstr.gt.0) call secstrp2dihc
1457 !      call geom_to_var(nvar,x)
1458 !      call etotal(energia(0))
1459 !      call enerprint(energia(0))
1460 !      call briefout(0,etot)
1461 !      stop
1462 !d    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1463 !d    write (iout,'(a)') 'Variable list:'
1464 !d    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1465 #ifdef MPI
1466       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1467         write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1468         'Processor',myrank,': end reading molecular data.'
1469 #endif
1470       return
1471       end subroutine molread
1472 !-----------------------------------------------------------------------------
1473 !-----------------------------------------------------------------------------
1474       end module io