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(.true.,ncont_ref,icont_ref,co)
 282
 283 !      do k=1,nres
 284 !       write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
 285 !      enddo
 286       close(icsa_native_int)
 287
 288   200 format (8f10.4)
 289
 290       return
 291       end subroutine from_int
 292 !-----------------------------------------------------------------------------
 293       subroutine dihang_to_c(aarray)
 294
 295       use geometry_data
 296       use csa_data
 297       use geometry, only:chainbuild
 298 !      implicit real*8 (a-h,o-z)
 299 !      include 'DIMENSIONS'
 300 !      include 'COMMON.CSA'
 301 !      include 'COMMON.BANK'
 302 !      include 'COMMON.CHAIN'
 303 !      include 'COMMON.GEO'
 304 !      include 'COMMON.VAR'
 305       integer :: i
 306       real(kind=8),dimension(mxang,nres,mxch) :: aarray !(mxang,maxres,mxch)
 307
 308 !     do i=4,nres
 309 !      phi(i)=dihang_in(1,i-2,1,1)
 310 !     enddo
 311       do i=2,nres-1
 312        theta(i+1)=aarray(1,i,1)
 313        phi(i+2)=aarray(2,i,1)
 314        alph(i)=aarray(3,i,1)
 315        omeg(i)=aarray(4,i,1)
 316       enddo
 317
 318       call chainbuild
 319
 320       return
 321       end subroutine dihang_to_c
 322 !-----------------------------------------------------------------------------
 323 ! geomout.F
 324 !-----------------------------------------------------------------------------
 325 #ifdef NOXDR
 326       subroutine cartout(time)
 327 #else
 328       subroutine cartoutx(time)
 329 #endif
 330       use geometry_data, only: c,nres
 331       use energy_data
 332       use MD_data, only: potE,t_bath
 333 !      implicit real*8 (a-h,o-z)
 334 !      include 'DIMENSIONS'
 335 !      include 'COMMON.CHAIN'
 336 !      include 'COMMON.INTERACT'
 337 !      include 'COMMON.NAMES'
 338 !      include 'COMMON.IOUNITS'
 339 !      include 'COMMON.HEADER'
 340 !      include 'COMMON.SBRIDGE'
 341 !      include 'COMMON.DISTFIT'
 342 !      include 'COMMON.MD'
 343       real(kind=8) :: time
 344 !el  local variables
 345       integer :: j,k,i
 346
 347 #if defined(AIX) || defined(PGI)
 348       open(icart,file=cartname,position="append")
 349 #else
 350       open(icart,file=cartname,access="append")
 351 #endif
 352       write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 353       if (dyn_ss) then
 354        write (icart,'(i4,$)') &
 355          nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)
 356       else
 357        write (icart,'(i4,$)') &
 358          nss,(ihpb(j),jhpb(j),j=1,nss)
 359        endif
 360        write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,&
 361        (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),&
 362        (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
 363       write (icart,'(8f10.5)') &
 364        ((c(k,j),k=1,3),j=1,nres),&
 365        ((c(k,j+nres),k=1,3),j=nnt,nct)
 366       close(icart)
 367       return
 368
 369 #ifdef NOXDR
 370       end subroutine cartout
 371 #else
 372       end subroutine cartoutx
 373 #endif
 374 !-----------------------------------------------------------------------------
 375 #ifndef NOXDR
 376       subroutine cartout(time)
 377 !      implicit real*8 (a-h,o-z)
 378 !      include 'DIMENSIONS'
 379       use geometry_data, only: c,nres
 380       use energy_data
 381       use MD_data, only: potE,t_bath
 382 #ifdef MPI
 383       use MPI_data
 384       include 'mpif.h'
 385 !      include 'COMMON.SETUP'
 386 #else
 387       integer,parameter :: me=0
 388 #endif
 389 !      include 'COMMON.CHAIN'
 390 !      include 'COMMON.INTERACT'
 391 !      include 'COMMON.NAMES'
 392 !      include 'COMMON.IOUNITS'
 393 !      include 'COMMON.HEADER'
 394 !      include 'COMMON.SBRIDGE'
 395 !      include 'COMMON.DISTFIT'
 396 !      include 'COMMON.MD'
 397       real(kind=8) :: time
 398       integer :: iret,itmp
 399       real(kind=4) :: prec
 400       real(kind=4),dimension(3,2*nres+2) :: xcoord      !(3,maxres2+2)  (maxres2=2*maxres
 401 !el  local variables
 402       integer :: j,i,ixdrf
 403
 404 #ifdef AIX
 405       call xdrfopen_(ixdrf,cartname, "a", iret)
 406       call xdrffloat_(ixdrf, real(time), iret)
 407       call xdrffloat_(ixdrf, real(potE), iret)
 408       call xdrffloat_(ixdrf, real(uconst), iret)
 409       call xdrffloat_(ixdrf, real(uconst_back), iret)
 410       call xdrffloat_(ixdrf, real(t_bath), iret)
 411       call xdrfint_(ixdrf, nss, iret)
 412       do j=1,nss
 413        if (dyn_ss) then
 414         call xdrfint_(ixdrf, idssb(j)+nres, iret)
 415         call xdrfint_(ixdrf, jdssb(j)+nres, iret)
 416        else
 417         call xdrfint_(ixdrf, ihpb(j), iret)
 418         call xdrfint_(ixdrf, jhpb(j), iret)
 419        endif
 420       enddo
 421       call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
 422       do i=1,nfrag
 423         call xdrffloat_(ixdrf, real(qfrag(i)), iret)
 424       enddo
 425       do i=1,npair
 426         call xdrffloat_(ixdrf, real(qpair(i)), iret)
 427       enddo
 428       do i=1,nfrag_back
 429         call xdrffloat_(ixdrf, real(utheta(i)), iret)
 430         call xdrffloat_(ixdrf, real(ugamma(i)), iret)
 431         call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
 432       enddo
 433 #else
 434       call xdrfopen(ixdrf,cartname, "a", iret)
 435       call xdrffloat(ixdrf, real(time), iret)
 436       call xdrffloat(ixdrf, real(potE), iret)
 437       call xdrffloat(ixdrf, real(uconst), iret)
 438       call xdrffloat(ixdrf, real(uconst_back), iret)
 439       call xdrffloat(ixdrf, real(t_bath), iret)
 440       call xdrfint(ixdrf, nss, iret)
 441       do j=1,nss
 442        if (dyn_ss) then
 443         call xdrfint(ixdrf, idssb(j)+nres, iret)
 444         call xdrfint(ixdrf, jdssb(j)+nres, iret)
 445        else
 446         call xdrfint(ixdrf, ihpb(j), iret)
 447         call xdrfint(ixdrf, jhpb(j), iret)
 448        endif
 449       enddo
 450       call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
 451       do i=1,nfrag
 452         call xdrffloat(ixdrf, real(qfrag(i)), iret)
 453       enddo
 454       do i=1,npair
 455         call xdrffloat(ixdrf, real(qpair(i)), iret)
 456       enddo
 457       do i=1,nfrag_back
 458         call xdrffloat(ixdrf, real(utheta(i)), iret)
 459         call xdrffloat(ixdrf, real(ugamma(i)), iret)
 460         call xdrffloat(ixdrf, real(uscdiff(i)), iret)
 461       enddo
 462 #endif
 463       prec=10000.0
 464       do i=1,nres
 465        do j=1,3
 466         xcoord(j,i)=c(j,i)
 467        enddo
 468       enddo
 469       do i=nnt,nct
 470        do j=1,3
 471         xcoord(j,nres+i-nnt+1)=c(j,i+nres)
 472        enddo
 473       enddo
 474
 475       itmp=nres+nct-nnt+1
 476 #ifdef AIX
 477       call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
 478       call xdrfclose_(ixdrf, iret)
 479 #else
 480       call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
 481       call xdrfclose(ixdrf, iret)
 482 #endif
 483       return
 484       end subroutine cartout
 485 #endif
 486 !-----------------------------------------------------------------------------
 487       subroutine statout(itime)
 488
 489       use energy_data
 490       use control_data
 491       use MD_data
 492       use MPI_data
 493       use compare, only:rms_nac_nnc
 494 !      implicit real*8 (a-h,o-z)
 495 !      include 'DIMENSIONS'
 496 !      include 'COMMON.CONTROL'
 497 !      include 'COMMON.CHAIN'
 498 !      include 'COMMON.INTERACT'
 499 !      include 'COMMON.NAMES'
 500 !      include 'COMMON.IOUNITS'
 501 !      include 'COMMON.HEADER'
 502 !      include 'COMMON.SBRIDGE'
 503 !      include 'COMMON.DISTFIT'
 504 !      include 'COMMON.MD'
 505 !      include 'COMMON.REMD'
 506 !      include 'COMMON.SETUP'
 507       integer :: itime
 508       real(kind=8),dimension(0:n_ene) :: energia
 509 !      double precision gyrate
 510 !el      external gyrate
 511 !el      common /gucio/ cm
 512       character(len=256) :: line1,line2
 513       character(len=4) :: format1,format2
 514       character(len=30) :: format
 515 !el  local variables
 516       integer :: i,ii1,ii2,j
 517       real(kind=8) :: rms,frac,frac_nn,co,distance
 518
 519 #ifdef AIX
 520       if(itime.eq.0) then
 521        open(istat,file=statname,position="append")
 522       endif
 523 #else
 524 #ifdef PGI
 525       open(istat,file=statname,position="append")
 526 #else
 527       open(istat,file=statname,access="append")
 528 #endif
 529 #endif
 530        if (AFMlog.gt.0) then
 531        if (refstr) then
 532          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
 533           write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')&
 534                itime,totT,EK,potE+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,i3,a,i3)')'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       call alloc_ener_arrays
1243 !--------------------------------
1244 ! 8/13/98 Set limits to generating the dihedral angles
1245       do i=1,nres
1246         phibound(1,i)=-pi
1247         phibound(2,i)=pi
1248       enddo
1249       read (inp,*) ndih_constr
1250       if (ndih_constr.gt.0) then
1251         raw_psipred=.false.
1252         allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1253         allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1254         allocate(ftors(ndih_constr)) !(maxdih_constr)
1255
1256 !        read (inp,*) ftors
1257         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), &
1258         i=1,ndih_constr)
1259         if(me.eq.king.or..not.out1file)then
1260          write (iout,*) &
1261          'There are',ndih_constr,' constraints on phi angles.'
1262          do i=1,ndih_constr
1263           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), &
1264           ftors(i)
1265          enddo
1266         endif
1267         do i=1,ndih_constr
1268           phi0(i)=deg2rad*phi0(i)
1269           drange(i)=deg2rad*drange(i)
1270         enddo
1271 !        if(me.eq.king.or..not.out1file) &
1272 !         write (iout,*) 'FTORS',ftors
1273         do i=1,ndih_constr
1274           ii = idih_constr(i)
1275           phibound(1,ii) = phi0(i)-drange(i)
1276           phibound(2,ii) = phi0(i)+drange(i)
1277         enddo
1278       else if (ndih_constr.lt.0) then
1279         raw_psipred=.true.
1280         allocate(idih_constr(nres))
1281         allocate(secprob(3,nres))
1282         allocate(vpsipred(3,nres))
1283         allocate(sdihed(2,nres))
1284         call card_concat(weightcard,.true.)
1285         call reada(weightcard,"PHIHEL",phihel,50.0D0)
1286         call reada(weightcard,"PHIBET",phibet,180.0D0)
1287         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
1288         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
1289         call reada(weightcard,"WDIHC",wdihc,0.591D0)
1290         write (iout,*) "Weight of dihedral angle restraints",wdihc
1291         read(inp,'(9x,3f7.3)') &
1292           (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
1293         write (iout,*) "The secprob array"
1294         do i=nnt,nct
1295           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
1296         enddo
1297         ndih_constr=0
1298         do i=nnt+3,nct
1299           if (itype(i-3,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1 &
1300           .and. itype(i-1,1).ne.ntyp1 .and. itype(i,1).ne.ntyp1) then
1301             ndih_constr=ndih_constr+1
1302             idih_constr(ndih_constr)=i
1303             sumv=0.0d0
1304             do j=1,3
1305               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
1306               sumv=sumv+vpsipred(j,ndih_constr)
1307             enddo
1308             do j=1,3
1309               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
1310             enddo
1311             phibound(1,ndih_constr)=phihel*deg2rad
1312             phibound(2,ndih_constr)=phibet*deg2rad
1313             sdihed(1,ndih_constr)=sigmahel*deg2rad
1314             sdihed(2,ndih_constr)=sigmabet*deg2rad
1315           endif
1316         enddo
1317
1318       endif
1319       if (with_theta_constr) then
1320 !C with_theta_constr is keyword allowing for occurance of theta constrains
1321       read (inp,*) ntheta_constr
1322 !C ntheta_constr is the number of theta constrains
1323       if (ntheta_constr.gt.0) then
1324 !C        read (inp,*) ftors
1325         allocate(itheta_constr(ntheta_constr))
1326         allocate(theta_constr0(ntheta_constr))
1327         allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr))
1328         read (inp,*) (itheta_constr(i),theta_constr0(i), &
1329        theta_drange(i),for_thet_constr(i), &
1330        i=1,ntheta_constr)
1331 !C the above code reads from 1 to ntheta_constr
1332 !C itheta_constr(i) residue i for which is theta_constr
1333 !C theta_constr0 the global minimum value
1334 !C theta_drange is range for which there is no energy penalty
1335 !C for_thet_constr is the force constant for quartic energy penalty
1336 !C E=k*x**4
1337         if(me.eq.king.or..not.out1file)then
1338          write (iout,*) &
1339         'There are',ntheta_constr,' constraints on phi angles.'
1340          do i=1,ntheta_constr
1341           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), &
1342          theta_drange(i), &
1343          for_thet_constr(i)
1344          enddo
1345         endif
1346         do i=1,ntheta_constr
1347           theta_constr0(i)=deg2rad*theta_constr0(i)
1348           theta_drange(i)=deg2rad*theta_drange(i)
1349         enddo
1350 !C        if(me.eq.king.or..not.out1file)
1351 !C     &   write (iout,*) 'FTORS',ftors
1352 !C        do i=1,ntheta_constr
1353 !C          ii = itheta_constr(i)
1354 !C          thetabound(1,ii) = phi0(i)-drange(i)
1355 !C          thetabound(2,ii) = phi0(i)+drange(i)
1356 !C        enddo
1357       endif ! ntheta_constr.gt.0
1358       endif! with_theta_constr
1359
1360       nnt=1
1361 #ifdef MPI
1362       if (me.eq.king) then
1363 #endif
1364        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1365        do i=1,nres
1366          write (iout,'(a3,i5,2f10.1)') &
1367          restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1368        enddo
1369 #ifdef MP
1370       endif
1371 #endif
1372       nct=nres
1373       allocate(ireschain(nres))
1374       ireschain=0
1375       write(iout,*),"before seq2chains",ireschain
1376       call seq2chains
1377       write(iout,*) "after seq2chains",nchain
1378       allocate ( chain_border1(2,nchain))
1379       chain_border1(1,1)=1
1380       chain_border1(2,1)=chain_border(2,1)+1
1381       do i=2,nchain-1
1382         chain_border1(1,i)=chain_border(1,i)-1
1383         chain_border1(2,i)=chain_border(2,i)+1
1384       enddo
1385       if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
1386       chain_border1(2,nchain)=nres
1387       write(iout,*) "nres",nres," nchain",nchain
1388       do i=1,nchain
1389         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),&
1390          chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
1391       enddo
1392       allocate(tabpermchain(50,5040))
1393       call chain_symmetry(npermchain,tabpermchain)
1394       print *,'NNT=',NNT,' NCT=',NCT
1395       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1396       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1397       if (pdbref) then
1398         if(me.eq.king.or..not.out1file) &
1399          write (iout,'(a,i3)') 'nsup=',nsup
1400         nstart_seq=nnt
1401         if (nsup.le.(nct-nnt+1)) then
1402           do i=0,nct-nnt+1-nsup
1403             if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
1404               nstart_seq=nnt+i
1405               goto 111
1406             endif
1407           enddo
1408           write (iout,'(a)') &
1409                   'Error - sequences to be superposed do not match.'
1410           stop
1411         else
1412           do i=0,nsup-(nct-nnt+1)
1413             if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1414             then
1415               nstart_sup=nstart_sup+i
1416               nsup=nct-nnt+1
1417               goto 111
1418             endif
1419           enddo
1420           write (iout,'(a)') &
1421                   'Error - sequences to be superposed do not match.'
1422         endif
1423   111   continue
1424         if (nsup.eq.0) nsup=nct-nnt+1
1425         if (nstart_sup.eq.0) nstart_sup=nnt
1426         if (nstart_seq.eq.0) nstart_seq=nnt
1427         if(me.eq.king.or..not.out1file) &
1428          write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1429                        ' nstart_seq=',nstart_seq !,"242343453254"
1430       endif
1431 !--- Zscore rms -------
1432       if (nz_start.eq.0) nz_start=nnt
1433       if (nz_end.eq.0 .and. nsup.gt.0) then
1434         nz_end=nnt+nsup-1
1435       else if (nz_end.eq.0) then
1436         nz_end=nct
1437       endif
1438       if(me.eq.king.or..not.out1file)then
1439        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1440        write (iout,*) 'IZ_SC=',iz_sc, 'NCT=',nct
1441       endif
1442 !----------------------
1443       call init_int_table
1444       if (refstr) then
1445         if (.not.pdbref) then
1446           call read_angles(inp,*38)
1447           goto 39
1448    38     write (iout,'(a)') 'Error reading reference structure.'
1449 #ifdef MPI
1450           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1451           stop 'Error reading reference structure'
1452 #endif
1453    39     call chainbuild
1454           call setup_var
1455 !zscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1456           nstart_sup=nnt
1457           nstart_seq=nnt
1458           nsup=nct-nnt+1
1459           kkk=1
1460           do i=1,2*nres
1461             do j=1,3
1462               cref(j,i,kkk)=c(j,i)
1463             enddo
1464           enddo
1465           call contact(.true.,ncont_ref,icont_ref,co)
1466         endif
1467 !        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1468 !        call flush(iout)
1469 !EL        if (constr_dist.gt.0) call read_dist_constr
1470 !EL        write (iout,*) "After read_dist_constr nhpb",nhpb
1471 !EL        if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1472 !EL        call hpb_partition
1473         if(me.eq.king.or..not.out1file) &
1474          write (iout,*) 'Contact order:',co
1475         if (pdbref) then
1476         if(me.eq.king.or..not.out1file) &
1477          write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1478         do i=1,ncont_ref
1479           do j=1,2
1480             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1481           enddo
1482           if(me.eq.king.or..not.out1file) &
1483            write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1484            icont_ref(1,i),' ',&
1485            restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1486         enddo
1487         endif
1488       if (constr_homology.gt.0) then
1489 !        write (iout,*) "Calling read_constr_homology"
1490 !        call flush(iout)
1491         call read_constr_homology
1492         if (indpdb.gt.0 .or. pdbref) then
1493           do i=1,2*nres
1494             do j=1,3
1495               c(j,i)=crefjlee(j,i)
1496               cref(j,i,1)=crefjlee(j,i)
1497             enddo
1498           enddo
1499         endif
1500 #define DEBUG
1501 #ifdef DEBUG
1502         write (iout,*) "sc_loc_geom: Array C"
1503         do i=1,nres
1504           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),&
1505            (c(j,i+nres),j=1,3)
1506         enddo
1507         write (iout,*) "Array Cref"
1508         do i=1,nres
1509           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),&
1510            (cref(j,i+nres,1),j=1,3)
1511         enddo
1512 #endif
1513        call int_from_cart1(.false.)
1514        call sc_loc_geom(.false.)
1515        do i=1,nres
1516          thetaref(i)=theta(i)
1517          phiref(i)=phi(i)
1518        enddo
1519        do i=1,nres-1
1520          do j=1,3
1521            dc(j,i)=c(j,i+1)-c(j,i)
1522            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1523          enddo
1524        enddo
1525        do i=2,nres-1
1526          do j=1,3
1527            dc(j,i+nres)=c(j,i+nres)-c(j,i)
1528            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1529          enddo
1530        enddo
1531       else
1532         homol_nset=0
1533         if (start_from_model) then
1534           nmodel_start=0
1535           do
1536             read(inp,'(a)',end=332,err=332) pdbfile
1537             if (me.eq.king .or. .not. out1file)&
1538              write (iout,'(a,5x,a)') 'Opening PDB file',&
1539              pdbfile(:ilen(pdbfile))
1540             open(ipdbin,file=pdbfile,status='old',err=336)
1541             goto 335
1542  336        write (iout,'(a,5x,a)') 'Error opening PDB file',&
1543            pdbfile(:ilen(pdbfile))
1544             call flush(iout)
1545             stop
1546  335        continue
1547             unres_pdb=.false.
1548             nres_temp=nres
1549 !            call readpdb
1550             call readpdb_template(nmodel_start+1)
1551             close(ipdbin)
1552             if (nres.ge.nres_temp) then
1553               nmodel_start=nmodel_start+1
1554               pdbfiles_chomo(nmodel_start)=pdbfile
1555               do i=1,2*nres
1556                 do j=1,3
1557                   chomo(j,i,nmodel_start)=c(j,i)
1558                 enddo
1559               enddo
1560             else
1561               if (me.eq.king .or. .not. out1file) &
1562                write (iout,'(a,2i7,1x,a)') &
1563                 "Different number of residues",nres_temp,nres, &
1564                 " model skipped."
1565             endif
1566             nres=nres_temp
1567           enddo
1568   332     continue
1569           if (nmodel_start.eq.0) then
1570             if (me.eq.king .or. .not. out1file) &
1571              write (iout,'(a)') &
1572              "No valid starting model found START_FROM_MODELS is OFF"
1573               start_from_model=.false.
1574           endif
1575           write (iout,*) "nmodel_start",nmodel_start
1576         endif
1577       endif
1578
1579       endif
1580         if (constr_dist.gt.0) call read_dist_constr
1581         write (iout,*) "After read_dist_constr nhpb",nhpb
1582         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1583         if (velnanoconst.ne.0) call read_afmnano
1584         call hpb_partition
1585
1586       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1587           .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1588           modecalc.ne.10) then
1589 ! If input structure hasn't been supplied from the PDB file read or generate
1590 ! initial geometry.
1591         if (iranconf.eq.0 .and. .not. extconf) then
1592           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1593            write (iout,'(a)') 'Initial geometry will be read in.'
1594           if (read_cart) then
1595             read(inp,'(8f10.5)',end=36,err=36) &
1596              ((c(l,k),l=1,3),k=1,nres),&
1597              ((c(l,k+nres),l=1,3),k=nnt,nct)
1598             write (iout,*) "Exit READ_CART"
1599             write (iout,'(8f10.5)') &
1600              ((c(l,k),l=1,3),k=1,nres)
1601             write (iout,'(8f10.5)') &
1602              ((c(l,k+nres),l=1,3),k=nnt,nct)
1603             call int_from_cart1(.true.)
1604             write (iout,*) "Finish INT_TO_CART"
1605             do i=1,nres-1
1606               do j=1,3
1607                 dc(j,i)=c(j,i+1)-c(j,i)
1608                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1609               enddo
1610             enddo
1611             do i=nnt,nct
1612               if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1613                 do j=1,3
1614                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1615                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1616                 enddo
1617               endif
1618             enddo
1619             return
1620           else
1621            write(iout,*) "read angles from input"
1622            call read_angles(inp,*36)
1623             call chainbuild
1624
1625           endif
1626           goto 37
1627    36     write (iout,'(a)') 'Error reading angle file.'
1628 #ifdef MPI
1629           call mpi_finalize( MPI_COMM_WORLD,IERR )
1630 #endif
1631           stop 'Error reading angle file.'
1632    37     continue
1633         else if (extconf) then
1634          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1635           write (iout,'(a)') 'Extended chain initial geometry.'
1636          do i=3,nres
1637           theta(i)=90d0*deg2rad
1638           if (molnum(i).eq.2) theta(i)=160d0*deg2rad
1639          enddo
1640          do i=4,nres
1641           phi(i)=180d0*deg2rad
1642           if (molnum(i).eq.2) phi(i)=30d0*deg2rad
1643          enddo
1644          do i=2,nres-1
1645           alph(i)=110d0*deg2rad
1646           if (molnum(i).eq.2) alph(i)=30d0*deg2rad
1647          enddo
1648          do i=2,nres-1
1649           omeg(i)=-120d0*deg2rad
1650           if (molnum(i).eq.2) omeg(i)=60d0*deg2rad
1651           if (itype(i,1).le.0) omeg(i)=-omeg(i)
1652          enddo
1653          call chainbuild
1654         else
1655           if(me.eq.king.or..not.out1file) &
1656            write (iout,'(a)') 'Random-generated initial geometry.'
1657
1658
1659 #ifdef MPI
1660           if (me.eq.king  .or. fg_rank.eq.0 .and. &
1661                  ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1662 #endif
1663             do itrial=1,100
1664               itmp=1
1665               call gen_rand_conf(itmp,*30)
1666               goto 40
1667    30         write (iout,*) 'Failed to generate random conformation',&
1668                 ', itrial=',itrial
1669               write (*,*) 'Processor:',me,&
1670                 ' Failed to generate random conformation',&
1671                 ' itrial=',itrial
1672               call intout
1673
1674 #ifdef AIX
1675               call flush_(iout)
1676 #else
1677               call flush(iout)
1678 #endif
1679             enddo
1680             write (iout,'(a,i3,a)') 'Processor:',me,&
1681               ' error in generating random conformation.'
1682             write (*,'(a,i3,a)') 'Processor:',me,&
1683               ' error in generating random conformation.'
1684             call flush(iout)
1685 #ifdef MPI
1686             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1687    40       continue
1688           endif
1689 #else
1690           do itrial=1,100
1691             itmp=1
1692             call gen_rand_conf(itmp,*335)
1693             goto 40
1694   335       write (iout,*) 'Failed to generate random conformation',&
1695               ', itrial=',itrial
1696             write (*,*) 'Failed to generate random conformation',&
1697               ', itrial=',itrial
1698           enddo
1699           write (iout,'(a,i3,a)') 'Processor:',me,&
1700             ' error in generating random conformation.'
1701           write (*,'(a,i3,a)') 'Processor:',me,&
1702             ' error in generating random conformation.'
1703           stop
1704    40     continue
1705 #endif
1706         endif
1707       elseif (modecalc.eq.4) then
1708         read (inp,'(a)') intinname
1709         open (intin,file=intinname,status='old',err=333)
1710         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1711         write (iout,'(a)') 'intinname',intinname
1712         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1713         goto 334
1714   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1715 #ifdef MPI
1716         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1717 #endif
1718         stop 'Error opening angle file.'
1719   334   continue
1720
1721       endif
1722 ! Generate distance constraints, if the PDB structure is to be regularized.
1723       if (nthread.gt.0) then
1724         call read_threadbase
1725       endif
1726       call setup_var
1727       if (me.eq.king .or. .not. out1file) &
1728        call intout
1729       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1730         write (iout,'(/a,i3,a)') &
1731         'The chain contains',ns,' disulfide-bridging cysteines.'
1732         write (iout,'(20i4)') (iss(i),i=1,ns)
1733        if (dyn_ss) then
1734           write(iout,*)"Running with dynamic disulfide-bond formation"
1735        else
1736         write (iout,'(/a/)') 'Pre-formed links are:'
1737         do i=1,nss
1738           i1=ihpb(i)-nres
1739           i2=jhpb(i)-nres
1740           it1=itype(i1,1)
1741           it2=itype(i2,1)
1742           if (me.eq.king.or..not.out1file) &
1743           write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1744           restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1745           ebr,forcon(i)
1746         enddo
1747         write (iout,'(a)')
1748        endif
1749       endif
1750       if (ns.gt.0.and.dyn_ss) then
1751           do i=nss+1,nhpb
1752             ihpb(i-nss)=ihpb(i)
1753             jhpb(i-nss)=jhpb(i)
1754             forcon(i-nss)=forcon(i)
1755             dhpb(i-nss)=dhpb(i)
1756           enddo
1757           nhpb=nhpb-nss
1758           nss=0
1759           call hpb_partition
1760           do i=1,ns
1761             dyn_ss_mask(iss(i))=.true.
1762           enddo
1763       endif
1764       if (i2ndstr.gt.0) call secstrp2dihc
1765       if (indpdb.gt.0) then
1766           write(iout,*) "WCHODZE TU!!"
1767           call int_from_cart1(.true.)
1768       endif
1769 !      call geom_to_var(nvar,x)
1770 !      call etotal(energia(0))
1771 !      call enerprint(energia(0))
1772 !      call briefout(0,etot)
1773 !      stop
1774 !d    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1775 !d    write (iout,'(a)') 'Variable list:'
1776 !d    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1777 #ifdef MPI
1778       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1779         write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1780         'Processor',myrank,': end reading molecular data.'
1781 #endif
1782       return
1783       end subroutine molread
1784 !-----------------------------------------------------------------------------
1785       subroutine read_constr_homology
1786       use control, only:init_int_table,homology_partition
1787       use MD_data, only:iset
1788 !      implicit none
1789 !      include 'DIMENSIONS'
1790 !#ifdef MPI
1791 !      include 'mpif.h'
1792 !#endif
1793 !      include 'COMMON.SETUP'
1794 !      include 'COMMON.CONTROL'
1795 !      include 'COMMON.HOMOLOGY'
1796 !      include 'COMMON.CHAIN'
1797 !      include 'COMMON.IOUNITS'
1798 !      include 'COMMON.MD'
1799 !      include 'COMMON.QRESTR'
1800 !      include 'COMMON.GEO'
1801 !      include 'COMMON.INTERACT'
1802 !      include 'COMMON.NAMES'
1803 !      include 'COMMON.VAR'
1804 !
1805
1806 !     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1807 !    &                 dist_cut
1808 !     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1809 !    &    sigma_odl_temp(maxres,maxres,max_template)
1810       character*2 kic2
1811       character*24 model_ki_dist, model_ki_angle
1812       character*500 controlcard
1813       integer :: ki,i,ii,j,k,l
1814       integer, dimension (:), allocatable :: ii_in_use
1815       integer :: i_tmp,idomain_tmp,&
1816       irec,ik,iistart,nres_temp
1817 !      integer :: iset
1818 !      external :: ilen
1819       logical :: liiflag,lfirst
1820       integer :: i01,i10
1821 !
1822 !     FP - Nov. 2014 Temporary specifications for new vars
1823 !
1824       real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
1825                        rescore3_tmp, dist_cut
1826       real(kind=8), dimension (:,:),allocatable :: rescore
1827       real(kind=8), dimension (:,:),allocatable :: rescore2
1828       real(kind=8), dimension (:,:),allocatable :: rescore3
1829       real(kind=8) :: distal
1830       character*24 tpl_k_rescore
1831       character*256 pdbfile
1832
1833 ! -----------------------------------------------------------------
1834 ! Reading multiple PDB ref structures and calculation of retraints
1835 ! not using pre-computed ones stored in files model_ki_{dist,angle}
1836 ! FP (Nov., 2014)
1837 ! -----------------------------------------------------------------
1838 !
1839 !
1840 ! Alternative: reading from input
1841       call card_concat(controlcard,.true.)
1842       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1843       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1844       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1845       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1846       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1847       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1848       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1849       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1850       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
1851       if(.not.read2sigma.and.start_from_model) then
1852           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
1853            write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
1854           start_from_model=.false.
1855           iranconf=(indpdb.le.0)
1856       endif
1857       if(start_from_model .and. (me.eq.king .or. .not. out1file))&
1858          write(iout,*) 'START_FROM_MODELS is ON'
1859 !      if(start_from_model .and. rest) then
1860 !        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1861 !         write(iout,*) 'START_FROM_MODELS is OFF'
1862 !         write(iout,*) 'remove restart keyword from input'
1863 !        endif
1864 !      endif
1865       if (start_from_model) nmodel_start=constr_homology
1866       if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
1867       if (homol_nset.gt.1)then
1868          call card_concat(controlcard,.true.)
1869          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1870          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1871 !          write(iout,*) "iset homology_weight "
1872           do i=1,homol_nset
1873            write(iout,*) i,waga_homology(i)
1874           enddo
1875          endif
1876          iset=mod(kolor,homol_nset)+1
1877       else
1878        iset=1
1879        waga_homology(1)=1.0
1880       endif
1881
1882 !d      write (iout,*) "nnt",nnt," nct",nct
1883 !d      call flush(iout)
1884
1885       if (read_homol_frag) then
1886         call read_klapaucjusz
1887       else
1888
1889       lim_odl=0
1890       lim_dih=0
1891 !
1892 !      write(iout,*) 'nnt=',nnt,'nct=',nct
1893 !
1894 !      do i = nnt,nct
1895 !        do k=1,constr_homology
1896 !         idomain(k,i)=0
1897 !        enddo
1898 !      enddo
1899        idomain=0
1900
1901 !      ii=0
1902 !      do i = nnt,nct-2
1903 !        do j=i+2,nct
1904 !        ii=ii+1
1905 !        ii_in_use(ii)=0
1906 !        enddo
1907 !      enddo
1908       ii_in_use=0
1909       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
1910       if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
1911
1912       do k=1,constr_homology
1913
1914         read(inp,'(a)') pdbfile
1915         pdbfiles_chomo(k)=pdbfile
1916         if(me.eq.king .or. .not. out1file) &
1917          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
1918         pdbfile(:ilen(pdbfile))
1919         open(ipdbin,file=pdbfile,status='old',err=33)
1920         goto 34
1921   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
1922         pdbfile(:ilen(pdbfile))
1923         stop
1924   34    continue
1925 !        print *,'Begin reading pdb data'
1926 !
1927 ! Files containing res sim or local scores (former containing sigmas)
1928 !
1929
1930         write(kic2,'(bz,i2.2)') k
1931
1932         tpl_k_rescore="template"//kic2//".sco"
1933         write(iout,*) "tpl_k_rescore=",tpl_k_rescore
1934         unres_pdb=.false.
1935         nres_temp=nres
1936         write(iout,*) "read2sigma",read2sigma
1937
1938         if (read2sigma) then
1939           call readpdb_template(k)
1940         else
1941           call readpdb
1942         endif
1943         write(iout,*) "after readpdb"
1944         if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
1945         nres_chomo(k)=nres
1946         nres=nres_temp
1947         if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
1948         if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
1949         if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
1950         if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
1951         if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
1952         if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
1953         if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
1954         if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
1955         if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
1956         if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
1957         if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
1958         if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
1959         if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1960         if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
1961 !        if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1962         if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
1963         if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
1964         if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
1965         if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
1966 !        if(.not.allocated(distance)) allocate(distance(constr_homology))
1967 !        if(.not.allocated(distancek)) allocate(distancek(constr_homology))
1968
1969
1970 !
1971 !     Distance restraints
1972 !
1973 !          ... --> odl(k,ii)
1974 ! Copy the coordinates from reference coordinates (?)
1975         do i=1,2*nres_chomo(k)
1976           do j=1,3
1977             c(j,i)=cref(j,i,1)
1978 !           write (iout,*) "c(",j,i,") =",c(j,i)
1979           enddo
1980         enddo
1981 !
1982 ! From read_dist_constr (commented out 25/11/2014 <-> res sim)
1983 !
1984 !         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1985           open (ientin,file=tpl_k_rescore,status='old')
1986           if (nnt.gt.1) rescore(k,1)=0.0d0
1987           do irec=nnt,nct ! loop for reading res sim
1988             if (read2sigma) then
1989              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
1990                                      rescore3_tmp,idomain_tmp
1991              i_tmp=i_tmp+nnt-1
1992              idomain(k,i_tmp)=idomain_tmp
1993              rescore(k,i_tmp)=rescore_tmp
1994              rescore2(k,i_tmp)=rescore2_tmp
1995              rescore3(k,i_tmp)=rescore3_tmp
1996              if (.not. out1file .or. me.eq.king)&
1997              write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
1998                            i_tmp,rescore2_tmp,rescore_tmp,&
1999                                      rescore3_tmp,idomain_tmp
2000             else
2001              idomain(k,irec)=1
2002              read (ientin,*,end=1401) rescore_tmp
2003
2004 !           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2005              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2006 !           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2007             endif
2008           enddo
2009  1401   continue
2010         close (ientin)
2011         if (waga_dist.ne.0.0d0) then
2012           ii=0
2013           do i = nnt,nct-2
2014             do j=i+2,nct
2015
2016               x12=c(1,i)-c(1,j)
2017               y12=c(2,i)-c(2,j)
2018               z12=c(3,i)-c(3,j)
2019               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2020 !              write (iout,*) k,i,j,distal,dist2_cut
2021
2022             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2023                  .and. distal.le.dist2_cut ) then
2024
2025               ii=ii+1
2026               ii_in_use(ii)=1
2027               l_homo(k,ii)=.true.
2028
2029 !             write (iout,*) "k",k
2030 !             write (iout,*) "i",i," j",j," constr_homology",
2031 !    &                       constr_homology
2032               ires_homo(ii)=i
2033               jres_homo(ii)=j
2034               odl(k,ii)=distal
2035               if (read2sigma) then
2036                 sigma_odl(k,ii)=0
2037                 do ik=i,j
2038                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2039                 enddo
2040                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2041                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2042               sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2043               else
2044                 if (odl(k,ii).le.dist_cut) then
2045                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2046                 else
2047 #ifdef OLDSIGMA
2048                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2049                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2050 #else
2051                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2052                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2053 #endif
2054                 endif
2055               endif
2056               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2057             else
2058 !              ii=ii+1
2059 !              l_homo(k,ii)=.false.
2060             endif
2061             enddo
2062           enddo
2063         lim_odl=ii
2064         endif
2065 !        write (iout,*) "Distance restraints set"
2066 !        call flush(iout)
2067 !
2068 !     Theta, dihedral and SC retraints
2069 !
2070         if (waga_angle.gt.0.0d0) then
2071 !         open (ientin,file=tpl_k_sigma_dih,status='old')
2072 !         do irec=1,maxres-3 ! loop for reading sigma_dih
2073 !            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2074 !            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2075 !            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2076 !    &                            sigma_dih(k,i+nnt-1)
2077 !         enddo
2078 !1402   continue
2079 !         close (ientin)
2080           do i = nnt+3,nct
2081             if (idomain(k,i).eq.0) then
2082                sigma_dih(k,i)=0.0
2083                cycle
2084             endif
2085             dih(k,i)=phiref(i) ! right?
2086 !           read (ientin,*) sigma_dih(k,i) ! original variant
2087 !             write (iout,*) "dih(",k,i,") =",dih(k,i)
2088 !             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2089 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2090 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
2091 !    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
2092
2093             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2094                           rescore(k,i-2)+rescore(k,i-3))/4.0
2095 !            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2096 !           write (iout,*) "Raw sigmas for dihedral angle restraints"
2097 !           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2098 !           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2099 !                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
2100 !   Instead of res sim other local measure of b/b str reliability possible
2101             if (sigma_dih(k,i).ne.0) &
2102             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2103 !           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2104           enddo
2105           lim_dih=nct-nnt-2
2106         endif
2107 !        write (iout,*) "Dihedral angle restraints set"
2108 !        call flush(iout)
2109
2110         if (waga_theta.gt.0.0d0) then
2111 !         open (ientin,file=tpl_k_sigma_theta,status='old')
2112 !         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2113 !            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2114 !            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2115 !    &                              sigma_theta(k,i+nnt-1)
2116 !         enddo
2117 !1403   continue
2118 !         close (ientin)
2119
2120           do i = nnt+2,nct ! right? without parallel.
2121 !         do i = i=1,nres ! alternative for bounds acc to readpdb?
2122 !         do i=ithet_start,ithet_end ! with FG parallel.
2123              if (idomain(k,i).eq.0) then
2124               sigma_theta(k,i)=0.0
2125               cycle
2126              endif
2127              thetatpl(k,i)=thetaref(i)
2128 !            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2129 !            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
2130 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2131 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
2132 !            read (ientin,*) sigma_theta(k,i) ! 1st variant
2133              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2134                              rescore(k,i-2))/3.0
2135 !             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2136              if (sigma_theta(k,i).ne.0) &
2137              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2138
2139 !            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2140 !                             rescore(k,i-2) !  right expression ?
2141 !            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2142           enddo
2143         endif
2144 !        write (iout,*) "Angle restraints set"
2145 !        call flush(iout)
2146
2147         if (waga_d.gt.0.0d0) then
2148 !       open (ientin,file=tpl_k_sigma_d,status='old')
2149 !         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2150 !            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2151 !            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2152 !    &                          sigma_d(k,i+nnt-1)
2153 !         enddo
2154 !1404   continue
2155
2156           do i = nnt,nct ! right? without parallel.
2157 !         do i=2,nres-1 ! alternative for bounds acc to readpdb?
2158 !         do i=loc_start,loc_end ! with FG parallel.
2159                if (itype(i,1).eq.10) cycle
2160                if (idomain(k,i).eq.0 ) then
2161                   sigma_d(k,i)=0.0
2162                   cycle
2163                endif
2164                xxtpl(k,i)=xxref(i)
2165                yytpl(k,i)=yyref(i)
2166                zztpl(k,i)=zzref(i)
2167 !              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2168 !              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2169 !              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2170 !              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
2171                sigma_d(k,i)=rescore3(k,i) !  right expression ?
2172                if (sigma_d(k,i).ne.0) &
2173                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2174
2175 !              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
2176 !              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2177 !              read (ientin,*) sigma_d(k,i) ! 1st variant
2178           enddo
2179         endif
2180       enddo
2181 !      write (iout,*) "SC restraints set"
2182 !      call flush(iout)
2183 !
2184 ! remove distance restraints not used in any model from the list
2185 ! shift data in all arrays
2186 !
2187 !      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
2188       if (waga_dist.ne.0.0d0) then
2189         ii=0
2190         liiflag=.true.
2191         lfirst=.true.
2192         do i=nnt,nct-2
2193          do j=i+2,nct
2194           ii=ii+1
2195 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2196 !     &            .and. distal.le.dist2_cut ) then
2197 !          write (iout,*) "i",i," j",j," ii",ii
2198 !          call flush(iout)
2199           if (ii_in_use(ii).eq.0.and.liiflag.or. &
2200           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2201             liiflag=.false.
2202             i10=ii
2203             if (lfirst) then
2204               lfirst=.false.
2205               iistart=ii
2206             else
2207               if(i10.eq.lim_odl) i10=i10+1
2208               do ki=0,i10-i01-1
2209                ires_homo(iistart+ki)=ires_homo(ki+i01)
2210                jres_homo(iistart+ki)=jres_homo(ki+i01)
2211                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2212                do k=1,constr_homology
2213                 odl(k,iistart+ki)=odl(k,ki+i01)
2214                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2215                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2216                enddo
2217               enddo
2218               iistart=iistart+i10-i01
2219             endif
2220           endif
2221           if (ii_in_use(ii).ne.0.and..not.liiflag) then
2222              i01=ii
2223              liiflag=.true.
2224           endif
2225          enddo
2226         enddo
2227         lim_odl=iistart-1
2228       endif
2229 !      write (iout,*) "Removing distances completed"
2230 !      call flush(iout)
2231       endif ! .not. klapaucjusz
2232
2233       if (constr_homology.gt.0) call homology_partition
2234       write (iout,*) "After homology_partition"
2235       call flush(iout)
2236       if (constr_homology.gt.0) call init_int_table
2237       write (iout,*) "After init_int_table"
2238       call flush(iout)
2239 !      endif ! .not. klapaucjusz
2240 !      endif
2241 !      if (constr_homology.gt.0) call homology_partition
2242 !      write (iout,*) "After homology_partition"
2243 !      call flush(iout)
2244 !      if (constr_homology.gt.0) call init_int_table
2245 !      write (iout,*) "After init_int_table"
2246 !      call flush(iout)
2247 !      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2248 !      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2249 !
2250 ! Print restraints
2251 !
2252       if (.not.out_template_restr) return
2253 !d      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2254       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2255        write (iout,*) "Distance restraints from templates"
2256        do ii=1,lim_odl
2257        write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
2258         ii,ires_homo(ii),jres_homo(ii),&
2259         (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
2260         ki=1,constr_homology)
2261        enddo
2262        write (iout,*) "Dihedral angle restraints from templates"
2263        do i=nnt+3,nct
2264         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2265             (rad2deg*dih(ki,i),&
2266             rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2267        enddo
2268        write (iout,*) "Virtual-bond angle restraints from templates"
2269        do i=nnt+2,nct
2270         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2271             (rad2deg*thetatpl(ki,i),&
2272             rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2273        enddo
2274        write (iout,*) "SC restraints from templates"
2275        do i=nnt,nct
2276         write(iout,'(i7,100(4f8.2,4x))') i,&
2277         (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
2278          1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2279        enddo
2280       endif
2281       return
2282       end subroutine read_constr_homology
2283 !-----------------------------------------------------------------------------
2284       subroutine read_klapaucjusz
2285       use energy_data
2286       implicit none
2287 !     include 'DIMENSIONS'
2288 !#ifdef MPI
2289 !     include 'mpif.h'
2290 !#endif
2291 !     include 'COMMON.SETUP'
2292 !     include 'COMMON.CONTROL'
2293 !     include 'COMMON.HOMOLOGY'
2294 !     include 'COMMON.CHAIN'
2295 !     include 'COMMON.IOUNITS'
2296 !     include 'COMMON.MD'
2297 !     include 'COMMON.GEO'
2298 !     include 'COMMON.INTERACT'
2299 !     include 'COMMON.NAMES'
2300       character(len=256) fragfile
2301       integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,&
2302                          ii_in_use
2303       integer, dimension(:,:), allocatable :: iresclust,inclust
2304       integer :: nclust
2305
2306       character(len=2) :: kic2
2307       character(len=24) :: model_ki_dist, model_ki_angle
2308       character(len=500) :: controlcard
2309       integer :: ki, i, j, jj,k, l, i_tmp,&
2310       idomain_tmp,&
2311       ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,&
2312       i01,i10,nnt_chain,nct_chain
2313       real(kind=8) :: distal
2314       logical :: lprn = .true.
2315       integer :: nres_temp
2316 !      integer :: ilen
2317 !      external :: ilen
2318       logical :: liiflag,lfirst
2319
2320       real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut
2321       real(kind=8), dimension (:,:), allocatable  :: rescore
2322       real(kind=8), dimension (:,:), allocatable :: rescore2
2323       character(len=24) :: tpl_k_rescore
2324       character(len=256) :: pdbfile
2325
2326 !
2327 ! For new homol impl
2328 !
2329 !     include 'COMMON.VAR'
2330 !
2331 !      write (iout,*) "READ_KLAPAUCJUSZ"
2332 !      print *,"READ_KLAPAUCJUSZ"
2333 !      call flush(iout)
2334       call getenv("FRAGFILE",fragfile)
2335       write (iout,*) "Opening", fragfile
2336       call flush(iout)
2337       open(ientin,file=fragfile,status="old",err=10)
2338 !      write (iout,*) " opened"
2339 !      call flush(iout)
2340
2341       sigma_theta=0.0
2342       sigma_d=0.0
2343       sigma_dih=0.0
2344       l_homo = .false.
2345
2346       nres_temp=nres
2347       itype_temp(:)=itype(:,1)
2348       ii=0
2349       lim_odl=0
2350
2351 !      write (iout,*) "Entering loop"
2352 !      call flush(iout)
2353
2354       DO IGR = 1,NCHAIN_GROUP
2355
2356 !      write (iout,*) "igr",igr
2357       call flush(iout)
2358       read(ientin,*) constr_homology,nclust
2359       if (start_from_model) then
2360         nmodel_start=constr_homology
2361       else
2362         nmodel_start=0
2363       endif
2364
2365       ii_old=lim_odl
2366
2367       ichain=iequiv(1,igr)
2368       nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
2369       nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
2370 !      write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
2371
2372 ! Read pdb files
2373       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
2374       if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
2375       do k=1,constr_homology
2376         read(ientin,'(a)') pdbfile
2377         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', &
2378         pdbfile(:ilen(pdbfile))
2379         pdbfiles_chomo(k)=pdbfile
2380         open(ipdbin,file=pdbfile,status='old',err=33)
2381         goto 34
2382   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
2383         pdbfile(:ilen(pdbfile))
2384         stop
2385   34    continue
2386         unres_pdb=.false.
2387 !        nres_temp=nres
2388         call readpdb_template(k)
2389         nres_chomo(k)=nres
2390 !        nres=nres_temp
2391         do i=1,nres
2392           rescore(k,i)=0.2d0
2393           rescore2(k,i)=1.0d0
2394         enddo
2395       enddo
2396 ! Read clusters
2397       do i=1,nclust
2398         read(ientin,*) ninclust(i),nresclust(i)
2399         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
2400         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
2401       enddo
2402 !
2403 ! Loop over clusters
2404 !
2405       do l=1,nclust
2406         do ll = 1,ninclust(l)
2407
2408         k = inclust(ll,l)
2409 !        write (iout,*) "l",l," ll",ll," k",k
2410         do i=1,nres
2411           idomain(k,i)=0
2412         enddo
2413         do i=1,nresclust(l)
2414           if (nnt.gt.1)  then
2415             idomain(k,iresclust(i,l)+1) = 1
2416           else
2417             idomain(k,iresclust(i,l)) = 1
2418           endif
2419         enddo
2420 !
2421 !     Distance restraints
2422 !
2423 !          ... --> odl(k,ii)
2424 ! Copy the coordinates from reference coordinates (?)
2425 !        nres_temp=nres
2426         nres=nres_chomo(k)
2427         do i=1,2*nres
2428           do j=1,3
2429             c(j,i)=chomo(j,i,k)
2430 !           write (iout,*) "c(",j,i,") =",c(j,i)
2431           enddo
2432         enddo
2433         call int_from_cart(.true.,.false.)
2434         call sc_loc_geom(.false.)
2435         do i=1,nres
2436           thetaref(i)=theta(i)
2437           phiref(i)=phi(i)
2438         enddo
2439 !        nres=nres_temp
2440         if (waga_dist.ne.0.0d0) then
2441           ii=ii_old
2442 !          do i = nnt,nct-2
2443           do i = nnt_chain,nct_chain-2
2444 !            do j=i+2,nct
2445             do j=i+2,nct_chain
2446
2447               x12=c(1,i)-c(1,j)
2448               y12=c(2,i)-c(2,j)
2449               z12=c(3,i)-c(3,j)
2450               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2451 !              write (iout,*) k,i,j,distal,dist2_cut
2452
2453             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2454                  .and. distal.le.dist2_cut ) then
2455
2456               ii=ii+1
2457               ii_in_use(ii)=1
2458               l_homo(k,ii)=.true.
2459
2460 !             write (iout,*) "k",k
2461 !             write (iout,*) "i",i," j",j," constr_homology",
2462 !     &                       constr_homology
2463               ires_homo(ii)=i+chain_border1(1,igr)-1
2464               jres_homo(ii)=j+chain_border1(1,igr)-1
2465               odl(k,ii)=distal
2466               if (read2sigma) then
2467                 sigma_odl(k,ii)=0
2468                 do ik=i,j
2469                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2470                 enddo
2471                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2472                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2473              sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2474               else
2475                 if (odl(k,ii).le.dist_cut) then
2476                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2477                 else
2478 #ifdef OLDSIGMA
2479                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2480                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2481 #else
2482                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2483                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2484 #endif
2485                 endif
2486               endif
2487               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2488             else
2489               ii=ii+1
2490 !              l_homo(k,ii)=.false.
2491             endif
2492             enddo
2493           enddo
2494         lim_odl=ii
2495         endif
2496 !
2497 !     Theta, dihedral and SC retraints
2498 !
2499         if (waga_angle.gt.0.0d0) then
2500           do i = nnt_chain+3,nct_chain
2501             iii=i+chain_border1(1,igr)-1
2502             if (idomain(k,i).eq.0) then
2503 !               sigma_dih(k,i)=0.0
2504                cycle
2505             endif
2506             dih(k,iii)=phiref(i)
2507             sigma_dih(k,iii)= &
2508                (rescore(k,i)+rescore(k,i-1)+ &
2509                            rescore(k,i-2)+rescore(k,i-3))/4.0
2510 !            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
2511 !     &       " sigma_dihed",sigma_dih(k,i)
2512             if (sigma_dih(k,iii).ne.0) &
2513              sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
2514           enddo
2515 !          lim_dih=nct-nnt-2
2516         endif
2517
2518         if (waga_theta.gt.0.0d0) then
2519           do i = nnt_chain+2,nct_chain
2520              iii=i+chain_border1(1,igr)-1
2521              if (idomain(k,i).eq.0) then
2522 !              sigma_theta(k,i)=0.0
2523               cycle
2524              endif
2525              thetatpl(k,iii)=thetaref(i)
2526              sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ &
2527                               rescore(k,i-2))/3.0
2528              if (sigma_theta(k,iii).ne.0) &
2529              sigma_theta(k,iii)=1.0d0/ &
2530              (sigma_theta(k,iii)*sigma_theta(k,iii))
2531           enddo
2532         endif
2533
2534         if (waga_d.gt.0.0d0) then
2535           do i = nnt_chain,nct_chain
2536              iii=i+chain_border1(1,igr)-1
2537                if (itype(i,1).eq.10) cycle
2538                if (idomain(k,i).eq.0 ) then
2539 !                  sigma_d(k,i)=0.0
2540                   cycle
2541                endif
2542                xxtpl(k,iii)=xxref(i)
2543                yytpl(k,iii)=yyref(i)
2544                zztpl(k,iii)=zzref(i)
2545                sigma_d(k,iii)=rescore(k,i)
2546                if (sigma_d(k,iii).ne.0) &
2547                 sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
2548 !               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
2549           enddo
2550         endif
2551       enddo ! l
2552       enddo ! ll
2553 !
2554 ! remove distance restraints not used in any model from the list
2555 ! shift data in all arrays
2556 !
2557 !      write (iout,*) "ii_old",ii_old
2558       if (waga_dist.ne.0.0d0) then
2559 #ifdef DEBUG
2560        write (iout,*) "Distance restraints from templates"
2561        do iii=1,lim_odl
2562        write(iout,'(4i5,100(2f8.2,1x,l1,4x))') &
2563         iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), &
2564         (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), &
2565         ki=1,constr_homology)
2566        enddo
2567 #endif
2568         ii=ii_old
2569         liiflag=.true.
2570         lfirst=.true.
2571         do i=nnt_chain,nct_chain-2
2572          do j=i+2,nct_chain
2573           ii=ii+1
2574 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2575 !     &            .and. distal.le.dist2_cut ) then
2576 !          write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
2577 !          call flush(iout)
2578           if (ii_in_use(ii).eq.0.and.liiflag.or. &
2579           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2580             liiflag=.false.
2581             i10=ii
2582             if (lfirst) then
2583               lfirst=.false.
2584               iistart=ii
2585             else
2586               if(i10.eq.lim_odl) i10=i10+1
2587               do ki=0,i10-i01-1
2588                ires_homo(iistart+ki)=ires_homo(ki+i01)
2589                jres_homo(iistart+ki)=jres_homo(ki+i01)
2590                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2591                do k=1,constr_homology
2592                 odl(k,iistart+ki)=odl(k,ki+i01)
2593                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2594                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2595                enddo
2596               enddo
2597               iistart=iistart+i10-i01
2598             endif
2599           endif
2600           if (ii_in_use(ii).ne.0.and..not.liiflag) then
2601              i01=ii
2602              liiflag=.true.
2603           endif
2604          enddo
2605         enddo
2606         lim_odl=iistart-1
2607       endif
2608
2609       lll=lim_odl-ii_old
2610
2611       do i=2,nequiv(igr)
2612
2613         ichain=iequiv(i,igr)
2614
2615         do j=nnt_chain,nct_chain
2616           jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
2617           do k=1,constr_homology
2618             dih(k,jj)=dih(k,j)
2619             sigma_dih(k,jj)=sigma_dih(k,j)
2620             thetatpl(k,jj)=thetatpl(k,j)
2621             sigma_theta(k,jj)=sigma_theta(k,j)
2622             xxtpl(k,jj)=xxtpl(k,j)
2623             yytpl(k,jj)=yytpl(k,j)
2624             zztpl(k,jj)=zztpl(k,j)
2625             sigma_d(k,jj)=sigma_d(k,j)
2626           enddo
2627         enddo
2628
2629         jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
2630 !        write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
2631         do j=ii_old+1,lim_odl
2632           ires_homo(j+lll)=ires_homo(j)+jj
2633           jres_homo(j+lll)=jres_homo(j)+jj
2634           do k=1,constr_homology
2635             odl(k,j+lll)=odl(k,j)
2636             sigma_odl(k,j+lll)=sigma_odl(k,j)
2637             l_homo(k,j+lll)=l_homo(k,j)
2638           enddo
2639         enddo
2640
2641         ii_old=ii_old+lll
2642         lim_odl=lim_odl+lll
2643
2644       enddo
2645
2646       ENDDO ! IGR
2647
2648       if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
2649       nres=nres_temp
2650       itype(:,1)=itype_temp(:)
2651
2652       return
2653    10 stop "Error in fragment file"
2654       end subroutine read_klapaucjusz
2655
2656 !-----------------------------------------------------------
2657       subroutine seq2chains
2658 !c
2659 !c Split the total UNRES sequence, which has dummy residues separating
2660 !c the chains, into separate chains. The length of  chain ichain is
2661 !c contained in chain_length(ichain), the first and last non-dummy
2662 !c residues are in chain_border(1,ichain) and chain_border(2,ichain),
2663 !c respectively. The lengths pertain to non-dummy residues only.
2664 !c
2665 !      implicit none
2666 !      include 'DIMENSIONS'
2667       use energy_data, only:molnum,nchain,chain_length,ireschain
2668       implicit none
2669 !      integer ireschain(nres)
2670       integer ii,ichain,i,j,mnum
2671       logical new_chain
2672       print *,"in seq2"
2673       ichain=1
2674       new_chain=.true.
2675       if (.not.allocated(chain_length)) allocate(chain_length(50))
2676       if (.not.allocated(chain_border)) allocate(chain_border(2,50))
2677
2678       chain_length(ichain)=0
2679       ii=1
2680       do while (ii.lt.nres)
2681         write(iout,*) "in seq2chains",ii,new_chain
2682         mnum=molnum(ii)
2683         if (itype(ii,mnum).eq.ntyp1_molec(mnum)) then
2684           if (.not.new_chain) then
2685             new_chain=.true.
2686             chain_border(2,ichain)=ii-1
2687             ichain=ichain+1
2688             chain_border(1,ichain)=ii+1
2689             chain_length(ichain)=0
2690           endif
2691         else
2692           if (new_chain) then
2693             chain_border(1,ichain)=ii
2694             new_chain=.false.
2695           endif
2696           chain_length(ichain)=chain_length(ichain)+1
2697         endif
2698         ii=ii+1
2699       enddo
2700       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) then
2701         ii=ii-1
2702       else
2703         chain_length(ichain)=chain_length(ichain)+1
2704       endif
2705       if (chain_length(ichain).gt.0) then
2706         chain_border(2,ichain)=ii
2707         nchain=ichain
2708       else
2709         nchain=ichain-1
2710       endif
2711       ireschain=0
2712       do i=1,nchain
2713         do j=chain_border(1,i),chain_border(2,i)
2714           ireschain(j)=i
2715         enddo
2716       enddo
2717       return
2718       end subroutine
2719 !---------------------------------------------------------------------
2720       subroutine chain_symmetry(npermchain,tabpermchain)
2721 !c
2722 !c Determine chain symmetry. nperm is the number of permutations and
2723 !c tabperchain contains the allowed permutations of the chains.
2724 !c
2725 !      implicit none
2726 !      include "DIMENSIONS"
2727 !      include "COMMON.IOUNITS"
2728       implicit none
2729       integer itemp(50),&
2730        npermchain,tabpermchain(50,5040),&
2731        tabperm(50,5040),mapchain(50),&
2732        iequiv(50,nres),iflag(nres)
2733       integer i,j,k,l,ii,nchain_group,nequiv(50),iieq,&
2734        nperm,npermc,ind,mnum
2735       if (nchain.eq.1) then
2736         npermchain=1
2737         tabpermchain(1,1)=1
2738 !c        print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
2739         return
2740       endif
2741 !c
2742 !c Look for equivalent chains
2743 #ifdef DEBUG
2744       write(iout,*) "nchain",nchain
2745       do i=1,nchain
2746         write(iout,*) "chain",i," from",chain_border(1,i),&
2747            " to",chain_border(2,i)
2748         write(iout,*)&
2749         "sequence ",(itype(j,molnum(j)),j=chain_border(1,i),chain_border(2,i))
2750       enddo
2751 #endif
2752       do i=1,nchain
2753         iflag(i)=0
2754       enddo
2755       nchain_group=0
2756       do i=1,nchain
2757         if (iflag(i).gt.0) cycle
2758         iflag(i)=1
2759         nchain_group=nchain_group+1
2760         iieq=1
2761         iequiv(iieq,nchain_group)=i
2762         do j=i+1,nchain
2763           if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
2764 !c          k=0
2765 !c          do while(k.lt.chain_length(i) .and.
2766 !c     &     itype(chain_border(1,i)+k).eq.itype(chain_border(1,j)+k))
2767           do k=0,chain_length(i)-1
2768 !c            k=k+1
2769             mnum=molnum(k)
2770             if (itype(chain_border(1,i)+k,mnum).ne.&
2771                itype(chain_border(1,j)+k,mnum)) exit
2772           enddo
2773           if (k.lt.chain_length(i)) cycle
2774           iflag(j)=1
2775           iieq=iieq+1
2776           iequiv(iieq,nchain_group)=j
2777         enddo
2778         nequiv(nchain_group)=iieq
2779       enddo
2780       write(iout,*) "Number of equivalent chain groups:",nchain_group
2781       write(iout,*) "Equivalent chain groups"
2782       do i=1,nchain_group
2783         write(iout,*) "group",i," #members",nequiv(i)," chains",&
2784            (iequiv(j,i),j=1,nequiv(i))
2785       enddo
2786       ind=0
2787       do i=1,nchain_group
2788         do j=1,nequiv(i)
2789           ind=ind+1
2790           mapchain(ind)=iequiv(j,i)
2791         enddo
2792       enddo
2793       write (iout,*) "mapchain"
2794       do i=1,nchain
2795         write (iout,*) i,mapchain(i)
2796       enddo
2797       ii=0
2798       do i=1,nchain_group
2799         call permut(nequiv(i),nperm,tabperm)
2800         if (ii.eq.0) then
2801           ii=nequiv(i)
2802           npermchain=nperm
2803           do j=1,nperm
2804             do k=1,ii
2805               tabpermchain(k,j)=iequiv(tabperm(k,j),i)
2806             enddo
2807           enddo
2808         else
2809           npermc=npermchain
2810           npermchain=npermchain*nperm
2811           ind=0
2812           do k=1,nperm
2813             do j=1,npermc
2814               ind=ind+1
2815               do l=1,ii
2816                 tabpermchain(l,ind)=tabpermchain(l,j)
2817               enddo
2818               do l=1,nequiv(i)
2819                 tabpermchain(ii+l,ind)=iequiv(tabperm(l,k),i)
2820               enddo
2821             enddo
2822           enddo
2823           ii=ii+nequiv(i)
2824         endif
2825       enddo
2826       do i=1,npermchain
2827         do j=1,nchain
2828           itemp(mapchain(j))=tabpermchain(j,i)
2829         enddo
2830         do j=1,nchain
2831           tabpermchain(j,i)=itemp(j)
2832         enddo
2833       enddo
2834       write(iout,*) "Number of chain permutations",npermchain
2835       write(iout,*) "Permutations"
2836       do i=1,npermchain
2837         write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
2838       enddo
2839       return
2840       end
2841 !c---------------------------------------------------------------------
2842       integer function tperm(i,iperm,tabpermchain)
2843 !      implicit none
2844 !      include 'DIMENSIONS'
2845       integer i,iperm
2846       integer tabpermchain(50,5040)
2847       if (i.eq.0) then
2848         tperm=0
2849       else
2850         tperm=tabpermchain(i,iperm)
2851       endif
2852       return
2853       end function
2854
2855 !-----------------------------------------------------------------------------
2856       end module io