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,i7,a,i7)')'nsup=',nsup,&
1019          ' nstart_sup=',nstart_sup !,"ergwergewrgae"
1020 !el        if(.not.allocated(itype_pdb))
1021         allocate(itype_pdb(nres))
1022         do i=1,nres
1023           itype_pdb(i)=itype(i,1)
1024         enddo
1025         close (ipdbin)
1026         nnt=nstart_sup
1027         nct=nstart_sup+nsup-1
1028 !el        if(.not.allocated(icont_ref))
1029         allocate(icont_ref(2,(nres/2)*nres)) ! maxcont=12*maxres
1030         call contact(.false.,ncont_ref,icont_ref,co)
1031
1032         if (sideadd) then
1033          if(me.eq.king.or..not.out1file) &
1034           write(iout,*)'Adding sidechains'
1035          maxsi=1000
1036          do i=2,nres-1
1037           iti=itype(i,1)
1038           if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then
1039             nsi=0
1040             fail=.true.
1041             do while (fail.and.nsi.le.maxsi)
1042               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail,molnum(i))
1043               nsi=nsi+1
1044             enddo
1045             if(fail) write(iout,*)'Adding sidechain failed for res ',&
1046                     i,' after ',nsi,' trials'
1047           endif
1048          enddo
1049         endif
1050       endif
1051       print *,"CZY TU DOCHODZE"
1052       if (indpdb.eq.0) then
1053       nres_molec(:)=0
1054         allocate(sequence(maxres,5))
1055 !      itype(:,:)=0
1056       itmp=0
1057       if (protein) then
1058 ! Read sequence if not taken from the pdb file.
1059         molec=1
1060         read (inp,*) nres_molec(molec)
1061         print *,'nres=',nres
1062         if (iscode.gt.0) then
1063           read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec))
1064         else
1065           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1066         endif
1067 !        read(inp,*) weightcard_t
1068 !        print *,"po seq" weightcard_t
1069 ! Convert sequence to numeric code
1070
1071         do i=1,nres_molec(molec)
1072           itmp=itmp+1
1073           itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
1074           print *,itype(i,1)
1075
1076         enddo
1077        endif
1078 !        read(inp,*) weightcard_t
1079 !        print *,"po seq", weightcard_t
1080
1081        if (nucleic) then
1082 ! Read sequence if not taken from the pdb file.
1083         molec=2
1084         read (inp,*) nres_molec(molec)
1085         print *,'nres=',nres_molec(molec)
1086 !        allocate(sequence(maxres,5))
1087 !        if (iscode.gt.0) then
1088           read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
1089           print *,"KUR**"
1090           print *,(sequence(i,molec),i=1,nres_molec(molec))
1091 ! Convert sequence to numeric code
1092
1093         do i=1,nres_molec(molec)
1094           itmp=itmp+1
1095           istype(itmp)=sugarcode(sequence(i,molec)(1:1),i)
1096           sequence(i,molec)=sequence(i,molec)(1:2)
1097           itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
1098           write(iout,*) "NUCLE=", itype(itmp,molec)
1099         enddo
1100        endif
1101
1102        if (ions) then
1103 ! Read sequence if not taken from the pdb file.
1104         molec=5
1105         read (inp,*) nres_molec(molec)
1106 !        print *,'nres=',nres
1107           read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1108 ! Convert sequence to numeric code
1109         print *,nres_molec(molec)
1110         do i=1,nres_molec(molec)
1111           itmp=itmp+1
1112           print *,itmp,"itmp"
1113           itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
1114         enddo
1115        endif
1116        nres=0
1117        do i=1,5
1118         nres=nres+nres_molec(i)
1119         print *,"nres_molec",nres,nres_molec(i)
1120        enddo
1121
1122 ! Assign initial virtual bond lengths
1123         if(.not.allocated(molnum)) then
1124          allocate(molnum(nres+1))
1125          itmp=0
1126         do i=1,5
1127                do j=1,nres_molec(i)
1128                itmp=itmp+1
1129               molnum(itmp)=i
1130                enddo
1131          enddo
1132 !        print *,nres_molec(i)
1133         endif
1134         print *,nres,"nres"
1135         if(.not.allocated(vbld)) then
1136            print *, "I DO ENTER"
1137            allocate(vbld(2*nres))
1138         endif
1139         if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
1140         do i=2,nres
1141            mnum=molnum(i)
1142           if (molnum(i).eq.1) then
1143           vbld(i)=vbl
1144           vbld_inv(i)=vblinv
1145
1146           else
1147           write(iout,*) "typy",itype(i,mnum),ntyp1_molec(mnum),i
1148           vbld(i)=6.0
1149           vbld_inv(i)=1.0/6.0
1150           if ((itype(i,mnum).eq.ntyp1_molec(mnum)).or.&
1151           (itype(i-1,mnum).eq.ntyp1_molec(mnum))) then
1152           vbld(i)=1.9
1153           vbld_inv(i)=1.0/1.9
1154           endif
1155           endif
1156         enddo
1157         do i=2,nres-1
1158            mnum=molnum(i)
1159            if (molnum(i).eq.1) then
1160 !          print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i
1161           vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
1162           vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
1163            else
1164           write(iout,*) "typy2",itype(i,mnum),ntyp1_molec(mnum),i
1165           if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
1166           vbld_inv(i+nres)=1.0
1167           vbld(i+nres)=0.0
1168           else
1169           vbld(i+nres)=vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
1170           vbld_inv(i+nres)=1.0/vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
1171           endif
1172            endif
1173 !          write (iout,*) "i",i," itype",itype(i,1),
1174 !     &      " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
1175         enddo
1176       endif
1177 !      print *,nres
1178 !      print '(20i4)',(itype(i,1),i=1,nres)
1179 !----------------------------
1180 !el reallocate tables
1181 !      do i=1,maxres2
1182 !        do j=1,3
1183 !          c_alloc(j,i)=c(j,i)
1184 !          dc_alloc(j,i)=dc(j,i)
1185 !        enddo
1186 !      enddo
1187 !      do i=1,maxres
1188 !elwrite(iout,*) "itype",i,itype(i,1)
1189 !        itype_alloc(i)=itype(i,1)
1190 !      enddo
1191
1192 !      deallocate(c)
1193 !      deallocate(dc)
1194 !      deallocate(itype)
1195 !      allocate(c(3,2*nres+4))
1196 !      allocate(dc(3,0:2*nres+2))
1197 !      allocate(itype(nres+2))
1198       allocate(itel(nres+2))
1199       itel(:)=0
1200
1201 !      do i=1,2*nres+2
1202 !        do j=1,3
1203 !          c(j,i)=c_alloc(j,i)
1204 !          dc(j,i)=dc_alloc(j,i)
1205 !        enddo
1206 !      enddo
1207 !      do i=1,nres+2
1208 !        itype(i,1)=itype_alloc(i)
1209 !        itel(i)=0
1210 !      enddo
1211 !--------------------------
1212       do i=1,nres
1213 #ifdef PROCOR
1214         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then
1215 #else
1216         if (itype(i,1).eq.ntyp1) then
1217 #endif
1218           itel(i)=0
1219 #ifdef PROCOR
1220         else if (iabs(itype(i+1,1)).ne.20) then
1221 #else
1222         else if (iabs(itype(i,1)).ne.20) then
1223 #endif
1224           itel(i)=1
1225         else
1226           itel(i)=2
1227         endif
1228       enddo
1229       if(me.eq.king.or..not.out1file)then
1230        write (iout,*) "ITEL"
1231        print *,nres,"nres"
1232        do i=1,nres-1
1233          write (iout,*) i,itype(i,1),itel(i)
1234        enddo
1235        print *,'Call Read_Bridge.'
1236       endif
1237       call read_bridge
1238 !--------------------------------
1239 !       print *,"tu dochodze"
1240 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1241       call alloc_geo_arrays
1242       print *,"after  alloc_geo_arrays"
1243       call alloc_ener_arrays
1244       print *,"after  alloc_ener_arrays"
1245 !--------------------------------
1246 ! 8/13/98 Set limits to generating the dihedral angles
1247       do i=1,nres
1248         phibound(1,i)=-pi
1249         phibound(2,i)=pi
1250       enddo
1251       read (inp,*) ndih_constr
1252       if (ndih_constr.gt.0) then
1253         raw_psipred=.false.
1254         allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1255         allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1256         allocate(ftors(ndih_constr)) !(maxdih_constr)
1257
1258 !        read (inp,*) ftors
1259         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), &
1260         i=1,ndih_constr)
1261         if(me.eq.king.or..not.out1file)then
1262          write (iout,*) &
1263          'There are',ndih_constr,' constraints on phi angles.'
1264          do i=1,ndih_constr
1265           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), &
1266           ftors(i)
1267          enddo
1268         endif
1269         do i=1,ndih_constr
1270           phi0(i)=deg2rad*phi0(i)
1271           drange(i)=deg2rad*drange(i)
1272         enddo
1273 !        if(me.eq.king.or..not.out1file) &
1274 !         write (iout,*) 'FTORS',ftors
1275         do i=1,ndih_constr
1276           ii = idih_constr(i)
1277           phibound(1,ii) = phi0(i)-drange(i)
1278           phibound(2,ii) = phi0(i)+drange(i)
1279         enddo
1280       else if (ndih_constr.lt.0) then
1281         raw_psipred=.true.
1282         allocate(idih_constr(nres))
1283         allocate(secprob(3,nres))
1284         allocate(vpsipred(3,nres))
1285         allocate(sdihed(2,nres))
1286         call card_concat(weightcard,.true.)
1287         call reada(weightcard,"PHIHEL",phihel,50.0D0)
1288         call reada(weightcard,"PHIBET",phibet,180.0D0)
1289         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
1290         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
1291         call reada(weightcard,"WDIHC",wdihc,0.591D0)
1292         write (iout,*) "Weight of dihedral angle restraints",wdihc
1293         read(inp,'(9x,3f7.3)') &
1294           (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
1295         write (iout,*) "The secprob array"
1296         do i=nnt,nct
1297           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
1298         enddo
1299         ndih_constr=0
1300         do i=nnt+3,nct
1301           if (itype(i-3,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1 &
1302           .and. itype(i-1,1).ne.ntyp1 .and. itype(i,1).ne.ntyp1) then
1303             ndih_constr=ndih_constr+1
1304             idih_constr(ndih_constr)=i
1305             sumv=0.0d0
1306             do j=1,3
1307               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
1308               sumv=sumv+vpsipred(j,ndih_constr)
1309             enddo
1310             do j=1,3
1311               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
1312             enddo
1313             phibound(1,ndih_constr)=phihel*deg2rad
1314             phibound(2,ndih_constr)=phibet*deg2rad
1315             sdihed(1,ndih_constr)=sigmahel*deg2rad
1316             sdihed(2,ndih_constr)=sigmabet*deg2rad
1317           endif
1318         enddo
1319
1320       endif
1321       if (with_theta_constr) then
1322 !C with_theta_constr is keyword allowing for occurance of theta constrains
1323       read (inp,*) ntheta_constr
1324 !C ntheta_constr is the number of theta constrains
1325       if (ntheta_constr.gt.0) then
1326 !C        read (inp,*) ftors
1327         allocate(itheta_constr(ntheta_constr))
1328         allocate(theta_constr0(ntheta_constr))
1329         allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr))
1330         read (inp,*) (itheta_constr(i),theta_constr0(i), &
1331        theta_drange(i),for_thet_constr(i), &
1332        i=1,ntheta_constr)
1333 !C the above code reads from 1 to ntheta_constr
1334 !C itheta_constr(i) residue i for which is theta_constr
1335 !C theta_constr0 the global minimum value
1336 !C theta_drange is range for which there is no energy penalty
1337 !C for_thet_constr is the force constant for quartic energy penalty
1338 !C E=k*x**4
1339         if(me.eq.king.or..not.out1file)then
1340          write (iout,*) &
1341         'There are',ntheta_constr,' constraints on phi angles.'
1342          do i=1,ntheta_constr
1343           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), &
1344          theta_drange(i), &
1345          for_thet_constr(i)
1346          enddo
1347         endif
1348         do i=1,ntheta_constr
1349           theta_constr0(i)=deg2rad*theta_constr0(i)
1350           theta_drange(i)=deg2rad*theta_drange(i)
1351         enddo
1352 !C        if(me.eq.king.or..not.out1file)
1353 !C     &   write (iout,*) 'FTORS',ftors
1354 !C        do i=1,ntheta_constr
1355 !C          ii = itheta_constr(i)
1356 !C          thetabound(1,ii) = phi0(i)-drange(i)
1357 !C          thetabound(2,ii) = phi0(i)+drange(i)
1358 !C        enddo
1359       endif ! ntheta_constr.gt.0
1360       endif! with_theta_constr
1361
1362       nnt=1
1363 #ifdef MPI
1364       if (me.eq.king) then
1365 #endif
1366        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1367        do i=1,nres
1368          write (iout,'(a3,i5,2f10.1)') &
1369          restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1370        enddo
1371 #ifdef MP
1372       endif
1373 #endif
1374       nct=nres
1375       allocate(ireschain(nres))
1376       ireschain=0
1377       write(iout,*),"before seq2chains",ireschain
1378       call seq2chains
1379       write(iout,*) "after seq2chains",nchain
1380       allocate ( chain_border1(2,nchain))
1381       chain_border1(1,1)=1
1382       chain_border1(2,1)=chain_border(2,1)+1
1383       do i=2,nchain-1
1384         chain_border1(1,i)=chain_border(1,i)-1
1385         chain_border1(2,i)=chain_border(2,i)+1
1386       enddo
1387       if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
1388       chain_border1(2,nchain)=nres
1389       write(iout,*) "nres",nres," nchain",nchain
1390       do i=1,nchain
1391         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),&
1392          chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
1393       enddo
1394       allocate(tabpermchain(50,5040))
1395       call chain_symmetry(npermchain,tabpermchain)
1396       print *,'NNT=',NNT,' NCT=',NCT
1397       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1398       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1399       if (pdbref) then
1400         if(me.eq.king.or..not.out1file) &
1401          write (iout,'(a,i7)') 'nsup=',nsup
1402         nstart_seq=nnt
1403         if (nsup.le.(nct-nnt+1)) then
1404           do i=0,nct-nnt+1-nsup
1405             if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
1406               nstart_seq=nnt+i
1407               goto 111
1408             endif
1409           enddo
1410           write (iout,'(a)') &
1411                   'Error - sequences to be superposed do not match.'
1412           stop
1413         else
1414           do i=0,nsup-(nct-nnt+1)
1415             if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1416             then
1417               nstart_sup=nstart_sup+i
1418               nsup=nct-nnt+1
1419               goto 111
1420             endif
1421           enddo
1422           write (iout,'(a)') &
1423                   'Error - sequences to be superposed do not match.'
1424         endif
1425   111   continue
1426         if (nsup.eq.0) nsup=nct-nnt+1
1427         if (nstart_sup.eq.0) nstart_sup=nnt
1428         if (nstart_seq.eq.0) nstart_seq=nnt
1429         if(me.eq.king.or..not.out1file) &
1430          write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1431                        ' nstart_seq=',nstart_seq !,"242343453254"
1432       endif
1433 !--- Zscore rms -------
1434       if (nz_start.eq.0) nz_start=nnt
1435       if (nz_end.eq.0 .and. nsup.gt.0) then
1436         nz_end=nnt+nsup-1
1437       else if (nz_end.eq.0) then
1438         nz_end=nct
1439       endif
1440       if(me.eq.king.or..not.out1file)then
1441        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1442        write (iout,*) 'IZ_SC=',iz_sc, 'NCT=',nct
1443       endif
1444 !----------------------
1445       call init_int_table
1446       if (refstr) then
1447         if (.not.pdbref) then
1448           call read_angles(inp,*38)
1449           goto 39
1450    38     write (iout,'(a)') 'Error reading reference structure.'
1451 #ifdef MPI
1452           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1453           stop 'Error reading reference structure'
1454 #endif
1455    39     call chainbuild
1456           call setup_var
1457 !zscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1458           nstart_sup=nnt
1459           nstart_seq=nnt
1460           nsup=nct-nnt+1
1461           kkk=1
1462           do i=1,2*nres
1463             do j=1,3
1464               cref(j,i,kkk)=c(j,i)
1465             enddo
1466           enddo
1467           call contact(.true.,ncont_ref,icont_ref,co)
1468         endif
1469 !        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1470 !        call flush(iout)
1471 !EL        if (constr_dist.gt.0) call read_dist_constr
1472 !EL        write (iout,*) "After read_dist_constr nhpb",nhpb
1473 !EL        if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1474 !EL        call hpb_partition
1475         if(me.eq.king.or..not.out1file) &
1476          write (iout,*) 'Contact order:',co
1477         if (pdbref) then
1478         if(me.eq.king.or..not.out1file) &
1479          write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1480         do i=1,ncont_ref
1481           do j=1,2
1482             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1483           enddo
1484           if(me.eq.king.or..not.out1file) &
1485            write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1486            icont_ref(1,i),' ',&
1487            restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1488         enddo
1489         endif
1490       if (constr_homology.gt.0) then
1491 !        write (iout,*) "Calling read_constr_homology"
1492 !        call flush(iout)
1493         call read_constr_homology
1494         if (indpdb.gt.0 .or. pdbref) then
1495           do i=1,2*nres
1496             do j=1,3
1497               c(j,i)=crefjlee(j,i)
1498               cref(j,i,1)=crefjlee(j,i)
1499             enddo
1500           enddo
1501         endif
1502 #define DEBUG
1503 #ifdef DEBUG
1504         write (iout,*) "sc_loc_geom: Array C"
1505         do i=1,nres
1506           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),&
1507            (c(j,i+nres),j=1,3)
1508         enddo
1509         write (iout,*) "Array Cref"
1510         do i=1,nres
1511           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),&
1512            (cref(j,i+nres,1),j=1,3)
1513         enddo
1514 #endif
1515        call int_from_cart1(.false.)
1516        call sc_loc_geom(.false.)
1517        do i=1,nres
1518          thetaref(i)=theta(i)
1519          phiref(i)=phi(i)
1520        enddo
1521        do i=1,nres-1
1522          do j=1,3
1523            dc(j,i)=c(j,i+1)-c(j,i)
1524            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1525          enddo
1526        enddo
1527        do i=2,nres-1
1528          do j=1,3
1529            dc(j,i+nres)=c(j,i+nres)-c(j,i)
1530            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1531          enddo
1532        enddo
1533       else
1534         homol_nset=0
1535         if (start_from_model) then
1536           nmodel_start=0
1537           do
1538             read(inp,'(a)',end=332,err=332) pdbfile
1539             if (me.eq.king .or. .not. out1file)&
1540              write (iout,'(a,5x,a)') 'Opening PDB file',&
1541              pdbfile(:ilen(pdbfile))
1542             open(ipdbin,file=pdbfile,status='old',err=336)
1543             goto 335
1544  336        write (iout,'(a,5x,a)') 'Error opening PDB file',&
1545            pdbfile(:ilen(pdbfile))
1546             call flush(iout)
1547             stop
1548  335        continue
1549             unres_pdb=.false.
1550             nres_temp=nres
1551 !            call readpdb
1552             call readpdb_template(nmodel_start+1)
1553             close(ipdbin)
1554             if (nres.ge.nres_temp) then
1555               nmodel_start=nmodel_start+1
1556               pdbfiles_chomo(nmodel_start)=pdbfile
1557               do i=1,2*nres
1558                 do j=1,3
1559                   chomo(j,i,nmodel_start)=c(j,i)
1560                 enddo
1561               enddo
1562             else
1563               if (me.eq.king .or. .not. out1file) &
1564                write (iout,'(a,2i7,1x,a)') &
1565                 "Different number of residues",nres_temp,nres, &
1566                 " model skipped."
1567             endif
1568             nres=nres_temp
1569           enddo
1570   332     continue
1571           if (nmodel_start.eq.0) then
1572             if (me.eq.king .or. .not. out1file) &
1573              write (iout,'(a)') &
1574              "No valid starting model found START_FROM_MODELS is OFF"
1575               start_from_model=.false.
1576           endif
1577           write (iout,*) "nmodel_start",nmodel_start
1578         endif
1579       endif
1580
1581       endif
1582         if (constr_dist.gt.0) call read_dist_constr
1583         write (iout,*) "After read_dist_constr nhpb",nhpb
1584         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1585         if (velnanoconst.ne.0) call read_afmnano
1586         call hpb_partition
1587
1588       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1589           .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1590           modecalc.ne.10) then
1591 ! If input structure hasn't been supplied from the PDB file read or generate
1592 ! initial geometry.
1593         if (iranconf.eq.0 .and. .not. extconf) then
1594           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1595            write (iout,'(a)') 'Initial geometry will be read in.'
1596           if (read_cart) then
1597             read(inp,'(8f10.5)',end=36,err=36) &
1598              ((c(l,k),l=1,3),k=1,nres),&
1599              ((c(l,k+nres),l=1,3),k=nnt,nct)
1600             write (iout,*) "Exit READ_CART"
1601             write (iout,'(8f10.5)') &
1602              ((c(l,k),l=1,3),k=1,nres)
1603             write (iout,'(8f10.5)') &
1604              ((c(l,k+nres),l=1,3),k=nnt,nct)
1605             call int_from_cart1(.true.)
1606             write (iout,*) "Finish INT_TO_CART"
1607             do i=1,nres-1
1608               do j=1,3
1609                 dc(j,i)=c(j,i+1)-c(j,i)
1610                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1611               enddo
1612             enddo
1613             do i=nnt,nct
1614               if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1615                 do j=1,3
1616                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1617                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1618                 enddo
1619               endif
1620             enddo
1621             return
1622           else
1623            write(iout,*) "read angles from input"
1624            call read_angles(inp,*36)
1625             call chainbuild
1626
1627           endif
1628           goto 37
1629    36     write (iout,'(a)') 'Error reading angle file.'
1630 #ifdef MPI
1631           call mpi_finalize( MPI_COMM_WORLD,IERR )
1632 #endif
1633           stop 'Error reading angle file.'
1634    37     continue
1635         else if (extconf) then
1636          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1637           write (iout,'(a)') 'Extended chain initial geometry.'
1638          do i=3,nres
1639           theta(i)=90d0*deg2rad
1640           if (molnum(i).eq.2) theta(i)=160d0*deg2rad
1641          enddo
1642          do i=4,nres
1643           phi(i)=180d0*deg2rad
1644           if (molnum(i).eq.2) phi(i)=30d0*deg2rad
1645          enddo
1646          do i=2,nres-1
1647           alph(i)=110d0*deg2rad
1648           if (molnum(i).eq.2) alph(i)=30d0*deg2rad
1649          enddo
1650          do i=2,nres-1
1651           omeg(i)=-120d0*deg2rad
1652           if (molnum(i).eq.2) omeg(i)=60d0*deg2rad
1653           if (itype(i,1).le.0) omeg(i)=-omeg(i)
1654          enddo
1655          call chainbuild
1656         else
1657           if(me.eq.king.or..not.out1file) &
1658            write (iout,'(a)') 'Random-generated initial geometry.'
1659
1660
1661 #ifdef MPI
1662           if (me.eq.king  .or. fg_rank.eq.0 .and. &
1663                  ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1664 #endif
1665             do itrial=1,100
1666               itmp=1
1667               call gen_rand_conf(itmp,*30)
1668               goto 40
1669    30         write (iout,*) 'Failed to generate random conformation',&
1670                 ', itrial=',itrial
1671               write (*,*) 'Processor:',me,&
1672                 ' Failed to generate random conformation',&
1673                 ' itrial=',itrial
1674               call intout
1675
1676 #ifdef AIX
1677               call flush_(iout)
1678 #else
1679               call flush(iout)
1680 #endif
1681             enddo
1682             write (iout,'(a,i3,a)') 'Processor:',me,&
1683               ' error in generating random conformation.'
1684             write (*,'(a,i3,a)') 'Processor:',me,&
1685               ' error in generating random conformation.'
1686             call flush(iout)
1687 #ifdef MPI
1688             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1689    40       continue
1690           endif
1691 #else
1692           do itrial=1,100
1693             itmp=1
1694             call gen_rand_conf(itmp,*335)
1695             goto 40
1696   335       write (iout,*) 'Failed to generate random conformation',&
1697               ', itrial=',itrial
1698             write (*,*) 'Failed to generate random conformation',&
1699               ', itrial=',itrial
1700           enddo
1701           write (iout,'(a,i3,a)') 'Processor:',me,&
1702             ' error in generating random conformation.'
1703           write (*,'(a,i3,a)') 'Processor:',me,&
1704             ' error in generating random conformation.'
1705           stop
1706    40     continue
1707 #endif
1708         endif
1709       elseif (modecalc.eq.4) then
1710         read (inp,'(a)') intinname
1711         open (intin,file=intinname,status='old',err=333)
1712         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1713         write (iout,'(a)') 'intinname',intinname
1714         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1715         goto 334
1716   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1717 #ifdef MPI
1718         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1719 #endif
1720         stop 'Error opening angle file.'
1721   334   continue
1722
1723       endif
1724 ! Generate distance constraints, if the PDB structure is to be regularized.
1725       if (nthread.gt.0) then
1726         call read_threadbase
1727       endif
1728       call setup_var
1729       if (me.eq.king .or. .not. out1file) &
1730        call intout
1731       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1732         write (iout,'(/a,i3,a)') &
1733         'The chain contains',ns,' disulfide-bridging cysteines.'
1734         write (iout,'(20i4)') (iss(i),i=1,ns)
1735        if (dyn_ss) then
1736           write(iout,*)"Running with dynamic disulfide-bond formation"
1737        else
1738         write (iout,'(/a/)') 'Pre-formed links are:'
1739         do i=1,nss
1740           i1=ihpb(i)-nres
1741           i2=jhpb(i)-nres
1742           it1=itype(i1,1)
1743           it2=itype(i2,1)
1744           if (me.eq.king.or..not.out1file) &
1745           write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1746           restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1747           ebr,forcon(i)
1748         enddo
1749         write (iout,'(a)')
1750        endif
1751       endif
1752       if (ns.gt.0.and.dyn_ss) then
1753           do i=nss+1,nhpb
1754             ihpb(i-nss)=ihpb(i)
1755             jhpb(i-nss)=jhpb(i)
1756             forcon(i-nss)=forcon(i)
1757             dhpb(i-nss)=dhpb(i)
1758           enddo
1759           nhpb=nhpb-nss
1760           nss=0
1761           call hpb_partition
1762           do i=1,ns
1763             dyn_ss_mask(iss(i))=.true.
1764           enddo
1765       endif
1766       if (i2ndstr.gt.0) call secstrp2dihc
1767       if (indpdb.gt.0) then
1768           write(iout,*) "WCHODZE TU!!"
1769           call int_from_cart1(.true.)
1770       endif
1771 !      call geom_to_var(nvar,x)
1772 !      call etotal(energia(0))
1773 !      call enerprint(energia(0))
1774 !      call briefout(0,etot)
1775 !      stop
1776 !d    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1777 !d    write (iout,'(a)') 'Variable list:'
1778 !d    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1779 #ifdef MPI
1780       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1781         write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1782         'Processor',myrank,': end reading molecular data.'
1783 #endif
1784       return
1785       end subroutine molread
1786 !-----------------------------------------------------------------------------
1787       subroutine read_constr_homology
1788       use control, only:init_int_table,homology_partition
1789       use MD_data, only:iset
1790 !      implicit none
1791 !      include 'DIMENSIONS'
1792 !#ifdef MPI
1793 !      include 'mpif.h'
1794 !#endif
1795 !      include 'COMMON.SETUP'
1796 !      include 'COMMON.CONTROL'
1797 !      include 'COMMON.HOMOLOGY'
1798 !      include 'COMMON.CHAIN'
1799 !      include 'COMMON.IOUNITS'
1800 !      include 'COMMON.MD'
1801 !      include 'COMMON.QRESTR'
1802 !      include 'COMMON.GEO'
1803 !      include 'COMMON.INTERACT'
1804 !      include 'COMMON.NAMES'
1805 !      include 'COMMON.VAR'
1806 !
1807
1808 !     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1809 !    &                 dist_cut
1810 !     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1811 !    &    sigma_odl_temp(maxres,maxres,max_template)
1812       character*2 kic2
1813       character*24 model_ki_dist, model_ki_angle
1814       character*500 controlcard
1815       integer :: ki,i,ii,j,k,l
1816       integer, dimension (:), allocatable :: ii_in_use
1817       integer :: i_tmp,idomain_tmp,&
1818       irec,ik,iistart,nres_temp
1819 !      integer :: iset
1820 !      external :: ilen
1821       logical :: liiflag,lfirst
1822       integer :: i01,i10
1823 !
1824 !     FP - Nov. 2014 Temporary specifications for new vars
1825 !
1826       real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
1827                        rescore3_tmp, dist_cut
1828       real(kind=8), dimension (:,:),allocatable :: rescore
1829       real(kind=8), dimension (:,:),allocatable :: rescore2
1830       real(kind=8), dimension (:,:),allocatable :: rescore3
1831       real(kind=8) :: distal
1832       character*24 tpl_k_rescore
1833       character*256 pdbfile
1834
1835 ! -----------------------------------------------------------------
1836 ! Reading multiple PDB ref structures and calculation of retraints
1837 ! not using pre-computed ones stored in files model_ki_{dist,angle}
1838 ! FP (Nov., 2014)
1839 ! -----------------------------------------------------------------
1840 !
1841 !
1842 ! Alternative: reading from input
1843       call card_concat(controlcard,.true.)
1844       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1845       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1846       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1847       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1848       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1849       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1850       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1851       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1852       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
1853       if(.not.read2sigma.and.start_from_model) then
1854           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
1855            write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
1856           start_from_model=.false.
1857           iranconf=(indpdb.le.0)
1858       endif
1859       if(start_from_model .and. (me.eq.king .or. .not. out1file))&
1860          write(iout,*) 'START_FROM_MODELS is ON'
1861 !      if(start_from_model .and. rest) then
1862 !        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1863 !         write(iout,*) 'START_FROM_MODELS is OFF'
1864 !         write(iout,*) 'remove restart keyword from input'
1865 !        endif
1866 !      endif
1867       if (start_from_model) nmodel_start=constr_homology
1868       if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
1869       if (homol_nset.gt.1)then
1870          call card_concat(controlcard,.true.)
1871          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1872          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1873 !          write(iout,*) "iset homology_weight "
1874           do i=1,homol_nset
1875            write(iout,*) i,waga_homology(i)
1876           enddo
1877          endif
1878          iset=mod(kolor,homol_nset)+1
1879       else
1880        iset=1
1881        waga_homology(1)=1.0
1882       endif
1883
1884 !d      write (iout,*) "nnt",nnt," nct",nct
1885 !d      call flush(iout)
1886
1887       if (read_homol_frag) then
1888         call read_klapaucjusz
1889       else
1890
1891       lim_odl=0
1892       lim_dih=0
1893 !
1894 !      write(iout,*) 'nnt=',nnt,'nct=',nct
1895 !
1896 !      do i = nnt,nct
1897 !        do k=1,constr_homology
1898 !         idomain(k,i)=0
1899 !        enddo
1900 !      enddo
1901        idomain=0
1902
1903 !      ii=0
1904 !      do i = nnt,nct-2
1905 !        do j=i+2,nct
1906 !        ii=ii+1
1907 !        ii_in_use(ii)=0
1908 !        enddo
1909 !      enddo
1910       ii_in_use=0
1911       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
1912       if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
1913
1914       do k=1,constr_homology
1915
1916         read(inp,'(a)') pdbfile
1917         pdbfiles_chomo(k)=pdbfile
1918         if(me.eq.king .or. .not. out1file) &
1919          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
1920         pdbfile(:ilen(pdbfile))
1921         open(ipdbin,file=pdbfile,status='old',err=33)
1922         goto 34
1923   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
1924         pdbfile(:ilen(pdbfile))
1925         stop
1926   34    continue
1927 !        print *,'Begin reading pdb data'
1928 !
1929 ! Files containing res sim or local scores (former containing sigmas)
1930 !
1931
1932         write(kic2,'(bz,i2.2)') k
1933
1934         tpl_k_rescore="template"//kic2//".sco"
1935         write(iout,*) "tpl_k_rescore=",tpl_k_rescore
1936         unres_pdb=.false.
1937         nres_temp=nres
1938         write(iout,*) "read2sigma",read2sigma
1939
1940         if (read2sigma) then
1941           call readpdb_template(k)
1942         else
1943           call readpdb
1944         endif
1945         write(iout,*) "after readpdb"
1946         if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
1947         nres_chomo(k)=nres
1948         nres=nres_temp
1949         if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
1950         if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
1951         if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
1952         if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
1953         if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
1954         if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
1955         if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
1956         if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
1957         if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
1958         if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
1959         if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
1960         if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
1961         if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1962         if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
1963 !        if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1964         if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
1965         if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
1966         if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
1967         if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
1968 !        if(.not.allocated(distance)) allocate(distance(constr_homology))
1969 !        if(.not.allocated(distancek)) allocate(distancek(constr_homology))
1970
1971
1972 !
1973 !     Distance restraints
1974 !
1975 !          ... --> odl(k,ii)
1976 ! Copy the coordinates from reference coordinates (?)
1977         do i=1,2*nres_chomo(k)
1978           do j=1,3
1979             c(j,i)=cref(j,i,1)
1980 !           write (iout,*) "c(",j,i,") =",c(j,i)
1981           enddo
1982         enddo
1983 !
1984 ! From read_dist_constr (commented out 25/11/2014 <-> res sim)
1985 !
1986 !         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1987           open (ientin,file=tpl_k_rescore,status='old')
1988           if (nnt.gt.1) rescore(k,1)=0.0d0
1989           do irec=nnt,nct ! loop for reading res sim
1990             if (read2sigma) then
1991              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
1992                                      rescore3_tmp,idomain_tmp
1993              i_tmp=i_tmp+nnt-1
1994              idomain(k,i_tmp)=idomain_tmp
1995              rescore(k,i_tmp)=rescore_tmp
1996              rescore2(k,i_tmp)=rescore2_tmp
1997              rescore3(k,i_tmp)=rescore3_tmp
1998              if (.not. out1file .or. me.eq.king)&
1999              write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
2000                            i_tmp,rescore2_tmp,rescore_tmp,&
2001                                      rescore3_tmp,idomain_tmp
2002             else
2003              idomain(k,irec)=1
2004              read (ientin,*,end=1401) rescore_tmp
2005
2006 !           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2007              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2008 !           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2009             endif
2010           enddo
2011  1401   continue
2012         close (ientin)
2013         if (waga_dist.ne.0.0d0) then
2014           ii=0
2015           do i = nnt,nct-2
2016             do j=i+2,nct
2017
2018               x12=c(1,i)-c(1,j)
2019               y12=c(2,i)-c(2,j)
2020               z12=c(3,i)-c(3,j)
2021               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2022 !              write (iout,*) k,i,j,distal,dist2_cut
2023
2024             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2025                  .and. distal.le.dist2_cut ) then
2026
2027               ii=ii+1
2028               ii_in_use(ii)=1
2029               l_homo(k,ii)=.true.
2030
2031 !             write (iout,*) "k",k
2032 !             write (iout,*) "i",i," j",j," constr_homology",
2033 !    &                       constr_homology
2034               ires_homo(ii)=i
2035               jres_homo(ii)=j
2036               odl(k,ii)=distal
2037               if (read2sigma) then
2038                 sigma_odl(k,ii)=0
2039                 do ik=i,j
2040                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2041                 enddo
2042                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2043                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2044               sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2045               else
2046                 if (odl(k,ii).le.dist_cut) then
2047                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2048                 else
2049 #ifdef OLDSIGMA
2050                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2051                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2052 #else
2053                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2054                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2055 #endif
2056                 endif
2057               endif
2058               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2059             else
2060 !              ii=ii+1
2061 !              l_homo(k,ii)=.false.
2062             endif
2063             enddo
2064           enddo
2065         lim_odl=ii
2066         endif
2067 !        write (iout,*) "Distance restraints set"
2068 !        call flush(iout)
2069 !
2070 !     Theta, dihedral and SC retraints
2071 !
2072         if (waga_angle.gt.0.0d0) then
2073 !         open (ientin,file=tpl_k_sigma_dih,status='old')
2074 !         do irec=1,maxres-3 ! loop for reading sigma_dih
2075 !            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2076 !            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2077 !            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2078 !    &                            sigma_dih(k,i+nnt-1)
2079 !         enddo
2080 !1402   continue
2081 !         close (ientin)
2082           do i = nnt+3,nct
2083             if (idomain(k,i).eq.0) then
2084                sigma_dih(k,i)=0.0
2085                cycle
2086             endif
2087             dih(k,i)=phiref(i) ! right?
2088 !           read (ientin,*) sigma_dih(k,i) ! original variant
2089 !             write (iout,*) "dih(",k,i,") =",dih(k,i)
2090 !             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2091 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2092 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
2093 !    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
2094
2095             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2096                           rescore(k,i-2)+rescore(k,i-3))/4.0
2097 !            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2098 !           write (iout,*) "Raw sigmas for dihedral angle restraints"
2099 !           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2100 !           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2101 !                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
2102 !   Instead of res sim other local measure of b/b str reliability possible
2103             if (sigma_dih(k,i).ne.0) &
2104             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2105 !           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2106           enddo
2107           lim_dih=nct-nnt-2
2108         endif
2109 !        write (iout,*) "Dihedral angle restraints set"
2110 !        call flush(iout)
2111
2112         if (waga_theta.gt.0.0d0) then
2113 !         open (ientin,file=tpl_k_sigma_theta,status='old')
2114 !         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2115 !            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2116 !            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2117 !    &                              sigma_theta(k,i+nnt-1)
2118 !         enddo
2119 !1403   continue
2120 !         close (ientin)
2121
2122           do i = nnt+2,nct ! right? without parallel.
2123 !         do i = i=1,nres ! alternative for bounds acc to readpdb?
2124 !         do i=ithet_start,ithet_end ! with FG parallel.
2125              if (idomain(k,i).eq.0) then
2126               sigma_theta(k,i)=0.0
2127               cycle
2128              endif
2129              thetatpl(k,i)=thetaref(i)
2130 !            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2131 !            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
2132 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2133 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
2134 !            read (ientin,*) sigma_theta(k,i) ! 1st variant
2135              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2136                              rescore(k,i-2))/3.0
2137 !             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2138              if (sigma_theta(k,i).ne.0) &
2139              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2140
2141 !            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2142 !                             rescore(k,i-2) !  right expression ?
2143 !            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2144           enddo
2145         endif
2146 !        write (iout,*) "Angle restraints set"
2147 !        call flush(iout)
2148
2149         if (waga_d.gt.0.0d0) then
2150 !       open (ientin,file=tpl_k_sigma_d,status='old')
2151 !         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2152 !            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2153 !            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2154 !    &                          sigma_d(k,i+nnt-1)
2155 !         enddo
2156 !1404   continue
2157
2158           do i = nnt,nct ! right? without parallel.
2159 !         do i=2,nres-1 ! alternative for bounds acc to readpdb?
2160 !         do i=loc_start,loc_end ! with FG parallel.
2161                if (itype(i,1).eq.10) cycle
2162                if (idomain(k,i).eq.0 ) then
2163                   sigma_d(k,i)=0.0
2164                   cycle
2165                endif
2166                xxtpl(k,i)=xxref(i)
2167                yytpl(k,i)=yyref(i)
2168                zztpl(k,i)=zzref(i)
2169 !              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2170 !              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2171 !              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2172 !              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
2173                sigma_d(k,i)=rescore3(k,i) !  right expression ?
2174                if (sigma_d(k,i).ne.0) &
2175                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2176
2177 !              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
2178 !              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2179 !              read (ientin,*) sigma_d(k,i) ! 1st variant
2180           enddo
2181         endif
2182       enddo
2183 !      write (iout,*) "SC restraints set"
2184 !      call flush(iout)
2185 !
2186 ! remove distance restraints not used in any model from the list
2187 ! shift data in all arrays
2188 !
2189 !      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
2190       if (waga_dist.ne.0.0d0) then
2191         ii=0
2192         liiflag=.true.
2193         lfirst=.true.
2194         do i=nnt,nct-2
2195          do j=i+2,nct
2196           ii=ii+1
2197 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2198 !     &            .and. distal.le.dist2_cut ) then
2199 !          write (iout,*) "i",i," j",j," ii",ii
2200 !          call flush(iout)
2201           if (ii_in_use(ii).eq.0.and.liiflag.or. &
2202           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2203             liiflag=.false.
2204             i10=ii
2205             if (lfirst) then
2206               lfirst=.false.
2207               iistart=ii
2208             else
2209               if(i10.eq.lim_odl) i10=i10+1
2210               do ki=0,i10-i01-1
2211                ires_homo(iistart+ki)=ires_homo(ki+i01)
2212                jres_homo(iistart+ki)=jres_homo(ki+i01)
2213                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2214                do k=1,constr_homology
2215                 odl(k,iistart+ki)=odl(k,ki+i01)
2216                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2217                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2218                enddo
2219               enddo
2220               iistart=iistart+i10-i01
2221             endif
2222           endif
2223           if (ii_in_use(ii).ne.0.and..not.liiflag) then
2224              i01=ii
2225              liiflag=.true.
2226           endif
2227          enddo
2228         enddo
2229         lim_odl=iistart-1
2230       endif
2231 !      write (iout,*) "Removing distances completed"
2232 !      call flush(iout)
2233       endif ! .not. klapaucjusz
2234
2235       if (constr_homology.gt.0) call homology_partition
2236       write (iout,*) "After homology_partition"
2237       call flush(iout)
2238       if (constr_homology.gt.0) call init_int_table
2239       write (iout,*) "After init_int_table"
2240       call flush(iout)
2241 !      endif ! .not. klapaucjusz
2242 !      endif
2243 !      if (constr_homology.gt.0) call homology_partition
2244 !      write (iout,*) "After homology_partition"
2245 !      call flush(iout)
2246 !      if (constr_homology.gt.0) call init_int_table
2247 !      write (iout,*) "After init_int_table"
2248 !      call flush(iout)
2249 !      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2250 !      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2251 !
2252 ! Print restraints
2253 !
2254       if (.not.out_template_restr) return
2255 !d      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2256       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2257        write (iout,*) "Distance restraints from templates"
2258        do ii=1,lim_odl
2259        write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
2260         ii,ires_homo(ii),jres_homo(ii),&
2261         (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
2262         ki=1,constr_homology)
2263        enddo
2264        write (iout,*) "Dihedral angle restraints from templates"
2265        do i=nnt+3,nct
2266         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2267             (rad2deg*dih(ki,i),&
2268             rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2269        enddo
2270        write (iout,*) "Virtual-bond angle restraints from templates"
2271        do i=nnt+2,nct
2272         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2273             (rad2deg*thetatpl(ki,i),&
2274             rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2275        enddo
2276        write (iout,*) "SC restraints from templates"
2277        do i=nnt,nct
2278         write(iout,'(i7,100(4f8.2,4x))') i,&
2279         (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
2280          1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2281        enddo
2282       endif
2283       return
2284       end subroutine read_constr_homology
2285 !-----------------------------------------------------------------------------
2286       subroutine read_klapaucjusz
2287       use energy_data
2288       implicit none
2289 !     include 'DIMENSIONS'
2290 !#ifdef MPI
2291 !     include 'mpif.h'
2292 !#endif
2293 !     include 'COMMON.SETUP'
2294 !     include 'COMMON.CONTROL'
2295 !     include 'COMMON.HOMOLOGY'
2296 !     include 'COMMON.CHAIN'
2297 !     include 'COMMON.IOUNITS'
2298 !     include 'COMMON.MD'
2299 !     include 'COMMON.GEO'
2300 !     include 'COMMON.INTERACT'
2301 !     include 'COMMON.NAMES'
2302       character(len=256) fragfile
2303       integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,&
2304                          ii_in_use
2305       integer, dimension(:,:), allocatable :: iresclust,inclust
2306       integer :: nclust
2307
2308       character(len=2) :: kic2
2309       character(len=24) :: model_ki_dist, model_ki_angle
2310       character(len=500) :: controlcard
2311       integer :: ki, i, j, jj,k, l, i_tmp,&
2312       idomain_tmp,&
2313       ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,&
2314       i01,i10,nnt_chain,nct_chain
2315       real(kind=8) :: distal
2316       logical :: lprn = .true.
2317       integer :: nres_temp
2318 !      integer :: ilen
2319 !      external :: ilen
2320       logical :: liiflag,lfirst
2321
2322       real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut
2323       real(kind=8), dimension (:,:), allocatable  :: rescore
2324       real(kind=8), dimension (:,:), allocatable :: rescore2
2325       character(len=24) :: tpl_k_rescore
2326       character(len=256) :: pdbfile
2327
2328 !
2329 ! For new homol impl
2330 !
2331 !     include 'COMMON.VAR'
2332 !
2333 !      write (iout,*) "READ_KLAPAUCJUSZ"
2334 !      print *,"READ_KLAPAUCJUSZ"
2335 !      call flush(iout)
2336       call getenv("FRAGFILE",fragfile)
2337       write (iout,*) "Opening", fragfile
2338       call flush(iout)
2339       open(ientin,file=fragfile,status="old",err=10)
2340 !      write (iout,*) " opened"
2341 !      call flush(iout)
2342
2343       sigma_theta=0.0
2344       sigma_d=0.0
2345       sigma_dih=0.0
2346       l_homo = .false.
2347
2348       nres_temp=nres
2349       itype_temp(:)=itype(:,1)
2350       ii=0
2351       lim_odl=0
2352
2353 !      write (iout,*) "Entering loop"
2354 !      call flush(iout)
2355
2356       DO IGR = 1,NCHAIN_GROUP
2357
2358 !      write (iout,*) "igr",igr
2359       call flush(iout)
2360       read(ientin,*) constr_homology,nclust
2361       if (start_from_model) then
2362         nmodel_start=constr_homology
2363       else
2364         nmodel_start=0
2365       endif
2366
2367       ii_old=lim_odl
2368
2369       ichain=iequiv(1,igr)
2370       nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
2371       nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
2372 !      write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
2373
2374 ! Read pdb files
2375       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
2376       if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
2377       do k=1,constr_homology
2378         read(ientin,'(a)') pdbfile
2379         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', &
2380         pdbfile(:ilen(pdbfile))
2381         pdbfiles_chomo(k)=pdbfile
2382         open(ipdbin,file=pdbfile,status='old',err=33)
2383         goto 34
2384   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
2385         pdbfile(:ilen(pdbfile))
2386         stop
2387   34    continue
2388         unres_pdb=.false.
2389 !        nres_temp=nres
2390         call readpdb_template(k)
2391         nres_chomo(k)=nres
2392 !        nres=nres_temp
2393         do i=1,nres
2394           rescore(k,i)=0.2d0
2395           rescore2(k,i)=1.0d0
2396         enddo
2397       enddo
2398 ! Read clusters
2399       do i=1,nclust
2400         read(ientin,*) ninclust(i),nresclust(i)
2401         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
2402         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
2403       enddo
2404 !
2405 ! Loop over clusters
2406 !
2407       do l=1,nclust
2408         do ll = 1,ninclust(l)
2409
2410         k = inclust(ll,l)
2411 !        write (iout,*) "l",l," ll",ll," k",k
2412         do i=1,nres
2413           idomain(k,i)=0
2414         enddo
2415         do i=1,nresclust(l)
2416           if (nnt.gt.1)  then
2417             idomain(k,iresclust(i,l)+1) = 1
2418           else
2419             idomain(k,iresclust(i,l)) = 1
2420           endif
2421         enddo
2422 !
2423 !     Distance restraints
2424 !
2425 !          ... --> odl(k,ii)
2426 ! Copy the coordinates from reference coordinates (?)
2427 !        nres_temp=nres
2428         nres=nres_chomo(k)
2429         do i=1,2*nres
2430           do j=1,3
2431             c(j,i)=chomo(j,i,k)
2432 !           write (iout,*) "c(",j,i,") =",c(j,i)
2433           enddo
2434         enddo
2435         call int_from_cart(.true.,.false.)
2436         call sc_loc_geom(.false.)
2437         do i=1,nres
2438           thetaref(i)=theta(i)
2439           phiref(i)=phi(i)
2440         enddo
2441 !        nres=nres_temp
2442         if (waga_dist.ne.0.0d0) then
2443           ii=ii_old
2444 !          do i = nnt,nct-2
2445           do i = nnt_chain,nct_chain-2
2446 !            do j=i+2,nct
2447             do j=i+2,nct_chain
2448
2449               x12=c(1,i)-c(1,j)
2450               y12=c(2,i)-c(2,j)
2451               z12=c(3,i)-c(3,j)
2452               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2453 !              write (iout,*) k,i,j,distal,dist2_cut
2454
2455             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2456                  .and. distal.le.dist2_cut ) then
2457
2458               ii=ii+1
2459               ii_in_use(ii)=1
2460               l_homo(k,ii)=.true.
2461
2462 !             write (iout,*) "k",k
2463 !             write (iout,*) "i",i," j",j," constr_homology",
2464 !     &                       constr_homology
2465               ires_homo(ii)=i+chain_border1(1,igr)-1
2466               jres_homo(ii)=j+chain_border1(1,igr)-1
2467               odl(k,ii)=distal
2468               if (read2sigma) then
2469                 sigma_odl(k,ii)=0
2470                 do ik=i,j
2471                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2472                 enddo
2473                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2474                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2475              sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2476               else
2477                 if (odl(k,ii).le.dist_cut) then
2478                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2479                 else
2480 #ifdef OLDSIGMA
2481                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2482                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2483 #else
2484                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2485                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2486 #endif
2487                 endif
2488               endif
2489               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2490             else
2491               ii=ii+1
2492 !              l_homo(k,ii)=.false.
2493             endif
2494             enddo
2495           enddo
2496         lim_odl=ii
2497         endif
2498 !
2499 !     Theta, dihedral and SC retraints
2500 !
2501         if (waga_angle.gt.0.0d0) then
2502           do i = nnt_chain+3,nct_chain
2503             iii=i+chain_border1(1,igr)-1
2504             if (idomain(k,i).eq.0) then
2505 !               sigma_dih(k,i)=0.0
2506                cycle
2507             endif
2508             dih(k,iii)=phiref(i)
2509             sigma_dih(k,iii)= &
2510                (rescore(k,i)+rescore(k,i-1)+ &
2511                            rescore(k,i-2)+rescore(k,i-3))/4.0
2512 !            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
2513 !     &       " sigma_dihed",sigma_dih(k,i)
2514             if (sigma_dih(k,iii).ne.0) &
2515              sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
2516           enddo
2517 !          lim_dih=nct-nnt-2
2518         endif
2519
2520         if (waga_theta.gt.0.0d0) then
2521           do i = nnt_chain+2,nct_chain
2522              iii=i+chain_border1(1,igr)-1
2523              if (idomain(k,i).eq.0) then
2524 !              sigma_theta(k,i)=0.0
2525               cycle
2526              endif
2527              thetatpl(k,iii)=thetaref(i)
2528              sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ &
2529                               rescore(k,i-2))/3.0
2530              if (sigma_theta(k,iii).ne.0) &
2531              sigma_theta(k,iii)=1.0d0/ &
2532              (sigma_theta(k,iii)*sigma_theta(k,iii))
2533           enddo
2534         endif
2535
2536         if (waga_d.gt.0.0d0) then
2537           do i = nnt_chain,nct_chain
2538              iii=i+chain_border1(1,igr)-1
2539                if (itype(i,1).eq.10) cycle
2540                if (idomain(k,i).eq.0 ) then
2541 !                  sigma_d(k,i)=0.0
2542                   cycle
2543                endif
2544                xxtpl(k,iii)=xxref(i)
2545                yytpl(k,iii)=yyref(i)
2546                zztpl(k,iii)=zzref(i)
2547                sigma_d(k,iii)=rescore(k,i)
2548                if (sigma_d(k,iii).ne.0) &
2549                 sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
2550 !               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
2551           enddo
2552         endif
2553       enddo ! l
2554       enddo ! ll
2555 !
2556 ! remove distance restraints not used in any model from the list
2557 ! shift data in all arrays
2558 !
2559 !      write (iout,*) "ii_old",ii_old
2560       if (waga_dist.ne.0.0d0) then
2561 #ifdef DEBUG
2562        write (iout,*) "Distance restraints from templates"
2563        do iii=1,lim_odl
2564        write(iout,'(4i5,100(2f8.2,1x,l1,4x))') &
2565         iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), &
2566         (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), &
2567         ki=1,constr_homology)
2568        enddo
2569 #endif
2570         ii=ii_old
2571         liiflag=.true.
2572         lfirst=.true.
2573         do i=nnt_chain,nct_chain-2
2574          do j=i+2,nct_chain
2575           ii=ii+1
2576 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2577 !     &            .and. distal.le.dist2_cut ) then
2578 !          write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
2579 !          call flush(iout)
2580           if (ii_in_use(ii).eq.0.and.liiflag.or. &
2581           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2582             liiflag=.false.
2583             i10=ii
2584             if (lfirst) then
2585               lfirst=.false.
2586               iistart=ii
2587             else
2588               if(i10.eq.lim_odl) i10=i10+1
2589               do ki=0,i10-i01-1
2590                ires_homo(iistart+ki)=ires_homo(ki+i01)
2591                jres_homo(iistart+ki)=jres_homo(ki+i01)
2592                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2593                do k=1,constr_homology
2594                 odl(k,iistart+ki)=odl(k,ki+i01)
2595                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2596                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2597                enddo
2598               enddo
2599               iistart=iistart+i10-i01
2600             endif
2601           endif
2602           if (ii_in_use(ii).ne.0.and..not.liiflag) then
2603              i01=ii
2604              liiflag=.true.
2605           endif
2606          enddo
2607         enddo
2608         lim_odl=iistart-1
2609       endif
2610
2611       lll=lim_odl-ii_old
2612
2613       do i=2,nequiv(igr)
2614
2615         ichain=iequiv(i,igr)
2616
2617         do j=nnt_chain,nct_chain
2618           jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
2619           do k=1,constr_homology
2620             dih(k,jj)=dih(k,j)
2621             sigma_dih(k,jj)=sigma_dih(k,j)
2622             thetatpl(k,jj)=thetatpl(k,j)
2623             sigma_theta(k,jj)=sigma_theta(k,j)
2624             xxtpl(k,jj)=xxtpl(k,j)
2625             yytpl(k,jj)=yytpl(k,j)
2626             zztpl(k,jj)=zztpl(k,j)
2627             sigma_d(k,jj)=sigma_d(k,j)
2628           enddo
2629         enddo
2630
2631         jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
2632 !        write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
2633         do j=ii_old+1,lim_odl
2634           ires_homo(j+lll)=ires_homo(j)+jj
2635           jres_homo(j+lll)=jres_homo(j)+jj
2636           do k=1,constr_homology
2637             odl(k,j+lll)=odl(k,j)
2638             sigma_odl(k,j+lll)=sigma_odl(k,j)
2639             l_homo(k,j+lll)=l_homo(k,j)
2640           enddo
2641         enddo
2642
2643         ii_old=ii_old+lll
2644         lim_odl=lim_odl+lll
2645
2646       enddo
2647
2648       ENDDO ! IGR
2649
2650       if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
2651       nres=nres_temp
2652       itype(:,1)=itype_temp(:)
2653
2654       return
2655    10 stop "Error in fragment file"
2656       end subroutine read_klapaucjusz
2657
2658 !-----------------------------------------------------------
2659       subroutine seq2chains
2660 !c
2661 !c Split the total UNRES sequence, which has dummy residues separating
2662 !c the chains, into separate chains. The length of  chain ichain is
2663 !c contained in chain_length(ichain), the first and last non-dummy
2664 !c residues are in chain_border(1,ichain) and chain_border(2,ichain),
2665 !c respectively. The lengths pertain to non-dummy residues only.
2666 !c
2667 !      implicit none
2668 !      include 'DIMENSIONS'
2669       use energy_data, only:molnum,nchain,chain_length,ireschain
2670       implicit none
2671 !      integer ireschain(nres)
2672       integer ii,ichain,i,j,mnum
2673       logical new_chain
2674       print *,"in seq2"
2675       ichain=1
2676       new_chain=.true.
2677       if (.not.allocated(chain_length)) allocate(chain_length(50))
2678       if (.not.allocated(chain_border)) allocate(chain_border(2,50))
2679
2680       chain_length(ichain)=0
2681       ii=1
2682       do while (ii.lt.nres)
2683         write(iout,*) "in seq2chains",ii,new_chain
2684         mnum=molnum(ii)
2685         if (itype(ii,mnum).eq.ntyp1_molec(mnum)) then
2686           if (.not.new_chain) then
2687             new_chain=.true.
2688             chain_border(2,ichain)=ii-1
2689             ichain=ichain+1
2690             chain_border(1,ichain)=ii+1
2691             chain_length(ichain)=0
2692           endif
2693         else
2694           if (new_chain) then
2695             chain_border(1,ichain)=ii
2696             new_chain=.false.
2697           endif
2698           chain_length(ichain)=chain_length(ichain)+1
2699         endif
2700         ii=ii+1
2701       enddo
2702       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) then
2703         ii=ii-1
2704       else
2705         chain_length(ichain)=chain_length(ichain)+1
2706       endif
2707       if (chain_length(ichain).gt.0) then
2708         chain_border(2,ichain)=ii
2709         nchain=ichain
2710       else
2711         nchain=ichain-1
2712       endif
2713       ireschain=0
2714       do i=1,nchain
2715         do j=chain_border(1,i),chain_border(2,i)
2716           ireschain(j)=i
2717         enddo
2718       enddo
2719       return
2720       end subroutine
2721 !---------------------------------------------------------------------
2722       subroutine chain_symmetry(npermchain,tabpermchain)
2723 !c
2724 !c Determine chain symmetry. nperm is the number of permutations and
2725 !c tabperchain contains the allowed permutations of the chains.
2726 !c
2727 !      implicit none
2728 !      include "DIMENSIONS"
2729 !      include "COMMON.IOUNITS"
2730       implicit none
2731       integer itemp(50),&
2732        npermchain,tabpermchain(50,5040),&
2733        tabperm(50,5040),mapchain(50),&
2734        iequiv(50,nres),iflag(nres)
2735       integer i,j,k,l,ii,nchain_group,nequiv(50),iieq,&
2736        nperm,npermc,ind,mnum
2737       if (nchain.eq.1) then
2738         npermchain=1
2739         tabpermchain(1,1)=1
2740 !c        print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
2741         return
2742       endif
2743 !c
2744 !c Look for equivalent chains
2745 #ifdef DEBUG
2746       write(iout,*) "nchain",nchain
2747       do i=1,nchain
2748         write(iout,*) "chain",i," from",chain_border(1,i),&
2749            " to",chain_border(2,i)
2750         write(iout,*)&
2751         "sequence ",(itype(j,molnum(j)),j=chain_border(1,i),chain_border(2,i))
2752       enddo
2753 #endif
2754       do i=1,nchain
2755         iflag(i)=0
2756       enddo
2757       nchain_group=0
2758       do i=1,nchain
2759         if (iflag(i).gt.0) cycle
2760         iflag(i)=1
2761         nchain_group=nchain_group+1
2762         iieq=1
2763         iequiv(iieq,nchain_group)=i
2764         do j=i+1,nchain
2765           if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
2766 !c          k=0
2767 !c          do while(k.lt.chain_length(i) .and.
2768 !c     &     itype(chain_border(1,i)+k).eq.itype(chain_border(1,j)+k))
2769           do k=0,chain_length(i)-1
2770 !c            k=k+1
2771             mnum=molnum(k)
2772             if (itype(chain_border(1,i)+k,mnum).ne.&
2773                itype(chain_border(1,j)+k,mnum)) exit
2774           enddo
2775           if (k.lt.chain_length(i)) cycle
2776           iflag(j)=1
2777           iieq=iieq+1
2778           iequiv(iieq,nchain_group)=j
2779         enddo
2780         nequiv(nchain_group)=iieq
2781       enddo
2782       write(iout,*) "Number of equivalent chain groups:",nchain_group
2783       write(iout,*) "Equivalent chain groups"
2784       do i=1,nchain_group
2785         write(iout,*) "group",i," #members",nequiv(i)," chains",&
2786            (iequiv(j,i),j=1,nequiv(i))
2787       enddo
2788       ind=0
2789       do i=1,nchain_group
2790         do j=1,nequiv(i)
2791           ind=ind+1
2792           mapchain(ind)=iequiv(j,i)
2793         enddo
2794       enddo
2795       write (iout,*) "mapchain"
2796       do i=1,nchain
2797         write (iout,*) i,mapchain(i)
2798       enddo
2799       ii=0
2800       do i=1,nchain_group
2801         call permut(nequiv(i),nperm,tabperm)
2802         if (ii.eq.0) then
2803           ii=nequiv(i)
2804           npermchain=nperm
2805           do j=1,nperm
2806             do k=1,ii
2807               tabpermchain(k,j)=iequiv(tabperm(k,j),i)
2808             enddo
2809           enddo
2810         else
2811           npermc=npermchain
2812           npermchain=npermchain*nperm
2813           ind=0
2814           do k=1,nperm
2815             do j=1,npermc
2816               ind=ind+1
2817               do l=1,ii
2818                 tabpermchain(l,ind)=tabpermchain(l,j)
2819               enddo
2820               do l=1,nequiv(i)
2821                 tabpermchain(ii+l,ind)=iequiv(tabperm(l,k),i)
2822               enddo
2823             enddo
2824           enddo
2825           ii=ii+nequiv(i)
2826         endif
2827       enddo
2828       do i=1,npermchain
2829         do j=1,nchain
2830           itemp(mapchain(j))=tabpermchain(j,i)
2831         enddo
2832         do j=1,nchain
2833           tabpermchain(j,i)=itemp(j)
2834         enddo
2835       enddo
2836       write(iout,*) "Number of chain permutations",npermchain
2837       write(iout,*) "Permutations"
2838       do i=1,npermchain
2839         write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
2840       enddo
2841       return
2842       end
2843 !c---------------------------------------------------------------------
2844       integer function tperm(i,iperm,tabpermchain)
2845 !      implicit none
2846 !      include 'DIMENSIONS'
2847       integer i,iperm
2848       integer tabpermchain(50,5040)
2849       if (i.eq.0) then
2850         tperm=0
2851       else
2852         tperm=tabpermchain(i,iperm)
2853       endif
2854       return
2855       end function
2856
2857 !-----------------------------------------------------------------------------
2858       end module io