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,molnum(i))
 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(.false.,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+potEcomp(51),totE+potEcomp(51),&
 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+potEcomp(51),totE+potEcomp(51),&
 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+potEcomp(51),totE+potEcomp(51),&
 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+potEcomp(51),totE+potEcomp(51), &
 565                 kinetic_T,t_bath,gyrate(),&
 566                 distance,potEcomp(23),me
 567           format1="a114"
 568         endif
 569        else if (velnanoconst.ne.0 ) then
 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,2f6.3,f12.3,f10.1,2f8.2, &
 573          f9.3,i5,$)') &
 574                itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
 575                rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
 576                vecsim,vectrue,me
 577           format1="a133"
 578 !C          print *,"CHUJOWO"
 579          else
 580 !C          print *,'A CHUJ',potEcomp(23)
 581           write (line1,'(i10,f15.2,8f12.3,i5,$)')&
 582                 itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51), &
 583                 kinetic_T,t_bath,gyrate(),&
 584                 vecsim,vectrue,me
 585           format1="a114"
 586         endif
 587        else
 588        if (refstr) then
 589          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 590           write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') &
 591                 itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
 592                 rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
 593           format1="a133"
 594         else
 595           write (line1,'(i10,f15.2,7f12.3,i5,$)') &
 596                  itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
 597                  amax,kinetic_T,t_bath,gyrate(),me
 598           format1="a114"
 599         endif
 600         ENDIF
 601         if(usampl.and.totT.gt.eq_time) then
 602            write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,&
 603             (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),&
 604             (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 605            write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair &
 606                    +21*nfrag_back
 607         else
 608            format2="a001"
 609            line2=' '
 610         endif
 611         if (print_compon) then
 612           if(itime.eq.0) then
 613            write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
 614                                                            ",20a12)"
 615            write (istat,format) "#","",&
 616             (ename(print_order(i)),i=1,nprint_ene)
 617           endif
 618           write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
 619                                                            ",20f12.3)"
 620           write (istat,format) line1,line2,&
 621             (potEcomp(print_order(i)),i=1,nprint_ene)
 622         else
 623           write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
 624           write (istat,format) line1,line2
 625         endif
 626 #if defined(AIX)
 627         call flush(istat)
 628 #else
 629         close(istat)
 630 #endif
 631       return
 632       end subroutine  statout
 633 !-----------------------------------------------------------------------------
 634 ! readrtns_CSA.F
 635 !-----------------------------------------------------------------------------
 636       subroutine readrtns
 637
 638       use control_data
 639       use energy_data
 640       use MPI_data
 641       use muca_md, only:read_muca
 642 !      implicit real*8 (a-h,o-z)
 643 !      include 'DIMENSIONS'
 644 #ifdef MPI
 645       include 'mpif.h'
 646 #endif
 647 !      include 'COMMON.SETUP'
 648 !      include 'COMMON.CONTROL'
 649 !      include 'COMMON.SBRIDGE'
 650 !      include 'COMMON.IOUNITS'
 651       logical :: file_exist
 652       integer :: i
 653 ! Read force-field parameters except weights
 654 !      call parmread
 655 ! Read job setup parameters
 656       call read_control
 657 ! Read force-field parameters except weights
 658       call parmread
 659
 660 ! Read control parameters for energy minimzation if required
 661       if (minim) call read_minim
 662 ! Read MCM control parameters if required
 663       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
 664 ! Read MD control parameters if reqjuired
 665       if (modecalc.eq.12) call read_MDpar
 666 ! Read MREMD control parameters if required
 667       if (modecalc.eq.14) then
 668          call read_MDpar
 669          call read_REMDpar
 670       endif
 671 ! Read MUCA control parameters if required
 672       if (lmuca) call read_muca
 673 ! Read CSA control parameters if required (from fort.40 if exists
 674 ! otherwise from general input file)
 675       if (modecalc.eq.8) then
 676        inquire (file="fort.40",exist=file_exist)
 677        if (.not.file_exist) call csaread
 678       endif
 679 !fmc      if (modecalc.eq.10) call mcmfread
 680 ! Read molecule information, molecule geometry, energy-term weights, and
 681 ! restraints if requested
 682       call molread
 683 ! Print restraint information
 684 #ifdef MPI
 685       if (.not. out1file .or. me.eq.king) then
 686 #endif
 687       if (nhpb.gt.nss) &
 688       write (iout,'(a,i5,a)') "The following",nhpb-nss,&
 689        " distance constraints have been imposed"
 690       do i=nss+1,nhpb
 691         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
 692       enddo
 693 #ifdef MPI
 694       endif
 695 #endif
 696 !      print *,"Processor",myrank," leaves READRTNS"
 697 !      write(iout,*) "end readrtns"
 698       return
 699       end subroutine readrtns
 700 !-----------------------------------------------------------------------------
 701       subroutine molread
 702 !
 703 ! Read molecular data.
 704 !
 705 !      use control, only: ilen
 706       use control_data
 707       use geometry_data
 708       use energy_data
 709       use energy
 710       use compare_data
 711       use MD_data, only: t_bath
 712       use MPI_data
 713       use compare, only:seq_comp,contact
 714       use control
 715       implicit none
 716 !      implicit real*8 (a-h,o-z)
 717 !      include 'DIMENSIONS'
 718 #ifdef MPI
 719       include 'mpif.h'
 720       integer :: error_msg,ierror,ierr,ierrcode
 721 #endif
 722 !      include 'COMMON.IOUNITS'
 723 !      include 'COMMON.GEO'
 724 !      include 'COMMON.VAR'
 725 !      include 'COMMON.INTERACT'
 726 !      include 'COMMON.LOCAL'
 727 !      include 'COMMON.NAMES'
 728 !      include 'COMMON.CHAIN'
 729 !      include 'COMMON.FFIELD'
 730 !      include 'COMMON.SBRIDGE'
 731 !      include 'COMMON.HEADER'
 732 !      include 'COMMON.CONTROL'
 733 !      include 'COMMON.DBASE'
 734 !      include 'COMMON.THREAD'
 735 !      include 'COMMON.CONTACTS'
 736 !      include 'COMMON.TORCNSTR'
 737 !      include 'COMMON.TIME1'
 738 !      include 'COMMON.BOUNDS'
 739 !      include 'COMMON.MD'
 740 !      include 'COMMON.SETUP'
 741       character(len=4),dimension(:,:),allocatable :: sequence   !(maxres,maxmolecules)
 742 !      integer :: rescode
 743 !      double precision x(maxvar)
 744       character(len=256) :: pdbfile
 745       character(len=800) :: weightcard
 746       character(len=80) :: weightcard_t!,ucase
 747 !      integer,dimension(:),allocatable :: itype_pdb    !(maxres)
 748 !      common /pizda/ itype_pdb
 749       logical :: fail   !seq_comp,
 750       real(kind=8) :: energia(0:n_ene)
 751 !      integer ilen
 752 !el      external ilen
 753 !el local varables
 754       integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2,mnum
 755
 756       real(kind=8),dimension(3,maxres2+2) :: c_alloc
 757       real(kind=8),dimension(3,0:maxres2) :: dc_alloc
 758       real(kind=8),dimension(:,:), allocatable :: secprob
 759       integer,dimension(maxres) :: itype_alloc
 760
 761       integer :: iti,nsi,maxsi,itrial,itmp
 762       real(kind=8) :: wlong,scalscp,co,ssscale,phihel,phibet,sigmahel,&
 763       sigmabet,sumv,nres_temp
 764       allocate(weights(n_ene))
 765 !-----------------------------
 766       allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
 767       allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
 768       allocate(itype(maxres,5)) !(maxres)
 769       allocate(istype(maxres))
 770 !
 771 ! Zero out tables.
 772 !
 773       c(:,:)=0.0D0
 774       dc(:,:)=0.0D0
 775       itype(:,:)=0
 776 !-----------------------------
 777 !
 778 ! Body
 779 !
 780 ! Read weights of the subsequent energy terms.
 781       call card_concat(weightcard,.true.)
 782       call reada(weightcard,'WLONG',wlong,1.0D0)
 783       call reada(weightcard,'WSC',wsc,wlong)
 784       call reada(weightcard,'WSCP',wscp,wlong)
 785       call reada(weightcard,'WELEC',welec,1.0D0)
 786       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 787       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 788       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 789       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 790       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 791       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 792       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 793       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 794       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 795       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 796       call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
 797       call reada(weightcard,'WELPP',welpp,0.0d0)
 798       call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
 799       call reada(weightcard,'WELPSB',welpsb,0.0D0)
 800       call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
 801       call reada(weightcard,'WELSB',welsb,0.0D0)
 802       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
 803       call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
 804       call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
 805       call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
 806 !      print *,"KUR..",wtor_nucl
 807       call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
 808       call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
 809       call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
 810       call reada(weightcard,'WBOND',wbond,1.0D0)
 811       call reada(weightcard,'WTOR',wtor,1.0D0)
 812       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 813       call reada(weightcard,'WSHIELD',wshield,0.05D0)
 814       call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
 815       call reada(weightcard,'WLT',wliptran,1.0D0)
 816       call reada(weightcard,'WTUBE',wtube,1.0d0)
 817       call reada(weightcard,'WANG',wang,1.0D0)
 818       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 819       call reada(weightcard,'SCAL14',scal14,0.4D0)
 820       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 821       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 822       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 823       call reada(weightcard,'TEMP0',temp0,300.0d0)
 824       call reada(weightcard,'WSCBASE',wscbase,0.0D0)
 825       if (index(weightcard,'SOFT').gt.0) ipot=6
 826       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
 827       call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
 828       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
 829       call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
 830       call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0)
 831       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
 832       call reada(weightcard,'WSCPHO',wscpho,0.0d0)
 833       call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
 834
 835 ! 12/1/95 Added weight for the multi-body term WCORR
 836       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 837       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 838       weights(1)=wsc
 839       weights(2)=wscp
 840       weights(3)=welec
 841       weights(4)=wcorr
 842       weights(5)=wcorr5
 843       weights(6)=wcorr6
 844       weights(7)=wel_loc
 845       weights(8)=wturn3
 846       weights(9)=wturn4
 847       weights(10)=wturn6
 848       weights(11)=wang
 849       weights(12)=wscloc
 850       weights(13)=wtor
 851       weights(14)=wtor_d
 852       weights(15)=wstrain
 853       weights(16)=wvdwpp
 854       weights(17)=wbond
 855       weights(18)=scal14
 856       weights(21)=wsccor
 857           weights(26)=wvdwpp_nucl
 858           weights(27)=welpp
 859           weights(28)=wvdwpsb
 860           weights(29)=welpsb
 861           weights(30)=wvdwsb
 862           weights(31)=welsb
 863           weights(32)=wbond_nucl
 864           weights(33)=wang_nucl
 865           weights(34)=wsbloc
 866           weights(35)=wtor_nucl
 867           weights(36)=wtor_d_nucl
 868           weights(37)=wcorr_nucl
 869           weights(38)=wcorr3_nucl
 870           weights(41)=wcatcat
 871           weights(42)=wcatprot
 872           weights(46)=wscbase
 873           weights(47)=wpepbase
 874           weights(48)=wscpho
 875           weights(49)=wpeppho
 876           weights(50)=wcatnucl
 877           weights(56)=wcat_tran
 878       if(me.eq.king.or..not.out1file) &
 879        write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 880         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 881         wturn4,wturn6
 882    10 format (/'Energy-term weights (unscaled):'// &
 883        'WSCC=   ',f10.6,' (SC-SC)'/ &
 884        'WSCP=   ',f10.6,' (SC-p)'/ &
 885        'WELEC=  ',f10.6,' (p-p electr)'/ &
 886        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 887        'WBOND=  ',f10.6,' (stretching)'/ &
 888        'WANG=   ',f10.6,' (bending)'/ &
 889        'WSCLOC= ',f10.6,' (SC local)'/ &
 890        'WTOR=   ',f10.6,' (torsional)'/ &
 891        'WTORD=  ',f10.6,' (double torsional)'/ &
 892        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 893        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 894        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 895        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 896        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 897        'WSCCOR= ',f10.6,' (back-scloc correlation)'/ &
 898        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 899        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 900        'WTURN6= ',f10.6,' (turns, 6th order)')
 901       if(me.eq.king.or..not.out1file)then
 902        if (wcorr4.gt.0.0d0) then
 903         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
 904          'between contact pairs of peptide groups'
 905         write (iout,'(2(a,f5.3/))') &
 906         'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
 907         'Range of quenching the correlation terms:',2*delt_corr
 908        else if (wcorr.gt.0.0d0) then
 909         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
 910          'between contact pairs of peptide groups'
 911        endif
 912        write (iout,'(a,f8.3)') &
 913         'Scaling factor of 1,4 SC-p interactions:',scal14
 914        write (iout,'(a,f8.3)') &
 915         'General scaling factor of SC-p interactions:',scalscp
 916       endif
 917       r0_corr=cutoff_corr-delt_corr
 918       do i=1,ntyp
 919         aad(i,1)=scalscp*aad(i,1)
 920         aad(i,2)=scalscp*aad(i,2)
 921         bad(i,1)=scalscp*bad(i,1)
 922         bad(i,2)=scalscp*bad(i,2)
 923       enddo
 924       call rescale_weights(t_bath)
 925       if(me.eq.king.or..not.out1file) &
 926        write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
 927         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
 928         wturn4,wturn6
 929    22 format (/'Energy-term weights (scaled):'// &
 930        'WSCC=   ',f10.6,' (SC-SC)'/ &
 931        'WSCP=   ',f10.6,' (SC-p)'/ &
 932        'WELEC=  ',f10.6,' (p-p electr)'/ &
 933        'WVDWPP= ',f10.6,' (p-p VDW)'/ &
 934        'WBOND=  ',f10.6,' (stretching)'/ &
 935        'WANG=   ',f10.6,' (bending)'/ &
 936        'WSCLOC= ',f10.6,' (SC local)'/ &
 937        'WTOR=   ',f10.6,' (torsional)'/ &
 938        'WTORD=  ',f10.6,' (double torsional)'/ &
 939        'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
 940        'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
 941        'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
 942        'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
 943        'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
 944        'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ &
 945        'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
 946        'WTURN4= ',f10.6,' (turns, 4th order)'/ &
 947        'WTURN6= ',f10.6,' (turns, 6th order)')
 948       if(me.eq.king.or..not.out1file) &
 949        write (iout,*) "Reference temperature for weights calculation:",&
 950         temp0
 951       call reada(weightcard,"D0CM",d0cm,3.78d0)
 952       call reada(weightcard,"AKCM",akcm,15.1d0)
 953       call reada(weightcard,"AKTH",akth,11.0d0)
 954       call reada(weightcard,"AKCT",akct,12.0d0)
 955       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 956       call reada(weightcard,"V2SS",v2ss,7.61d0)
 957       call reada(weightcard,"V3SS",v3ss,13.7d0)
 958       call reada(weightcard,"EBR",ebr,-5.50D0)
 959       call reada(weightcard,"ATRISS",atriss,0.301D0)
 960       call reada(weightcard,"BTRISS",btriss,0.021D0)
 961       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 962       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 963       call reada(weightcard,"SSSCALE",ssscale,1.0D0)
 964       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 965
 966       call reada(weightcard,"HT",Ht,0.0D0)
 967       if (dyn_ss) then
 968        ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
 969         Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
 970         akcm=akcm/wsc*ssscale
 971         akth=akth/wsc*ssscale
 972         akct=akct/wsc*ssscale
 973         v1ss=v1ss/wsc*ssscale
 974         v2ss=v2ss/wsc*ssscale
 975         v3ss=v3ss/wsc*ssscale
 976       else
 977         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 978       endif
 979
 980       if(me.eq.king.or..not.out1file) then
 981        write (iout,*) "Parameters of the SS-bond potential:"
 982        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,&
 983        " AKCT",akct
 984        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 985        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 986        write (iout,*)" HT",Ht
 987        print *,'indpdb=',indpdb,' pdbref=',pdbref
 988       endif
 989       if (indpdb.gt.0 .or. pdbref) then
 990         read(inp,'(a)') pdbfile
 991         if(me.eq.king.or..not.out1file) &
 992          write (iout,'(2a)') 'PDB data will be read from file ',&
 993          pdbfile(:ilen(pdbfile))
 994         open(ipdbin,file=pdbfile,status='old',err=33)
 995         goto 34
 996   33    write (iout,'(a)') 'Error opening PDB file.'
 997         stop
 998   34    continue
 999 !        print *,'Begin reading pdb data'
