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).ne.10) then
 224 !d             print *,'i=',i,' itype=',itype(i),' theta=',theta(i+1)
 225                fail=.true.
 226                ii=0
 227                do while (fail .and. ii .le. maxsi)
 228                  call gen_side(itype(i),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, only:refstr
 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
 517       real(kind=8) :: rms,frac,frac_nn,co
 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 (refstr) then
 531          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 532           write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') &
 533                 itime,totT,EK,potE,totE,&
 534                 rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
 535           format1="a133"
 536         else
 537           write (line1,'(i10,f15.2,7f12.3,i5,$)') &
 538                  itime,totT,EK,potE,totE,&
 539                  amax,kinetic_T,t_bath,gyrate(),me
 540           format1="a114"
 541         endif
 542         if(usampl.and.totT.gt.eq_time) then
 543            write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,&
 544             (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),&
 545             (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 546            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair &
 547                    +21*nfrag_back
 548         else
 549            format2="a001"
 550            line2=' '
 551         endif
 552         if (print_compon) then
 553           if(itime.eq.0) then
 554            write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
 555                                                            ",20a12)"
 556            write (istat,format) "#","",&
 557             (ename(print_order(i)),i=1,nprint_ene)
 558           endif
 559           write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
 560                                                            ",20f12.3)"
 561           write (istat,format) line1,line2,&
 562             (potEcomp(print_order(i)),i=1,nprint_ene)
 563         else
 564           write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
 565           write (istat,format) line1,line2
 566         endif
 567 #if defined(AIX)
 568         call flush(istat)
 569 #else
 570         close(istat)
 571 #endif
 572       return
 573       end subroutine  statout
 574 !-----------------------------------------------------------------------------
 575 ! readrtns_CSA.F
 576 !-----------------------------------------------------------------------------
 577       subroutine readrtns
 578
 579       use control_data
 580       use energy_data
 581       use MPI_data
 582       use muca_md, only:read_muca
 583 !      implicit real*8 (a-h,o-z)
 584 !      include 'DIMENSIONS'
 585 #ifdef MPI
 586       include 'mpif.h'
 587 #endif
 588 !      include 'COMMON.SETUP'
 589 !      include 'COMMON.CONTROL'
 590 !      include 'COMMON.SBRIDGE'
 591 !      include 'COMMON.IOUNITS'
 592       logical :: file_exist
 593       integer :: i
 594 ! Read force-field parameters except weights
 595       call parmread
 596 ! Read job setup parameters
 597       call read_control
 598 ! Read control parameters for energy minimzation if required
 599       if (minim) call read_minim
 600 ! Read MCM control parameters if required
 601       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
 602 ! Read MD control parameters if reqjuired
 603       if (modecalc.eq.12) call read_MDpar
 604 ! Read MREMD control parameters if required
 605       if (modecalc.eq.14) then
 606          call read_MDpar
 607          call read_REMDpar
 608       endif
 609 ! Read MUCA control parameters if required
 610       if (lmuca) call read_muca
 611 ! Read CSA control parameters if required (from fort.40 if exists
 612 ! otherwise from general input file)
 613       if (modecalc.eq.8) then
 614        inquire (file="fort.40",exist=file_exist)
 615        if (.not.file_exist) call csaread
 616       endif
 617 !fmc      if (modecalc.eq.10) call mcmfread
 618 ! Read molecule information, molecule geometry, energy-term weights, and
 619 ! restraints if requested
 620       call molread
 621 ! Print restraint information
 622 #ifdef MPI
 623       if (.not. out1file .or. me.eq.king) then
 624 #endif
 625       if (nhpb.gt.nss) &
 626       write (iout,'(a,i5,a)') "The following",nhpb-nss,&
 627        " distance constraints have been imposed"
 628       do i=nss+1,nhpb
 629         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
 630       enddo
 631 #ifdef MPI
 632       endif
 633 #endif
 634 !      print *,"Processor",myrank," leaves READRTNS"
 635 !      write(iout,*) "end readrtns"
 636       return
 637       end subroutine readrtns
 638 !-----------------------------------------------------------------------------
 639       subroutine molread
 640 !
 641 ! Read molecular data.
 642 !
 643 !      use control, only: ilen
 644       use control_data
 645       use geometry_data
 646       use energy_data
 647       use energy
 648       use compare_data
 649       use MD_data, only: t_bath
 650       use MPI_data
 651       use compare, only:seq_comp,contact
 652       use control
 653 !      implicit real*8 (a-h,o-z)
 654 !      include 'DIMENSIONS'
 655 #ifdef MPI
 656       include 'mpif.h'
 657       integer :: error_msg,ierror,ierr,ierrcode
 658 #endif
 659 !      include 'COMMON.IOUNITS'
 660 !      include 'COMMON.GEO'
 661 !      include 'COMMON.VAR'
 662 !      include 'COMMON.INTERACT'
 663 !      include 'COMMON.LOCAL'
 664 !      include 'COMMON.NAMES'
 665 !      include 'COMMON.CHAIN'
 666 !      include 'COMMON.FFIELD'
 667 !      include 'COMMON.SBRIDGE'
 668 !      include 'COMMON.HEADER'
 669 !      include 'COMMON.CONTROL'
 670 !      include 'COMMON.DBASE'
 671 !      include 'COMMON.THREAD'
 672 !      include 'COMMON.CONTACTS'
 673 !      include 'COMMON.TORCNSTR'
 674 !      include 'COMMON.TIME1'
 675 !      include 'COMMON.BOUNDS'
 676 !      include 'COMMON.MD'
 677 !      include 'COMMON.SETUP'
 678       character(len=4),dimension(:),allocatable :: sequence     !(maxres)
 679 !      integer :: rescode
 680 !      double precision x(maxvar)
 681       character(len=256) :: pdbfile
 682       character(len=320) :: weightcard
 683       character(len=80) :: weightcard_t!,ucase
 684 !      integer,dimension(:),allocatable :: itype_pdb    !(maxres)
 685 !      common /pizda/ itype_pdb
 686       logical :: fail   !seq_comp,
 687       real(kind=8) :: energia(0:n_ene)
 688 !      integer ilen
 689 !el      external ilen
 690 !el local varables
 691       integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2
 692
 693       real(kind=8),dimension(3,maxres2+2) :: c_alloc
 694       real(kind=8),dimension(3,0:maxres2) :: dc_alloc
 695       integer,dimension(maxres) :: itype_alloc
 696
 697       integer :: iti,nsi,maxsi,itrial,itmp
 698       real(kind=8) :: wlong,scalscp,co
 699       allocate(weights(n_ene))
 700 !-----------------------------
 701       allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
 702       allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
 703       allocate(itype(maxres)) !(maxres)
 704 !
 705 ! Zero out tables.
 706 !
 707       c(:,:)=0.0D0
 708       dc(:,:)=0.0D0
 709       itype(:)=0
 710 !-----------------------------
 711 !
 712 ! Body
 713 !
 714 ! Read weights of the subsequent energy terms.
 715       call card_concat(weightcard,.true.)
 716       call reada(weightcard,'WLONG',wlong,1.0D0)
 717       call reada(weightcard,'WSC',wsc,wlong)
 718       call reada(weightcard,'WSCP',wscp,wlong)
 719       call reada(weightcard,'WELEC',welec,1.0D0)
 720       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 721       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 722       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 723       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 724       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 725       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 726       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 727       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 728       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 729       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 730       call reada(weightcard,'WBOND',wbond,1.0D0)
 731       call reada(weightcard,'WTOR',wtor,1.0D0)
 732       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 733       call reada(weightcard,'WSHIELD',wshield,0.05D0)
 734       call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
 735       call reada(weightcard,'WLT',wliptran,1.0D0)
 736       call reada(weightcard,'WTUBE',wtube,1.0d0)
 737       call reada(weightcard,'WANG',wang,1.0D0)
 738       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 739       call reada(weightcard,'SCAL14',scal14,0.4D0)
 740       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 741       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 742       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 743       call reada(weightcard,'TEMP0',temp0,300.0d0)
 744       if (index(weightcard,'SOFT').gt.0) ipot=6
 745 ! 12/1/95 Added weight for the multi-body term WCORR
 746       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 747       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 748       weights(1)=wsc
 749       weights(2)=wscp
 750       weights(3)=welec
 751       weights(4)=wcorr
 752       weights(5)=wcorr5
 753       weights(6)=wcorr6
 754       weights(7)=wel_loc
 755       weights(8)=wturn3
 756       weights(9)=wturn4
 757       weights(10)=wturn6
 758       weights(11)=wang
 759       weights(12)=wscloc
 760       weights(13)=wtor
 761       weights(14)=wtor_d
 762       weights(15)=wstrain
 763       weights(16)=wvdwpp
 764       weights(17)=wbond
 765       weights(18)=scal14
 766       weights(21)=wsccor
 767       if(me.eq.king.or..not.out1file) &
 768        write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 769         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 770         wturn4,wturn6
 771    10 format (/'Energy-term weights (unscaled):'// &
 772        'WSCC=   ',f10.6,' (SC-SC)'/ &
 773        'WSCP=   ',f10.6,' (SC-p)'/ &
 774        'WELEC=  ',f10.6,' (p-p electr)'/ &
 775        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 776        'WBOND=  ',f10.6,' (stretching)'/ &
 777        'WANG=   ',f10.6,' (bending)'/ &
 778        'WSCLOC= ',f10.6,' (SC local)'/ &
 779        'WTOR=   ',f10.6,' (torsional)'/ &
 780        'WTORD=  ',f10.6,' (double torsional)'/ &
 781        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 782        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 783        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 784        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 785        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 786        'WSCCOR= ',f10.6,' (back-scloc correlation)'/ &
 787        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 788        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 789        'WTURN6= ',f10.6,' (turns, 6th order)')
 790       if(me.eq.king.or..not.out1file)then
 791        if (wcorr4.gt.0.0d0) then
 792         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
 793          'between contact pairs of peptide groups'
 794         write (iout,'(2(a,f5.3/))') &
 795         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
 796         'Range of quenching the correlation terms:',2*delt_corr
 797        else if (wcorr.gt.0.0d0) then
 798         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
 799          'between contact pairs of peptide groups'
 800        endif
 801        write (iout,'(a,f8.3)') &
 802         'Scaling factor of 1,4 SC-p interactions:',scal14
 803        write (iout,'(a,f8.3)') &
 804         'General scaling factor of SC-p interactions:',scalscp
 805       endif
 806       r0_corr=cutoff_corr-delt_corr
 807       do i=1,ntyp
 808         aad(i,1)=scalscp*aad(i,1)
 809         aad(i,2)=scalscp*aad(i,2)
 810         bad(i,1)=scalscp*bad(i,1)
 811         bad(i,2)=scalscp*bad(i,2)
 812       enddo
 813       call rescale_weights(t_bath)
 814       if(me.eq.king.or..not.out1file) &
 815        write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 816         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 817         wturn4,wturn6
 818    22 format (/'Energy-term weights (scaled):'// &
 819        'WSCC=   ',f10.6,' (SC-SC)'/ &
 820        'WSCP=   ',f10.6,' (SC-p)'/ &
 821        'WELEC=  ',f10.6,' (p-p electr)'/ &
 822        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 823        'WBOND=  ',f10.6,' (stretching)'/ &
 824        'WANG=   ',f10.6,' (bending)'/ &
 825        'WSCLOC= ',f10.6,' (SC local)'/ &
 826        'WTOR=   ',f10.6,' (torsional)'/ &
 827        'WTORD=  ',f10.6,' (double torsional)'/ &
 828        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 829        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 830        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 831        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 832        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 833        'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ &
 834        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 835        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 836        'WTURN6= ',f10.6,' (turns, 6th order)')
 837       if(me.eq.king.or..not.out1file) &
 838        write (iout,*) "Reference temperature for weights calculation:",&
 839         temp0
 840       call reada(weightcard,"D0CM",d0cm,3.78d0)
 841       call reada(weightcard,"AKCM",akcm,15.1d0)
 842       call reada(weightcard,"AKTH",akth,11.0d0)
 843       call reada(weightcard,"AKCT",akct,12.0d0)
 844       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 845       call reada(weightcard,"V2SS",v2ss,7.61d0)
 846       call reada(weightcard,"V3SS",v3ss,13.7d0)
 847       call reada(weightcard,"EBR",ebr,-5.50D0)
 848       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 849
 850       call reada(weightcard,"HT",Ht,0.0D0)
 851       if (dyn_ss) then
 852         ss_depth=ebr/wsc-0.25*eps(1,1)
 853         Ht=Ht/wsc-0.25*eps(1,1)
 854         akcm=akcm*wstrain/wsc
 855         akth=akth*wstrain/wsc
 856         akct=akct*wstrain/wsc
 857         v1ss=v1ss*wstrain/wsc
 858         v2ss=v2ss*wstrain/wsc
 859         v3ss=v3ss*wstrain/wsc
 860       else
 861         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 862       endif
 863
 864       if(me.eq.king.or..not.out1file) then
 865        write (iout,*) "Parameters of the SS-bond potential:"
 866        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,&
 867        " AKCT",akct
 868        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 869        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 870        write (iout,*)" HT",Ht
 871        print *,'indpdb=',indpdb,' pdbref=',pdbref
 872       endif
 873       if (indpdb.gt.0 .or. pdbref) then
 874         read(inp,'(a)') pdbfile
 875         if(me.eq.king.or..not.out1file) &
 876          write (iout,'(2a)') 'PDB data will be read from file ',&
 877          pdbfile(:ilen(pdbfile))
 878         open(ipdbin,file=pdbfile,status='old',err=33)
 879         goto 34
 880   33    write (iout,'(a)') 'Error opening PDB file.'
 881         stop
 882   34    continue
 883 !        print *,'Begin reading pdb data'
 884         call readpdb
 885 !        print *,'Finished reading pdb data'
 886         if(me.eq.king.or..not.out1file) &
 887          write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
 888          ' nstart_sup=',nstart_sup !,"ergwergewrgae"
 889 !el        if(.not.allocated(itype_pdb))
 890         allocate(itype_pdb(nres))
 891         do i=1,nres
 892           itype_pdb(i)=itype(i)
 893         enddo
 894         close (ipdbin)
 895         nnt=nstart_sup
 896         nct=nstart_sup+nsup-1
 897 !el        if(.not.allocated(icont_ref))
 898         allocate(icont_ref(2,12*nres)) ! maxcont=12*maxres
 899         call contact(.false.,ncont_ref,icont_ref,co)
 900
 901         if (sideadd) then
 902          if(me.eq.king.or..not.out1file) &
 903           write(iout,*)'Adding sidechains'
 904          maxsi=1000
 905          do i=2,nres-1
 906           iti=itype(i)
 907           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 908             nsi=0
 909             fail=.true.
 910             do while (fail.and.nsi.le.maxsi)
 911               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 912               nsi=nsi+1
 913             enddo
 914             if(fail) write(iout,*)'Adding sidechain failed for res ',&
 915                     i,' after ',nsi,' trials'
 916           endif
 917          enddo
 918         endif
 919       endif
 920
 921       if (indpdb.eq.0) then
 922 ! Read sequence if not taken from the pdb file.
 923         read (inp,*) nres
 924 !        print *,'nres=',nres
 925         allocate(sequence(nres))
 926         if (iscode.gt.0) then
 927           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 928         else
 929           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 930         endif
 931 ! Convert sequence to numeric code
 932         do i=1,nres
 933           itype(i)=rescode(i,sequence(i),iscode)
 934         enddo
 935 ! Assign initial virtual bond lengths
 936 !elwrite(iout,*) "test_alloc"
 937         if(.not.allocated(vbld)) allocate(vbld(2*nres))
 938 !elwrite(iout,*) "test_alloc"
 939         if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
 940 !elwrite(iout,*) "test_alloc"
 941         do i=2,nres
 942           vbld(i)=vbl
 943           vbld_inv(i)=vblinv
 944         enddo
 945         do i=2,nres-1
 946           vbld(i+nres)=dsc(iabs(itype(i)))
 947           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 948 !          write (iout,*) "i",i," itype",itype(i),
 949 !     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 950         enddo
 951       endif
 952 !      print *,nres
 953 !      print '(20i4)',(itype(i),i=1,nres)
 954 !----------------------------
 955 !el reallocate tables
 956 !      do i=1,maxres2
 957 !        do j=1,3
 958 !          c_alloc(j,i)=c(j,i)
 959 !          dc_alloc(j,i)=dc(j,i)
 960 !        enddo
 961 !      enddo
 962 !      do i=1,maxres
 963 !elwrite(iout,*) "itype",i,itype(i)
 964 !        itype_alloc(i)=itype(i)
 965 !      enddo
 966
 967 !      deallocate(c)
 968 !      deallocate(dc)
 969 !      deallocate(itype)
 970 !      allocate(c(3,2*nres+4))
 971 !      allocate(dc(3,0:2*nres+2))
 972 !      allocate(itype(nres+2))
 973       allocate(itel(nres+2))
 974       itel(:)=0
 975
 976 !      do i=1,2*nres+2
 977 !        do j=1,3
 978 !          c(j,i)=c_alloc(j,i)
 979 !          dc(j,i)=dc_alloc(j,i)
 980 !        enddo
 981 !      enddo
 982 !      do i=1,nres+2
 983 !        itype(i)=itype_alloc(i)
 984 !        itel(i)=0
 985 !      enddo
 986 !--------------------------
 987       do i=1,nres
 988 #ifdef PROCOR
 989         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 990 #else
 991         if (itype(i).eq.ntyp1) then
 992 #endif
 993           itel(i)=0
 994 #ifdef PROCOR
 995         else if (iabs(itype(i+1)).ne.20) then
 996 #else
 997         else if (iabs(itype(i)).ne.20) then
 998 #endif
 999           itel(i)=1
1000         else
1001           itel(i)=2
1002         endif
1003       enddo
1004       if(me.eq.king.or..not.out1file)then
1005        write (iout,*) "ITEL"
1006        do i=1,nres-1
1007          write (iout,*) i,itype(i),itel(i)
1008        enddo
1009        print *,'Call Read_Bridge.'
1010       endif
1011       call read_bridge
1012 !--------------------------------
1013 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1014       call alloc_geo_arrays
1015       call alloc_ener_arrays
1016 !--------------------------------
1017 ! 8/13/98 Set limits to generating the dihedral angles
1018       do i=1,nres
1019         phibound(1,i)=-pi
1020         phibound(2,i)=pi
1021       enddo
1022       read (inp,*) ndih_constr
1023       if (ndih_constr.gt.0) then
1024         allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1025         allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1026         read (inp,*) ftors
1027         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1028         if(me.eq.king.or..not.out1file)then
1029          write (iout,*) &
1030          'There are',ndih_constr,' constraints on phi angles.'
1031          do i=1,ndih_constr
1032           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1033          enddo
1034         endif
1035         do i=1,ndih_constr
1036           phi0(i)=deg2rad*phi0(i)
1037           drange(i)=deg2rad*drange(i)
1038         enddo
1039         if(me.eq.king.or..not.out1file) &
1040          write (iout,*) 'FTORS',ftors
1041         do i=1,ndih_constr
1042           ii = idih_constr(i)
1043           phibound(1,ii) = phi0(i)-drange(i)
1044           phibound(2,ii) = phi0(i)+drange(i)
1045         enddo
1046       endif
1047       nnt=1
1048 #ifdef MPI
1049       if (me.eq.king) then
1050 #endif
1051        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1052        do i=1,nres
1053          write (iout,'(a3,i5,2f10.1)') &
1054          restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1055        enddo
1056 #ifdef MP
1057       endif
1058 #endif
1059       nct=nres
1060 !d      print *,'NNT=',NNT,' NCT=',NCT
1061       if (itype(1).eq.ntyp1) nnt=2
1062       if (itype(nres).eq.ntyp1) nct=nct-1
1063       if (pdbref) then
1064         if(me.eq.king.or..not.out1file) &
1065          write (iout,'(a,i3)') 'nsup=',nsup
1066         nstart_seq=nnt
1067         if (nsup.le.(nct-nnt+1)) then
1068           do i=0,nct-nnt+1-nsup
1069             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1070               nstart_seq=nnt+i
1071               goto 111
1072             endif
1073           enddo
1074           write (iout,'(a)') &
1075                   'Error - sequences to be superposed do not match.'
1076           stop
1077         else
1078           do i=0,nsup-(nct-nnt+1)
1079             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1080             then
1081               nstart_sup=nstart_sup+i
1082               nsup=nct-nnt+1
1083               goto 111
1084             endif
1085           enddo
1086           write (iout,'(a)') &
1087                   'Error - sequences to be superposed do not match.'
1088         endif
1089   111   continue
1090         if (nsup.eq.0) nsup=nct-nnt
1091         if (nstart_sup.eq.0) nstart_sup=nnt
1092         if (nstart_seq.eq.0) nstart_seq=nnt
1093         if(me.eq.king.or..not.out1file) &
1094          write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1095                        ' nstart_seq=',nstart_seq !,"242343453254"
1096       endif
1097 !--- Zscore rms -------
1098       if (nz_start.eq.0) nz_start=nnt
1099       if (nz_end.eq.0 .and. nsup.gt.0) then
1100         nz_end=nnt+nsup-1
1101       else if (nz_end.eq.0) then
1102         nz_end=nct
1103       endif
1104       if(me.eq.king.or..not.out1file)then
1105        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1106        write (iout,*) 'IZ_SC=',iz_sc
1107       endif
1108 !----------------------
1109       call init_int_table
1110       if (refstr) then
1111         if (.not.pdbref) then
1112           call read_angles(inp,*38)
1113           goto 39
1114    38     write (iout,'(a)') 'Error reading reference structure.'
1115 #ifdef MPI
1116           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1117           stop 'Error reading reference structure'
1118 #endif
1119    39     call chainbuild
1120           call setup_var
1121 !zscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1122           nstart_sup=nnt
1123           nstart_seq=nnt
1124           nsup=nct-nnt+1
1125           kkk=1
1126           do i=1,2*nres
1127             do j=1,3
1128               cref(j,i,kkk)=c(j,i)
1129             enddo
1130           enddo
1131           call contact(.true.,ncont_ref,icont_ref,co)
1132         endif
1133 !        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1134         call flush(iout)
1135         if (constr_dist.gt.0) call read_dist_constr
1136         write (iout,*) "After read_dist_constr nhpb",nhpb
1137         call hpb_partition
1138         if(me.eq.king.or..not.out1file) &
1139          write (iout,*) 'Contact order:',co
1140         if (pdbref) then
1141         if(me.eq.king.or..not.out1file) &
1142          write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1143         do i=1,ncont_ref
1144           do j=1,2
1145             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1146           enddo
1147           if(me.eq.king.or..not.out1file) &
1148            write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',&
1149            icont_ref(1,i),' ',&
1150            restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1151         enddo
1152         endif
1153       endif
1154       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1155           .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1156           modecalc.ne.10) then
1157 ! If input structure hasn't been supplied from the PDB file read or generate
1158 ! initial geometry.
1159         if (iranconf.eq.0 .and. .not. extconf) then
1160           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1161            write (iout,'(a)') 'Initial geometry will be read in.'
1162           if (read_cart) then
1163             read(inp,'(8f10.5)',end=36,err=36) &
1164              ((c(l,k),l=1,3),k=1,nres),&
1165              ((c(l,k+nres),l=1,3),k=nnt,nct)
1166             write (iout,*) "Exit READ_CART"
1167             write (iout,'(8f10.5)') &
1168              ((c(l,k),l=1,3),k=1,nres),&
1169              ((c(l,k+nres),l=1,3),k=nnt,nct)
1170             call int_from_cart1(.true.)
1171             write (iout,*) "Finish INT_TO_CART"
1172             do i=1,nres-1
1173               do j=1,3
1174                 dc(j,i)=c(j,i+1)-c(j,i)
1175                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1176               enddo
1177             enddo
1178             do i=nnt,nct
1179               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1180                 do j=1,3
1181                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1182                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1183                 enddo
1184               endif
1185             enddo
1186             return
1187           else
1188             call read_angles(inp,*36)
1189           endif
1190           goto 37
1191    36     write (iout,'(a)') 'Error reading angle file.'
1192 #ifdef MPI
1193           call mpi_finalize( MPI_COMM_WORLD,IERR )
1194 #endif
1195           stop 'Error reading angle file.'
1196    37     continue
1197         else if (extconf) then
1198          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1199           write (iout,'(a)') 'Extended chain initial geometry.'
1200          do i=3,nres
1201           theta(i)=90d0*deg2rad
1202          enddo
1203          do i=4,nres
1204           phi(i)=180d0*deg2rad
1205          enddo
1206          do i=2,nres-1
1207           alph(i)=110d0*deg2rad
1208          enddo
1209 !elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
1210          do i=2,nres-1
1211           omeg(i)=-120d0*deg2rad
1212           if (itype(i).le.0) omeg(i)=-omeg(i)
1213          enddo
1214         else
1215           if(me.eq.king.or..not.out1file) &
1216            write (iout,'(a)') 'Random-generated initial geometry.'
1217
1218
1219 #ifdef MPI
1220           if (me.eq.king  .or. fg_rank.eq.0 .and. &
1221                  ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1222 #endif
1223             do itrial=1,100
1224               itmp=1
1225               call gen_rand_conf(itmp,*30)
1226               goto 40
1227    30         write (iout,*) 'Failed to generate random conformation',&
1228                 ', itrial=',itrial
1229               write (*,*) 'Processor:',me,&
1230                 ' Failed to generate random conformation',&
1231                 ' itrial=',itrial
1232               call intout
1233
1234 #ifdef AIX
1235               call flush_(iout)
1236 #else
1237               call flush(iout)
1238 #endif
1239             enddo
1240             write (iout,'(a,i3,a)') 'Processor:',me,&
1241               ' error in generating random conformation.'
1242             write (*,'(a,i3,a)') 'Processor:',me,&
1243               ' error in generating random conformation.'
1244             call flush(iout)
1245 #ifdef MPI
1246             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1247    40       continue
1248           endif
1249 #else
1250           do itrial=1,100
1251             itmp=1
1252             call gen_rand_conf(itmp,*335)
1253             goto 40
1254   335       write (iout,*) 'Failed to generate random conformation',&
1255               ', itrial=',itrial
1256             write (*,*) 'Failed to generate random conformation',&
1257               ', itrial=',itrial
1258           enddo
1259           write (iout,'(a,i3,a)') 'Processor:',me,&
1260             ' error in generating random conformation.'
1261           write (*,'(a,i3,a)') 'Processor:',me,&
1262             ' error in generating random conformation.'
1263           stop
1264    40     continue
1265 #endif
1266         endif
1267       elseif (modecalc.eq.4) then
1268         read (inp,'(a)') intinname
1269         open (intin,file=intinname,status='old',err=333)
1270         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1271         write (iout,'(a)') 'intinname',intinname
1272         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1273         goto 334
1274   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1275 #ifdef MPI
1276         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1277 #endif
1278         stop 'Error opening angle file.'
1279   334   continue
1280
1281       endif
1282 ! Generate distance constraints, if the PDB structure is to be regularized.
1283       if (nthread.gt.0) then
1284         call read_threadbase
1285       endif
1286       call setup_var
1287 !elwrite (iout,*)"alph(i)*deg2rad",(alph(i), i=1,nres)
1288       if (me.eq.king .or. .not. out1file) &
1289        call intout
1290 !elwrite (iout,*)"alph(i)*rad2deg",(alph(i)*rad2deg, i=1,nres)
1291       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1292         write (iout,'(/a,i3,a)') &
1293         'The chain contains',ns,' disulfide-bridging cysteines.'
1294         write (iout,'(20i4)') (iss(i),i=1,ns)
1295        if (dyn_ss) then
1296           write(iout,*)"Running with dynamic disulfide-bond formation"
1297        else
1298         write (iout,'(/a/)') 'Pre-formed links are:'
1299         do i=1,nss
1300           i1=ihpb(i)-nres
1301           i2=jhpb(i)-nres
1302           it1=itype(i1)
1303           it2=itype(i2)
1304           if (me.eq.king.or..not.out1file) &
1305           write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1306           restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),&
1307           ebr,forcon(i)
1308         enddo
1309         write (iout,'(a)')
1310        endif
1311       endif
1312       if (ns.gt.0.and.dyn_ss) then
1313           do i=nss+1,nhpb
1314             ihpb(i-nss)=ihpb(i)
1315             jhpb(i-nss)=jhpb(i)
1316             forcon(i-nss)=forcon(i)
1317             dhpb(i-nss)=dhpb(i)
1318           enddo
1319           nhpb=nhpb-nss
1320           nss=0
1321           call hpb_partition
1322           do i=1,ns
1323             dyn_ss_mask(iss(i))=.true.
1324           enddo
1325       endif
1326       if (i2ndstr.gt.0) call secstrp2dihc
1327 !      call geom_to_var(nvar,x)
1328 !      call etotal(energia(0))
1329 !      call enerprint(energia(0))
1330 !      call briefout(0,etot)
1331 !      stop
1332 !d    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1333 !d    write (iout,'(a)') 'Variable list:'
1334 !d    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1335 #ifdef MPI
1336       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1337         write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1338         'Processor',myrank,': end reading molecular data.'
1339 #endif
1340       return
1341       end subroutine molread
1342 !-----------------------------------------------------------------------------
1343 !-----------------------------------------------------------------------------
1344       end module io