1000         call readpdb
1001         if (.not.allocated(crefjlee)) allocate (crefjlee(3,2*nres+2))
1002         do i=1,2*nres
1003           do j=1,3
1004             crefjlee(j,i)=c(j,i)
1005           enddo
1006         enddo
1007 #ifdef DEBUG
1008         do i=1,nres
1009           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1010      &      (crefjlee(j,i+nres),j=1,3)
1011         enddo
1012 #endif
1013
1014 !        call int_from_cart1(.true.)
1015
1016 !        print *,'Finished reading pdb data'
1017         if(me.eq.king.or..not.out1file) &
1018          write (iout,'(a,i7,a,i7)')'nsup=',nsup,&
1019          ' nstart_sup=',nstart_sup !,"ergwergewrgae"
1020 !el        if(.not.allocated(itype_pdb))
1021         allocate(itype_pdb(nres))
1022         do i=1,nres
1023           itype_pdb(i)=itype(i,1)
1024         enddo
1025         close (ipdbin)
1026         nnt=nstart_sup
1027         nct=nstart_sup+nsup-1
1028 !el        if(.not.allocated(icont_ref))
1029         allocate(icont_ref(2,(nres/2)*nres)) ! maxcont=12*maxres
1030         call contact(.false.,ncont_ref,icont_ref,co)
1031
1032         if (sideadd) then
1033          if(me.eq.king.or..not.out1file) &
1034           write(iout,*)'Adding sidechains'
1035          maxsi=1000
1036          do i=2,nres-1
1037           iti=itype(i,1)
1038           if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then
1039             nsi=0
1040             fail=.true.
1041             do while (fail.and.nsi.le.maxsi)
1042               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail,molnum(i))
1043               nsi=nsi+1
1044             enddo
1045             if(fail) write(iout,*)'Adding sidechain failed for res ',&
1046                     i,' after ',nsi,' trials'
1047           endif
1048          enddo
1049         endif
1050       endif
1051       print *,"CZY TU DOCHODZE"
1052       if (indpdb.eq.0) then
1053       nres_molec(:)=0
1054         allocate(sequence(maxres,5))
1055 !      itype(:,:)=0
1056       itmp=0
1057       if (protein) then
1058 ! Read sequence if not taken from the pdb file.
1059         molec=1
1060         read (inp,*) nres_molec(molec)
1061         print *,'nres=',nres
1062         if (iscode.gt.0) then
1063           read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec))
1064         else
1065           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1066         endif
1067 !        read(inp,*) weightcard_t
1068 !        print *,"po seq" weightcard_t
1069 ! Convert sequence to numeric code
1070
1071         do i=1,nres_molec(molec)
1072           itmp=itmp+1
1073           itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
1074           print *,itype(i,1)
1075
1076         enddo
1077        endif
1078 !        read(inp,*) weightcard_t
1079 !        print *,"po seq", weightcard_t
1080
1081        if (nucleic) then
1082 ! Read sequence if not taken from the pdb file.
1083         molec=2
1084         read (inp,*) nres_molec(molec)
1085         print *,'nres=',nres_molec(molec)
1086 !        allocate(sequence(maxres,5))
1087 !        if (iscode.gt.0) then
1088           read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
1089           print *,"KUR**"
1090           print *,(sequence(i,molec),i=1,nres_molec(molec))
1091 ! Convert sequence to numeric code
1092
1093         do i=1,nres_molec(molec)
1094           itmp=itmp+1
1095           istype(itmp)=sugarcode(sequence(i,molec)(1:1),i)
1096           sequence(i,molec)=sequence(i,molec)(1:2)
1097           itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
1098           write(iout,*) "NUCLE=", itype(itmp,molec)
1099         enddo
1100        endif
1101
1102        if (ions) then
1103 ! Read sequence if not taken from the pdb file.
1104         molec=5
1105         read (inp,*) nres_molec(molec)
1106 !        print *,'nres=',nres
1107           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1108 ! Convert sequence to numeric code
1109         print *,nres_molec(molec)
1110         do i=1,nres_molec(molec)
1111           itmp=itmp+1
1112           print *,itmp,"itmp"
1113           itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
1114         enddo
1115        endif
1116        nres=0
1117        do i=1,5
1118         nres=nres+nres_molec(i)
1119         print *,"nres_molec",nres,nres_molec(i)
1120        enddo
1121
1122 ! Assign initial virtual bond lengths
1123         if(.not.allocated(molnum)) then
1124          allocate(molnum(nres+1))
1125          itmp=0
1126         do i=1,5
1127                do j=1,nres_molec(i)
1128                itmp=itmp+1
1129               molnum(itmp)=i
1130                enddo
1131          enddo
1132 !        print *,nres_molec(i)
1133         endif
1134         print *,nres,"nres"
1135         if(.not.allocated(vbld)) then
1136            print *, "I DO ENTER"
1137            allocate(vbld(2*nres))
1138         endif
1139         if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
1140         do i=2,nres
1141            mnum=molnum(i)
1142           if (molnum(i).eq.1) then
1143           vbld(i)=vbl
1144           vbld_inv(i)=vblinv
1145
1146           else
1147           write(iout,*) "typy",itype(i,mnum),ntyp1_molec(mnum),i
1148           vbld(i)=6.0
1149           vbld_inv(i)=1.0/6.0
1150           if ((itype(i,mnum).eq.ntyp1_molec(mnum)).or.&
1151           (itype(i-1,mnum).eq.ntyp1_molec(mnum))) then
1152           vbld(i)=1.9
1153           vbld_inv(i)=1.0/1.9
1154           endif
1155           endif
1156         enddo
1157         do i=2,nres-1
1158            mnum=molnum(i)
1159            if (molnum(i).eq.1) then
1160 !          print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i
1161           vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
1162           vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
1163            else
1164           write(iout,*) "typy2",itype(i,mnum),ntyp1_molec(mnum),i
1165           if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
1166           vbld_inv(i+nres)=1.0
1167           vbld(i+nres)=0.0
1168           else
1169           vbld(i+nres)=vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
1170           vbld_inv(i+nres)=1.0/vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
1171           endif
1172            endif
1173 !          write (iout,*) "i",i," itype",itype(i,1),
1174 !     &      " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
1175         enddo
1176       endif
1177 !      print *,nres
1178 !      print '(20i4)',(itype(i,1),i=1,nres)
1179 !----------------------------
1180 !el reallocate tables
1181 !      do i=1,maxres2
1182 !        do j=1,3
1183 !          c_alloc(j,i)=c(j,i)
1184 !          dc_alloc(j,i)=dc(j,i)
1185 !        enddo
1186 !      enddo
1187 !      do i=1,maxres
1188 !elwrite(iout,*) "itype",i,itype(i,1)
1189 !        itype_alloc(i)=itype(i,1)
1190 !      enddo
1191
1192 !      deallocate(c)
1193 !      deallocate(dc)
1194 !      deallocate(itype)
1195 !      allocate(c(3,2*nres+4))
1196 !      allocate(dc(3,0:2*nres+2))
1197 !      allocate(itype(nres+2))
1198       allocate(itel(nres+2))
1199       itel(:)=0
1200
1201 !      do i=1,2*nres+2
1202 !        do j=1,3
1203 !          c(j,i)=c_alloc(j,i)
1204 !          dc(j,i)=dc_alloc(j,i)
1205 !        enddo
1206 !      enddo
1207 !      do i=1,nres+2
1208 !        itype(i,1)=itype_alloc(i)
1209 !        itel(i)=0
1210 !      enddo
1211 !--------------------------
1212       do i=1,nres
1213 #ifdef PROCOR
1214         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then
1215 #else
1216         if (itype(i,1).eq.ntyp1) then
1217 #endif
1218           itel(i)=0
1219 #ifdef PROCOR
1220         else if (iabs(itype(i+1,1)).ne.20) then
1221 #else
1222         else if (iabs(itype(i,1)).ne.20) then
1223 #endif
1224           itel(i)=1
1225         else
1226           itel(i)=2
1227         endif
1228       enddo
1229       if(me.eq.king.or..not.out1file)then
1230        write (iout,*) "ITEL"
1231        print *,nres,"nres"
1232        do i=1,nres-1
1233          write (iout,*) i,itype(i,1),itel(i)
1234        enddo
1235        print *,'Call Read_Bridge.'
1236       endif
1237       call read_bridge
1238 !--------------------------------
1239 !       print *,"tu dochodze"
1240 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1241       call alloc_geo_arrays
1242       print *,"after  alloc_geo_arrays"
1243       call alloc_ener_arrays
1244       print *,"after  alloc_ener_arrays"
1245 !--------------------------------
1246 ! 8/13/98 Set limits to generating the dihedral angles
1247       do i=1,nres
1248         phibound(1,i)=-pi
1249         phibound(2,i)=pi
1250       enddo
1251       read (inp,*) ndih_constr
1252       if (ndih_constr.gt.0) then
1253         raw_psipred=.false.
1254         allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1255         allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1256         allocate(ftors(ndih_constr)) !(maxdih_constr)
1257
1258 !        read (inp,*) ftors
1259         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), &
1260         i=1,ndih_constr)
1261         if(me.eq.king.or..not.out1file)then
1262          write (iout,*) &
1263          'There are',ndih_constr,' constraints on phi angles.'
1264          do i=1,ndih_constr
1265           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), &
1266           ftors(i)
1267          enddo
1268         endif
1269         do i=1,ndih_constr
1270           phi0(i)=deg2rad*phi0(i)
1271           drange(i)=deg2rad*drange(i)
1272         enddo
1273 !        if(me.eq.king.or..not.out1file) &
1274 !         write (iout,*) 'FTORS',ftors
1275         do i=1,ndih_constr
1276           ii = idih_constr(i)
1277           phibound(1,ii) = phi0(i)-drange(i)
1278           phibound(2,ii) = phi0(i)+drange(i)
1279         enddo
1280       else if (ndih_constr.lt.0) then
1281         raw_psipred=.true.
1282         allocate(idih_constr(nres))
1283         allocate(secprob(3,nres))
1284         allocate(vpsipred(3,nres))
1285         allocate(sdihed(2,nres))
1286         call card_concat(weightcard,.true.)
1287         call reada(weightcard,"PHIHEL",phihel,50.0D0)
1288         call reada(weightcard,"PHIBET",phibet,180.0D0)
1289         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
1290         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
1291         call reada(weightcard,"WDIHC",wdihc,0.591D0)
1292         write (iout,*) "Weight of dihedral angle restraints",wdihc
1293         read(inp,'(9x,3f7.3)') &
1294           (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
1295         write (iout,*) "The secprob array"
1296         do i=nnt,nct
1297           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
1298         enddo
1299         ndih_constr=0
1300         do i=nnt+3,nct
1301           if (itype(i-3,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1 &
1302           .and. itype(i-1,1).ne.ntyp1 .and. itype(i,1).ne.ntyp1) then
1303             ndih_constr=ndih_constr+1
1304             idih_constr(ndih_constr)=i
1305             sumv=0.0d0
1306             do j=1,3
1307               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
1308               sumv=sumv+vpsipred(j,ndih_constr)
1309             enddo
1310             do j=1,3
1311               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
1312             enddo
1313             phibound(1,ndih_constr)=phihel*deg2rad
1314             phibound(2,ndih_constr)=phibet*deg2rad
1315             sdihed(1,ndih_constr)=sigmahel*deg2rad
1316             sdihed(2,ndih_constr)=sigmabet*deg2rad
1317           endif
1318         enddo
1319
1320       endif
1321       if (with_theta_constr) then
1322 !C with_theta_constr is keyword allowing for occurance of theta constrains
1323       read (inp,*) ntheta_constr
1324 !C ntheta_constr is the number of theta constrains
1325       if (ntheta_constr.gt.0) then
1326 !C        read (inp,*) ftors
1327         allocate(itheta_constr(ntheta_constr))
1328         allocate(theta_constr0(ntheta_constr))
1329         allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr))
1330         read (inp,*) (itheta_constr(i),theta_constr0(i), &
1331        theta_drange(i),for_thet_constr(i), &
1332        i=1,ntheta_constr)
1333 !C the above code reads from 1 to ntheta_constr
1334 !C itheta_constr(i) residue i for which is theta_constr
1335 !C theta_constr0 the global minimum value
1336 !C theta_drange is range for which there is no energy penalty
1337 !C for_thet_constr is the force constant for quartic energy penalty
1338 !C E=k*x**4
1339         if(me.eq.king.or..not.out1file)then
1340          write (iout,*) &
1341         'There are',ntheta_constr,' constraints on phi angles.'
1342          do i=1,ntheta_constr
1343           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), &
1344          theta_drange(i), &
1345          for_thet_constr(i)
1346          enddo
1347         endif
1348         do i=1,ntheta_constr
1349           theta_constr0(i)=deg2rad*theta_constr0(i)
1350           theta_drange(i)=deg2rad*theta_drange(i)
1351         enddo
1352 !C        if(me.eq.king.or..not.out1file)
1353 !C     &   write (iout,*) 'FTORS',ftors
1354 !C        do i=1,ntheta_constr
1355 !C          ii = itheta_constr(i)
1356 !C          thetabound(1,ii) = phi0(i)-drange(i)
1357 !C          thetabound(2,ii) = phi0(i)+drange(i)
1358 !C        enddo
1359       endif ! ntheta_constr.gt.0
1360       endif! with_theta_constr
1361
1362       nnt=1
1363 #ifdef MPI
1364       if (me.eq.king) then
1365 #endif
1366        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1367        do i=1,nres
1368          write (iout,'(a3,i5,2f10.1)') &
1369          restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1370        enddo
1371 #ifdef MP
1372       endif
1373 #endif
1374       nct=nres
1375       allocate(ireschain(nres))
1376       ireschain=0
1377       write(iout,*),"before seq2chains",ireschain
1378       call seq2chains
1379       write(iout,*) "after seq2chains",nchain
1380       allocate ( chain_border1(2,nchain))
1381       chain_border1(1,1)=1
1382       chain_border1(2,1)=chain_border(2,1)+1
1383       do i=2,nchain-1
1384         chain_border1(1,i)=chain_border(1,i)-1
1385         chain_border1(2,i)=chain_border(2,i)+1
1386       enddo
1387       if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
1388       chain_border1(2,nchain)=nres
1389       write(iout,*) "nres",nres," nchain",nchain
1390       do i=1,nchain
1391         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),&
1392          chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
1393       enddo
1394       allocate(tabpermchain(50,5040))
1395       call chain_symmetry(npermchain,tabpermchain)
1396       print *,'NNT=',NNT,' NCT=',NCT
1397       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1398       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1399       if (pdbref) then
1400         if(me.eq.king.or..not.out1file) &
1401          write (iout,'(a,i7)') 'nsup=',nsup
1402         nstart_seq=nnt
1403         if (nsup.le.(nct-nnt+1)) then
1404           do i=0,nct-nnt+1-nsup
1405             if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
1406               nstart_seq=nnt+i
1407               goto 111
1408             endif
1409           enddo
1410           write (iout,'(a)') &
1411                   'Error - sequences to be superposed do not match.'
1412           stop
1413         else
1414           do i=0,nsup-(nct-nnt+1)
1415             if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1416             then
1417               nstart_sup=nstart_sup+i
1418               nsup=nct-nnt+1
1419               goto 111
1420             endif
1421           enddo
1422           write (iout,'(a)') &
1423                   'Error - sequences to be superposed do not match.'
1424         endif
1425   111   continue
1426         if (nsup.eq.0) nsup=nct-nnt+1
1427         if (nstart_sup.eq.0) nstart_sup=nnt
1428         if (nstart_seq.eq.0) nstart_seq=nnt
1429         if(me.eq.king.or..not.out1file) &
1430          write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1431                        ' nstart_seq=',nstart_seq !,"242343453254"
1432       endif
1433 !--- Zscore rms -------
1434       if (nz_start.eq.0) nz_start=nnt
1435       if (nz_end.eq.0 .and. nsup.gt.0) then
1436         nz_end=nnt+nsup-1
1437       else if (nz_end.eq.0) then
1438         nz_end=nct
1439       endif
1440       if(me.eq.king.or..not.out1file)then
1441        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1442        write (iout,*) 'IZ_SC=',iz_sc, 'NCT=',nct
1443       endif
1444 !----------------------
1445       call init_int_table
1446       if (refstr) then
1447         if (.not.pdbref) then
1448           call read_angles(inp,*38)
1449           goto 39
1450    38     write (iout,'(a)') 'Error reading reference structure.'
1451 #ifdef MPI
1452           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1453           stop 'Error reading reference structure'
1454 #endif
1455    39     call chainbuild
1456           call setup_var
1457 !zscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1458           nstart_sup=nnt
1459           nstart_seq=nnt
1460           nsup=nct-nnt+1
1461           kkk=1
1462           do i=1,2*nres
1463             do j=1,3
1464               cref(j,i,kkk)=c(j,i)
1465             enddo
1466           enddo
1467           call contact(.false.,ncont_ref,icont_ref,co)
1468         endif
1469 !        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1470 !        call flush(iout)
1471 !EL        if (constr_dist.gt.0) call read_dist_constr
1472 !EL        write (iout,*) "After read_dist_constr nhpb",nhpb
1473 !EL        if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1474 !EL        call hpb_partition
1475         if(me.eq.king.or..not.out1file) &
1476          write (iout,*) 'Contact order:',co
1477         if (pdbref) then
1478         if(me.eq.king.or..not.out1file) &
1479          write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1480          if (.false.) then
1481         do i=1,ncont_ref
1482           do j=1,2
1483             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1484           enddo
1485           if(me.eq.king.or..not.out1file) &
1486            write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1487            icont_ref(1,i),' ',&
1488            restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1489         enddo
1490         endif
1491         endif
1492       if (constr_homology.gt.0) then
1493 !        write (iout,*) "Calling read_constr_homology"
1494 !        call flush(iout)
1495         call read_constr_homology
1496         if (indpdb.gt.0 .or. pdbref) then
1497           do i=1,2*nres
1498             do j=1,3
1499               c(j,i)=crefjlee(j,i)
1500               cref(j,i,1)=crefjlee(j,i)
1501             enddo
1502           enddo
1503         endif
1504 #define DEBUG
1505 #ifdef DEBUG
1506         write (iout,*) "sc_loc_geom: Array C"
1507         do i=1,nres
1508           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),&
1509            (c(j,i+nres),j=1,3)
1510         enddo
1511         write (iout,*) "Array Cref"
1512         do i=1,nres
1513           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),&
1514            (cref(j,i+nres,1),j=1,3)
1515         enddo
1516 #endif
1517        call int_from_cart1(.false.)
1518        call sc_loc_geom(.false.)
1519        do i=1,nres
1520          thetaref(i)=theta(i)
1521          phiref(i)=phi(i)
1522        enddo
1523        do i=1,nres-1
1524          do j=1,3
1525            dc(j,i)=c(j,i+1)-c(j,i)
1526            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1527          enddo
1528        enddo
1529        do i=2,nres-1
1530          do j=1,3
1531            dc(j,i+nres)=c(j,i+nres)-c(j,i)
1532            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1533          enddo
1534        enddo
1535       else
1536         homol_nset=0
1537         if (start_from_model) then
1538           nmodel_start=0
1539           do
1540             read(inp,'(a)',end=332,err=332) pdbfile
1541             if (me.eq.king .or. .not. out1file)&
1542              write (iout,'(a,5x,a)') 'Opening PDB file',&
1543              pdbfile(:ilen(pdbfile))
1544             open(ipdbin,file=pdbfile,status='old',err=336)
1545             goto 335
1546  336        write (iout,'(a,5x,a)') 'Error opening PDB file',&
1547            pdbfile(:ilen(pdbfile))
1548             call flush(iout)
1549             stop
1550  335        continue
1551             unres_pdb=.false.
1552             nres_temp=nres
1553 !            call readpdb
1554             call readpdb_template(nmodel_start+1)
1555             close(ipdbin)
1556             if (nres.ge.nres_temp) then
1557               nmodel_start=nmodel_start+1
1558               pdbfiles_chomo(nmodel_start)=pdbfile
1559               do i=1,2*nres
1560                 do j=1,3
1561                   chomo(j,i,nmodel_start)=c(j,i)
1562                 enddo
1563               enddo
1564             else
1565               if (me.eq.king .or. .not. out1file) &
1566                write (iout,'(a,2i7,1x,a)') &
1567                 "Different number of residues",nres_temp,nres, &
1568                 " model skipped."
1569             endif
1570             nres=nres_temp
1571           enddo
1572   332     continue
1573           if (nmodel_start.eq.0) then
1574             if (me.eq.king .or. .not. out1file) &
1575              write (iout,'(a)') &
1576              "No valid starting model found START_FROM_MODELS is OFF"
1577               start_from_model=.false.
1578           endif
1579           write (iout,*) "nmodel_start",nmodel_start
1580         endif
1581       endif
1582
1583       endif
1584         if (constr_dist.gt.0) call read_dist_constr
1585         write (iout,*) "After read_dist_constr nhpb",nhpb
1586         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1587         if (velnanoconst.ne.0) call read_afmnano
1588         call hpb_partition
1589
1590       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1591           .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1592           modecalc.ne.10) then
1593 ! If input structure hasn't been supplied from the PDB file read or generate
1594 ! initial geometry.
1595         if (iranconf.eq.0 .and. .not. extconf) then
1596           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1597            write (iout,'(a)') 'Initial geometry will be read in.'
1598           if (read_cart) then
1599             read(inp,'(8f10.5)',end=36,err=36) &
1600              ((c(l,k),l=1,3),k=1,nres),&
1601              ((c(l,k+nres),l=1,3),k=nnt,nct)
1602             write (iout,*) "Exit READ_CART"
1603             write (iout,'(8f10.5)') &
1604              ((c(l,k),l=1,3),k=1,nres)
1605             write (iout,'(8f10.5)') &
1606              ((c(l,k+nres),l=1,3),k=nnt,nct)
1607             call int_from_cart1(.true.)
1608             write (iout,*) "Finish INT_TO_CART"
1609             do i=1,nres-1
1610               do j=1,3
1611                 dc(j,i)=c(j,i+1)-c(j,i)
1612                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1613               enddo
1614             enddo
1615             do i=nnt,nct
1616               if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1617                 do j=1,3
1618                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1619                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1620                 enddo
1621               endif
1622             enddo
1623             return
1624           else
1625            write(iout,*) "read angles from input"
1626            call read_angles(inp,*36)
1627             call chainbuild
1628
1629           endif
1630           goto 37
1631    36     write (iout,'(a)') 'Error reading angle file.'
1632 #ifdef MPI
1633           call mpi_finalize( MPI_COMM_WORLD,IERR )
1634 #endif
1635           stop 'Error reading angle file.'
1636    37     continue
1637         else if (extconf) then
1638          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1639           write (iout,'(a)') 'Extended chain initial geometry.'
1640          do i=3,nres
1641           theta(i)=90d0*deg2rad
1642           if (molnum(i).eq.2) theta(i)=160d0*deg2rad
1643          enddo
1644          do i=4,nres
1645           phi(i)=180d0*deg2rad
1646           if (molnum(i).eq.2) phi(i)=30d0*deg2rad
1647          enddo
1648          do i=2,nres-1
1649           alph(i)=110d0*deg2rad
1650           if (molnum(i).eq.2) alph(i)=30d0*deg2rad
1651          enddo
1652          do i=2,nres-1
1653           omeg(i)=-120d0*deg2rad
1654           if (molnum(i).eq.2) omeg(i)=60d0*deg2rad
1655           if (itype(i,1).le.0) omeg(i)=-omeg(i)
1656          enddo
1657          call chainbuild
1658         else
1659           if(me.eq.king.or..not.out1file) &
1660            write (iout,'(a)') 'Random-generated initial geometry.'
1661
1662
1663 #ifdef MPI
1664           if (me.eq.king  .or. fg_rank.eq.0 .and. &
1665                  ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1666 #endif
1667             do itrial=1,100
1668               itmp=1
1669               call gen_rand_conf(itmp,*30)
1670               goto 40
1671    30         write (iout,*) 'Failed to generate random conformation',&
1672                 ', itrial=',itrial
1673               write (*,*) 'Processor:',me,&
1674                 ' Failed to generate random conformation',&
1675                 ' itrial=',itrial
1676               call intout
1677
1678 #ifdef AIX
1679               call flush_(iout)
1680 #else
1681               call flush(iout)
1682 #endif
1683             enddo
1684             write (iout,'(a,i3,a)') 'Processor:',me,&
1685               ' error in generating random conformation.'
1686             write (*,'(a,i3,a)') 'Processor:',me,&
1687               ' error in generating random conformation.'
1688             call flush(iout)
1689 #ifdef MPI
1690             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1691    40       continue
1692           endif
1693 #else
1694           do itrial=1,100
1695             itmp=1
1696             call gen_rand_conf(itmp,*335)
1697             goto 40
1698   335       write (iout,*) 'Failed to generate random conformation',&
1699               ', itrial=',itrial
1700             write (*,*) 'Failed to generate random conformation',&
1701               ', itrial=',itrial
1702           enddo
1703           write (iout,'(a,i3,a)') 'Processor:',me,&
1704             ' error in generating random conformation.'
1705           write (*,'(a,i3,a)') 'Processor:',me,&
1706             ' error in generating random conformation.'
1707           stop
1708    40     continue
1709 #endif
1710         endif
1711       elseif (modecalc.eq.4) then
1712         read (inp,'(a)') intinname
1713         open (intin,file=intinname,status='old',err=333)
1714         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1715         write (iout,'(a)') 'intinname',intinname
1716         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1717         goto 334
1718   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1719 #ifdef MPI
1720         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1721 #endif
1722         stop 'Error opening angle file.'
1723   334   continue
1724
1725       endif
1726 ! Generate distance constraints, if the PDB structure is to be regularized.
1727       if (nthread.gt.0) then
1728         call read_threadbase
1729       endif
1730       call setup_var
1731       if (me.eq.king .or. .not. out1file) &
1732        call intout
1733       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1734         write (iout,'(/a,i3,a)') &
1735         'The chain contains',ns,' disulfide-bridging cysteines.'
1736         write (iout,'(20i4)') (iss(i),i=1,ns)
1737        if (dyn_ss) then
1738           write(iout,*)"Running with dynamic disulfide-bond formation"
1739        else
1740         write (iout,'(/a/)') 'Pre-formed links are:'
1741         do i=1,nss
1742           i1=ihpb(i)-nres
1743           i2=jhpb(i)-nres
1744           it1=itype(i1,1)
1745           it2=itype(i2,1)
1746           if (me.eq.king.or..not.out1file) &
1747           write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1748           restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1749           ebr,forcon(i)
1750         enddo
1751         write (iout,'(a)')
1752        endif
1753       endif
1754       if (ns.gt.0.and.dyn_ss) then
1755           do i=nss+1,nhpb
1756             ihpb(i-nss)=ihpb(i)
1757             jhpb(i-nss)=jhpb(i)
1758             forcon(i-nss)=forcon(i)
1759             dhpb(i-nss)=dhpb(i)
1760           enddo
1761           nhpb=nhpb-nss
1762           nss=0
1763           call hpb_partition
1764           do i=1,ns
1765             dyn_ss_mask(iss(i))=.true.
1766           enddo
1767       endif
1768       if (i2ndstr.gt.0) call secstrp2dihc
1769       if (indpdb.gt.0) then
1770           write(iout,*) "WCHODZE TU!!"
1771           call int_from_cart1(.true.)
1772       endif
1773 !      call geom_to_var(nvar,x)
1774 !      call etotal(energia(0))
1775 !      call enerprint(energia(0))
1776 !      call briefout(0,etot)
1777 !      stop
1778 !d    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1779 !d    write (iout,'(a)') 'Variable list:'
1780 !d    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1781 #ifdef MPI
1782       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1783         write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1784         'Processor',myrank,': end reading molecular data.'
1785 #endif
1786       return
1787       end subroutine molread
1788 !-----------------------------------------------------------------------------
1789       subroutine read_constr_homology
1790       use control, only:init_int_table,homology_partition
1791       use MD_data, only:iset
1792 !      implicit none
1793 !      include 'DIMENSIONS'
1794 !#ifdef MPI
1795 !      include 'mpif.h'
1796 !#endif
1797 !      include 'COMMON.SETUP'
1798 !      include 'COMMON.CONTROL'
1799 !      include 'COMMON.HOMOLOGY'
1800 !      include 'COMMON.CHAIN'
1801 !      include 'COMMON.IOUNITS'
1802 !      include 'COMMON.MD'
1803 !      include 'COMMON.QRESTR'
1804 !      include 'COMMON.GEO'
1805 !      include 'COMMON.INTERACT'
1806 !      include 'COMMON.NAMES'
1807 !      include 'COMMON.VAR'
1808 !
1809
1810 !     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1811 !    &                 dist_cut
1812 !     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1813 !    &    sigma_odl_temp(maxres,maxres,max_template)
1814       character*2 kic2
1815       character*24 model_ki_dist, model_ki_angle
1816       character*500 controlcard
1817       integer :: ki,i,ii,j,k,l
1818       integer, dimension (:), allocatable :: ii_in_use
1819       integer :: i_tmp,idomain_tmp,&
1820       irec,ik,iistart,nres_temp
1821 !      integer :: iset
1822 !      external :: ilen
1823       logical :: liiflag,lfirst
1824       integer :: i01,i10
1825 !
1826 !     FP - Nov. 2014 Temporary specifications for new vars
1827 !
1828       real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
1829                        rescore3_tmp, dist_cut
1830       real(kind=8), dimension (:,:),allocatable :: rescore
1831       real(kind=8), dimension (:,:),allocatable :: rescore2
1832       real(kind=8), dimension (:,:),allocatable :: rescore3
1833       real(kind=8) :: distal
1834       character*24 tpl_k_rescore
1835       character*256 pdbfile
1836
1837 ! -----------------------------------------------------------------
1838 ! Reading multiple PDB ref structures and calculation of retraints
1839 ! not using pre-computed ones stored in files model_ki_{dist,angle}
1840 ! FP (Nov., 2014)
1841 ! -----------------------------------------------------------------
1842 !
1843 !
1844 ! Alternative: reading from input
1845       call card_concat(controlcard,.true.)
1846       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1847       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1848       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1849       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1850       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1851       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1852       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1853       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1854       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
1855       if(.not.read2sigma.and.start_from_model) then
1856           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
1857            write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
1858           start_from_model=.false.
1859           iranconf=(indpdb.le.0)
1860       endif
1861       if(start_from_model .and. (me.eq.king .or. .not. out1file))&
1862          write(iout,*) 'START_FROM_MODELS is ON'
1863 !      if(start_from_model .and. rest) then
1864 !        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1865 !         write(iout,*) 'START_FROM_MODELS is OFF'
1866 !         write(iout,*) 'remove restart keyword from input'
1867 !        endif
1868 !      endif
1869       if (start_from_model) nmodel_start=constr_homology
1870       if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
1871       if (homol_nset.gt.1)then
1872          call card_concat(controlcard,.true.)
1873          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1874          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1875 !          write(iout,*) "iset homology_weight "
1876           do i=1,homol_nset
1877            write(iout,*) i,waga_homology(i)
1878           enddo
1879          endif
1880          iset=mod(kolor,homol_nset)+1
1881       else
1882        iset=1
1883        waga_homology(1)=1.0
1884       endif
1885
1886 !d      write (iout,*) "nnt",nnt," nct",nct
1887 !d      call flush(iout)
1888
1889       if (read_homol_frag) then
1890         call read_klapaucjusz
1891       else
1892
1893       lim_odl=0
1894       lim_dih=0
1895 !
1896 !      write(iout,*) 'nnt=',nnt,'nct=',nct
1897 !
1898 !      do i = nnt,nct
1899 !        do k=1,constr_homology
1900 !         idomain(k,i)=0
1901 !        enddo
1902 !      enddo
1903        idomain=0
1904
1905 !      ii=0
1906 !      do i = nnt,nct-2
1907 !        do j=i+2,nct
1908 !        ii=ii+1
1909 !        ii_in_use(ii)=0
1910 !        enddo
1911 !      enddo
1912       ii_in_use=0
1913       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
1914       if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
1915
1916       do k=1,constr_homology
1917
1918         read(inp,'(a)') pdbfile
1919         pdbfiles_chomo(k)=pdbfile
1920         if(me.eq.king .or. .not. out1file) &
1921          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
1922         pdbfile(:ilen(pdbfile))
1923         open(ipdbin,file=pdbfile,status='old',err=33)
1924         goto 34
1925   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
1926         pdbfile(:ilen(pdbfile))
1927         stop
1928   34    continue
1929 !        print *,'Begin reading pdb data'
1930 !
1931 ! Files containing res sim or local scores (former containing sigmas)
1932 !
1933
1934         write(kic2,'(bz,i2.2)') k
1935
1936         tpl_k_rescore="template"//kic2//".sco"
1937         write(iout,*) "tpl_k_rescore=",tpl_k_rescore
1938         unres_pdb=.false.
1939         nres_temp=nres
1940         write(iout,*) "read2sigma",read2sigma
1941
1942         if (read2sigma) then
1943           call readpdb_template(k)
1944         else
1945           call readpdb
1946         endif
1947         write(iout,*) "after readpdb"
1948         if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
1949         nres_chomo(k)=nres
1950         nres=nres_temp
1951         if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
1952         if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
1953         if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
1954         if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
1955         if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
1956         if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
1957         if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
1958         if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
1959         if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
1960         if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
1961         if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
1962         if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
1963         if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1964         if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
1965 !        if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1966         if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
1967         if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
1968         if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
1969         if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
1970 !        if(.not.allocated(distance)) allocate(distance(constr_homology))
1971 !        if(.not.allocated(distancek)) allocate(distancek(constr_homology))
1972
1973
1974 !
1975 !     Distance restraints
1976 !
1977 !          ... --> odl(k,ii)
1978 ! Copy the coordinates from reference coordinates (?)
1979         do i=1,2*nres_chomo(k)
1980           do j=1,3
1981             c(j,i)=cref(j,i,1)
1982 !           write (iout,*) "c(",j,i,") =",c(j,i)
1983           enddo
1984         enddo
1985 !
1986 ! From read_dist_constr (commented out 25/11/2014 <-> res sim)
1987 !
1988 !         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1989           open (ientin,file=tpl_k_rescore,status='old')
1990           if (nnt.gt.1) rescore(k,1)=0.0d0
1991           do irec=nnt,nct ! loop for reading res sim
1992             if (read2sigma) then
1993              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
1994                                      rescore3_tmp,idomain_tmp
1995              i_tmp=i_tmp+nnt-1
1996              idomain(k,i_tmp)=idomain_tmp
1997              rescore(k,i_tmp)=rescore_tmp
1998              rescore2(k,i_tmp)=rescore2_tmp
1999              rescore3(k,i_tmp)=rescore3_tmp
2000              if (.not. out1file .or. me.eq.king)&
2001              write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
2002                            i_tmp,rescore2_tmp,rescore_tmp,&
2003                                      rescore3_tmp,idomain_tmp
2004             else
2005              idomain(k,irec)=1
2006              read (ientin,*,end=1401) rescore_tmp
2007
2008 !           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2009              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2010 !           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2011             endif
2012           enddo
2013  1401   continue
2014         close (ientin)
2015         if (waga_dist.ne.0.0d0) then
2016           ii=0
2017           do i = nnt,nct-2
2018             do j=i+2,nct
2019
2020               x12=c(1,i)-c(1,j)
2021               y12=c(2,i)-c(2,j)
2022               z12=c(3,i)-c(3,j)
2023               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2024 !              write (iout,*) k,i,j,distal,dist2_cut
2025
2026             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2027                  .and. distal.le.dist2_cut ) then
2028
2029               ii=ii+1
2030               ii_in_use(ii)=1
2031               l_homo(k,ii)=.true.
2032
2033 !             write (iout,*) "k",k
2034 !             write (iout,*) "i",i," j",j," constr_homology",
2035 !    &                       constr_homology
2036               ires_homo(ii)=i
2037               jres_homo(ii)=j
2038               odl(k,ii)=distal
2039               if (read2sigma) then
2040                 sigma_odl(k,ii)=0
2041                 do ik=i,j
2042                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2043                 enddo
2044                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2045                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2046               sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2047               else
2048                 if (odl(k,ii).le.dist_cut) then
2049                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2050                 else
2051 #ifdef OLDSIGMA
2052                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2053                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2054 #else
2055                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2056                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2057 #endif
2058                 endif
2059               endif
2060               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2061             else
2062 !              ii=ii+1
2063 !              l_homo(k,ii)=.false.
2064             endif
2065             enddo
2066           enddo
2067         lim_odl=ii
2068         endif
2069 !        write (iout,*) "Distance restraints set"
2070 !        call flush(iout)
2071 !
2072 !     Theta, dihedral and SC retraints
2073 !
2074         if (waga_angle.gt.0.0d0) then
2075 !         open (ientin,file=tpl_k_sigma_dih,status='old')
2076 !         do irec=1,maxres-3 ! loop for reading sigma_dih
2077 !            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2078 !            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2079 !            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2080 !    &                            sigma_dih(k,i+nnt-1)
2081 !         enddo
2082 !1402   continue
2083 !         close (ientin)
2084           do i = nnt+3,nct
2085             if (idomain(k,i).eq.0) then
2086                sigma_dih(k,i)=0.0
2087                cycle
2088             endif
2089             dih(k,i)=phiref(i) ! right?
2090 !           read (ientin,*) sigma_dih(k,i) ! original variant
2091 !             write (iout,*) "dih(",k,i,") =",dih(k,i)
2092 !             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2093 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2094 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
2095 !    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
2096
2097             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2098                           rescore(k,i-2)+rescore(k,i-3))/4.0
2099 !            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2100 !           write (iout,*) "Raw sigmas for dihedral angle restraints"
2101 !           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2102 !           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2103 !                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
2104 !   Instead of res sim other local measure of b/b str reliability possible
2105             if (sigma_dih(k,i).ne.0) &
2106             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2107 !           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2108           enddo
2109           lim_dih=nct-nnt-2
2110         endif
2111 !        write (iout,*) "Dihedral angle restraints set"
2112 !        call flush(iout)
2113
2114         if (waga_theta.gt.0.0d0) then
2115 !         open (ientin,file=tpl_k_sigma_theta,status='old')
2116 !         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2117 !            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2118 !            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2119 !    &                              sigma_theta(k,i+nnt-1)
2120 !         enddo
2121 !1403   continue
2122 !         close (ientin)
2123
2124           do i = nnt+2,nct ! right? without parallel.
2125 !         do i = i=1,nres ! alternative for bounds acc to readpdb?
2126 !         do i=ithet_start,ithet_end ! with FG parallel.
2127              if (idomain(k,i).eq.0) then
2128               sigma_theta(k,i)=0.0
2129               cycle
2130              endif
2131              thetatpl(k,i)=thetaref(i)
2132 !            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2133 !            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
2134 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2135 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
2136 !            read (ientin,*) sigma_theta(k,i) ! 1st variant
2137              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2138                              rescore(k,i-2))/3.0
2139 !             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2140              if (sigma_theta(k,i).ne.0) &
2141              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2142
2143 !            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2144 !                             rescore(k,i-2) !  right expression ?
2145 !            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2146           enddo
2147         endif
2148 !        write (iout,*) "Angle restraints set"
2149 !        call flush(iout)
2150
2151         if (waga_d.gt.0.0d0) then
2152 !       open (ientin,file=tpl_k_sigma_d,status='old')
2153 !         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2154 !            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2155 !            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2156 !    &                          sigma_d(k,i+nnt-1)
2157 !         enddo
2158 !1404   continue
2159
2160           do i = nnt,nct ! right? without parallel.
2161 !         do i=2,nres-1 ! alternative for bounds acc to readpdb?
2162 !         do i=loc_start,loc_end ! with FG parallel.
2163                if (itype(i,1).eq.10) cycle
2164                if (idomain(k,i).eq.0 ) then
2165                   sigma_d(k,i)=0.0
2166                   cycle
2167                endif
2168                xxtpl(k,i)=xxref(i)
2169                yytpl(k,i)=yyref(i)
2170                zztpl(k,i)=zzref(i)
2171 !              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2172 !              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2173 !              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2174 !              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
2175                sigma_d(k,i)=rescore3(k,i) !  right expression ?
2176                if (sigma_d(k,i).ne.0) &
2177                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2178
2179 !              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
2180 !              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2181 !              read (ientin,*) sigma_d(k,i) ! 1st variant
2182           enddo
2183         endif
2184       enddo
2185 !      write (iout,*) "SC restraints set"
2186 !      call flush(iout)
2187 !
2188 ! remove distance restraints not used in any model from the list
2189 ! shift data in all arrays
2190 !
2191 !      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
2192       if (waga_dist.ne.0.0d0) then
2193         ii=0
2194         liiflag=.true.
2195         lfirst=.true.
2196         do i=nnt,nct-2
2197          do j=i+2,nct
2198           ii=ii+1
2199 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2200 !     &            .and. distal.le.dist2_cut ) then
2201 !          write (iout,*) "i",i," j",j," ii",ii
2202 !          call flush(iout)
2203           if (ii_in_use(ii).eq.0.and.liiflag.or. &
2204           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2205             liiflag=.false.
2206             i10=ii
2207             if (lfirst) then
2208               lfirst=.false.
2209               iistart=ii
2210             else
2211               if(i10.eq.lim_odl) i10=i10+1
2212               do ki=0,i10-i01-1
2213                ires_homo(iistart+ki)=ires_homo(ki+i01)
2214                jres_homo(iistart+ki)=jres_homo(ki+i01)
2215                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2216                do k=1,constr_homology
2217                 odl(k,iistart+ki)=odl(k,ki+i01)
2218                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2219                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2220                enddo
2221               enddo
2222               iistart=iistart+i10-i01
2223             endif
2224           endif
2225           if (ii_in_use(ii).ne.0.and..not.liiflag) then
2226              i01=ii
2227              liiflag=.true.
2228           endif
2229          enddo
2230         enddo
2231         lim_odl=iistart-1
2232       endif
2233 !      write (iout,*) "Removing distances completed"
2234 !      call flush(iout)
2235       endif ! .not. klapaucjusz
2236
2237       if (constr_homology.gt.0) call homology_partition
2238       write (iout,*) "After homology_partition"
2239       call flush(iout)
2240       if (constr_homology.gt.0) call init_int_table
2241       write (iout,*) "After init_int_table"
2242       call flush(iout)
2243 !      endif ! .not. klapaucjusz
2244 !      endif
2245 !      if (constr_homology.gt.0) call homology_partition
2246 !      write (iout,*) "After homology_partition"
2247 !      call flush(iout)
2248 !      if (constr_homology.gt.0) call init_int_table
2249 !      write (iout,*) "After init_int_table"
2250 !      call flush(iout)
2251 !      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2252 !      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2253 !
2254 ! Print restraints
2255 !
2256       if (.not.out_template_restr) return
2257 !d      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2258       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2259        write (iout,*) "Distance restraints from templates"
2260        do ii=1,lim_odl
2261        write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
2262         ii,ires_homo(ii),jres_homo(ii),&
2263         (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
2264         ki=1,constr_homology)
2265        enddo
2266        write (iout,*) "Dihedral angle restraints from templates"
2267        do i=nnt+3,nct
2268         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2269             (rad2deg*dih(ki,i),&
2270             rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2271        enddo
2272        write (iout,*) "Virtual-bond angle restraints from templates"
2273        do i=nnt+2,nct
2274         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2275             (rad2deg*thetatpl(ki,i),&
2276             rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2277        enddo
2278        write (iout,*) "SC restraints from templates"
2279        do i=nnt,nct
2280         write(iout,'(i7,100(4f8.2,4x))') i,&
2281         (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
2282          1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2283        enddo
2284       endif
2285       return
2286       end subroutine read_constr_homology
2287 !-----------------------------------------------------------------------------
2288       subroutine read_klapaucjusz
2289       use energy_data
2290       implicit none
2291 !     include 'DIMENSIONS'
2292 !#ifdef MPI
2293 !     include 'mpif.h'
2294 !#endif
2295 !     include 'COMMON.SETUP'
2296 !     include 'COMMON.CONTROL'
2297 !     include 'COMMON.HOMOLOGY'
2298 !     include 'COMMON.CHAIN'
2299 !     include 'COMMON.IOUNITS'
2300 !     include 'COMMON.MD'
2301 !     include 'COMMON.GEO'
2302 !     include 'COMMON.INTERACT'
2303 !     include 'COMMON.NAMES'
2304       character(len=256) fragfile
2305       integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,&
2306                          ii_in_use
2307       integer, dimension(:,:), allocatable :: iresclust,inclust
2308       integer :: nclust
2309
2310       character(len=2) :: kic2
2311       character(len=24) :: model_ki_dist, model_ki_angle
2312       character(len=500) :: controlcard
2313       integer :: ki, i, j, jj,k, l, i_tmp,&
2314       idomain_tmp,&
2315       ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,&
2316       i01,i10,nnt_chain,nct_chain
2317       real(kind=8) :: distal
2318       logical :: lprn = .true.
2319       integer :: nres_temp
2320 !      integer :: ilen
2321 !      external :: ilen
2322       logical :: liiflag,lfirst
2323
2324       real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut
2325       real(kind=8), dimension (:,:), allocatable  :: rescore
2326       real(kind=8), dimension (:,:), allocatable :: rescore2
2327       character(len=24) :: tpl_k_rescore
2328       character(len=256) :: pdbfile
2329
2330 !
2331 ! For new homol impl
2332 !
2333 !     include 'COMMON.VAR'
2334 !
2335 !      write (iout,*) "READ_KLAPAUCJUSZ"
2336 !      print *,"READ_KLAPAUCJUSZ"
2337 !      call flush(iout)
2338       call getenv("FRAGFILE",fragfile)
2339       write (iout,*) "Opening", fragfile
2340       call flush(iout)
2341       open(ientin,file=fragfile,status="old",err=10)
2342 !      write (iout,*) " opened"
2343 !      call flush(iout)
2344
2345       sigma_theta=0.0
2346       sigma_d=0.0
2347       sigma_dih=0.0
2348       l_homo = .false.
2349
2350       nres_temp=nres
2351       itype_temp(:)=itype(:,1)
2352       ii=0
2353       lim_odl=0
2354
2355 !      write (iout,*) "Entering loop"
2356 !      call flush(iout)
2357
2358       DO IGR = 1,NCHAIN_GROUP
2359
2360 !      write (iout,*) "igr",igr
2361       call flush(iout)
2362       read(ientin,*) constr_homology,nclust
2363       if (start_from_model) then
2364         nmodel_start=constr_homology
2365       else
2366         nmodel_start=0
2367       endif
2368
2369       ii_old=lim_odl
2370
2371       ichain=iequiv(1,igr)
2372       nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
2373       nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
2374 !      write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
2375
2376 ! Read pdb files
2377       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
2378       if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
2379       do k=1,constr_homology
2380         read(ientin,'(a)') pdbfile
2381         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', &
2382         pdbfile(:ilen(pdbfile))
2383         pdbfiles_chomo(k)=pdbfile
2384         open(ipdbin,file=pdbfile,status='old',err=33)
2385         goto 34
2386   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
2387         pdbfile(:ilen(pdbfile))
2388         stop
2389   34    continue
2390         unres_pdb=.false.
2391 !        nres_temp=nres
2392         call readpdb_template(k)
2393         nres_chomo(k)=nres
2394 !        nres=nres_temp
2395         do i=1,nres
2396           rescore(k,i)=0.2d0
2397           rescore2(k,i)=1.0d0
2398         enddo
2399       enddo
2400 ! Read clusters
2401       do i=1,nclust
2402         read(ientin,*) ninclust(i),nresclust(i)
2403         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
2404         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
2405       enddo
2406 !
2407 ! Loop over clusters
2408 !
2409       do l=1,nclust
2410         do ll = 1,ninclust(l)
2411
2412         k = inclust(ll,l)
2413 !        write (iout,*) "l",l," ll",ll," k",k
2414         do i=1,nres
2415           idomain(k,i)=0
2416         enddo
2417         do i=1,nresclust(l)
2418           if (nnt.gt.1)  then
2419             idomain(k,iresclust(i,l)+1) = 1
2420           else
2421             idomain(k,iresclust(i,l)) = 1
2422           endif
2423         enddo
2424 !
2425 !     Distance restraints
2426 !
2427 !          ... --> odl(k,ii)
2428 ! Copy the coordinates from reference coordinates (?)
2429 !        nres_temp=nres
2430         nres=nres_chomo(k)
2431         do i=1,2*nres
2432           do j=1,3
2433             c(j,i)=chomo(j,i,k)
2434 !           write (iout,*) "c(",j,i,") =",c(j,i)
2435           enddo
2436         enddo
2437         call int_from_cart(.true.,.false.)
2438         call sc_loc_geom(.false.)
2439         do i=1,nres
2440           thetaref(i)=theta(i)
2441           phiref(i)=phi(i)
2442         enddo
2443 !        nres=nres_temp
2444         if (waga_dist.ne.0.0d0) then
2445           ii=ii_old
2446 !          do i = nnt,nct-2
2447           do i = nnt_chain,nct_chain-2
2448 !            do j=i+2,nct
2449             do j=i+2,nct_chain
2450
2451               x12=c(1,i)-c(1,j)
2452               y12=c(2,i)-c(2,j)
2453               z12=c(3,i)-c(3,j)
2454               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2455 !              write (iout,*) k,i,j,distal,dist2_cut
2456
2457             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2458                  .and. distal.le.dist2_cut ) then
2459
2460               ii=ii+1
2461               ii_in_use(ii)=1
2462               l_homo(k,ii)=.true.
2463
2464 !             write (iout,*) "k",k
2465 !             write (iout,*) "i",i," j",j," constr_homology",
2466 !     &                       constr_homology
2467               ires_homo(ii)=i+chain_border1(1,igr)-1
2468               jres_homo(ii)=j+chain_border1(1,igr)-1
2469               odl(k,ii)=distal
2470               if (read2sigma) then
2471                 sigma_odl(k,ii)=0
2472                 do ik=i,j
2473                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2474                 enddo
2475                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2476                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2477              sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2478               else
2479                 if (odl(k,ii).le.dist_cut) then
2480                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2481                 else
2482 #ifdef OLDSIGMA
2483                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2484                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2485 #else
2486                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2487                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2488 #endif
2489                 endif
2490               endif
2491               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2492             else
2493               ii=ii+1
2494 !              l_homo(k,ii)=.false.
2495             endif
2496             enddo
2497           enddo
2498         lim_odl=ii
2499         endif
2500 !
2501 !     Theta, dihedral and SC retraints
2502 !
2503         if (waga_angle.gt.0.0d0) then
2504           do i = nnt_chain+3,nct_chain
2505             iii=i+chain_border1(1,igr)-1
2506             if (idomain(k,i).eq.0) then
2507 !               sigma_dih(k,i)=0.0
2508                cycle
2509             endif
2510             dih(k,iii)=phiref(i)
2511             sigma_dih(k,iii)= &
2512                (rescore(k,i)+rescore(k,i-1)+ &
2513                            rescore(k,i-2)+rescore(k,i-3))/4.0
2514 !            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
2515 !     &       " sigma_dihed",sigma_dih(k,i)
2516             if (sigma_dih(k,iii).ne.0) &
2517              sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
2518           enddo
2519 !          lim_dih=nct-nnt-2
2520         endif
2521
2522         if (waga_theta.gt.0.0d0) then
2523           do i = nnt_chain+2,nct_chain
2524              iii=i+chain_border1(1,igr)-1
2525              if (idomain(k,i).eq.0) then
2526 !              sigma_theta(k,i)=0.0
2527               cycle
2528              endif
2529              thetatpl(k,iii)=thetaref(i)
2530              sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ &
2531                               rescore(k,i-2))/3.0
2532              if (sigma_theta(k,iii).ne.0) &
2533              sigma_theta(k,iii)=1.0d0/ &
2534              (sigma_theta(k,iii)*sigma_theta(k,iii))
2535           enddo
2536         endif
2537
2538         if (waga_d.gt.0.0d0) then
2539           do i = nnt_chain,nct_chain
2540              iii=i+chain_border1(1,igr)-1
2541                if (itype(i,1).eq.10) cycle
2542                if (idomain(k,i).eq.0 ) then
2543 !                  sigma_d(k,i)=0.0
2544                   cycle
2545                endif
2546                xxtpl(k,iii)=xxref(i)
2547                yytpl(k,iii)=yyref(i)
2548                zztpl(k,iii)=zzref(i)
2549                sigma_d(k,iii)=rescore(k,i)
2550                if (sigma_d(k,iii).ne.0) &
2551                 sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
2552 !               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
2553           enddo
2554         endif
2555       enddo ! l
2556       enddo ! ll
2557 !
2558 ! remove distance restraints not used in any model from the list
2559 ! shift data in all arrays
2560 !
2561 !      write (iout,*) "ii_old",ii_old
2562       if (waga_dist.ne.0.0d0) then
2563 #ifdef DEBUG
2564        write (iout,*) "Distance restraints from templates"
2565        do iii=1,lim_odl
2566        write(iout,'(4i5,100(2f8.2,1x,l1,4x))') &
2567         iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), &
2568         (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), &
2569         ki=1,constr_homology)
2570        enddo
2571 #endif
2572         ii=ii_old
2573         liiflag=.true.
2574         lfirst=.true.
2575         do i=nnt_chain,nct_chain-2
2576          do j=i+2,nct_chain
2577           ii=ii+1
2578 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2579 !     &            .and. distal.le.dist2_cut ) then
2580 !          write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
2581 !          call flush(iout)
2582           if (ii_in_use(ii).eq.0.and.liiflag.or. &
2583           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2584             liiflag=.false.
2585             i10=ii
2586             if (lfirst) then
2587               lfirst=.false.
2588               iistart=ii
2589             else
2590               if(i10.eq.lim_odl) i10=i10+1
2591               do ki=0,i10-i01-1
2592                ires_homo(iistart+ki)=ires_homo(ki+i01)
2593                jres_homo(iistart+ki)=jres_homo(ki+i01)
2594                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2595                do k=1,constr_homology
2596                 odl(k,iistart+ki)=odl(k,ki+i01)
2597                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2598                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2599                enddo
2600               enddo
2601               iistart=iistart+i10-i01
2602             endif
2603           endif
2604           if (ii_in_use(ii).ne.0.and..not.liiflag) then
2605              i01=ii
2606              liiflag=.true.
2607           endif
2608          enddo
2609         enddo
2610         lim_odl=iistart-1
2611       endif
2612
2613       lll=lim_odl-ii_old
2614
2615       do i=2,nequiv(igr)
2616
2617         ichain=iequiv(i,igr)
2618
2619         do j=nnt_chain,nct_chain
2620           jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
2621           do k=1,constr_homology
2622             dih(k,jj)=dih(k,j)
2623             sigma_dih(k,jj)=sigma_dih(k,j)
2624             thetatpl(k,jj)=thetatpl(k,j)
2625             sigma_theta(k,jj)=sigma_theta(k,j)
2626             xxtpl(k,jj)=xxtpl(k,j)
2627             yytpl(k,jj)=yytpl(k,j)
2628             zztpl(k,jj)=zztpl(k,j)
2629             sigma_d(k,jj)=sigma_d(k,j)
2630           enddo
2631         enddo
2632
2633         jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
2634 !        write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
2635         do j=ii_old+1,lim_odl
2636           ires_homo(j+lll)=ires_homo(j)+jj
2637           jres_homo(j+lll)=jres_homo(j)+jj
2638           do k=1,constr_homology
2639             odl(k,j+lll)=odl(k,j)
2640             sigma_odl(k,j+lll)=sigma_odl(k,j)
2641             l_homo(k,j+lll)=l_homo(k,j)
2642           enddo
2643         enddo
2644
2645         ii_old=ii_old+lll
2646         lim_odl=lim_odl+lll
2647
2648       enddo
2649
2650       ENDDO ! IGR
2651
2652       if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
2653       nres=nres_temp
2654       itype(:,1)=itype_temp(:)
2655
2656       return
2657    10 stop "Error in fragment file"
2658       end subroutine read_klapaucjusz
2659
2660 !-----------------------------------------------------------
2661       subroutine seq2chains
2662 !c
2663 !c Split the total UNRES sequence, which has dummy residues separating
2664 !c the chains, into separate chains. The length of  chain ichain is
2665 !c contained in chain_length(ichain), the first and last non-dummy
2666 !c residues are in chain_border(1,ichain) and chain_border(2,ichain),
2667 !c respectively. The lengths pertain to non-dummy residues only.
2668 !c
2669 !      implicit none
2670 !      include 'DIMENSIONS'
2671       use energy_data, only:molnum,nchain,chain_length,ireschain
2672       implicit none
2673 !      integer ireschain(nres)
2674       integer ii,ichain,i,j,mnum
2675       logical new_chain
2676       print *,"in seq2"
2677       ichain=1
2678       new_chain=.true.
2679       if (.not.allocated(chain_length)) allocate(chain_length(50))
2680       if (.not.allocated(chain_border)) allocate(chain_border(2,50))
2681
2682       chain_length(ichain)=0
2683       ii=1
2684       do while (ii.lt.nres)
2685         write(iout,*) "in seq2chains",ii,new_chain
2686         mnum=molnum(ii)
2687         if (itype(ii,mnum).eq.ntyp1_molec(mnum)) then
2688           if (.not.new_chain) then
2689             new_chain=.true.
2690             chain_border(2,ichain)=ii-1
2691             ichain=ichain+1
2692             chain_border(1,ichain)=ii+1
2693             chain_length(ichain)=0
2694           endif
2695         else
2696           if (new_chain) then
2697             chain_border(1,ichain)=ii
2698             new_chain=.false.
2699           endif
2700           chain_length(ichain)=chain_length(ichain)+1
2701         endif
2702         ii=ii+1
2703       enddo
2704       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) then
2705         ii=ii-1
2706       else
2707         chain_length(ichain)=chain_length(ichain)+1
2708       endif
2709       if (chain_length(ichain).gt.0) then
2710         chain_border(2,ichain)=ii
2711         nchain=ichain
2712       else
2713         nchain=ichain-1
2714       endif
2715       ireschain=0
2716       do i=1,nchain
2717         do j=chain_border(1,i),chain_border(2,i)
2718           ireschain(j)=i
2719         enddo
2720       enddo
2721       return
2722       end subroutine
2723 !---------------------------------------------------------------------
2724       subroutine chain_symmetry(npermchain,tabpermchain)
2725 !c
2726 !c Determine chain symmetry. nperm is the number of permutations and
2727 !c tabperchain contains the allowed permutations of the chains.
2728 !c
2729 !      implicit none
2730 !      include "DIMENSIONS"
2731 !      include "COMMON.IOUNITS"
2732       implicit none
2733       integer itemp(50),&
2734        npermchain,tabpermchain(50,5040),&
2735        tabperm(50,5040),mapchain(50),&
2736        iequiv(50,nres),iflag(nres)
2737       integer i,j,k,l,ii,nchain_group,nequiv(50),iieq,&
2738        nperm,npermc,ind,mnum
2739       if (nchain.eq.1) then
2740         npermchain=1
2741         tabpermchain(1,1)=1
2742 !c        print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
2743         return
2744       endif
2745 !c
2746 !c Look for equivalent chains
2747 #ifdef DEBUG
2748       write(iout,*) "nchain",nchain
2749       do i=1,nchain
2750         write(iout,*) "chain",i," from",chain_border(1,i),&
2751            " to",chain_border(2,i)
2752         write(iout,*)&
2753         "sequence ",(itype(j,molnum(j)),j=chain_border(1,i),chain_border(2,i))
2754       enddo
2755 #endif
2756       do i=1,nchain
2757         iflag(i)=0
2758       enddo
2759       nchain_group=0
2760       do i=1,nchain
2761         if (iflag(i).gt.0) cycle
2762         iflag(i)=1
2763         nchain_group=nchain_group+1
2764         iieq=1
2765         iequiv(iieq,nchain_group)=i
2766         do j=i+1,nchain
2767           if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
2768 !c          k=0
2769 !c          do while(k.lt.chain_length(i) .and.
2770 !c     &     itype(chain_border(1,i)+k).eq.itype(chain_border(1,j)+k))
2771           do k=0,chain_length(i)-1
2772 !c            k=k+1
2773             mnum=molnum(k+1)
2774             if (itype(chain_border(1,i)+k,mnum).ne.&
2775                itype(chain_border(1,j)+k,mnum)) exit
2776           enddo
2777           if (k.lt.chain_length(i)) cycle
2778           iflag(j)=1
2779           iieq=iieq+1
2780           iequiv(iieq,nchain_group)=j
2781         enddo
2782         nequiv(nchain_group)=iieq
2783       enddo
2784       write(iout,*) "Number of equivalent chain groups:",nchain_group
2785       write(iout,*) "Equivalent chain groups"
2786       do i=1,nchain_group
2787         write(iout,*) "group",i," #members",nequiv(i)," chains",&
2788            (iequiv(j,i),j=1,nequiv(i))
2789       enddo
2790       ind=0
2791       do i=1,nchain_group
2792         do j=1,nequiv(i)
2793           ind=ind+1
2794           mapchain(ind)=iequiv(j,i)
2795         enddo
2796       enddo
2797       write (iout,*) "mapchain"
2798       do i=1,nchain
2799         write (iout,*) i,mapchain(i)
2800       enddo
2801       ii=0
2802       do i=1,nchain_group
2803         call permut(nequiv(i),nperm,tabperm)
2804         if (ii.eq.0) then
2805           ii=nequiv(i)
2806           npermchain=nperm
2807           do j=1,nperm
2808             do k=1,ii
2809               tabpermchain(k,j)=iequiv(tabperm(k,j),i)
2810             enddo
2811           enddo
2812         else
2813           npermc=npermchain
2814           npermchain=npermchain*nperm
2815           ind=0
2816           do k=1,nperm
2817             do j=1,npermc
2818               ind=ind+1
2819               do l=1,ii
2820                 tabpermchain(l,ind)=tabpermchain(l,j)
2821               enddo
2822               do l=1,nequiv(i)
2823                 tabpermchain(ii+l,ind)=iequiv(tabperm(l,k),i)
2824               enddo
2825             enddo
2826           enddo
2827           ii=ii+nequiv(i)
2828         endif
2829       enddo
2830       do i=1,npermchain
2831         do j=1,nchain
2832           itemp(mapchain(j))=tabpermchain(j,i)
2833         enddo
2834         do j=1,nchain
2835           tabpermchain(j,i)=itemp(j)
2836         enddo
2837       enddo
2838       write(iout,*) "Number of chain permutations",npermchain
2839       write(iout,*) "Permutations"
2840       do i=1,npermchain
2841         write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
2842       enddo
2843       return
2844       end subroutine
2845 !c---------------------------------------------------------------------
2846       integer function tperm(i,iperm,tabpermchain)
2847 !      implicit none
2848 !      include 'DIMENSIONS'
2849       integer i,iperm
2850       integer tabpermchain(50,5040)
2851       if (i.eq.0) then
2852         tperm=0
2853       else
2854         tperm=tabpermchain(i,iperm)
2855       endif
2856       return
2857       end function
2858
2859 !-----------------------------------------------------------------------------
2860       end module io